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
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
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
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
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
00076 UgliGeomHandle ugh(vc);
00077
00078 if (! ugh.IsValid()) {
00079 MSG("NueDisplayModule",Msg::kWarning) << "Got invalid Ugli\n";
00080
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
00092
00093 double tmin=0, tmax=0;
00094 bool first = true;
00095 NtpSRSlice *slice = SntpHelpers::GetSlice(event->slc,st);
00096
00097
00098 if (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++){
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
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
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
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
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
00239
00240 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
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
00260 if(strip->planeview==PlaneView::kU){
00261
00262 }
00263 else if(strip->planeview==PlaneView::kV){
00264
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
00283
00284 st0.push_back(strip->time1-fiberDist/PropagationVelocity::Velocity());
00285
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
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
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
00320 if(strip->planeview==PlaneView::kU){
00321
00322 }
00323 else if(strip->planeview==PlaneView::kV){
00324
00325 }
00326 }
00327 }
00328 }
00329
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
00340 Float_t ShwEventPulseHeight = TimeHstShw->Integral();
00341 Float_t TrkEventPulseHeight = TimeHstTrk->Integral();
00342 Float_t TotalEventPulseHeight = TotalHst->Integral();
00343
00344
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
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
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
00429
00430
00431
00432
00433
00434 }
00435