loon: TimingVarsAna.cxx Source File

00001 #include "TH1F.h"
00002 #include "TF1.h"
00003 #include "TCanvas.h"
00004 #include "TClonesArray.h"
00005 #include "TMinuit.h"
00006 #include "StandardNtuple/NtpStRecord.h"
00007 #include "CandNtupleSR/NtpSRRecord.h"
00008 #include "CandNtupleSR/NtpSREvent.h"
00009 #include "CandNtupleSR/NtpSRTrack.h"
00010 #include "CandNtupleSR/NtpSRStrip.h"
00011 #include "CandNtupleSR/NtpSRShower.h"
00012 #include "CandNtupleSR/NtpSRSlice.h"
00013 #include "MessageService/MsgService.h"
00014 #include "NueAna/NueAnaTools/SntpHelpers.h"
00015 #include "NueAna/TimingVarsAna.h"
00016 #include "AnalysisNtuples/ANtpDefaultValue.h"
00017 #include "Calibrator/CalMIPCalibration.h"
00018 #include "Plex/PlexStripEndId.h"
00019 #include "RecoBase/PropagationVelocity.h"
00020 #include "UgliGeometry/UgliStripHandle.h"
00021 #include "UgliGeometry/UgliGeomHandle.h"
00022 //#include "Midad/Base/Mint.h"
00023 
00024 CVSID("$Id: TimingVarsAna.cxx,v 1.9 2009/07/03 14:45:34 vahle Exp $");
00025 
00026 #include "DatabaseInterface/DbiResultPtr.tpl"
00027 
00028 TimingVarsAna::TimingVarsAna(TimingVars &tv):
00029    fTimingVars(tv)
00030 {
00031 
00032 //    TimeHst = new THStack("TimeHist","Digit Times(ns) Red - Track, Blue - Shower.");
00033     TimeHstTrk = new TH1F();
00034     TimeHstTrk->SetName("TimeHstTrk");
00035     TimeHstShw = new TH1F();
00036     TimeHstShw->SetName("TimeHstShw");
00037 
00038     TotalHst = new TH1F();
00039     TotalHst->SetName("TotalHst");
00040 
00041 //    plot = new TCanvas("plot","plot",200,10,900,600);
00042 
00043 }
00044 
00045 
00046 TimingVarsAna::~TimingVarsAna()
00047 {
00048   if(TimeHstTrk) delete TimeHstTrk;
00049   if(TimeHstShw) delete TimeHstShw;
00050   if(TotalHst) delete TotalHst;
00051 }
00052 
00053 void TimingVarsAna::Analyze(int evtn, RecRecordImp<RecCandHeader> *srobj){
00054   if(srobj==0){
00055     return;
00056   }
00057   if(((dynamic_cast<NtpStRecord *>(srobj))==0)&&
00058      ((dynamic_cast<NtpSRRecord *>(srobj))==0)){
00059     return;
00060   }
00061 
00062   //Reset(srobj->GetHeader().GetSnarl(),evtn);
00063   NtpStRecord *st = dynamic_cast<NtpStRecord *>(srobj);
00064   NtpSREvent *event = SntpHelpers::GetEvent(evtn,srobj);
00065 
00066  
00067   if(!event){
00068       MSG("ShwfitAna",Msg::kError)<<"Couldn't get event "<<evtn
00069                                    <<" from Snarl "<<srobj->GetHeader().GetSnarl()<<endl;
00070       return;
00071    }
00072 
00073    VldContext vc=st->GetHeader().GetVldContext();
00074 
00075    //VldTimeStamp vldts = vc.GetTimeStamp();
00076    UgliGeomHandle ugh(vc);
00077    //UgliGeomHandle ugh = gMint->GetUgliGeomHandle();
00078    if (! ugh.IsValid()) {
00079      MSG("NueDisplayModule",Msg::kWarning) << "Got invalid Ugli\n";
00080      //return;
00081    }
00082 
00083   fDetectorType = vc.GetDetector();
00084   ReleaseType::Release_t fRel;
00085   string relName = st->GetTitle();
00086   string reco = relName.substr(0,relName.find_first_of("("));
00087   if(reco == "CEDAR") fRel = ReleaseType::kCedar;
00088   if(reco == "DOGWOOD") fRel = ReleaseType::kDogwood;
00089   else fRel = ReleaseType::kBirch;    
00090   int ntrks = event->ntrack;
00091       //timing histogram
00092     //find range
00093     double tmin=0, tmax=0;
00094     bool first = true;
00095     NtpSRSlice *slice = SntpHelpers::GetSlice(event->slc,st);
00096 
00097     //shamelessly stolen from Mad
00098     if (slice){//slice
00099       float highest_plane = 0;
00100       float lowest_plane = 500;
00101       float highest_strip0 = 0;
00102       float lowest_strip0 = 192;
00103       float highest_strip1 = 0;
00104       float lowest_strip1 = 192;
00105 
00106       float highest_z = 0.;
00107       float lowest_z = 30.;
00108       float highest_t0 = -4.0;
00109       float lowest_t0 = 4.0;
00110       float highest_t1 = -4.0;
00111       float lowest_t1 = 4.0;
00112 
00113       for (int i = 0; i<slice->nstrip; i++){//loop over strips
00114         NtpSRStrip *strip = SntpHelpers::GetStrip(slice->stp[i],srobj);
00115         if(strip == 0) continue;
00116         double t = strip->time0;
00117         if (t>-100){
00118           if (first){
00119             tmin = tmax = t;
00120             first = false;
00121           }
00122           if (t < tmin) tmin = t;
00123           if (t > tmax) tmax = t;
00124         }
00125         t = strip->time1;
00126         if (t>-100){
00127           if (first){
00128             tmin = tmax = t;
00129             first = false;
00130           }
00131           if (t < tmin) tmin = t;
00132           if (t > tmax) tmax = t;
00133         }
00134         int tempo_pln = strip->plane;
00135         int tempo_stp = strip->strip;
00136         float tempo_tpos = strip->tpos;
00137         if(tempo_pln<lowest_plane) {
00138           lowest_plane=tempo_pln;
00139           lowest_z=strip->z;
00140         }
00141         if(tempo_pln>highest_plane) {
00142           highest_plane=tempo_pln;
00143           highest_z=strip->z;
00144         }
00145 
00146         if (strip->planeview==PlaneView::kU){
00147          // fSlcUZ->Fill(strip->z,strip->tpos,(strip->ph0.sigcor + strip->ph1.sigcor)/SIGMAPMEU);
00148           if(tempo_tpos<lowest_t0) {
00149             lowest_strip0=tempo_stp;
00150             lowest_t0=tempo_tpos;
00151           }
00152           if(tempo_tpos>highest_t0) {
00153             highest_strip0=tempo_stp;
00154             highest_t0=tempo_tpos;
00155           }
00156         }
00157         else if (strip->planeview==PlaneView::kV){
00158           //fSlcVZ->Fill(strip->z,strip->tpos,(strip->ph0.sigcor + strip->ph1.sigcor)/SIGMAPMEU);
00159           if(tempo_tpos<lowest_t1) {
00160             lowest_strip1=tempo_stp;
00161             lowest_t1=tempo_tpos;
00162           }
00163           if(tempo_tpos>highest_t1) {
00164             highest_strip1=tempo_stp;
00165             highest_t1=tempo_tpos;
00166           }
00167         }
00168       }
00169       if(lowest_plane-10>=0) {
00170         lowest_plane-=10;
00171         lowest_z-=10.*0.06;
00172       }
00173       else {
00174         lowest_plane=0;
00175         lowest_z=0.;
00176       }
00177 
00178       if(lowest_strip0-5>=0) {
00179         lowest_strip0-=5;
00180         lowest_t0-=5.*0.041;
00181       }
00182       else {
00183         lowest_strip0=0;
00184         lowest_t0=-4.0;
00185       }
00186 
00187       if(lowest_strip1-5>=0) {
00188         lowest_strip1-=5;
00189         lowest_t1-=5.*0.041;
00190       }
00191       else {
00192         lowest_strip1=0;
00193         lowest_t1=-4.0;
00194       }
00195 
00196       if(highest_plane+10<=485) {
00197         highest_plane+=10;
00198         highest_z+=10.*0.06;
00199       }
00200       else {
00201         highest_plane=485;
00202         highest_z=30.;
00203       }
00204 
00205       if(highest_strip0+5<=191) {
00206         highest_strip0+=5;
00207         highest_t0+=5.*0.041;
00208       }
00209       else {
00210         highest_strip0=191;
00211         highest_t0=4.0;
00212       }
00213 
00214       if(highest_strip1+5<=191) {
00215         highest_strip1+=5;
00216         highest_t1+=5.*0.041;
00217       }
00218       else {
00219         highest_strip1=191;
00220         highest_t1=4.0;
00221       }
00222     }
00223     // give some buffer at either end...
00224     tmin -= 50e-9;
00225     tmax += 20e-9;
00226 
00227     double eps = 1.0e-8;
00228     if (tmin == tmax) { tmin -= eps; tmax += eps; }
00229 
00230     // for far detector, suppress display of pre-trigger time interval
00231     if (fDetectorType == Detector::kFar){
00232       if ((tmax - tmin)*1e9>1500) tmin = tmax - 1500e-9;
00233     }
00234 
00235     TimeHstTrk->SetBins(50,0,100);
00236     TimeHstShw->SetBins(50,0,100);
00237     TotalHst->SetBins(50,0,100);
00238     //record track hits
00239     //int ntrks = event->ntrack;
00240     if (ntrks){//if(ntrks)
00241       int trkidx = -1;
00242       int trkplanes = -1;
00243       for (int i = 0; i<ntrks; i++){
00244         int index = SntpHelpers::GetTrackIndex(i,event);
00245         NtpSRTrack *track = SntpHelpers::GetTrack(index,st);
00246         if(track == 0) continue;
00247         if (track->plane.n>trkplanes){
00248           trkplanes = track->plane.n;
00249           trkidx = index;
00250         }
00251         for (int istp = 0; istp<track->nstrip; istp++){
00252           //ftrkshw->Fill(track->stpx[istp],track->stpy[istp],1);
00253           NtpSRStrip *strip = SntpHelpers::GetStrip(track->stp[istp],st);
00254           if(strip == 0) continue;
00255           TimeHstTrk->Fill((track->stpt0[istp]-tmin-(track->ds-track->stpds[istp])/3e8)/1e-9,strip->ph0.pe);
00256           TimeHstTrk->Fill((track->stpt1[istp]-tmin-(track->ds-track->stpds[istp])/3e8)/1e-9,strip->ph1.pe);
00257           TotalHst->Fill((track->stpt0[istp]-tmin-(track->ds-track->stpds[istp])/3e8)/1e-9,strip->ph0.pe);
00258           TotalHst->Fill((track->stpt1[istp]-tmin-(track->ds-track->stpds[istp])/3e8)/1e-9,strip->ph1.pe);
00259           //record track strip information
00260           if(strip->planeview==PlaneView::kU){
00261             //fTrkUZ->Fill(strip->z,strip->tpos,100000);//paranoia
00262           }
00263           else if(strip->planeview==PlaneView::kV){
00264             //fTrkVZ->Fill(strip->z,strip->tpos,100000);
00265           }
00266         }
00267       }
00268       vector<double> spathLength;
00269       vector<double> st0;
00270       if (trkidx != -1 && fDetectorType == Detector::kNear){
00271         NtpSRTrack *track = SntpHelpers::GetTrack(trkidx,st);
00272         for (int istp = 0; istp<track->nstrip; istp++){
00273           NtpSRStrip *strip = SntpHelpers::GetStrip(track->stp[istp],st);
00274           if(strip == 0) continue;
00275           spathLength.push_back(track->ds-track->stpds[istp]);
00276           PlexStripEndId seid(Detector::kNear,strip->plane,strip->strip,StripEnd::kWest);
00277           UgliStripHandle stripHandle = ugh.GetStripHandle(seid);
00278           float halfLength = stripHandle.GetHalfLength();
00279           const TVector3 ghitxyz(track->stpx[istp],track->stpy[istp],track->stpz[istp]);
00280           TVector3 lhitxyz = stripHandle.GlobalToLocal(ghitxyz);
00281           float fiberDist = (halfLength - lhitxyz.x() + stripHandle.ClearFiber(StripEnd::kWest) + stripHandle.WlsPigtail(StripEnd::kWest));
00282           //using strip time, I don't understand track time TJ
00283           //st0.push_back(track->stpt1[istp]-fiberDist/PropagationVelocity::Velocity());
00284           st0.push_back(strip->time1-fiberDist/PropagationVelocity::Velocity());
00285           //cout<<strip->plane<<" "<<strip->strip<<" "<<track->ds-track->stpds[istp]<<Form(" %.9f %.9f %.9f",track->stpt1[istp],strip->time1,strip->time1-fiberDist/PropagationVelocity::Velocity())<<endl;
00286         }
00287         float trms1 = 0;
00288         float trms2 = 0;
00289 
00290         trms1/=st0.size();
00291         trms2/=st0.size();
00292         trms1=sqrt(trms1)*1e9;
00293         trms2=sqrt(trms2)*1e9;
00294       }
00295     }
00296 
00297     //record shower hits
00298     int nshws = event->nshower;
00299     if (nshws){
00300       for (int i = 0; i<nshws; i++){
00301         int index = SntpHelpers::GetShowerIndex(i,event);
00302         NtpSRShower *shower = SntpHelpers::GetShower(index,st);
00303         if(shower == 0) continue;
00304         for (int istp = 0; istp<shower->nstrip; istp++){
00305           NtpSRStrip *strip = SntpHelpers::GetStrip(shower->stp[istp],st);
00306           if(strip == 0) continue;
00307           if(ReleaseType::IsCedar(fRel)||ReleaseType::IsDogwood(fRel)){
00308             TimeHstShw->Fill((shower->stpt0[istp]-tmin)/1e-9,strip->ph0.pe);
00309             TimeHstShw->Fill((shower->stpt1[istp]-tmin)/1e-9,strip->ph1.pe);
00310             TotalHst->Fill((shower->stpt0[istp]-tmin)/1e-9,strip->ph0.pe);
00311             TotalHst->Fill((shower->stpt1[istp]-tmin)/1e-9,strip->ph1.pe);
00312             //cout<<"Bin: "<<(shower->stpt1[istp]-tmin)/1e-9<<" Weight: "<<strip->ph0.pe+strip->ph1.pe<<endl;
00313           }else {
00314             TimeHstShw->Fill((strip->time0-tmin)/1e-9,strip->ph0.pe);
00315             TimeHstShw->Fill((strip->time1-tmin)/1e-9,strip->ph1.pe);
00316             TotalHst->Fill((strip->time0-tmin)/1e-9,strip->ph0.pe);
00317             TotalHst->Fill((strip->time1-tmin)/1e-9,strip->ph1.pe);
00318           }
00319           //record shower strip information
00320           if(strip->planeview==PlaneView::kU){
00321           //  fShwUZ->Fill(strip->z,strip->tpos,100000);//paranoia
00322           }
00323           else if(strip->planeview==PlaneView::kV){
00324            // fShwVZ->Fill(strip->z,strip->tpos,100000);
00325           }
00326         }
00327       }
00328     }
00329 //Max bins and contents
00330     Int_t ShwMaxBin = TimeHstShw->GetMaximumBin();
00331     Float_t ShwMaxBinCont = TimeHstShw->GetBinContent(ShwMaxBin);
00332 
00333     Int_t TrkMaxBin = TimeHstTrk->GetMaximumBin();
00334     Float_t TrkMaxBinCont = TimeHstTrk->GetBinContent(TrkMaxBin);
00335 
00336     Int_t TotalMaxBin = TotalHst->GetMaximumBin();
00337     Float_t TotalMaxBinCont = TotalHst->GetBinContent(TotalMaxBin);
00338 
00339 //Total Pulse height
00340     Float_t ShwEventPulseHeight = TimeHstShw->Integral();
00341     Float_t TrkEventPulseHeight = TimeHstTrk->Integral();
00342     Float_t TotalEventPulseHeight = TotalHst->Integral();
00343 
00344 //Second max bin contents
00345     Float_t ShwSecondMaxBin;
00346     if(TimeHstShw->GetBinContent(ShwMaxBin+1)>TimeHstShw->GetBinContent(ShwMaxBin-1)){
00347       ShwSecondMaxBin = TimeHstShw->GetBinContent(ShwMaxBin+1);
00348     }else{
00349       ShwSecondMaxBin = TimeHstShw->GetBinContent(ShwMaxBin-1);
00350     }
00351     Float_t TrkSecondMaxBin;
00352     if(TimeHstTrk->GetBinContent(TrkMaxBin+1)>TimeHstTrk->GetBinContent(TrkMaxBin-1)){
00353       TrkSecondMaxBin = TimeHstTrk->GetBinContent(TrkMaxBin+1);
00354     }else{
00355       TrkSecondMaxBin = TimeHstTrk->GetBinContent(TrkMaxBin-1);
00356     }
00357     Float_t TotalSecondMaxBin;
00358     if(TotalHst->GetBinContent(TotalMaxBin+1)>TotalHst->GetBinContent(TotalMaxBin-1)){
00359       TotalSecondMaxBin = TotalHst->GetBinContent(TotalMaxBin+1);
00360     }else{
00361       TotalSecondMaxBin = TotalHst->GetBinContent(TotalMaxBin-1);
00362     }
00363 
00364 //Percent of pulse height in Max 2 bins
00365     Float_t TrkPercentofTotal = -9999.99;
00366     Float_t ShwPercentofTotal = -9999.99;
00367     Float_t TotalPercentofTotal = -9999.99;
00368 
00369     if (ShwEventPulseHeight>0){
00370       ShwPercentofTotal = (ShwSecondMaxBin + ShwMaxBinCont)/(ShwEventPulseHeight);
00371     }
00372     if (TrkEventPulseHeight>0){
00373       TrkPercentofTotal = (TrkSecondMaxBin + TrkMaxBinCont)/(TrkEventPulseHeight);
00374     }
00375     if (TotalEventPulseHeight>0){
00376       TotalPercentofTotal = (TotalSecondMaxBin + TotalMaxBinCont)/(TotalEventPulseHeight);
00377     }
00378 
00379 //how many of the two bins prior to the max bin are at least 20% of the max bin ph
00380     Int_t ShwbinsPrior = 0;
00381     if(TimeHstShw->GetBinContent(ShwMaxBin-1)>(0.2*ShwMaxBinCont)){ShwbinsPrior++;}
00382     if(TimeHstShw->GetBinContent(ShwMaxBin-2)>(0.2*ShwMaxBinCont)){ShwbinsPrior++;}
00383     Int_t TrkbinsPrior = 0;
00384     if(TimeHstTrk->GetBinContent(TrkMaxBin-1)>(0.2*TrkMaxBinCont)){TrkbinsPrior++;}
00385     if(TimeHstTrk->GetBinContent(TrkMaxBin-2)>(0.2*TrkMaxBinCont)){TrkbinsPrior++;}
00386     Int_t TotalbinsPrior = 0;
00387     if(TotalHst->GetBinContent(TotalMaxBin-1)>(0.2*TotalMaxBinCont)){TotalbinsPrior++;}
00388     if(TotalHst->GetBinContent(TotalMaxBin-2)>(0.2*TotalMaxBinCont)){TotalbinsPrior++;}
00389     
00391     Int_t ShwBinsOver40 = 0;
00392     Int_t TrkBinsOver40 = 0;
00393     Int_t TotalBinsOver40 = 0;
00394     for (int k=0;k<50;k++){
00395       if(TimeHstShw->GetBinContent(k)>20)ShwBinsOver40++;
00396       if(TimeHstTrk->GetBinContent(k)>20)TrkBinsOver40++;
00397       if(TotalHst->GetBinContent(k)>20)TotalBinsOver40++;
00398     }
00399 
00400     fTimingVars.ShwMaxBin = ShwMaxBin;
00401     fTimingVars.TrkMaxBin = TrkMaxBin;
00402     fTimingVars.TotalMaxBin = TotalMaxBin;
00403 
00404     fTimingVars.ShwMaxBinCont = ShwMaxBinCont;
00405     fTimingVars.TrkMaxBinCont = TrkMaxBinCont;
00406     fTimingVars.TotalMaxBinCont = TotalMaxBinCont;
00407 
00408     fTimingVars.ShwSecondMaxBin = ShwSecondMaxBin;
00409     fTimingVars.TrkSecondMaxBin = TrkSecondMaxBin;
00410     fTimingVars.TotalSecondMaxBin = TotalSecondMaxBin;
00411 
00412     fTimingVars.ShwEventPulseHeight = ShwEventPulseHeight;
00413     fTimingVars.TrkEventPulseHeight = TrkEventPulseHeight;
00414     fTimingVars.TotalEventPulseHeight = TotalEventPulseHeight;
00415 
00416     fTimingVars.ShwPercentofTotal = ShwPercentofTotal;
00417     fTimingVars.TrkPercentofTotal = TrkPercentofTotal;
00418     fTimingVars.TotalPercentofTotal = TotalPercentofTotal;
00419 
00420     fTimingVars.ShwbinsPrior = ShwbinsPrior;
00421     fTimingVars.TrkbinsPrior = TrkbinsPrior;
00422     fTimingVars.TotalbinsPrior = TotalbinsPrior;
00423 
00424     fTimingVars.ShwBinsOver40 = ShwBinsOver40;
00425     fTimingVars.TrkBinsOver40 = TrkBinsOver40;
00426     fTimingVars.TotalBinsOver40 = TotalBinsOver40;
00427     /*
00428     char name[100];
00429     sprintf(name,"/home/danche/ViewScan5/plot/ShwHst_snarl_%d.png",st->GetHeader().GetSnarl());
00430 
00431     TimeHstShw->Draw();
00432     plot->Print(name);
00433     */
00434 }
00435