loon: AlgShieldPlank.cxx Source File

00001 
00002 //
00003 // $Id: AlgShieldPlank.cxx,v 1.2 2005/05/27 21:03:53 rhatcher Exp $ 
00004 //
00005 // AlgShieldPlank
00006 //
00007 // Package: ShieldPlank
00008 //
00009 // First go at alg for making shield strips
00010 //
00011 // Contact: Caius Howcroft<howcroft@hep.phy.cam.ac.uk>
00012 //
00013 //
00015 
00016 
00017 #include <cassert>
00018 #include <vector>
00019 
00020 
00021 #include "Algorithm/AlgFactory.h"
00022 #include "Algorithm/AlgHandle.h"
00023 #include "Algorithm/AlgConfig.h"
00024 #include "Candidate/CandContext.h"
00025 
00026 #include "ShieldPlank/AlgShieldPlank.h"
00027 #include "ShieldPlank/CandShieldPlank.h"
00028 #include "ShieldPlank/CandShieldPlankHandle.h"
00029 #include "CandDigit/CandDigitHandle.h"
00030 #include "CandDigit/CandDigitListHandle.h"
00031 
00032 #include "Conventions/Mphysical.h"
00033 #include "Conventions/Munits.h"
00034 #include "Conventions/PlaneView.h"
00035 #include "MessageService/MsgService.h"
00036 #include "MinosObjectMap/MomNavigator.h"
00037 #include "Navigation/NavKey.h"
00038 #include "Navigation/NavSet.h"
00039 
00040 #include "UgliGeometry/UgliGeomHandle.h"
00041 #include "Validity/VldContext.h"
00042 
00043 #include "Conventions/Munits.h"
00044 
00045 #define SUBTRACT_ABSTIME_NOT
00046 
00047 
00048 CVSID("$Id: AlgShieldPlank.cxx,v 1.2 2005/05/27 21:03:53 rhatcher Exp $");
00049 //______________________________________________________________________
00050 AlgShieldPlank::AlgShieldPlank()
00051 {
00052 }
00053 
00054 //______________________________________________________________________
00055 AlgShieldPlank::~AlgShieldPlank()
00056 {
00057 }
00058 
00059 //______________________________________________________________________
00060 void AlgShieldPlank::RunAlg(AlgConfig& /* ac */, CandHandle& ch,
00061                             CandContext& cx)
00062 {
00063   MSG("ShieldPlank",Msg::kDebug) << "Start ALGSHIELDSTRIP"<<endl;
00064   //======================================================
00065   //Set up all pointer junk
00066   //
00067   assert(cx.GetDataIn());
00068   if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) {
00069     MSG("ShieldPlank",Msg::kError) << "AlgShieldPlank wasnt passed a CandDigit"<<endl;
00070     return;
00071   }
00072   
00073   const TObjArray *tary =
00074     dynamic_cast<const TObjArray*>(cx.GetDataIn());
00075   
00076 
00077   CandShieldPlankHandle& cssh = (CandShieldPlankHandle&) ch;
00078   const CandRecord *candrec = cx.GetCandRecord();
00079   assert(candrec);
00080   const VldContext *vldcptr = candrec->GetVldContext();
00081   assert(vldcptr);
00082   VldContext vldc = *vldcptr;
00083   UgliGeomHandle ugh(vldc);
00084 
00085   //some constants
00086   //
00087   const int minend = min(StripEnd::kNegative,StripEnd::kPositive);
00088     
00089   //======================================================
00090   //
00091   //
00092   double meanX[3]={0.0, 0.0, 0.0};
00093   double meanY[3]={0.0, 0.0, 0.0};
00094   double meanZ[3]={0.0, 0.0, 0.0};
00095   int    ndigit  = 0;
00096   double chargePE[3]={0.0, 0.0, 0.0};
00097   double meantime[3]= {0.,0.,0.};
00098   int    nstrip[3] ={0, 0, 0};
00099   double mean_corrected_time[3]={0., 0., 0.};
00100 
00101   int section       = -1;
00102   int subsection       = -1;
00103 
00104   int uglierrors = 0;
00105   double  dtoff = 0;
00106 
00107   int plane =-1;
00108   
00109   for (Int_t i=0; i<=tary->GetLast(); i++) {
00110     TObject *tobj = tary->At(i);
00111     assert(tobj->InheritsFrom("CandDigitHandle"));
00112     CandDigitHandle *cdh = dynamic_cast<CandDigitHandle*>(tobj);
00113     cssh.AddDaughterLink(*cdh);
00114     ndigit++;
00115     MSG("ShieldPlank",Msg::kDebug) << "Adding digit to strip"<<endl; 
00116     if (!i) {
00117       const CandDigitHandle *constcdh = dynamic_cast<const CandDigitHandle*>(tobj);
00118       const CandDigitListHandle *cdlh = dynamic_cast<const CandDigitListHandle*>(constcdh->GetMother());
00119       assert(cdlh);
00120 #ifdef SUBTRACT_ABSTIME
00121       dtoff = cdlh->GetAbsTime();
00122 #endif //SUBTRACT_ABSTIME
00123 
00124       MSG("ShieldPlank",Msg::kDebug) << "Setting time offset to:" <<dtoff<<endl;
00125     }
00126     const PlexSEIdAltL& psealt =  cdh->GetPlexSEIdAltL();
00127     plane = psealt.GetPlane();
00128     double meanfiberlength = 0;
00129     const int stpend =  (int)(psealt.GetEnd())-minend;
00130     double digittime = cdh->GetTime(CalTimeType::kNone)-dtoff;   
00131     double qpe = cdh->GetCharge(CalDigitType::kPE);
00132     const double qraw = cdh->GetCharge(CalDigitType::kNone);
00133     if(qraw>0 && qpe/qraw > 0.1)
00134       {
00135         static int msgrpt = 0;
00136         if(msgrpt<10){
00137           MSG("ShieldPlank",Msg::kError) << "Uncalibrated shield charge ADC:"<< qraw << "PE:"<<qpe<<" use empirical clibration of 90 pe/adc"<<endl;
00138           ++msgrpt;
00139         }else if(msgrpt==10){
00140           MSG("ShieldPlank",Msg::kError) << "....Last warning about uncalibrated charge"<<endl;  
00141           ++msgrpt;
00142         }
00143           
00144         qpe = qraw / 90.;
00145       }
00146 
00147     MSG("ShieldPlank",Msg::kDebug) << "Find mean strip co-ords"<<endl;
00148     for(PlexSEIdAltL::const_iterator iter = psealt.begin();
00149         iter != psealt.end(); ++iter)
00150       {
00151         const PlexStripEndId pseid( (*iter).GetSEId() );
00152         const UgliStripHandle ush(ugh.GetStripHandle(pseid));
00153         if(ush.IsValid())
00154           {
00155             
00156             const double halflength = ush.GetHalfLength();          
00157             const TVector3 gpos0(ush.GlobalPos(0.0));
00158             const TVector3 gpos1(ush.GlobalPos(halflength));
00159             const TVector3 gpos2(ush.GlobalPos(-halflength));
00160 
00161             const TVector3 &gposN =((gpos1.Z()>gpos2.Z())?gpos1:gpos2);
00162             const TVector3 &gposS =((gpos1.Z()>gpos2.Z())?gpos2:gpos1);
00163             MSG("ShieldPlank",Msg::kDebug) << "  Strip Z co-ord 0, +, -:"<<gpos0.Z()<<", "<< gpos1.Z()<< ", "<<gpos2.Z()<<endl;
00164             MSG("ShieldPlank",Msg::kDebug) << "  Strip Z co-ord N,S:"<<gposN.Z()<<", "<<gposS.Z()<<endl;
00165             MSG("ShieldPlank",Msg::kDebug) << "  Strip co-ord:"<<gpos0.X()<<", "<< gpos0.Y()<< ", "<<gposN.Z()<<"/"<<gposS.Z()<<endl;
00166             
00167             meanfiberlength +=  (ush.ClearFiber(psealt.GetEnd())+ ush.WlsPigtail(psealt.GetEnd()));
00168             
00169             if(psealt.GetEnd()==StripEnd::kNorth){
00170               meanX[stpend] += gposN.X();
00171               meanY[stpend] += gposN.Y();
00172             }else{
00173               meanX[stpend] += gposS.X();
00174               meanY[stpend] += gposS.Y();
00175             }
00176             meanX[2] += gpos0.X();
00177             meanY[2] += gpos0.Y();
00178             
00179             meanZ[StripEnd::kNorth - minend] += gposN.Z();
00180             meanZ[StripEnd::kSouth - minend] += gposS.Z();          
00181             meanZ[2] += gpos0.Z();
00182             
00183             ++nstrip[stpend];
00184             ++nstrip[2];
00185             
00186           }else{
00187             MSG("ShieldPlank",Msg::kError)<< "No Valid UgliStripHandle for plane "<< (*iter).GetSEId().GetPlane()<<endl;
00188             ++uglierrors;
00189           }
00190       }//iterate over possible strip ends
00191     if(psealt.size())
00192       {
00193         meanfiberlength = meanfiberlength/((double)psealt.size());
00194       }
00195     MSG("ShieldPlank",Msg::kDebug)<< "Mean fiber length: "<< meanfiberlength<<endl;
00196     //
00197     //Calc the time corrected for delay of the fiber
00198     //
00199     double corrected_digit_time = digittime - (meanfiberlength / (Munits::c_light/1.77)) ;
00200     MSG("ShieldPlank",Msg::kDebug)<< "RawTime: "<< digittime<< " Corrected time: " << corrected_digit_time<<endl;
00201 
00202     MSG("ShieldPlank",Msg::kDebug) << "Adding in charge: "<<qpe<<" at time "<< digittime<< endl;
00203     chargePE[stpend]+= qpe;
00204     chargePE[2]+= qpe;
00205     meantime[stpend]+= qpe*digittime;
00206     meantime[2]+= qpe*digittime;
00207     mean_corrected_time[stpend] += qpe*corrected_digit_time;
00208     mean_corrected_time[2] += qpe*corrected_digit_time;
00209     
00210   }//iterator over digits
00211 
00212   for(unsigned int i=0; i<3; ++i)
00213     {
00214       if(nstrip[i])
00215         {
00216           meanX[i]=meanX[i] / (double)nstrip[i];
00217           meanY[i]=meanY[i] / (double)nstrip[i];
00218         }
00219       if(nstrip[2]) meanZ[i]=meanZ[i] / (double)nstrip[2];
00220       if(chargePE[i]){
00221         meantime[i] = meantime[i]/chargePE[i];
00222         mean_corrected_time[i] = mean_corrected_time[i] / chargePE[i];
00223       }
00224       MSG("ShieldPlank",Msg::kDebug) << "Final XYZT:" <<         meanX[i]<< ", "<< meanY[i]<<", "<<meanZ[i]<<", "<< meantime[i]<< ", "<< mean_corrected_time[i]<<endl;
00225     }
00226 
00227   //decide on what section it is and what bit
00228   if(plane>=528 && plane<=575){
00229     section=1;
00230     if(plane>=543 && plane<=557) subsection = 0;
00231     if(plane>=528 && plane<=542) subsection = 0;
00232     if(plane>=572 && plane<=575) subsection = 1;
00233     if(plane>=563 && plane<=566) subsection =-1;
00234     if(plane>=567 && plane<=568) subsection = 2;
00235     if(plane>=558 && plane<=559) subsection =-2;
00236   }
00237   if(plane>=592 && plane<=639){
00238     section=2;
00239     if(plane>=607 && plane<=621) subsection = 0; //top upper later
00240     if(plane>=592 && plane<=606) subsection = 0; //top bottom layer
00241     if(plane>=636 && plane<=639) subsection = 1;
00242     if(plane>=627 && plane<=630) subsection =-1;
00243     if(plane>=631 && plane<=632) subsection = 2;
00244     if(plane>=622 && plane<=623) subsection =-2;
00245   }
00246   if(plane>=656 && plane<=703){
00247     section=3;
00248     if(plane>=671 && plane<=685) subsection = 0;
00249     if(plane>=656 && plane<=670) subsection = 0;
00250     if(plane>=700 && plane<=703) subsection =-1;
00251     if(plane>=691 && plane<=694) subsection = 1;
00252     if(plane>=695 && plane<=696) subsection = 2;
00253     if(plane>=686 && plane<=687) subsection =-2;
00254   }
00255   if(plane>=720 && plane<=767){
00256     section=4;
00257     if(plane>=735 && plane<=749) subsection = 0;
00258     if(plane>=720 && plane<=734) subsection = 0;
00259     if(plane>=764 && plane<=797) subsection = 1;
00260     if(plane>=755 && plane<=758) subsection =-1;
00261     if(plane>=759 && plane<=760) subsection = 2;
00262     if(plane>=750 && plane<=751) subsection =-2;
00263     
00264   }
00265   MSG("ShieldPlank",Msg::kDebug) << "Setting section to:"<<section<<", "<<subsection<<endl;    
00266   cssh.SetSection(section);
00267   cssh.SetSubSection(subsection);
00268   cssh.SetMeanX(meanX[2]);
00269   cssh.SetMeanY(meanY[2]);
00270   cssh.SetMeanZ(meanZ);
00271   cssh.SetMeanT(meantime);
00272   cssh.SetMeanCorrectedT(mean_corrected_time);
00273   MSG("ShieldPlank",Msg::kDebug) << " Set Corrected T: "<< cssh.GetMeanCorrectedTime(StripEnd::kNorth)<< ", "<<cssh.GetMeanCorrectedTime(StripEnd::kSouth)<<", "<<cssh.GetMeanCorrectedTime(StripEnd::kWhole)<<endl;
00274   cssh.SetNErrors(uglierrors);
00275   MSG("ShieldPlank",Msg::kDebug) << " setting charge:"<< chargePE[0]<<", "<< chargePE[1]<<", "<<chargePE[2]<<endl;
00276   cssh.SetChargePE(chargePE);
00277   MSG("ShieldPlank",Msg::kDebug) << "Done Shield Strip"<<endl;
00278 
00279  
00280 }
00281 
00282 
00283 //______________________________________________________________________
00284 void AlgShieldPlank::Trace(const char* /*c*/) const
00285 {
00286   //yet more pointless code.
00287 
00288 }
00289 
00290 ClassImp(AlgShieldPlank)
00291