00001
00002
00003
00004
00005
00007
00008 #include <set>
00009 #include <cmath>
00010
00011 #include "TClonesArray.h"
00012 #include "TProfile.h"
00013 #include "TVector3.h"
00014
00015 #include "AtNuEvent/AtmosEvent.h"
00016 #include "AtNuEvent/AtmosStrip.h"
00017 #include "AtNuEvent/AtmosTrack.h"
00018 #include "AtNuEvent/AtmosScintHit.h"
00019 #include "BeamDataUtil/BeamMonSpill.h"
00020 #include "BeamDataUtil/BMSpillAna.h"
00021 #include "MessageService/MsgService.h"
00022 #include "Conventions/Munits.h"
00023 #include "BeamDataNtuple/NtpBDLiteRecord.h"
00024 #include "MCNtuple/NtpMCDigiScintHit.h"
00025 #include "CandNtupleSR/NtpSREvent.h"
00026 #include "CandNtupleSR/NtpSRShower.h"
00027 #include "CandNtupleSR/NtpSRSlice.h"
00028 #include "CandNtupleSR/NtpSRStrip.h"
00029 #include "CandNtupleSR/NtpSRTrack.h"
00030 #include "StandardNtuple/NtpStRecord.h"
00031 #include "Conventions/PlaneView.h"
00032 #include "Conventions/Detector.h"
00033 #include "UgliGeometry/UgliGeomHandle.h"
00034
00035 #include "MeuCal/MeuCuts.h"
00036 #include "MeuCal/MeuHitInfo.h"
00037 #include "MeuCal/MeuReco.h"
00038 #include "MeuCal/MeuSummary.h"
00039
00040 using std::endl;
00041 using std::cout;
00042 using std::map;
00043 using std::vector;
00044
00045
00046
00047 CVSID("$Id: MeuCuts.cxx,v 1.18 2007/11/11 07:03:42 rhatcher Exp $");
00048
00049
00050
00051 MeuCuts::MeuCuts()
00052 {
00053 MSG("MeuCuts",Msg::kDebug)
00054 <<"Running MeuCuts Constructor..."<<endl;
00055
00056
00057 MSG("MeuCuts",Msg::kDebug)
00058 <<"Finished MeuCuts Constructor"<<endl;
00059 }
00060
00061
00062
00063 MeuCuts::~MeuCuts()
00064 {
00065 MSG("MeuCuts",Msg::kDebug)
00066 <<"Running MeuCuts Destructor..."<<endl;
00067
00068
00069 MSG("MeuCuts",Msg::kDebug)
00070 <<"Finished MeuCuts Destructor"<<endl;
00071 }
00072
00073
00074
00075 Bool_t MeuCuts::IsInPittFidVol(Float_t x,Float_t y,Float_t z,
00076 Float_t u,Float_t v) const
00077 {
00078
00079
00080
00081 if( !( z>0.6 && z<3.56) ) return kFALSE;
00082 if( !( u>0.3 && u<1.8) ) return kFALSE;
00083 if( !( v>-1.8 && v< -0.3) ) return kFALSE;
00084 if( x>=2.4) return kFALSE;
00085 const Float_t coil_cut=0.8*0.8;
00086 if( (x*x + y*y) <= coil_cut) return kFALSE;
00087 return kTRUE;
00088 }
00089
00090
00091
00092 Bool_t MeuCuts::CheckTrkViews(const MeuSummary& s,Int_t* diffAtVtx,
00093 Int_t* diffAtEnd) const
00094 {
00095 Bool_t noLargeDifference=true;
00096
00097 if (s.MinPlane3!=999 && s.MaxPlane3!=-1 &&
00098 s.MinPlane2!=999 && s.MaxPlane2!=-1){
00099
00100 Int_t diffAtLow=abs(s.MinPlane2-s.MinPlane3);
00101 Int_t diffAtHigh=abs(s.MaxPlane2-s.MaxPlane3);
00102 MAXMSG("MeuCuts",Msg::kDebug,100)
00103 <<"diffAtLow="<<diffAtLow
00104 <<", fMinPlane2="<<s.MinPlane2
00105 <<", fMinPlane3="<<s.MinPlane3<<endl
00106 <<"diffAtHigh="<<diffAtHigh
00107 <<", fMaxPlane2="<<s.MaxPlane2
00108 <<", fMaxPlane3="<<s.MaxPlane3<<endl;
00109
00110 Bool_t posDir=s.VtxPlane<s.EndPlane;
00111 Int_t diffE=diffAtHigh;
00112 Int_t diffV=diffAtLow;
00113 if (!posDir) {
00114 diffE=diffAtLow;
00115 diffV=diffAtHigh;
00116 }
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130 if (s.Detector==Detector::kFar){
00131 if (diffAtLow>=7 || diffAtHigh>=7){
00132 MAXMSG("MeuCuts",Msg::kDebug,100)
00133 <<"Found large gap:"<<endl
00134 <<"diffAtLow="<<diffAtLow
00135 <<", fMinPlane2="<<s.MinPlane2
00136 <<", fMinPlane3="<<s.MinPlane3<<endl
00137 <<"diffAtHigh="<<diffAtHigh
00138 <<", fMaxPlane2="<<s.MaxPlane2
00139 <<", fMaxPlane3="<<s.MaxPlane3
00140 <<endl;
00141 noLargeDifference=false;
00142 }
00143 }
00144 else if (s.Detector==Detector::kNear){
00145 if (diffE>=7 || diffV>=11){
00146 MAXMSG("MeuCuts",Msg::kDebug,100)
00147 <<"Found large gap:"<<endl
00148 <<"diffLow="<<diffAtLow
00149 <<", MinPl2="<<s.MinPlane2
00150 <<", MinPl3="<<s.MinPlane3<<endl
00151 <<"diffHigh="<<diffAtHigh
00152 <<", MaxPl2="<<s.MaxPlane2
00153 <<", MaxPl3="<<s.MaxPlane3
00154 <<", dE="<<diffE
00155 <<", dV="<<diffV
00156 <<endl;
00157 noLargeDifference=false;
00158 }
00159 }
00160 else cout<<"No detector type!"<<endl;
00161
00162
00163 if (diffAtVtx!=0 && diffAtEnd!=0) {
00164 *diffAtVtx=diffV;
00165 *diffAtEnd=diffE;
00166 }
00167 }
00168 else{
00169 MAXMSG("MeuCuts",Msg::kInfo,50)
00170 <<"Can't analyse track views because variables not set"<<endl;
00171 noLargeDifference=false;
00172 }
00173
00174 return noLargeDifference;
00175 }
00176
00177
00178
00179 Bool_t MeuCuts::FilterBadXY(const AtmosEvent& ev,
00180 const MeuSummary& s) const
00181 {
00182 const AtmosTrack* track=dynamic_cast<const AtmosTrack*>
00183 (ev.TrackList->At(0));
00184
00185 if (!track) {
00186
00187 return true;
00188 }
00189
00190 TClonesArray& strips=(*ev.StripList);
00191 Int_t numStrips=strips.GetEntries();
00192 if (numStrips<=0) {
00193 MSG("MeuCuts",Msg::kWarning)
00194 <<"numStrips="<<numStrips<<endl;
00195
00196 return true;
00197 }
00198
00199 if (s.VtxPlane<0 || s.EndPlane<0){
00200 MSG("MeuCuts",Msg::kWarning)
00201 <<"No vtx and/or end plane for this event!"<<endl;
00202
00203 return true;
00204 }
00205
00207
00209 for (Int_t hit=0;hit<numStrips;hit++){
00210 const AtmosStrip* strip=dynamic_cast<AtmosStrip*>(strips[hit]);
00211 const AtmosStrip& st=(*strip);
00212
00213
00214 if (!st.Trk) continue;
00215
00216
00217 Float_t x=-99999;
00218 Float_t y=-99999;
00219 Float_t z=-99999;
00220 this->CalcXYZ(ev,st,x,y,z);
00221 if (fabs(x)>4.1 || fabs(y)>4.1) {
00222 MSG("MeuCuts",Msg::kDebug)
00223 <<"silly x or y; x="<<x<<", y="<<y<<endl;
00224 return true;
00225 }
00226 }
00227
00228 return false;
00229 }
00230
00231
00232
00233 Bool_t MeuCuts::FilterBadEvtPerSlc(const NtpStRecord& ntp,Int_t slc,
00234 Int_t entry,Int_t* evts) const
00235 {
00236 static Int_t lastEntry=-1;
00237 static map<Int_t,Int_t> evtPerSlc;
00238
00239
00240 if (entry!=lastEntry){
00241 MAXMSG("MeuCuts",Msg::kDebug,30)
00242 <<"Building map..."<<endl;
00243
00244 evtPerSlc.clear();
00245
00246
00247 lastEntry=entry;
00248
00249 TClonesArray& evtTca=(*ntp.evt);
00250 const Int_t numEvts=evtTca.GetEntriesFast();
00251
00252
00253 for (Int_t ntpevt=0;ntpevt<numEvts;++ntpevt) {
00254 const NtpSREvent* pevt=
00255 dynamic_cast<NtpSREvent*>(evtTca[ntpevt]);
00256 const NtpSREvent& evt=(*pevt);
00257
00258
00259 evtPerSlc[evt.slc]++;
00260 }
00261
00262 Bool_t print=false;
00263 if (print){
00264
00265
00266 typedef map<Int_t,Int_t>::const_iterator slcIt;
00267 MAXMSG("MeuCuts",Msg::kDebug,30)
00268 <<"Events per slice for entry="<<entry<<endl;
00269 for (slcIt it=evtPerSlc.begin();it!=evtPerSlc.end();++it) {
00270 MAXMSG("MeuCuts",Msg::kDebug,200)
00271 <<" slice "<<it->first<<" has "
00272 <<it->second<<" event(s)"<<endl;
00273 }
00274 }
00275 }
00276
00277 if (evts!=0) {
00278 if (evtPerSlc.find(slc)!=evtPerSlc.end()) {
00279 *evts=evtPerSlc[slc];
00280 }
00281 else *evts=-1;
00282 }
00283
00284 return evtPerSlc[slc]!=1;
00285 }
00286
00287
00288
00289 Bool_t MeuCuts::FilterBadTrackEnd
00290 (const MeuSummary& s,const std::map<Int_t,MeuHitInfo>& plInfo,
00291 Float_t* sigBeyondEnd) const
00292 {
00293 Bool_t posDir=(s.EndPlane>s.VtxPlane);
00294 Bool_t badTracking=false;
00295
00296 Float_t sumSigCorBeyondEnd=0;
00297 Float_t sumAdcBeyondEnd=0;
00298
00299 typedef map<Int_t,MeuHitInfo>::const_iterator iter;
00300 for (iter it=plInfo.begin();it!=plInfo.end();++it){
00301 const MeuHitInfo& cp=it->second;
00302 Int_t pl=cp.Plane;
00303
00304
00305
00306 Bool_t beyondTrkEnd=pl<s.EndPlane;
00307 if (posDir) beyondTrkEnd=pl>s.EndPlane;
00308
00309 MAXMSG("MeuCuts",Msg::kVerbose,1000)
00310 <<"pl="<<pl<<", track:"<<s.VtxPlane<<"->"<<s.EndPlane<<endl;
00311
00312 if (!beyondTrkEnd) continue;
00313
00314 sumSigCorBeyondEnd+=cp.SigCor;
00315 sumAdcBeyondEnd+=cp.Adc;
00316
00317 MAXMSG("MeuCuts",Msg::kVerbose,1000)
00318 <<"Charge beyond trk end: pl="<<pl
00319 <<", track: "<<s.VtxPlane<<"->"<<s.EndPlane
00320 <<", sum="<<sumSigCorBeyondEnd<<endl;
00321
00322
00323
00324
00325
00326
00327
00328 if (sumSigCorBeyondEnd>500 || sumAdcBeyondEnd>500) badTracking=true;
00329 }
00330
00331
00332 if (sigBeyondEnd!=0) *sigBeyondEnd=sumAdcBeyondEnd;
00333
00334 return badTracking;
00335 }
00336
00337
00338
00339 Bool_t MeuCuts::FilterBadEndGaps
00340 (const MeuSummary& s,const std::map<Int_t,MeuHitInfo>& plInfo,
00341 Float_t* adcInGap) const
00342 {
00343 Bool_t posDir=(s.EndPlane>s.VtxPlane);
00344 Bool_t badGap=false;
00345 Float_t sumAdcInGap=0;
00346
00347
00348
00349
00350
00351
00352
00353
00354 Int_t startPl=s.EndPlane-6;
00355 if (!posDir) startPl=s.EndPlane+6;
00356
00357 MAXMSG("MeuCuts",Msg::kDebug,10)
00358 <<"FilterBadEndGaps: startPlane="<<startPl
00359 <<", track:"<<s.VtxPlane<<"->"<<s.EndPlane<<endl;
00360
00361 Int_t pl=startPl;
00362 while (pl!=s.EndPlane){
00363
00364 MAXMSG("MeuCuts",Msg::kDebug,200)
00365 <<"Looping: pl="<<pl<<endl;
00366
00367 map<Int_t,MeuHitInfo>::const_iterator it=plInfo.find(pl);
00368 if (it==plInfo.end()) {
00369 MSG("MeuCuts",Msg::kFatal)
00370 <<"Ahhh = .end(), pl="<<pl<<", startPlane="<<startPl
00371 <<", track:"<<s.VtxPlane<<"->"<<s.EndPlane<<endl
00372 <<endl<<"Will exit here..."<<endl;
00373 exit(1);
00374 }
00375
00376 const MeuHitInfo& cp=it->second;
00377
00378
00379 if (cp.Strip>=0) {
00380 if (posDir) pl++;
00381 else pl--;
00382 continue;
00383 }
00384
00385 MAXMSG("MeuCuts",Msg::kDebug,100)
00386 <<"Found untracked plane="<<pl
00387 <<", adc="<<cp.Adc<<", sumAdc="<<sumAdcInGap<<endl;
00388
00389 sumAdcInGap+=cp.Adc;
00390
00391 if (posDir) pl++;
00392 else pl--;
00393 }
00394
00395 if (sumAdcInGap>500) {
00396 MAXMSG("MeuCuts",Msg::kDebug,100)
00397 <<"Found bad track end gap, run="<<s.Run<<", snarl="<<s.Snarl
00398 <<", sumAdc="<<sumAdcInGap<<endl;
00399 badGap=true;
00400 }
00401
00402
00403 if (adcInGap) *adcInGap=sumAdcInGap;
00404
00405 return badGap;
00406 }
00407
00408
00409
00410 Bool_t MeuCuts::FilterBadDistEndStrip
00411 (const MeuSummary& s,const std::map<Int_t,MeuHitInfo>& plInfo,
00412 Float_t* minDistEndStrip) const
00413 {
00414 Bool_t filter=false;
00415
00416
00417
00418
00419 Float_t minDist=999;
00420 Int_t lowPl=s.WinVtxSidePl;
00421 Int_t highPl=s.WinStopSidePl;
00422 if (lowPl>highPl){
00423 lowPl=s.WinStopSidePl;
00424 highPl=s.WinVtxSidePl;
00425 }
00426
00427
00428 if (plInfo.find(lowPl)==plInfo.end() ||
00429 plInfo.find(highPl)==plInfo.end()) {
00430 MAXMSG("MeuCuts",Msg::kWarning,200)
00431 <<"FilterBadDistEndStrip: no window, can't do it!"<<endl;
00432 filter=true;
00433 return filter;
00434 }
00435
00436
00437 for (map<Int_t,MeuHitInfo>::const_iterator it=plInfo.find(lowPl);
00438 it!=++plInfo.find(highPl);++it){
00439
00440 const MeuHitInfo& cp=it->second;
00441
00442 Float_t distToStripEnd=cp.StripHalfLength-
00443 fabs(cp.DistToStripCentre);
00444 if (distToStripEnd<minDist) minDist=distToStripEnd;
00445
00446 if (distToStripEnd<0.2){
00447 filter=true;
00448 }
00449
00450 MAXMSG("MeuCuts",Msg::kDebug,500)
00451 <<"pl="<<cp.Plane<<", distToStripEnd="<<distToStripEnd
00452 <<", HL="<<cp.StripHalfLength
00453 <<", distToCent="<<cp.DistToStripCentre
00454 <<", filt="<<filter<<endl;
00455 }
00456
00457
00458 if (minDistEndStrip) *minDistEndStrip=minDist;
00459
00460 return filter;
00461 }
00462
00463
00464
00465 Bool_t MeuCuts::FilterLowMatTraversed(const MeuSummary& s) const
00466 {
00467 const Float_t material01Pl=(5.94)*Munits::cm;
00468 const Float_t material16Pl=16*material01Pl;
00469 const Float_t material14Pl=14*material01Pl;
00470
00471 Bool_t lowMatTrav=true;
00472
00473
00474 if (s.TotalMatTraversed>material16Pl+material14Pl){
00475 lowMatTrav=false;
00476 }
00477
00478 return lowMatTrav;
00479 }
00480
00481
00482
00483 void MeuCuts::SetSpecificCuts(MeuSummary& s) const
00484 {
00485 if (s.Detector==Detector::kNear){
00486 s.RFid=1.0;
00487 s.DistToEdgeFid=0.4;
00488 }
00489 else if (s.Detector==Detector::kFar) {
00490 s.RFid=3.0*Munits::m;
00491 s.RFidCoil=0.5*Munits::m;
00492 s.DistToEdgeFid=1.0*Munits::m;
00493 }
00494 else {
00495 cout<<"Ahhh detector not defined"<<endl;
00496 }
00497 }
00498
00499
00500
00501 void MeuCuts::FillSTSumDetails(const NtpStRecord& ntp,
00502 MeuSummary& s,Int_t evt,Int_t slc) const
00503 {
00504 const RecCandHeader& rec=ntp.GetHeader();
00505 TClonesArray& slcTca=(*ntp.slc);
00506
00507 s.Run=rec.GetRun();
00508 s.SubRun=rec.GetSubRun();
00509 s.Snarl=rec.GetSnarl();
00510 s.TimeSec=ntp.GetVldContext()->GetTimeStamp().GetSec();
00511 s.TimeNanoSec=rec.GetVldContext().GetTimeStamp().GetNanoSec();
00512
00513 s.NStrip=dynamic_cast<NtpSRSlice*>(slcTca[slc])->nstrip;
00514 s.Detector=rec.GetVldContext().GetDetector();
00515 s.SimFlag=rec.GetVldContext().GetSimFlag();
00516 s.TrigSrc=rec.GetTrigSrc();
00517 s.RunType=rec.GetRunType();
00518
00519 s.Evt=evt;
00520 s.Slc=slc;
00521
00522
00523 if (s.SimFlag==SimFlag::kMC && s.Detector==Detector::kNear) {
00524
00525 TClonesArray& mcTca=*ntp.mc;
00526 MAXMSG("MeuAnalysis",Msg::kInfo,1)
00527 <<"Number of entries in NtpMCTruth="<<mcTca.GetEntries()
00528 <<", if >0 then spill MC"<<endl;
00529 if (mcTca.GetEntries()==0){
00530 MAXMSG("MeuAnalysis",Msg::kWarning,5)
00531 <<"Performing hack to correct the fTrigSrc of"
00532 <<" ND cosmic MC, was="
00533 <<s.TrigSrc<<", setting to 4"<<endl;
00534 s.TrigSrc=4;
00535 }
00536 }
00537
00538
00539 if (this->IsGoodDataQuality(ntp)) s.GoodDataQuality = true;
00540 else s.GoodDataQuality = false;
00541
00542 NtpSRDataQuality dataqua = ntp.dataquality;
00543 s.CrateMask = dataqua.cratemask;
00544 s.BusyChips = dataqua.busychips;
00545 s.ColdChips = dataqua.coldchips;
00546 }
00547
00548
00549
00550 void MeuCuts::FillSTSumDetails(const AtmosEvent& ev,
00551 MeuSummary& s) const
00552 {
00553 s.Run=ev.Run;
00554 s.SubRun=ev.SubRun;
00555 s.Snarl=ev.Snarl;
00556 s.TimeSec=ev.UnixTime;
00557 s.TimeNanoSec=ev.NanoSec;
00558
00559 s.NStrip=(*ev.StripList).GetEntries();
00560 s.Detector=Detector::kFar;
00561
00562 if (ev.NScintHits<=0) s.SimFlag=SimFlag::kData;
00563 else s.SimFlag=SimFlag::kMC;
00564 s.TrigSrc=ev.TrigSrc;
00565
00566
00567 s.Evt=-1;
00568 s.Slc=-1;
00569 }
00570
00571
00572
00573 void MeuCuts::FillBeamMonDetails(const NtpBDLiteRecord& bd,
00574 MeuSummary& s) const
00575 {
00576 s.BDtortgt=bd.tortgt;
00577 s.BDtor101=bd.tor101;
00578 }
00579
00580
00581
00582 void MeuCuts::FillSTSumRecoDetails
00583 (MeuSummary& s,const std::map<Int_t,MeuHitInfo>& plInfo) const
00584 {
00585 map<Int_t,MeuHitInfo>::const_iterator
00586 it=plInfo.find(s.VtxPlane);
00587 if (it!=plInfo.end()){
00588 const MeuHitInfo& cp=it->second;
00589 s.VtxX=cp.X;
00590 s.VtxY=cp.Y;
00591 s.VtxZ=cp.Z;
00592 }
00593 it=plInfo.find(s.EndPlane);
00594 if (it!=plInfo.end()){
00595 const MeuHitInfo& cp=it->second;
00596 s.EndX=cp.X;
00597 s.EndY=cp.Y;
00598 s.EndZ=cp.Z;
00599 }
00600
00601 it=plInfo.find(s.WinVtxSidePl);
00602 if (it!=plInfo.end()){
00603 const MeuHitInfo& cp=it->second;
00604 s.WinVtxSideView=cp.View;
00605 s.WinVtxSideX=cp.X;
00606 s.WinVtxSideY=cp.Y;
00607 s.WinVtxSideZ=cp.Z;
00608 s.WinVtxSideTPos=cp.TPos;
00609 s.WinVtxSideLPos=cp.LPos;
00610 s.WinVtxSideStrip=cp.Strip;
00611 }
00612
00613 it=plInfo.find(s.WinStopSidePl);
00614 if (it!=plInfo.end()){
00615 const MeuHitInfo& cp=it->second;
00616 s.WinStopSideView=cp.View;
00617 s.WinStopSideX=cp.X;
00618 s.WinStopSideY=cp.Y;
00619 s.WinStopSideZ=cp.Z;
00620 s.WinStopSideTPos=cp.TPos;
00621 s.WinStopSideLPos=cp.LPos;
00622 s.WinStopSideStrip=cp.Strip;
00623 }
00624 }
00625
00626
00627
00628 void MeuCuts::FillTimeHistos(const NtpStRecord& ntp,
00629 const NtpSREvent& evt,
00630 const MeuSummary& s) const
00631 {
00632
00633 static TH1F* hTimeLength=0;
00634 static TH1F* hTimeBefore=0;
00635 static TH1F* hTimeAfter=0;
00636
00637 static TH1F* hSlcTimeLength=0;
00638 static TH1F* hSlcTimeBefore=0;
00639 static TH1F* hSlcTimeAfter=0;
00640
00641 static TH1F* hTrkTimeLength=0;
00642 static TH1F* hTrkTimeBefore=0;
00643 static TH1F* hTrkTimeAfter=0;
00644
00645 static TH1F* hMedianTrkSlcDiff=0;
00646
00647 static TProfile* pAdcVsTimeCut=0;
00648 static TProfile* pSigCorVsTimeCut=0;
00649
00650
00651
00652 static TH1F* hSigMapWinStp=0;
00653 static TH1F* hSigMapAllStp=0;
00654 static TH1F* hSigCorWinStp=0;
00655 static TH1F* hSigCorAllStp=0;
00656 static TH1F* hSigLinWinStp=0;
00657 static TH1F* hSigLinAllStp=0;
00658 static TH1F* hAdcWinStp=0;
00659 static TH1F* hAdcAllStp=0;
00660 static TH1F* hPeWinStp=0;
00661 static TH1F* hPeAllStp=0;
00662 static TH1F* hDigitWinStp=0;
00663 static TH1F* hDigitAllStp=0;
00664
00665 static TH1F* hSigMapWinStp1=0;
00666 static TH1F* hSigMapAllStp1=0;
00667 static TH1F* hSigCorWinStp1=0;
00668 static TH1F* hSigCorAllStp1=0;
00669 static TH1F* hSigLinWinStp1=0;
00670 static TH1F* hSigLinAllStp1=0;
00671 static TH1F* hAdcWinStp1=0;
00672 static TH1F* hAdcAllStp1=0;
00673 static TH1F* hPeWinStp1=0;
00674 static TH1F* hPeAllStp1=0;
00675 static TH1F* hDigitWinStp1=0;
00676 static TH1F* hDigitAllStp1=0;
00677
00678 static TH1F* hSigMapWinStp2=0;
00679 static TH1F* hSigMapAllStp2=0;
00680 static TH1F* hSigCorWinStp2=0;
00681 static TH1F* hSigCorAllStp2=0;
00682 static TH1F* hSigLinWinStp2=0;
00683 static TH1F* hSigLinAllStp2=0;
00684 static TH1F* hAdcWinStp2=0;
00685 static TH1F* hAdcAllStp2=0;
00686 static TH1F* hPeWinStp2=0;
00687 static TH1F* hPeAllStp2=0;
00688 static TH1F* hDigitWinStp2=0;
00689 static TH1F* hDigitAllStp2=0;
00690
00691 if (hTimeLength==0) {
00692 MSG("MeuCuts",Msg::kDebug)
00693 <<"Creating time histos..."<<endl;
00694
00695 hTimeLength=new TH1F("hTimeLength","hTimeLength",10000,-100,50000);
00696 hTimeLength->SetFillColor(0);
00697 hTimeLength->SetTitle("Snarl Length in Time");
00698 hTimeLength->GetXaxis()->SetTitle("Time (ns)");
00699 hTimeLength->GetXaxis()->CenterTitle();
00700
00701 hTimeBefore=new TH1F("hTimeBefore","hTimeBefore",10000,-100,50000);
00702 hTimeBefore->SetFillColor(0);
00703 hTimeBefore->SetTitle("Time Between Event and First Hit");
00704 hTimeBefore->GetXaxis()->SetTitle("Time (ns)");
00705 hTimeBefore->GetXaxis()->CenterTitle();
00706
00707 hTimeAfter=new TH1F("hTimeAfter","hTimeAfter",10000,-100,50000);
00708 hTimeAfter->SetFillColor(0);
00709 hTimeAfter->SetTitle("Time Between Event and Last Hit");
00710 hTimeAfter->GetXaxis()->SetTitle("Time (ns)");
00711 hTimeAfter->GetXaxis()->CenterTitle();
00712
00713 hSlcTimeLength=new TH1F
00714 ("hSlcTimeLength","hSlcTimeLength",10000,-100,50000);
00715 hSlcTimeLength->SetFillColor(0);
00716 hSlcTimeLength->SetTitle("Slice Length in Time");
00717 hSlcTimeLength->GetXaxis()->SetTitle("Time (ns)");
00718 hSlcTimeLength->GetXaxis()->CenterTitle();
00719
00720 hSlcTimeBefore=new TH1F
00721 ("hSlcTimeBefore","hSlcTimeBefore",10000,-100,50000);
00722 hSlcTimeBefore->SetFillColor(0);
00723 hSlcTimeBefore->SetTitle("Time Between Event and First Hit in Slice");
00724 hSlcTimeBefore->GetXaxis()->SetTitle("Time (ns)");
00725 hSlcTimeBefore->GetXaxis()->CenterTitle();
00726
00727 hSlcTimeAfter=new TH1F
00728 ("hSlcTimeAfter","hSlcTimeAfter",10000,-100,50000);
00729 hSlcTimeAfter->SetFillColor(0);
00730 hSlcTimeAfter->SetTitle("Time Between Event and Last Hit in Slice");
00731 hSlcTimeAfter->GetXaxis()->SetTitle("Time (ns)");
00732 hSlcTimeAfter->GetXaxis()->CenterTitle();
00733
00734 hTrkTimeLength=new TH1F
00735 ("hTrkTimeLength","hTrkTimeLength",10000,-100,50000);
00736 hTrkTimeLength->SetFillColor(0);
00737 hTrkTimeLength->SetTitle("Track Length in Time");
00738 hTrkTimeLength->GetXaxis()->SetTitle("Time (ns)");
00739 hTrkTimeLength->GetXaxis()->CenterTitle();
00740
00741 hTrkTimeBefore=new TH1F
00742 ("hTrkTimeBefore","hTrkTimeBefore",10000,-100,50000);
00743 hTrkTimeBefore->SetFillColor(0);
00744 hTrkTimeBefore->SetTitle("Time Between Event and First Tracked Hit");
00745 hTrkTimeBefore->GetXaxis()->SetTitle("Time (ns)");
00746 hTrkTimeBefore->GetXaxis()->CenterTitle();
00747
00748 hTrkTimeAfter=new TH1F
00749 ("hTrkTimeAfter","hTrkTimeAfter",10000,-100,50000);
00750 hTrkTimeAfter->SetFillColor(0);
00751 hTrkTimeAfter->SetTitle("Time Between Event and Last Tracked Hit");
00752 hTrkTimeAfter->GetXaxis()->SetTitle("Time (ns)");
00753 hTrkTimeAfter->GetXaxis()->CenterTitle();
00754
00755 hMedianTrkSlcDiff=new TH1F
00756 ("hMedianTrkSlcDiff","hMedianTrkSlcDiff",10000,-3000,3000);
00757 hMedianTrkSlcDiff->SetFillColor(0);
00758 hMedianTrkSlcDiff->SetTitle("Diff. between Median Slc and Median Trk Time");
00759 hMedianTrkSlcDiff->GetXaxis()->SetTitle("Time (ns)");
00760 hMedianTrkSlcDiff->GetXaxis()->CenterTitle();
00761
00762
00763
00764 pAdcVsTimeCut=new TProfile("pAdcVsTimeCut",
00765 "pAdcVsTimeCut",12,-100,11*100);
00766 pAdcVsTimeCut->GetXaxis()->SetTitle("Time-Evt Time");
00767 pAdcVsTimeCut->GetXaxis()->CenterTitle();
00768 pAdcVsTimeCut->GetYaxis()->SetTitle("Average Adc");
00769 pAdcVsTimeCut->GetYaxis()->CenterTitle();
00770
00771 pSigCorVsTimeCut=new TProfile("pSigCorVsTimeCut",
00772 "pSigCorVsTimeCut",12,-100,11*100);
00773 pSigCorVsTimeCut->GetXaxis()->SetTitle("Time-Evt Time");
00774 pSigCorVsTimeCut->GetXaxis()->CenterTitle();
00775 pSigCorVsTimeCut->GetYaxis()->SetTitle("Average SigCor");
00776 pSigCorVsTimeCut->GetYaxis()->CenterTitle();
00777
00778
00779
00780
00781 hSigMapWinStp=new TH1F("hSigMapWinStp","hSigMapWinStp",
00782 20320,-320,20000);
00783 hSigMapWinStp->SetFillColor(0);
00784 hSigMapWinStp->SetTitle("SigMap of Strips in Window (in Slice)");
00785 hSigMapWinStp->GetXaxis()->SetTitle("SigMap");
00786 hSigMapWinStp->GetXaxis()->CenterTitle();
00787
00788
00789 hSigMapAllStp=new TH1F("hSigMapAllStp","hSigMapAllStp",
00790 20320,-320,20000);
00791 hSigMapAllStp->SetFillColor(0);
00792 hSigMapAllStp->SetTitle("SigMap of All Strips in Slice");
00793 hSigMapAllStp->GetXaxis()->SetTitle("SigMap");
00794 hSigMapAllStp->GetXaxis()->CenterTitle();
00795
00796
00797
00798 hSigCorWinStp=new TH1F("hSigCorWinStp","hSigCorWinStp",
00799 20320,-320,20000);
00800 hSigCorWinStp->SetFillColor(0);
00801 hSigCorWinStp->SetTitle("SigCor of Strips in Window (in Slice)");
00802 hSigCorWinStp->GetXaxis()->SetTitle("SigCor");
00803 hSigCorWinStp->GetXaxis()->CenterTitle();
00804
00805
00806 hSigCorAllStp=new TH1F("hSigCorAllStp","hSigCorAllStp",
00807 20320,-320,20000);
00808 hSigCorAllStp->SetFillColor(0);
00809 hSigCorAllStp->SetTitle("SigCor of All Strips in Slice");
00810 hSigCorAllStp->GetXaxis()->SetTitle("SigCor");
00811 hSigCorAllStp->GetXaxis()->CenterTitle();
00812
00813
00814
00815 hSigLinWinStp=new TH1F("hSigLinWinStp","hSigLinWinStp",
00816 20320,-320,20000);
00817 hSigLinWinStp->SetFillColor(0);
00818 hSigLinWinStp->SetTitle("SigLin of Strips in Window (in Slice)");
00819 hSigLinWinStp->GetXaxis()->SetTitle("SigLin");
00820 hSigLinWinStp->GetXaxis()->CenterTitle();
00821
00822
00823 hSigLinAllStp=new TH1F("hSigLinAllStp","hSigLinAllStp",
00824 20320,-320,20000);
00825 hSigLinAllStp->SetFillColor(0);
00826 hSigLinAllStp->SetTitle("SigLin of All Strips in Slice");
00827 hSigLinAllStp->GetXaxis()->SetTitle("SigLin");
00828 hSigLinAllStp->GetXaxis()->CenterTitle();
00829
00830
00831
00832 hAdcWinStp=new TH1F("hAdcWinStp","hAdcWinStp",
00833 20320,-320,20000);
00834 hAdcWinStp->SetFillColor(0);
00835 hAdcWinStp->SetTitle("Adc of Strips in Window (in Slice)");
00836 hAdcWinStp->GetXaxis()->SetTitle("Adc");
00837 hAdcWinStp->GetXaxis()->CenterTitle();
00838
00839
00840 hAdcAllStp=new TH1F("hAdcAllStp","hAdcAllStp",
00841 20320,-320,20000);
00842 hAdcAllStp->SetFillColor(0);
00843 hAdcAllStp->SetTitle("Adc of All Strips in Slice");
00844 hAdcAllStp->GetXaxis()->SetTitle("Adc");
00845 hAdcAllStp->GetXaxis()->CenterTitle();
00846
00847
00848
00849 hPeWinStp=new TH1F("hPeWinStp","hPeWinStp",
00850 10160,-4,250);
00851 hPeWinStp->SetFillColor(0);
00852 hPeWinStp->SetTitle("Pe of Strips in Window (in Slice)");
00853 hPeWinStp->GetXaxis()->SetTitle("Pe");
00854 hPeWinStp->GetXaxis()->CenterTitle();
00855
00856
00857 hPeAllStp=new TH1F("hPeAllStp","hPeAllStp",
00858 10160,-4,250);
00859 hPeAllStp->SetFillColor(0);
00860 hPeAllStp->SetTitle("Pe of All Strips in Slice");
00861 hPeAllStp->GetXaxis()->SetTitle("Pe");
00862 hPeAllStp->GetXaxis()->CenterTitle();
00863
00864
00865
00866 hDigitWinStp=new TH1F("hDigitWinStp","hDigitWinStp",
00867 4020,-2,400);
00868 hDigitWinStp->SetFillColor(0);
00869 hDigitWinStp->SetTitle("Digits of Strips in Window (in Slice)");
00870 hDigitWinStp->GetXaxis()->SetTitle("Digits");
00871 hDigitWinStp->GetXaxis()->CenterTitle();
00872 hDigitWinStp->SetBit(TH1::kCanRebin);
00873
00874 hDigitAllStp=new TH1F("hDigitAllStp","hDigitAllStp",
00875 4020,-2,400);
00876 hDigitAllStp->SetFillColor(0);
00877 hDigitAllStp->SetTitle("Digits of All Strips in Slice");
00878 hDigitAllStp->GetXaxis()->SetTitle("Digits");
00879 hDigitAllStp->GetXaxis()->CenterTitle();
00880
00881
00882
00883
00884 hSigMapWinStp1=new TH1F("hSigMapWinStp1","hSigMapWinStp1",
00885 20320,-320,20000);
00886 hSigMapWinStp1->SetFillColor(0);
00887 hSigMapWinStp1->SetTitle("SigMap of Strips in Window (in Slice)");
00888 hSigMapWinStp1->GetXaxis()->SetTitle("SigMap");
00889 hSigMapWinStp1->GetXaxis()->CenterTitle();
00890
00891
00892 hSigMapAllStp1=new TH1F("hSigMapAllStp1","hSigMapAllStp1",
00893 20320,-320,20000);
00894 hSigMapAllStp1->SetFillColor(0);
00895 hSigMapAllStp1->SetTitle("SigMap of All Strips in Slice");
00896 hSigMapAllStp1->GetXaxis()->SetTitle("SigMap");
00897 hSigMapAllStp1->GetXaxis()->CenterTitle();
00898
00899
00900
00901 hSigCorWinStp1=new TH1F("hSigCorWinStp1","hSigCorWinStp1",
00902 20320,-320,20000);
00903 hSigCorWinStp1->SetFillColor(0);
00904 hSigCorWinStp1->SetTitle("SigCor of Strips in Window (in Slice)");
00905 hSigCorWinStp1->GetXaxis()->SetTitle("SigCor");
00906 hSigCorWinStp1->GetXaxis()->CenterTitle();
00907
00908
00909 hSigCorAllStp1=new TH1F("hSigCorAllStp1","hSigCorAllStp1",
00910 20320,-320,20000);
00911 hSigCorAllStp1->SetFillColor(0);
00912 hSigCorAllStp1->SetTitle("SigCor of All Strips in Slice");
00913 hSigCorAllStp1->GetXaxis()->SetTitle("SigCor");
00914 hSigCorAllStp1->GetXaxis()->CenterTitle();
00915
00916
00917
00918 hSigLinWinStp1=new TH1F("hSigLinWinStp1","hSigLinWinStp1",
00919 20320,-320,20000);
00920 hSigLinWinStp1->SetFillColor(0);
00921 hSigLinWinStp1->SetTitle("SigLin of Strips in Window (in Slice)");
00922 hSigLinWinStp1->GetXaxis()->SetTitle("SigLin");
00923 hSigLinWinStp1->GetXaxis()->CenterTitle();
00924
00925
00926 hSigLinAllStp1=new TH1F("hSigLinAllStp1","hSigLinAllStp1",
00927 20320,-320,20000);
00928 hSigLinAllStp1->SetFillColor(0);
00929 hSigLinAllStp1->SetTitle("SigLin of All Strips in Slice");
00930 hSigLinAllStp1->GetXaxis()->SetTitle("SigLin");
00931 hSigLinAllStp1->GetXaxis()->CenterTitle();
00932
00933
00934
00935 hAdcWinStp1=new TH1F("hAdcWinStp1","hAdcWinStp1",
00936 20320,-320,20000);
00937 hAdcWinStp1->SetFillColor(0);
00938 hAdcWinStp1->SetTitle("Adc of Strips in Window (in Slice)");
00939 hAdcWinStp1->GetXaxis()->SetTitle("Adc");
00940 hAdcWinStp1->GetXaxis()->CenterTitle();
00941
00942
00943 hAdcAllStp1=new TH1F("hAdcAllStp1","hAdcAllStp1",
00944 20320,-320,20000);
00945 hAdcAllStp1->SetFillColor(0);
00946 hAdcAllStp1->SetTitle("Adc of All Strips in Slice");
00947 hAdcAllStp1->GetXaxis()->SetTitle("Adc");
00948 hAdcAllStp1->GetXaxis()->CenterTitle();
00949
00950
00951
00952 hPeWinStp1=new TH1F("hPeWinStp1","hPeWinStp1",
00953 10160,-4,250);
00954 hPeWinStp1->SetFillColor(0);
00955 hPeWinStp1->SetTitle("Pe of Strips in Window (in Slice)");
00956 hPeWinStp1->GetXaxis()->SetTitle("Pe");
00957 hPeWinStp1->GetXaxis()->CenterTitle();
00958
00959
00960 hPeAllStp1=new TH1F("hPeAllStp1","hPeAllStp1",
00961 10160,-4,250);
00962 hPeAllStp1->SetFillColor(0);
00963 hPeAllStp1->SetTitle("Pe of All Strips in Slice");
00964 hPeAllStp1->GetXaxis()->SetTitle("Pe");
00965 hPeAllStp1->GetXaxis()->CenterTitle();
00966
00967
00968
00969 hDigitWinStp1=new TH1F("hDigitWinStp1","hDigitWinStp1",
00970 4020,-2,400);
00971 hDigitWinStp1->SetFillColor(0);
00972 hDigitWinStp1->SetTitle("Digits of Strips in Window (in Slice)");
00973 hDigitWinStp1->GetXaxis()->SetTitle("Digits");
00974 hDigitWinStp1->GetXaxis()->CenterTitle();
00975
00976
00977 hDigitAllStp1=new TH1F("hDigitAllStp1","hDigitAllStp1",
00978 4020,-2,400);
00979 hDigitAllStp1->SetFillColor(0);
00980 hDigitAllStp1->SetTitle("Digits of All Strips in Slice");
00981 hDigitAllStp1->GetXaxis()->SetTitle("Digits");
00982 hDigitAllStp1->GetXaxis()->CenterTitle();
00983
00984
00985
00986
00987 hSigMapWinStp2=new TH1F("hSigMapWinStp2","hSigMapWinStp2",
00988 20320,-320,20000);
00989 hSigMapWinStp2->SetFillColor(0);
00990 hSigMapWinStp2->SetTitle("SigMap of Strips in Window (in Slice)");
00991 hSigMapWinStp2->GetXaxis()->SetTitle("SigMap");
00992 hSigMapWinStp2->GetXaxis()->CenterTitle();
00993
00994
00995 hSigMapAllStp2=new TH1F("hSigMapAllStp2","hSigMapAllStp2",
00996 20320,-320,20000);
00997 hSigMapAllStp2->SetFillColor(0);
00998 hSigMapAllStp2->SetTitle("SigMap of All Strips in Slice");
00999 hSigMapAllStp2->GetXaxis()->SetTitle("SigMap");
01000 hSigMapAllStp2->GetXaxis()->CenterTitle();
01001
01002
01003
01004 hSigCorWinStp2=new TH1F("hSigCorWinStp2","hSigCorWinStp2",
01005 20320,-320,20000);
01006 hSigCorWinStp2->SetFillColor(0);
01007 hSigCorWinStp2->SetTitle("SigCor of Strips in Window (in Slice)");
01008 hSigCorWinStp2->GetXaxis()->SetTitle("SigCor");
01009 hSigCorWinStp2->GetXaxis()->CenterTitle();
01010
01011
01012 hSigCorAllStp2=new TH1F("hSigCorAllStp2","hSigCorAllStp2",
01013 20320,-320,20000);
01014 hSigCorAllStp2->SetFillColor(0);
01015 hSigCorAllStp2->SetTitle("SigCor of All Strips in Slice");
01016 hSigCorAllStp2->GetXaxis()->SetTitle("SigCor");
01017 hSigCorAllStp2->GetXaxis()->CenterTitle();
01018
01019
01020
01021 hSigLinWinStp2=new TH1F("hSigLinWinStp2","hSigLinWinStp2",
01022 20320,-320,20000);
01023 hSigLinWinStp2->SetFillColor(0);
01024 hSigLinWinStp2->SetTitle("SigLin of Strips in Window (in Slice)");
01025 hSigLinWinStp2->GetXaxis()->SetTitle("SigLin");
01026 hSigLinWinStp2->GetXaxis()->CenterTitle();
01027
01028
01029 hSigLinAllStp2=new TH1F("hSigLinAllStp2","hSigLinAllStp2",
01030 20320,-320,20000);
01031 hSigLinAllStp2->SetFillColor(0);
01032 hSigLinAllStp2->SetTitle("SigLin of All Strips in Slice");
01033 hSigLinAllStp2->GetXaxis()->SetTitle("SigLin");
01034 hSigLinAllStp2->GetXaxis()->CenterTitle();
01035
01036
01037
01038 hAdcWinStp2=new TH1F("hAdcWinStp2","hAdcWinStp2",
01039 20320,-320,20000);
01040 hAdcWinStp2->SetFillColor(0);
01041 hAdcWinStp2->SetTitle("Adc of Strips in Window (in Slice)");
01042 hAdcWinStp2->GetXaxis()->SetTitle("Adc");
01043 hAdcWinStp2->GetXaxis()->CenterTitle();
01044
01045
01046 hAdcAllStp2=new TH1F("hAdcAllStp2","hAdcAllStp2",
01047 20320,-320,20000);
01048 hAdcAllStp2->SetFillColor(0);
01049 hAdcAllStp2->SetTitle("Adc of All Strips in Slice");
01050 hAdcAllStp2->GetXaxis()->SetTitle("Adc");
01051 hAdcAllStp2->GetXaxis()->CenterTitle();
01052
01053
01054
01055 hPeWinStp2=new TH1F("hPeWinStp2","hPeWinStp2",
01056 10160,-4,250);
01057 hPeWinStp2->SetFillColor(0);
01058 hPeWinStp2->SetTitle("Pe of Strips in Window (in Slice)");
01059 hPeWinStp2->GetXaxis()->SetTitle("Pe");
01060 hPeWinStp2->GetXaxis()->CenterTitle();
01061
01062
01063 hPeAllStp2=new TH1F("hPeAllStp2","hPeAllStp2",
01064 10160,-4,250);
01065 hPeAllStp2->SetFillColor(0);
01066 hPeAllStp2->SetTitle("Pe of All Strips in Slice");
01067 hPeAllStp2->GetXaxis()->SetTitle("Pe");
01068 hPeAllStp2->GetXaxis()->CenterTitle();
01069
01070
01071
01072 hDigitWinStp2=new TH1F("hDigitWinStp2","hDigitWinStp2",
01073 4020,-2,400);
01074 hDigitWinStp2->SetFillColor(0);
01075 hDigitWinStp2->SetTitle("Digits of Strips in Window (in Slice)");
01076 hDigitWinStp2->GetXaxis()->SetTitle("Digits");
01077 hDigitWinStp2->GetXaxis()->CenterTitle();
01078
01079
01080 hDigitAllStp2=new TH1F("hDigitAllStp2","hDigitAllStp2",
01081 4020,-2,400);
01082 hDigitAllStp2->SetFillColor(0);
01083 hDigitAllStp2->SetTitle("Digits of All Strips in Slice");
01084 hDigitAllStp2->GetXaxis()->SetTitle("Digits");
01085 hDigitAllStp2->GetXaxis()->CenterTitle();
01086
01087 }
01088
01089 TClonesArray& stpTca=(*ntp.stp);
01090 Int_t numStps=stpTca.GetEntries();
01091
01092 Double_t earliestTime=999;
01093 Double_t latestTime=-999;
01094 multiset<Double_t> times;
01095
01096
01097 for (Int_t i=0;i<numStps;i++){
01098
01099 const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(stpTca[i]);
01100 const NtpSRStrip& stp=(*pstp);
01101
01102
01103
01104
01105 if (stp.time1>=0 && stp.time1>latestTime) latestTime=stp.time1;
01106 if (stp.time0>=0 && stp.time0>latestTime) latestTime=stp.time0;
01107 if (stp.time1>=0 && stp.time1<earliestTime) earliestTime=stp.time1;
01108 if (stp.time0>=0 && stp.time0<earliestTime) earliestTime=stp.time0;
01109
01110 Double_t time=stp.time0;
01111 if (time<0 || time>1) time=stp.time1;
01112
01113 if (time>0 && time<=1) {
01114 times.insert(time);
01115 }
01116
01117 MAXMSG("MeuCuts",Msg::kDebug,500)
01118 <<"t0="<<stp.time0<<", t1="<<stp.time1
01119 <<", lT="<<latestTime<<", eT="<<earliestTime
01120 <<", dT="<<(latestTime-earliestTime)/Munits::ns<<endl;
01121 }
01122
01123
01124 multiset<Double_t>::const_iterator it=times.begin();
01125 advance(it,times.size()/2);
01126 Double_t medianTime=*it;
01127 MAXMSG("MeuCuts",Msg::kDebug,200)
01128 <<"Median time="<<medianTime
01129 <<", dBefore="<<(medianTime-earliestTime)/Munits::ns
01130 <<", dAfter="<<(latestTime-medianTime)/Munits::ns<<endl;
01131
01132
01133
01134
01135 hTimeAfter->Fill((latestTime-medianTime)/Munits::ns);
01136 hTimeBefore->Fill((medianTime-earliestTime)/Munits::ns);
01137 hTimeLength->Fill((latestTime-earliestTime)/Munits::ns);
01138
01140
01142
01143
01144 earliestTime=999;
01145 latestTime=-999;
01146 times.clear();
01147
01148 TClonesArray& slcTca=(*ntp.slc);
01149
01150
01151 const NtpSRSlice* pslc=
01152 dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]);
01153 const NtpSRSlice& slc=(*pslc);
01154
01155
01156 for (Int_t i=0;i<slc.nstrip;++i) {
01157 const NtpSRStrip* pstp=
01158 dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
01159 const NtpSRStrip& stp=(*pstp);
01160
01161
01162 if (stp.time1>=0 && stp.time1>latestTime) latestTime=stp.time1;
01163 if (stp.time0>=0 && stp.time0>latestTime) latestTime=stp.time0;
01164 if (stp.time1>=0 && stp.time1<earliestTime) earliestTime=stp.time1;
01165 if (stp.time0>=0 && stp.time0<earliestTime) earliestTime=stp.time0;
01166
01167 Double_t time=stp.time0;
01168 if (time<0 || time>1) time=stp.time1;
01169
01170 if (time>0 && time<=1) {
01171 times.insert(time);
01172 }
01173
01174
01175
01176 if ((stp.plane>=s.WinVtxSidePl && stp.plane<=s.WinStopSidePl) ||
01177 (stp.plane<=s.WinVtxSidePl && stp.plane>=s.WinStopSidePl)) {
01178
01179 Double_t timeDiff=time-s.MedianTime;
01180 MAXMSG("MeuCuts",Msg::kDebug,50)
01181 <<"timeDiff="<<timeDiff/(Munits::ns)<<", time="<<time
01182 <<", Median="<<s.MedianTime<<endl;
01183 for (Int_t n=1;n<11;n++) {
01184 if (timeDiff<n*100.*(Munits::ns)) {
01185 MAXMSG("MeuCuts",Msg::kDebug,300)
01186 <<"n="<<n<<", filling..."<<endl;
01187 pAdcVsTimeCut->Fill(n*100,stp.ph1.raw);
01188 pSigCorVsTimeCut->Fill(n*100,stp.ph1.sigcor);
01189 }
01190 else if (timeDiff>10*100.*(Munits::ns) && n==10){
01191 MAXMSG("MeuCuts",Msg::kDebug,300)
01192 <<"n="<<n<<", filling last bin"<<endl;
01193 pAdcVsTimeCut->Fill(10*100,stp.ph1.raw);
01194 pSigCorVsTimeCut->Fill(10*100,stp.ph1.sigcor);
01195 }
01196 else {
01197 MAXMSG("MeuCuts",Msg::kDebug,300)
01198 <<"n="<<n<<", not filling bin"<<endl;
01199 }
01200 }
01201 }
01202
01204
01205
01206 hSigCorAllStp->Fill(stp.ph0.sigcor+stp.ph1.sigcor);
01207 hSigLinAllStp->Fill(stp.ph0.siglin+stp.ph1.siglin);
01208 hAdcAllStp->Fill(stp.ph0.raw+stp.ph1.raw);
01209 hPeAllStp->Fill(stp.ph0.pe+stp.ph1.pe);
01210 hDigitAllStp->Fill(stp.ndigit);
01211
01212
01213 hSigCorAllStp1->Fill(stp.ph0.sigcor);
01214 hSigLinAllStp1->Fill(stp.ph0.siglin);
01215 hAdcAllStp1->Fill(stp.ph0.raw);
01216 hPeAllStp1->Fill(stp.ph0.pe);
01217
01218
01219
01220 hSigCorAllStp2->Fill(stp.ph1.sigcor);
01221 hSigLinAllStp2->Fill(stp.ph1.siglin);
01222 hAdcAllStp2->Fill(stp.ph1.raw);
01223 hPeAllStp2->Fill(stp.ph1.pe);
01224
01225
01226
01227
01228 if ((stp.plane>=s.WinVtxSidePl && stp.plane<=s.WinStopSidePl) ||
01229 (stp.plane<=s.WinVtxSidePl && stp.plane>=s.WinStopSidePl)) {
01230
01231 hSigCorWinStp->Fill(stp.ph0.sigcor+stp.ph1.sigcor);
01232 hSigLinWinStp->Fill(stp.ph0.siglin+stp.ph1.siglin);
01233 hAdcWinStp->Fill(stp.ph0.raw+stp.ph1.raw);
01234 hPeWinStp->Fill(stp.ph0.pe+stp.ph1.pe);
01235 hDigitWinStp->Fill(stp.ndigit);
01236
01237
01238 hSigCorWinStp1->Fill(stp.ph0.sigcor);
01239 hSigLinWinStp1->Fill(stp.ph0.siglin);
01240 hAdcWinStp1->Fill(stp.ph0.raw);
01241 hPeWinStp1->Fill(stp.ph0.pe);
01242
01243
01244
01245 hSigCorWinStp2->Fill(stp.ph1.sigcor);
01246 hSigLinWinStp2->Fill(stp.ph1.siglin);
01247 hAdcWinStp2->Fill(stp.ph1.raw);
01248 hPeWinStp2->Fill(stp.ph1.pe);
01249
01250 }
01251 }
01252
01253
01254 it=times.begin();
01255 advance(it,times.size()/2);
01256 medianTime=*it;
01257
01258 hSlcTimeAfter->Fill((latestTime-medianTime)/Munits::ns);
01259 hSlcTimeBefore->Fill((medianTime-earliestTime)/Munits::ns);
01260 hSlcTimeLength->Fill((latestTime-earliestTime)/Munits::ns);
01261
01263
01265
01266
01267 earliestTime=999;
01268 latestTime=-999;
01269 times.clear();
01270
01271 TClonesArray& trkTca=(*ntp.trk);
01272 Int_t numTrks=evt.ntrack;
01273
01274
01275 for (Int_t itrk=0;itrk<numTrks;itrk++){
01276 const NtpSRTrack* ptrk=
01277 dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
01278 const NtpSRTrack& trk=(*ptrk);
01279
01280 TClonesArray& stpTca=(*ntp.stp);
01281
01282
01283 for (Int_t i=0;i<trk.nstrip;i++){
01284
01285
01286 if (trk.stp[i]<0) {
01287 MAXMSG("MeuCuts",Msg::kInfo,50)
01288 <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
01289 continue;
01290 }
01291 const NtpSRStrip* pstp=
01292 dynamic_cast<NtpSRStrip*>(stpTca[trk.stp[i]]);
01293 const NtpSRStrip& stp=(*pstp);
01294
01295
01296 if (stp.time1>=0 && stp.time1>latestTime) latestTime=stp.time1;
01297 if (stp.time0>=0 && stp.time0>latestTime) latestTime=stp.time0;
01298 if (stp.time1>=0 && stp.time1<earliestTime) earliestTime=
01299 stp.time1;
01300 if (stp.time0>=0 && stp.time0<earliestTime) earliestTime=
01301 stp.time0;
01302
01303 Double_t time=stp.time0;
01304 if (time<0 || time>1) time=stp.time1;
01305
01306 if (time>0 && time<=1) {
01307 times.insert(time);
01308 }
01309 }
01310 }
01311
01312
01313 it=times.begin();
01314 advance(it,times.size()/2);
01315 Double_t medianTimeTrk=*it;
01316
01317 MAXMSG("MeuCuts",Msg::kDebug,500)
01318 <<"Median time for slc and trk difference="
01319 <<(medianTime-medianTimeTrk)/Munits::ns<<endl;
01320
01321 hMedianTrkSlcDiff->Fill((medianTime-medianTimeTrk)/Munits::ns);
01322
01323 hTrkTimeAfter->Fill((latestTime-medianTimeTrk)/Munits::ns);
01324 hTrkTimeBefore->Fill((medianTimeTrk-earliestTime)/Munits::ns);
01325 hTrkTimeLength->Fill((latestTime-earliestTime)/Munits::ns);
01326 }
01327
01328
01329
01330 Bool_t MeuCuts::FilterBadTrkTimes(const NtpStRecord& ntp,MeuSummary& s,
01331 const NtpSREvent& evt) const
01332 {
01333 Bool_t badTrkTimes=false;
01334
01335
01336 Double_t earliestTime=999;
01337 Double_t latestTime=-999;
01338 multiset<Double_t> times;
01339
01340 TClonesArray& trkTca=(*ntp.trk);
01341 Int_t numTrks=evt.ntrack;
01342
01343 for (Int_t itrk=0;itrk<numTrks;itrk++){
01344 const NtpSRTrack* ptrk=
01345 dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
01346 const NtpSRTrack& trk=(*ptrk);
01347
01348 TClonesArray& stpTca=(*ntp.stp);
01349
01350
01351 for (Int_t i=0;i<trk.nstrip;i++){
01352
01353
01354 if (trk.stp[i]<0) {
01355 MAXMSG("MeuCuts",Msg::kInfo,50)
01356 <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
01357 continue;
01358 }
01359
01360 const NtpSRStrip* pstp=
01361 dynamic_cast<NtpSRStrip*>(stpTca[trk.stp[i]]);
01362 const NtpSRStrip& stp=(*pstp);
01363
01364
01365 if (stp.time1>=0 && stp.time1>latestTime) latestTime=stp.time1;
01366 if (stp.time0>=0 && stp.time0>latestTime) latestTime=stp.time0;
01367 if (stp.time1>=0 && stp.time1<earliestTime) earliestTime=
01368 stp.time1;
01369 if (stp.time0>=0 && stp.time0<earliestTime) earliestTime=
01370 stp.time0;
01371
01372 Double_t time=stp.time0;
01373 if (time<0 || time>1) time=stp.time1;
01374
01375 if (time>0 && time<=1) {
01376 times.insert(time);
01377 }
01378 }
01379 }
01380
01381
01382 multiset<Double_t>::const_iterator it=times.begin();
01383 advance(it,times.size()/2);
01384 const Double_t medianTimeTrk=*it;
01385
01386 if (latestTime>medianTimeTrk+(400*Munits::ns) ||
01387 earliestTime<medianTimeTrk-(400*Munits::ns)) {
01388 MAXMSG("MeuCuts",Msg::kDebug,200)
01389 <<"Median time="<<medianTimeTrk
01390 <<", dBefore="<<(medianTimeTrk-earliestTime)/Munits::ns
01391 <<", dAfter="<<(latestTime-medianTimeTrk)/Munits::ns<<endl;
01392 badTrkTimes=true;
01393 }
01394
01395
01396 s.MedianTime=medianTimeTrk;
01397
01398 return badTrkTimes;
01399 }
01400
01401
01402
01403 Bool_t MeuCuts::FilterBadShwDist(const NtpStRecord& ntp,MeuSummary& s,
01404 const NtpSREvent& evt,
01405 Int_t* shwDist) const
01406 {
01407
01408 if (s.Detector!=Detector::kNear) {
01409 MAXMSG("MeuCuts",Msg::kInfo,1)
01410 <<"Not running FilterBadShwDist for detector="<<s.Detector
01411 <<", it is only implemented for ND spill data"<<endl;
01412 return false;
01413 }
01414
01415
01416
01417 if (s.TrigSrc<100) {
01418 MAXMSG("MeuCuts",Msg::kInfo,1)
01419 <<"Not doing FilterBadShwDist, trigSrc="<<s.TrigSrc
01420 <<", simFlag="<<s.SimFlag<<endl;
01421 return false;
01422 }
01423
01424 Bool_t badShwDist=false;
01425
01426 Int_t shwVtxPl=999;
01427 Int_t shwEndPl=-999;
01428
01429
01430 TClonesArray& shwTca=(*ntp.shw);
01431 for(int i=0;i<evt.nshower;i++) {
01432 const NtpSRShower* pshw=
01433 dynamic_cast<NtpSRShower*>(shwTca[evt.shw[i]]);
01434 const NtpSRShower& shw=(*pshw);
01435
01436
01437 Bool_t vtxShw=abs(s.VtxPlane-shw.plane.beg)<=5;
01438
01439 if (vtxShw){
01440 if (shw.plane.beg<shwVtxPl) shwVtxPl=shw.plane.beg;
01441 if (shw.plane.end>shwEndPl) shwEndPl=shw.plane.end;
01442 }
01443 }
01444
01445 Bool_t vtxShw=false;
01446 if (shwEndPl!=-999 && shwVtxPl!=999) vtxShw=true;
01447
01448 Int_t winPlFromShwEnd=999;
01449
01450
01451 if (vtxShw){
01452 winPlFromShwEnd=s.WinVtxSidePl-shwEndPl;
01453
01454 if (winPlFromShwEnd<10) badShwDist=true;
01455 }
01456
01457
01458 if (shwDist) *shwDist=winPlFromShwEnd;
01459
01460 return badShwDist;
01461 }
01462
01463
01464
01465 void MeuCuts::CalcXYZ(const AtmosEvent& ,
01466 const AtmosStrip& st,
01467 Float_t& x,Float_t& y,Float_t& z) const
01468 {
01469 VldTimeStamp vldts;
01470
01471
01472 VldContext vc(Detector::kFar,SimFlag::kData,vldts);
01473
01474 UgliGeomHandle ugh(vc);
01475 TVector3 uvz;
01476 uvz.SetZ(st.Z);
01477 if ((st.View+2)==PlaneView::kU){
01478 uvz.SetX(st.T);
01479 uvz.SetY(st.L);
01480 }
01481 else if ((st.View+2)==PlaneView::kV){
01482 uvz.SetX(st.L);
01483 uvz.SetY(st.T);
01484 }
01485 else cout<<"Ahhhhhh"<<endl;
01486 TVector3 xyz = ugh.uvz2xyz(uvz);
01487
01488 x=xyz.X();
01489 y=xyz.Y();
01490 z=xyz.Z();
01491 }
01492
01493
01494
01495 Bool_t MeuCuts::AnalyseCoilProximity(AtmosEvent& ev,const MeuSummary& s,
01496 Float_t* r) const
01497 {
01498 Bool_t awayFromCoil=false;
01499 std::string sLogLevel="kDebug";
01500
01501 Msg::LogLevel_t logLevel=Msg::kInfo;
01502 if (sLogLevel=="kDebug") logLevel=Msg::kDebug;
01503 if (sLogLevel=="kVerbose") logLevel=Msg::kVerbose;
01504
01505 if (s.VtxPlane<0 || s.EndPlane<0){
01506 MSG("MeuCuts",Msg::kWarning)
01507 <<"No vtx and/or end plane for this event!"<<endl;
01508 return awayFromCoil;
01509 }
01510
01511 map<Int_t,Float_t> plEnDep;
01512 map<Int_t,TVector3> plPosition;
01513 Float_t smallestRadius=99999;
01514
01515 TClonesArray& strips=(*ev.StripList);
01516 Int_t numStrips=strips.GetEntries();
01517
01519
01521 for (Int_t hit=0;hit<numStrips;hit++){
01522 const AtmosStrip* strip=dynamic_cast<AtmosStrip*>(strips[hit]);
01523 const AtmosStrip& st=(*strip);
01524
01525
01526 if (!st.Trk) continue;
01527
01528
01529 Float_t x=-99999;
01530 Float_t y=-99999;
01531 Float_t z=-99999;
01532 this->CalcXYZ(ev,st,x,y,z);
01533 if (fabs(x)>6 || fabs(y)>6) MSG("MeuCuts",Msg::kWarning)
01534 <<"silly x or y; x="<<x<<", y="<<y<<endl;
01535
01536 if (!(plEnDep.count(st.Plane))) plEnDep[st.Plane]+=(st.QPEcorr[0]+st.QPEcorr[1]);
01537 else {
01538
01539 MSG("MeuCuts",logLevel)
01540 <<"Already filled en dep for this plane="<<st.Plane<<endl;
01541 }
01542
01543
01544
01545 if (plPosition.count(st.Plane)) {
01546
01547 MSG("MeuCuts",logLevel)<<"already filled plane"<<endl;
01548 }
01549 else plPosition[st.Plane]=TVector3(x,y,z);
01550
01551
01552 Float_t radius=sqrt(pow(x,2)+pow(y,2));
01553 if (radius<smallestRadius) smallestRadius=radius;
01554
01555 }
01556
01557
01558 if (smallestRadius>s.RFidCoil && smallestRadius<5) awayFromCoil=true;
01559
01560 MAXMSG("MeuCuts",logLevel,200)
01561 <<"Coil: smallest radius="<<smallestRadius
01562 <<", fRFidCoil="<<s.RFidCoil
01563 <<", awayFromCoil="<<awayFromCoil
01564 <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane<<endl;
01565
01566
01567 if (r) *r=smallestRadius;
01568
01569
01570 Bool_t printDiagnostics=false;
01571 if (!awayFromCoil && printDiagnostics){
01572 MAXMSG("MeuCuts",Msg::kInfo,5000)
01573 <<"Potentially hit coil: smallest radius="<<smallestRadius
01574 <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane<<endl;
01575
01576
01577 Bool_t tmpPosDir=true;
01578 if (s.VtxPlane-s.EndPlane>0) tmpPosDir=false;
01579 const Bool_t posDir=tmpPosDir;
01580
01582
01584 Int_t pl=s.VtxPlane;
01585 Int_t endPlane=s.EndPlane-1;
01586 if (posDir) endPlane=s.EndPlane+1;
01587 while (pl!=endPlane){
01588
01589 Float_t x=-99999;
01590 Float_t y=-99999;
01591 Float_t z=-99999;
01592 Float_t radius=-99999;
01593
01594
01595 map<Int_t,TVector3>::iterator plPosIt=plPosition.find(pl);
01596 map<Int_t,TVector3>::iterator plPosEndIt=plPosition.end();
01597 if (plPosIt!=plPosEndIt){
01598 TVector3& plPos=plPosIt->second;
01599 x=plPos.X();
01600 y=plPos.Y();
01601 z=plPos.Z();
01602
01603
01604 radius=sqrt(pow(x,2)+pow(y,2));
01605 }
01606
01607 MAXMSG("MeuCuts",Msg::kInfo,5000)
01608 <<"pl="<<pl<<": r="<<radius
01609 <<", x="<<x<<", y="<<y<<", z="<<z<<endl;
01610
01611
01612 if (posDir) pl++;
01613 else pl--;
01614 }
01615 }
01616
01617 return awayFromCoil;
01618 }
01619
01620
01621
01622 Bool_t MeuCuts::IsCorrectTrigSrc(const MeuSummary& s) const
01623 {
01624
01625 if (s.Detector==Detector::kFar && s.TrigSrc>100) {
01626 MAXMSG("MeuCuts",Msg::kInfo,10)
01627 <<"Not using this event, detector="<<s.Detector
01628 <<", trigSrc="<<s.TrigSrc<<endl;
01629 return false;
01630 }
01631 else return true;
01632 }
01633
01634
01635
01636 Bool_t MeuCuts::IsAwayFromCoil
01637 (const MeuSummary& s,const std::map<Int_t,MeuHitInfo>& plInfo,
01638 Float_t* r) const
01639 {
01640 if (s.Detector!=Detector::kFar) {
01641 MAXMSG("MeuCuts",Msg::kInfo,1)
01642 <<"Not running IsAwayFromCoil for detector="<<s.Detector
01643 <<", only needed for FD"<<endl;
01644 return true;
01645 }
01646
01647 Bool_t awayFromCoil=false;
01648 Msg::LogLevel_t logLevel=Msg::kDebug;
01649
01650 if (s.VtxPlane<0 || s.EndPlane<0){
01651 MSG("MeuCuts",Msg::kWarning)
01652 <<"No vtx and/or end plane for this event!"<<endl;
01653 return awayFromCoil;
01654 }
01655
01656 Float_t smallestRadius=99999;
01657
01659
01661 typedef map<Int_t,MeuHitInfo>::const_iterator plIter;
01662 for (plIter it=plInfo.begin();it!=plInfo.end();++it){
01663
01664 const MeuHitInfo& cp=it->second;
01665
01666
01667 Float_t radius=sqrt(pow(cp.X,2)+pow(cp.Y,2));
01668 if (radius<smallestRadius) smallestRadius=radius;
01669
01670 }
01671
01672
01673 if (smallestRadius>s.RFidCoil && smallestRadius<5) awayFromCoil=true;
01674
01675 MAXMSG("MeuCuts",logLevel,200)
01676 <<"Coil: smallest radius="<<smallestRadius
01677 <<", fRFidCoil="<<s.RFidCoil
01678 <<", awayFromCoil="<<awayFromCoil
01679 <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane<<endl;
01680
01681
01682 if (r) *r=smallestRadius;
01683
01684 return awayFromCoil;
01685 }
01686
01687
01688
01689 void MeuCuts::GetBDSelectSpillInfo(const NtpBDLiteRecord& ntpBD,
01690 MeuSummary& s,Int_t entry,
01691 Int_t& goodSpillCounter) const
01692 {
01693 static BeamMonSpill spill;
01694 static BMSpillAna bmsa;
01695
01696
01697 static Registry r;
01698 r.UnLockValues();
01699 r.Set("TargetIn",1);
01700 r.Set("BeamType",-1);
01701 r.LockValues();
01702 bmsa.Config(r);
01703
01704 bmsa.SetSpill(ntpBD,spill);
01705 s.BDSelectSpill=bmsa.SelectSpill();
01706
01707 MAXMSG("MeuAnalysis",Msg::kInfo,1)
01708 <<"Setting first fBDSelectSpill... ="<<s.BDSelectSpill
01709 <<", tortgt="<<ntpBD.tortgt
01710 <<", timeDiff="<<ntpBD.GetHeader().GetTimeDiffStreamSpill()
01711 <<endl;
01712
01713 static Int_t lastEntry=-1;
01714 if (!s.BDSelectSpill) {
01715 if (entry!=lastEntry) {
01716 MAXMSG("MeuAnalysis",Msg::kInfo,10)
01717 <<"bmsa.SelectSpill()="<<s.BDSelectSpill
01718 <<", tortgt="<<ntpBD.tortgt<<endl;
01719 }
01720 }
01721 else {
01722
01723 if (entry!=lastEntry) goodSpillCounter++;
01724 }
01725 lastEntry=entry;
01726 }
01727
01728
01729
01730 void MeuCuts::ConvertToXY
01731 (const MeuSummary& s,std::map<Int_t,MeuHitInfo>& plInfo) const
01732 {
01733
01734 Msg::LogLevel_t logLevel=Msg::kDebug;
01735
01736 MAXMSG("MeuCuts",logLevel,1000)
01737 <<endl<<"ConvertToXY: Turn TPos, LPos and z -> xyz"<<endl;
01738
01739 Float_t initValue=-999;
01740 Bool_t posDir=(s.EndPlane>s.VtxPlane);
01741
01742
01743 Int_t pl=s.VtxPlane;
01744 Int_t tmpEndPlane=s.EndPlane;
01745 if (posDir) tmpEndPlane++;
01746 else tmpEndPlane--;
01747
01748 VldTimeStamp vldts;
01749
01750 VldContext vc(static_cast<Detector::Detector_t>(s.Detector),
01751 SimFlag::kData,vldts);
01752
01753 UgliGeomHandle ugh(vc);
01754 TVector3 uvz;
01755
01756
01757
01758 typedef map<Int_t,MeuHitInfo>::iterator plIt;
01759 for (plIt it=plInfo.begin();it!=plInfo.end();++it){
01760
01761
01762 MeuHitInfo& cp=it->second;
01763
01764 if (cp.LPos==-999) continue;
01765
01766
01767 uvz.SetZ(cp.Z);
01768
01769
01770 if (cp.View==PlaneView::kU){
01771 uvz.SetX(cp.TPos);
01772 uvz.SetY(cp.LPos);
01773 }
01774 else if (cp.View==PlaneView::kV){
01775 uvz.SetY(cp.TPos);
01776 uvz.SetX(cp.LPos);
01777 }
01778 else cout<<"ahhhhhhh no plane view"<<endl;
01779
01780
01781 TVector3 xyz=ugh.uvz2xyz(uvz);
01782
01783
01784 cp.X=xyz.X();
01785 cp.Y=xyz.Y();
01786
01787 MAXMSG("MeuCuts",logLevel,1000)
01788 <<"pl="<<pl<<", cp: X="<<cp.X<<", Y="<<cp.Y<<", Z="<<cp.Z
01789 <<", LPos="<<cp.LPos<<endl;
01790
01791
01792
01793 cp.LPos=initValue;
01794 MAXMSG("MeuCuts",Msg::kInfo,1)
01795 <<"Doing nasty hack for AtmosEvent X and Y pos"<<endl;
01796 }
01797 }
01798
01799
01800
01801 void MeuCuts::ExtractPlInfo(const AtmosEvent& ev,MeuSummary& s,
01802 std::map<Int_t,MeuHitInfo>& plInfo) const
01803 {
01804
01805
01806
01807 map<Int_t,MeuHitInfo> plSigCorLast;
01808 map<Int_t,Double_t> time0;
01809 map<Int_t,Double_t> time1;
01810
01811 TClonesArray& strips=(*ev.StripList);
01812 Int_t numStrips=strips.GetEntries();
01814
01816 for (Int_t hit=0;hit<numStrips;hit++){
01817 const AtmosStrip* strip=dynamic_cast<AtmosStrip*>(strips[hit]);
01818 const AtmosStrip& st=(*strip);
01819
01820 Int_t pl=st.Plane;
01821 MeuHitInfo& cp=plInfo[pl];
01822
01823
01824 cp.Plane=pl;
01825 cp.View=(st.View+2);
01826 cp.Z=st.Z;
01827
01828
01829 cp.LPos=st.L;
01830
01831
01832 cp.Adc+=(st.Qadc[1]+st.Qadc[0]);
01833 cp.Adc1+=st.Qadc[0];
01834 cp.Adc2+=st.Qadc[1];
01835 cp.Pe+=st.QPE[0]+st.QPE[1];
01836 cp.Pe1+=st.QPE[0];
01837 cp.Pe2+=st.QPE[1];
01838 cp.SigLin+=(st.Qadc[1]+st.Qadc[0]);
01839 cp.SigLin1+=st.Qadc[0];
01840 cp.SigLin2+=st.Qadc[1];
01841 cp.SigCor+=(st.QPEcorr[0]+st.QPEcorr[1]);
01842 cp.SigCor1+=st.QPEcorr[0];
01843 cp.SigCor2+=st.QPEcorr[1];
01844
01845
01846 cp.NumDigits+=st.Ndigits;
01847 cp.NumStrips++;
01848
01849
01850 if (!st.Trk) continue;
01851
01852
01853 if (st.View==0){
01854 if (pl>s.MaxPlane2) s.MaxPlane2=pl;
01855 if (pl<s.MinPlane2) s.MinPlane2=pl;
01856 }
01857 else if (st.View==1){
01858 if (pl>s.MaxPlane3) s.MaxPlane3=pl;
01859 if (pl<s.MinPlane3) s.MinPlane3=pl;
01860 }
01861 else cout<<"No view"<<endl;
01862
01863
01864 MeuHitInfo& tmpcp=plSigCorLast[pl];
01865
01866 if (cp.Strip<0){
01867
01868 cp.Strip=st.Strip;
01869 cp.TPos=st.T;
01870
01871
01872
01873 cp.SigCorTrk1=st.QPEcorr[0];
01874 cp.SigCorTrk2=st.QPEcorr[1];
01875
01876
01877 tmpcp.Strip=st.Strip;
01878 tmpcp.TPos=st.T;
01879 tmpcp.SigCor=(st.QPEcorr[0]+st.QPEcorr[1]);
01880 time0[pl]=st.Tcal[1];
01881 time1[pl]=st.Tcal[0];
01882 }
01883 else{
01884 MAXMSG("MeuCuts",Msg::kWarning,10)
01885 <<"Choosing between strips... run="<<ev.Run
01886 <<", sn="<<ev.Snarl<<endl
01887 <<" prev: pl="<<pl<<", st="<<tmpcp.Strip
01888 <<", sigcor="<<tmpcp.SigCor
01889 <<", t0="<<time0[pl]*1e9<<", t1="<<time1[pl]*1e9
01890 <<endl
01891 <<" next: pl="<<pl<<", st="<<st.Strip
01892 <<", sigcor="<<(st.QPEcorr[0]+st.QPEcorr[1])
01893 <<", t0="<<st.Tcal[1]*1e9<<", t1="<<st.Tcal[0]*1e9<<endl;
01894
01895
01896 if ((st.QPEcorr[0]+st.QPEcorr[1])>tmpcp.SigCor){
01897 cp.Strip=st.Strip;
01898 cp.TPos=st.T;
01899
01900
01901
01902 cp.SigCorTrk1=st.QPEcorr[0];
01903 cp.SigCorTrk2=st.QPEcorr[1];
01904
01905
01906 tmpcp.Strip=st.Strip;
01907 tmpcp.TPos=st.T;
01908 tmpcp.SigCor=(st.QPEcorr[0]+st.QPEcorr[1]);
01909 time0[pl]=st.Tcal[1];
01910 time1[pl]=st.Tcal[0];
01911 MAXMSG("MeuCuts",Msg::kDebug,1000)
01912 <<"...using new strip"<<endl;
01913 }
01914 else{
01915 MAXMSG("MeuCuts",Msg::kDebug,1000)
01916 <<"...using previous strip"<<endl;
01917 }
01918 }
01919 }
01920
01921
01922 this->ConvertToXY(s,plInfo);
01923 }
01924
01925
01926
01927 void MeuCuts::ExtractPlInfo(const NtpStRecord& ntp,MeuSummary& s,
01928 std::map<Int_t,MeuHitInfo>& plInfo,
01929 const NtpSREvent& evt) const
01930 {
01931
01932 Msg::LogLevel_t logLevel=Msg::kDebug;
01933
01934
01935 static TH1F* hAdcPerSigLinNtp1=0;
01936 static TH1F* hSigLinPerSigCorNtp1=0;
01937 static TH1F* hSigCorPerSigMapNtp1=0;
01938 static TH1F* hSigMapPerMipNtp1=0;
01939 static TH1F* hMipPerGeVNtp1=0;
01940 static TH1F* hAdcPerPeNtp1=0;
01941
01942 static TH1F* hAdcPerSigLinNtp2=0;
01943 static TH1F* hSigLinPerSigCorNtp2=0;
01944 static TH1F* hSigCorPerSigMapNtp2=0;
01945 static TH1F* hSigMapPerMipNtp2=0;
01946 static TH1F* hMipPerGeVNtp2=0;
01947 static TH1F* hAdcPerPeNtp2=0;
01948
01949 if (hAdcPerSigLinNtp1==0) {
01950
01951 hAdcPerSigLinNtp1=new TH1F("hAdcPerSigLinNtp1","hAdcPerSigLinNtp1",
01952 1100,-1,10);
01953 hAdcPerSigLinNtp1->SetFillColor(0);
01954 hAdcPerSigLinNtp1->SetTitle("AdcPerSigLin");
01955 hAdcPerSigLinNtp1->GetXaxis()->SetTitle("AdcPerSigLin");
01956 hAdcPerSigLinNtp1->GetXaxis()->CenterTitle();
01957
01958
01959 hSigLinPerSigCorNtp1=new TH1F("hSigLinPerSigCorNtp1",
01960 "hSigLinPerSigCorNtp1",
01961 1100,-1,10);
01962 hSigLinPerSigCorNtp1->SetFillColor(0);
01963 hSigLinPerSigCorNtp1->SetTitle("SigLinPerSigCor");
01964 hSigLinPerSigCorNtp1->GetXaxis()->SetTitle("SigLinPerSigCor");
01965 hSigLinPerSigCorNtp1->GetXaxis()->CenterTitle();
01966
01967
01968 hSigCorPerSigMapNtp1=new TH1F("hSigCorPerSigMapNtp1",
01969 "hSigCorPerSigMapNtp1",
01970 10000,-1,100);
01971 hSigCorPerSigMapNtp1->SetFillColor(0);
01972 hSigCorPerSigMapNtp1->SetTitle("SigCorPerSigMap");
01973 hSigCorPerSigMapNtp1->GetXaxis()->SetTitle("SigCorPerSigMap");
01974 hSigCorPerSigMapNtp1->GetXaxis()->CenterTitle();
01975
01976
01977 hSigMapPerMipNtp1=new TH1F("hSigMapPerMipNtp1",
01978 "hSigMapPerMipNtp1",
01979 10100,-10,1000);
01980 hSigMapPerMipNtp1->SetFillColor(0);
01981 hSigMapPerMipNtp1->SetTitle("SigMapPerMip");
01982 hSigMapPerMipNtp1->GetXaxis()->SetTitle("SigMapPerMip");
01983 hSigMapPerMipNtp1->GetXaxis()->CenterTitle();
01984
01985
01986 hMipPerGeVNtp1=new TH1F("hMipPerGeVNtp1","hMipPerGeVNtp1",
01987 10100,-10,1000);
01988 hMipPerGeVNtp1->SetFillColor(0);
01989 hMipPerGeVNtp1->SetTitle("MipPerGeV");
01990 hMipPerGeVNtp1->GetXaxis()->SetTitle("MipPerGeV");
01991 hMipPerGeVNtp1->GetXaxis()->CenterTitle();
01992
01993
01994 hAdcPerPeNtp1=new TH1F("hAdcPerPeNtp1","hAdcPerPeNtp1",
01995 8080,-20,2000);
01996 hAdcPerPeNtp1->SetFillColor(0);
01997 hAdcPerPeNtp1->SetTitle("AdcPerPe");
01998 hAdcPerPeNtp1->GetXaxis()->SetTitle("AdcPerPe");
01999 hAdcPerPeNtp1->GetXaxis()->CenterTitle();
02000
02001
02002
02003
02004 hAdcPerSigLinNtp2=new TH1F("hAdcPerSigLinNtp2","hAdcPerSigLinNtp2",
02005 1100,-1,10);
02006 hAdcPerSigLinNtp2->SetFillColor(0);
02007 hAdcPerSigLinNtp2->SetTitle("AdcPerSigLin");
02008 hAdcPerSigLinNtp2->GetXaxis()->SetTitle("AdcPerSigLin");
02009 hAdcPerSigLinNtp2->GetXaxis()->CenterTitle();
02010
02011
02012 hSigLinPerSigCorNtp2=new TH1F("hSigLinPerSigCorNtp2",
02013 "hSigLinPerSigCorNtp2",
02014 1100,-1,10);
02015 hSigLinPerSigCorNtp2->SetFillColor(0);
02016 hSigLinPerSigCorNtp2->SetTitle("SigLinPerSigCor");
02017 hSigLinPerSigCorNtp2->GetXaxis()->SetTitle("SigLinPerSigCor");
02018 hSigLinPerSigCorNtp2->GetXaxis()->CenterTitle();
02019
02020
02021 hSigCorPerSigMapNtp2=new TH1F("hSigCorPerSigMapNtp2",
02022 "hSigCorPerSigMapNtp2",
02023 10000,-1,100);
02024 hSigCorPerSigMapNtp2->SetFillColor(0);
02025 hSigCorPerSigMapNtp2->SetTitle("SigCorPerSigMap");
02026 hSigCorPerSigMapNtp2->GetXaxis()->SetTitle("SigCorPerSigMap");
02027 hSigCorPerSigMapNtp2->GetXaxis()->CenterTitle();
02028
02029
02030 hSigMapPerMipNtp2=new TH1F("hSigMapPerMipNtp2",
02031 "hSigMapPerMipNtp2",
02032 10100,-10,1000);
02033 hSigMapPerMipNtp2->SetFillColor(0);
02034 hSigMapPerMipNtp2->SetTitle("SigMapPerMip");
02035 hSigMapPerMipNtp2->GetXaxis()->SetTitle("SigMapPerMip");
02036 hSigMapPerMipNtp2->GetXaxis()->CenterTitle();
02037
02038
02039 hMipPerGeVNtp2=new TH1F("hMipPerGeVNtp2","hMipPerGeVNtp2",
02040 10100,-10,1000);
02041 hMipPerGeVNtp2->SetFillColor(0);
02042 hMipPerGeVNtp2->SetTitle("MipPerGeV");
02043 hMipPerGeVNtp2->GetXaxis()->SetTitle("MipPerGeV");
02044 hMipPerGeVNtp2->GetXaxis()->CenterTitle();
02045
02046
02047 hAdcPerPeNtp2=new TH1F("hAdcPerPeNtp2","hAdcPerPeNtp2",
02048 8080,-20,2000);
02049 hAdcPerPeNtp2->SetFillColor(0);
02050 hAdcPerPeNtp2->SetTitle("AdcPerPe");
02051 hAdcPerPeNtp2->GetXaxis()->SetTitle("AdcPerPe");
02052 hAdcPerPeNtp2->GetXaxis()->CenterTitle();
02053
02054 }
02055
02056 map<Int_t,MeuHitInfo> plSigCorLast;
02057 map<Int_t,Double_t> time0;
02058 map<Int_t,Double_t> time1;
02059
02060 TClonesArray& stpTca=(*ntp.stp);
02061
02062
02063 TClonesArray& slcTca=(*ntp.slc);
02064 const Int_t numSlcs=slcTca.GetEntries();
02065 if (evt.slc>=numSlcs) cout<<"Ahhhhhh, num slcs"<<endl;
02066
02067
02068 const NtpSRSlice* pslc=
02069 dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]);
02070 const NtpSRSlice& slc=(*pslc);
02071
02072
02073 for (Int_t i=0;i<slc.nstrip;++i) {
02074 const NtpSRStrip* pstp=
02075 dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
02076 const NtpSRStrip& stp=(*pstp);
02077
02078
02079 Double_t time=stp.time0;
02080 if (time<0 || time>1) time=stp.time1;
02081 if (time<0 || time>1) cout<<"Ahhhhh, time"<<endl;;
02082
02083 if (time>s.MedianTime+(400*Munits::ns) ||
02084 time<s.MedianTime-(400*Munits::ns)) {
02085 MAXMSG("MeuCuts",Msg::kInfo,10)
02086 <<"Cut strip: time="<<time
02087 <<", pl="<<stp.plane<<", st="<<stp.strip
02088 <<", sigCor="<<stp.ph0.sigcor+stp.ph1.sigcor
02089 <<", dT="<<(time-s.MedianTime)/Munits::ns<<endl;
02090 continue;
02091 }
02092
02093 Int_t pl=stp.plane;
02094 MeuHitInfo& cp=plInfo[pl];
02095
02096
02097 cp.Plane=pl;
02098 cp.View=stp.planeview;
02099
02100
02101
02102 cp.Adc+=stp.ph0.raw+stp.ph1.raw;
02103 cp.Adc1+=stp.ph0.raw;
02104 cp.Adc2+=stp.ph1.raw;
02105 cp.Pe+=stp.ph0.pe+stp.ph1.pe;
02106 cp.Pe1+=stp.ph0.pe;
02107 cp.Pe2+=stp.ph1.pe;
02108 cp.SigLin+=stp.ph0.siglin+stp.ph1.siglin;
02109 cp.SigLin1+=stp.ph0.siglin;
02110 cp.SigLin2+=stp.ph1.siglin;
02111 cp.SigCor+=stp.ph0.sigcor+stp.ph1.sigcor;
02112 cp.SigCor1+=stp.ph0.sigcor;
02113 cp.SigCor2+=stp.ph1.sigcor;
02114
02115
02116 if (stp.ph0.siglin) hAdcPerSigLinNtp1->
02117 Fill(stp.ph0.raw/stp.ph0.siglin);
02118 if (stp.ph0.sigcor) hSigLinPerSigCorNtp1->
02119 Fill(stp.ph0.siglin/stp.ph0.sigcor);
02120 if (stp.ph0.pe) hAdcPerPeNtp1->
02121 Fill(stp.ph0.raw/stp.ph0.pe);
02122
02123 if (stp.ph1.siglin) hAdcPerSigLinNtp2->
02124 Fill(stp.ph1.raw/stp.ph1.siglin);
02125 if (stp.ph1.sigcor) hSigLinPerSigCorNtp2->
02126 Fill(stp.ph1.siglin/stp.ph1.sigcor);
02127 if (stp.ph1.pe) hAdcPerPeNtp2->
02128 Fill(stp.ph1.raw/stp.ph1.pe);
02129
02130
02131 cp.NumDigits+=stp.ndigit;
02132
02133 cp.NumStrips++;
02134 }
02135
02136 TClonesArray& trkTca=(*ntp.trk);
02137
02138 Int_t numTrks=evt.ntrack;
02139
02140
02141
02142
02143
02144
02145 for (Int_t itrk=0;itrk<numTrks;itrk++){
02146 const NtpSRTrack* ptrk=
02147 dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
02148 const NtpSRTrack& trk=(*ptrk);
02149
02150 TClonesArray& stpTca=(*ntp.stp);
02151
02152
02153 for (Int_t i=0;i<trk.nstrip;i++){
02154
02155
02156 if (trk.stp[i]<0) {
02157 MAXMSG("MeuCuts",Msg::kInfo,50)
02158 <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
02159 continue;
02160 }
02161
02162 const NtpSRStrip* pstp=
02163 dynamic_cast<NtpSRStrip*>(stpTca[trk.stp[i]]);
02164 const NtpSRStrip& stp=(*pstp);
02165
02166
02167
02168
02169
02170
02171
02172
02173 if (trk.stpph0sigmap[i]) hSigCorPerSigMapNtp1->
02174 Fill(stp.ph0.sigcor/
02175 trk.stpph0sigmap[i]);
02176 if (trk.stpph0mip[i]) hSigMapPerMipNtp1->
02177 Fill(trk.stpph0sigmap[i]/
02178 trk.stpph0mip[i]);
02179 if (trk.stpph0gev[i]) hMipPerGeVNtp1->
02180 Fill(trk.stpph0mip[i]/trk.stpph0gev[i]);
02181
02182 if (trk.stpph1sigmap[i]) hSigCorPerSigMapNtp2->
02183 Fill(stp.ph1.sigcor/
02184 trk.stpph1sigmap[i]);
02185 if (trk.stpph1mip[i]) hSigMapPerMipNtp2->
02186 Fill(trk.stpph1sigmap[i]/
02187 trk.stpph1mip[i]);
02188 if (trk.stpph1gev[i]) hMipPerGeVNtp2->
02189 Fill(trk.stpph1mip[i]/trk.stpph1gev[i]);
02190
02191
02192 Double_t time=stp.time0;
02193 if (time<0 || time>1) time=stp.time1;
02194 if (time<0 || time>1) cout<<"Ahhhhh, time"<<endl;;
02195
02196
02197
02198 if (time>s.MedianTime+(400*Munits::ns) ||
02199 time<s.MedianTime-(400*Munits::ns)) {
02200 MAXMSG("MeuCuts",Msg::kWarning,200)
02201 <<"Ahhhhhhh trk: time="<<time
02202 <<", pl="<<stp.plane<<", st="<<stp.strip
02203 <<", time-medTime="<<(time-s.MedianTime)/Munits::ns<<endl;
02204 continue;
02205 }
02206
02207 Int_t pl=stp.plane;
02208
02209
02210 if (s.Detector!=Detector::kNear ||
02211 (s.Detector==Detector::kNear && pl<=120)) {
02212 if (stp.planeview==PlaneView::kU){
02213 if (pl>s.MaxPlane2) s.MaxPlane2=pl;
02214 if (pl<s.MinPlane2) s.MinPlane2=pl;
02215 }
02216 else if (stp.planeview==PlaneView::kV){
02217 if (pl>s.MaxPlane3) s.MaxPlane3=pl;
02218 if (pl<s.MinPlane3) s.MinPlane3=pl;
02219 }
02220 else cout<<"No view"<<endl;
02221 }
02222
02223
02224 MeuHitInfo& cp=plInfo[pl];
02225
02226
02227 MeuHitInfo& tmpcp=plSigCorLast[pl];
02228
02229 if (cp.Strip<0){
02230
02231 cp.Strip=stp.strip;
02232 cp.TPos=stp.tpos;
02233
02234
02235
02236 cp.SigCorTrk1=stp.ph0.sigcor;
02237 cp.SigCorTrk2=stp.ph1.sigcor;
02238
02239 cp.X=trk.stpx[i];
02240 cp.Y=trk.stpy[i];
02241 cp.Z=trk.stpz[i];
02242
02243
02244 tmpcp.Strip=stp.strip;
02245 tmpcp.TPos=stp.tpos;
02246 tmpcp.SigCor=stp.ph0.sigcor+stp.ph1.sigcor;
02247 tmpcp.X=trk.stpx[i];
02248 tmpcp.Y=trk.stpy[i];
02249 tmpcp.Z=trk.stpz[i];
02250
02251
02252
02253
02254 time0[pl]=stp.time0;
02255 time1[pl]=stp.time1;
02256 }
02257 else{
02258 MAXMSG("MeuCuts",logLevel,100)
02259 <<"Choosing between strips... run="
02260
02261 <<" prev: pl="<<pl<<", st="<<tmpcp.Strip
02262 <<", sigcor="<<tmpcp.SigCor
02263 <<", t0="<<time0[pl]*1e9<<", t1="<<time1[pl]*1e9
02264 <<endl
02265 <<" next: pl="<<pl<<", st="<<stp.strip
02266 <<", sigcor="<<stp.ph0.sigcor+stp.ph1.sigcor
02267 <<", t0="<<stp.time0*1e9<<", t1="<<stp.time1*1e9
02268 <<endl;
02269
02270
02271 if (stp.ph0.sigcor+stp.ph1.sigcor>tmpcp.SigCor){
02272 cp.Strip=stp.strip;
02273 cp.TPos=stp.tpos;
02274
02275
02276
02277 cp.SigCorTrk1=stp.ph0.sigcor;
02278 cp.SigCorTrk2=stp.ph1.sigcor;
02279
02280 cp.X=trk.stpx[i];
02281 cp.Y=trk.stpy[i];
02282 cp.Z=trk.stpz[i];
02283
02284
02285 tmpcp.Strip=stp.strip;
02286 tmpcp.TPos=stp.tpos;
02287 tmpcp.SigCor=stp.ph0.sigcor+stp.ph1.sigcor;
02288 tmpcp.X=trk.stpx[i];
02289 tmpcp.Y=trk.stpy[i];
02290 tmpcp.Z=trk.stpz[i];
02291 time0[pl]=stp.time0;
02292 time1[pl]=stp.time1;
02293 MAXMSG("MeuCuts",logLevel,100)
02294 <<"...using new strip="<<cp.Strip<<endl;
02295 }
02296 else{
02297 MAXMSG("MeuCuts",logLevel,100)
02298 <<"...using previous strip="<<cp.Strip<<endl;
02299 }
02300 }
02301 }
02302 }
02303 }
02304
02305
02306
02307 void MeuCuts::PrintNtpSt(const NtpStRecord& ntp) const
02308 {
02309 TClonesArray& tcaTk=(*ntp.trk);
02310 Int_t numTrks=tcaTk.GetEntries();
02311 TClonesArray& stpTca=(*ntp.stp);
02312 Int_t numStps=stpTca.GetEntries();
02313
02314 MAXMSG("MeuCuts",Msg::kInfo,200)
02315 <<"PrintNtpSt: numTrks="<<numTrks
02316 <<", total numStps="<<numStps<<endl;
02317
02318
02319 for (Int_t itrk=0;itrk<numTrks;itrk++){
02320 const NtpSRTrack* ptrk=
02321 dynamic_cast<NtpSRTrack*>(tcaTk[itrk]);
02322 const NtpSRTrack& trk=(*ptrk);
02323
02324 TClonesArray& stpTca=(*ntp.stp);
02325
02326
02327 for (Int_t i=0;i<trk.nstrip;i++){
02328
02329
02330 if (trk.stp[i]<0) {
02331 MAXMSG("MeuCuts",Msg::kInfo,50)
02332 <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
02333 continue;
02334 }
02335
02336 const NtpSRStrip* pstp=
02337 dynamic_cast<NtpSRStrip*>(stpTca[trk.stp[i]]);
02338 const NtpSRStrip& stp=(*pstp);
02339 MAXMSG("NDMeuCuts",Msg::kInfo,20000)
02340 <<"Trk: st="<<stp.strip<<", pl="<<stp.plane
02341 <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
02342 }
02343 }
02344
02345
02346 for (Int_t i=0;i<numStps;i++){
02347 const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(stpTca[i]);
02348 const NtpSRStrip& stp=(*pstp);
02349
02350 MAXMSG("NDMeuCuts",Msg::kInfo,20000)
02351 <<"All strips: st="<<stp.strip<<", pl="<<stp.plane
02352 <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
02353 }
02354 }
02355
02356
02357
02358 void MeuCuts::ExtractTruthInfo
02359 (const AtmosEvent& ev,MeuSummary& s,
02360 std::map<Int_t,MeuHitInfo>& plInfo) const
02361 {
02362 if (s.SimFlag!=SimFlag::kMC) {
02363 MAXMSG("MeuCuts",Msg::kInfo,1)
02364 <<endl<<"Not extracting truth info, SimFlag="<<s.SimFlag<<endl;
02365 return;
02366 }
02367
02368
02369 Float_t lowEn=999999;
02370 Float_t highEn=-1;
02371 Int_t lowEnPl=-1;
02372 Int_t highEnPl=-1;
02373 Int_t mainParticle=0;
02374 Int_t muon=13;
02375 Double_t B=0.133;
02376
02377
02378 vector<Float_t> vEn(1000);
02379 vector<Float_t> vEnSupp(1000);
02380
02381
02382 TClonesArray& scintHits=(*ev.ScintHitList);
02383 Int_t numDshs=scintHits.GetEntries();
02384
02385 for (Int_t hit=0;hit<numDshs;hit++){
02386 const AtmosScintHit* scintHit=
02387 dynamic_cast<AtmosScintHit*>(scintHits[hit]);
02388 const AtmosScintHit& dsh=(*scintHit);
02389
02390
02391 if (abs(dsh.Id)==muon) {
02392 mainParticle=dsh.Id;
02393 if (dsh.ParticleE>highEn) {
02394 highEn=dsh.ParticleE;
02395 highEnPl=dsh.Plane;
02396 }
02397 if (dsh.ParticleE<lowEn) {
02398 lowEn=dsh.ParticleE;
02399 lowEnPl=dsh.Plane;
02400 }
02401 }
02402
02403 Double_t birksSupp=1;
02404 if (dsh.DS!=0) birksSupp=1-(1./(1+B*(dsh.DE/dsh.DS)));
02405 else {
02406 MAXMSG("MeuCuts",Msg::kInfo,5)
02407 <<"Zero pathlength in DSH: dE="<<dsh.DE
02408 <<", id="<<dsh.Id<<endl;
02409 }
02410 Double_t suppEnDep=(1-birksSupp)*dsh.DE;
02411 vEnSupp[dsh.Plane]+=suppEnDep;
02412 vEn[dsh.Plane]+=dsh.DE;
02413 }
02414
02415 if (numDshs>0) {
02416
02417 if (lowEn==999999) lowEn=-1;
02418
02419
02420 s.MCHighEn=highEn;
02421 s.MCLowEn=lowEn;
02422 s.MCParticleId=mainParticle;
02423 s.MCVtxPlane=highEnPl;
02424 s.MCEndPlane=lowEnPl;
02425
02426
02427 typedef map<Int_t,MeuHitInfo>::iterator iter;
02428 for (iter it=plInfo.begin();it!=plInfo.end();++it){
02429 MeuHitInfo& cp=it->second;
02430 cp.MCEnDep+=vEn[cp.Plane];
02431 cp.MCSuppEnDep+=vEnSupp[cp.Plane];
02432 }
02433 }
02434 else {
02435 MAXMSG("MeuCuts",Msg::kInfo,1)
02436 <<endl<<"Not extracting truth info, numDshs="<<numDshs
02437 <<", SimFlag="<<s.SimFlag<<endl;
02438 }
02439 }
02440
02441
02442
02443 void MeuCuts::ExtractTruthInfo
02444 (const NtpStRecord& ntp,MeuSummary& s,
02445 std::map<Int_t,MeuHitInfo>& plInfo) const
02446 {
02447 if (s.SimFlag!=SimFlag::kMC) {
02448 MAXMSG("MeuCuts",Msg::kInfo,1)
02449 <<endl<<"Not extracting truth info, SimFlag="<<s.SimFlag<<endl;
02450 return;
02451 }
02452
02453
02454 Float_t lowEn=999999;
02455 Float_t highEn=-1;
02456 Int_t lowEnPl=-1;
02457 Int_t highEnPl=-1;
02458 Int_t mainParticle=0;
02459 Int_t muon=13;
02460 Double_t B=0.133;
02461
02462
02463 vector<Float_t> vEn(1000);
02464 vector<Float_t> vEnSupp(1000);
02465
02466 TClonesArray& tcaDsh=(*ntp.digihit);
02467 Int_t numDshs=tcaDsh.GetEntries();
02468
02469
02470 for (Int_t i=0;i<numDshs;i++){
02471
02472 const NtpMCDigiScintHit* pdsh=
02473 dynamic_cast<NtpMCDigiScintHit*>(tcaDsh[i]);
02474 const NtpMCDigiScintHit& dsh=(*pdsh);
02475
02476 MAXMSG("MeuCuts",Msg::kDebug,1000)
02477 <<"DSH: pl="<<dsh.plane<<", st="<<dsh.strip
02478 <<", pE="<<dsh.pE<<", dE="<<dsh.dE<<endl;
02479
02480
02481 if (abs(dsh.pId)==muon) {
02482 mainParticle=dsh.pId;
02483 if (dsh.pE>highEn) {
02484 highEn=dsh.pE;
02485 highEnPl=dsh.plane;
02486 }
02487 if (dsh.pE<lowEn) {
02488 lowEn=dsh.pE;
02489 lowEnPl=dsh.plane;
02490 }
02491 }
02492
02493 Double_t birksSupp=1;
02494 if (dsh.dS!=0) birksSupp=1-(1./(1+B*(dsh.dE/dsh.dS)));
02495 else {
02496 MAXMSG("MeuCuts",Msg::kInfo,5)
02497 <<"Zero pathlength in DSH: dE="<<dsh.dE
02498 <<", id="<<dsh.pId<<endl;
02499 }
02500 Double_t suppEnDep=(1-birksSupp)*dsh.dE;
02501 vEnSupp[dsh.plane]+=suppEnDep;
02502 vEn[dsh.plane]+=dsh.dE;
02503 }
02504
02505 if (numDshs>0) {
02506
02507 if (lowEn==999999) lowEn=-1;
02508
02509
02510 s.MCHighEn=highEn;
02511 s.MCLowEn=lowEn;
02512 s.MCParticleId=mainParticle;
02513 s.MCVtxPlane=highEnPl;
02514 s.MCEndPlane=lowEnPl;
02515
02516
02517 typedef map<Int_t,MeuHitInfo>::iterator iter;
02518 for (iter it=plInfo.begin();it!=plInfo.end();++it){
02519 MeuHitInfo& cp=it->second;
02520 cp.MCEnDep+=vEn[cp.Plane];
02521 cp.MCSuppEnDep+=vEnSupp[cp.Plane];
02522 }
02523 }
02524 else {
02525 MAXMSG("MeuCuts",Msg::kInfo,1)
02526 <<endl<<"Not extracting truth info, numDshs="<<numDshs
02527 <<", SimFlag="<<s.SimFlag<<endl;
02528 }
02529 }
02530
02531
02532
02533 void MeuCuts::FillEnVsDist
02534 (const MeuSummary& s,const std::map<Int_t,MeuHitInfo>& plInfo) const
02535 {
02536
02537 Msg::LogLevel_t logLevel=Msg::kDebug;
02538
02539 Bool_t posDir=(s.EndPlane>s.VtxPlane);
02540 Float_t materialFromTrkEnd=0;
02541
02542 const Float_t material01Pl=(5.94)*Munits::cm;
02543
02544 static TProfile* pEnVsDist=0;
02545 static TProfile* pEnVsDist0=0;
02546 static TProfile* pEnVsDist0Cor=0;
02547 static TProfile* pEnVsDistNoEnd=0;
02548
02549
02550
02551 static TH1F* hSigMapWinPln=0;
02552 static TH1F* hSigMapAllPln=0;
02553 static TH1F* hSigCorWinPln=0;
02554 static TH1F* hSigCorAllPln=0;
02555 static TH1F* hSigLinWinPln=0;
02556 static TH1F* hSigLinAllPln=0;
02557 static TH1F* hAdcWinPln=0;
02558 static TH1F* hAdcAllPln=0;
02559 static TH1F* hPeWinPln=0;
02560 static TH1F* hPeAllPln=0;
02561 static TH1F* hDigitWinPln=0;
02562 static TH1F* hDigitAllPln=0;
02563 static TH1F* hStripWinPln=0;
02564 static TH1F* hStripAllPln=0;
02565
02566 static TH1F* hSigMapWinPln1=0;
02567 static TH1F* hSigMapAllPln1=0;
02568 static TH1F* hSigCorWinPln1=0;
02569 static TH1F* hSigCorAllPln1=0;
02570 static TH1F* hSigLinWinPln1=0;
02571 static TH1F* hSigLinAllPln1=0;
02572 static TH1F* hAdcWinPln1=0;
02573 static TH1F* hAdcAllPln1=0;
02574 static TH1F* hPeWinPln1=0;
02575 static TH1F* hPeAllPln1=0;
02576 static TH1F* hDigitWinPln1=0;
02577 static TH1F* hDigitAllPln1=0;
02578 static TH1F* hStripWinPln1=0;
02579 static TH1F* hStripAllPln1=0;
02580
02581 static TH1F* hSigMapWinPln2=0;
02582 static TH1F* hSigMapAllPln2=0;
02583 static TH1F* hSigCorWinPln2=0;
02584 static TH1F* hSigCorAllPln2=0;
02585 static TH1F* hSigLinWinPln2=0;
02586 static TH1F* hSigLinAllPln2=0;
02587 static TH1F* hAdcWinPln2=0;
02588 static TH1F* hAdcAllPln2=0;
02589 static TH1F* hPeWinPln2=0;
02590 static TH1F* hPeAllPln2=0;
02591 static TH1F* hDigitWinPln2=0;
02592 static TH1F* hDigitAllPln2=0;
02593 static TH1F* hStripWinPln2=0;
02594 static TH1F* hStripAllPln2=0;
02595
02596
02597 if (pEnVsDist==0) {
02598 pEnVsDist=new TProfile("pEnVsDist","pEnVsDist",60,0,
02599 (60*material01Pl)/Munits::cm);
02600 pEnVsDistNoEnd=new TProfile("pEnVsDistNoEnd","pEnVsDistNoEnd",60,0,
02601 (60*material01Pl)/Munits::cm);
02602 pEnVsDist0=new TProfile("pEnVsDist0","pEnVsDist0",60,0,
02603 (60*material01Pl)/Munits::cm);
02604 pEnVsDist0Cor=new TProfile("pEnVsDist0Cor","pEnVsDist0Cor",60,0,
02605 (60*material01Pl)/Munits::cm);
02606 MSG("MeuCuts",logLevel)
02607 <<"Creating TProfile, pEnVsDist, max bin="
02608 <<(60*material01Pl)/Munits::cm<<endl;
02609
02610
02611 hSigMapWinPln=new TH1F("hSigMapWinPln","hSigMapWinPln",
02612 20320,-320,20000);
02613 hSigMapWinPln->SetFillColor(0);
02614 hSigMapWinPln->SetTitle("SigMap of Planes in Window (in Slice)");
02615 hSigMapWinPln->GetXaxis()->SetTitle("SigMap");
02616 hSigMapWinPln->GetXaxis()->CenterTitle();
02617
02618
02619 hSigMapAllPln=new TH1F("hSigMapAllPln","hSigMapAllPln",
02620 20320,-320,20000);
02621 hSigMapAllPln->SetFillColor(0);
02622 hSigMapAllPln->SetTitle("SigMap of All Planes in Slice");
02623 hSigMapAllPln->GetXaxis()->SetTitle("SigMap");
02624 hSigMapAllPln->GetXaxis()->CenterTitle();
02625
02626
02627
02628 hSigCorWinPln=new TH1F("hSigCorWinPln","hSigCorWinPln",
02629 20320,-320,20000);
02630 hSigCorWinPln->SetFillColor(0);
02631 hSigCorWinPln->SetTitle("SigCor of Planes in Window (in Slice)");
02632 hSigCorWinPln->GetXaxis()->SetTitle("SigCor");
02633 hSigCorWinPln->GetXaxis()->CenterTitle();
02634
02635
02636 hSigCorAllPln=new TH1F("hSigCorAllPln","hSigCorAllPln",
02637 20320,-320,20000);
02638 hSigCorAllPln->SetFillColor(0);
02639 hSigCorAllPln->SetTitle("SigCor of All Planes in Slice");
02640 hSigCorAllPln->GetXaxis()->SetTitle("SigCor");
02641 hSigCorAllPln->GetXaxis()->CenterTitle();
02642
02643
02644
02645 hSigLinWinPln=new TH1F("hSigLinWinPln","hSigLinWinPln",
02646 20320,-320,20000);
02647 hSigLinWinPln->SetFillColor(0);
02648 hSigLinWinPln->SetTitle("SigLin of Planes in Window (in Slice)");
02649 hSigLinWinPln->GetXaxis()->SetTitle("SigLin");
02650 hSigLinWinPln->GetXaxis()->CenterTitle();
02651
02652
02653 hSigLinAllPln=new TH1F("hSigLinAllPln","hSigLinAllPln",
02654 20320,-320,20000);
02655 hSigLinAllPln->SetFillColor(0);
02656 hSigLinAllPln->SetTitle("SigLin of All Planes in Slice");
02657 hSigLinAllPln->GetXaxis()->SetTitle("SigLin");
02658 hSigLinAllPln->GetXaxis()->CenterTitle();
02659
02660
02661
02662 hAdcWinPln=new TH1F("hAdcWinPln","hAdcWinPln",
02663 20320,-320,20000);
02664 hAdcWinPln->SetFillColor(0);
02665 hAdcWinPln->SetTitle("Adc of Planes in Window (in Slice)");
02666 hAdcWinPln->GetXaxis()->SetTitle("Adc");
02667 hAdcWinPln->GetXaxis()->CenterTitle();
02668
02669
02670 hAdcAllPln=new TH1F("hAdcAllPln","hAdcAllPln",
02671 20320,-320,20000);
02672 hAdcAllPln->SetFillColor(0);
02673 hAdcAllPln->SetTitle("Adc of All Planes in Slice");
02674 hAdcAllPln->GetXaxis()->SetTitle("Adc");
02675 hAdcAllPln->GetXaxis()->CenterTitle();
02676
02677
02678
02679 hPeWinPln=new TH1F("hPeWinPln","hPeWinPln",
02680 10160,-4,250);
02681 hPeWinPln->SetFillColor(0);
02682 hPeWinPln->SetTitle("Pe of Planes in Window (in Slice)");
02683 hPeWinPln->GetXaxis()->SetTitle("Pe");
02684 hPeWinPln->GetXaxis()->CenterTitle();
02685
02686
02687 hPeAllPln=new TH1F("hPeAllPln","hPeAllPln",
02688 10160,-4,250);
02689 hPeAllPln->SetFillColor(0);
02690 hPeAllPln->SetTitle("Pe of All Planes in Slice");
02691 hPeAllPln->GetXaxis()->SetTitle("Pe");
02692 hPeAllPln->GetXaxis()->CenterTitle();
02693
02694
02695
02696 hDigitWinPln=new TH1F("hDigitWinPln","hDigitWinPln",
02697 4020,-2,400);
02698 hDigitWinPln->SetFillColor(0);
02699 hDigitWinPln->SetTitle("Digits of Planes in Window (in Slice)");
02700 hDigitWinPln->GetXaxis()->SetTitle("Digits");
02701 hDigitWinPln->GetXaxis()->CenterTitle();
02702 hDigitWinPln->SetBit(TH1::kCanRebin);
02703
02704 hDigitAllPln=new TH1F("hDigitAllPln","hDigitAllPln",
02705 4020,-2,400);
02706 hDigitAllPln->SetFillColor(0);
02707 hDigitAllPln->SetTitle("Digits of All Planes in Slice");
02708 hDigitAllPln->GetXaxis()->SetTitle("Digits");
02709 hDigitAllPln->GetXaxis()->CenterTitle();
02710
02711
02712 hStripWinPln=new TH1F("hStripWinPln","hStripWinPln",
02713 4020,-2,400);
02714 hStripWinPln->SetFillColor(0);
02715 hStripWinPln->SetTitle("Stp/Pln of Plns in Window");
02716 hStripWinPln->GetXaxis()->SetTitle("Strips");
02717 hStripWinPln->GetXaxis()->CenterTitle();
02718 hStripWinPln->SetBit(TH1::kCanRebin);
02719
02720 hStripAllPln=new TH1F("hStripAllPln","hStripAllPln",
02721 4020,-2,400);
02722 hStripAllPln->SetFillColor(0);
02723 hStripAllPln->SetTitle("Stp/Pln of All Plns in Trk");
02724 hStripAllPln->GetXaxis()->SetTitle("Strips");
02725 hStripAllPln->GetXaxis()->CenterTitle();
02726
02727
02728
02729
02730 hSigMapWinPln1=new TH1F("hSigMapWinPln1","hSigMapWinPln1",
02731 20320,-320,20000);
02732 hSigMapWinPln1->SetFillColor(0);
02733 hSigMapWinPln1->SetTitle("SigMap of Planes in Window (in Slice)");
02734 hSigMapWinPln1->GetXaxis()->SetTitle("SigMap");
02735 hSigMapWinPln1->GetXaxis()->CenterTitle();
02736
02737
02738 hSigMapAllPln1=new TH1F("hSigMapAllPln1","hSigMapAllPln1",
02739 20320,-320,20000);
02740 hSigMapAllPln1->SetFillColor(0);
02741 hSigMapAllPln1->SetTitle("SigMap of All Planes in Slice");
02742 hSigMapAllPln1->GetXaxis()->SetTitle("SigMap");
02743 hSigMapAllPln1->GetXaxis()->CenterTitle();
02744
02745
02746
02747 hSigCorWinPln1=new TH1F("hSigCorWinPln1","hSigCorWinPln1",
02748 20320,-320,20000);
02749 hSigCorWinPln1->SetFillColor(0);
02750 hSigCorWinPln1->SetTitle("SigCor of Planes in Window (in Slice)");
02751 hSigCorWinPln1->GetXaxis()->SetTitle("SigCor");
02752 hSigCorWinPln1->GetXaxis()->CenterTitle();
02753
02754
02755 hSigCorAllPln1=new TH1F("hSigCorAllPln1","hSigCorAllPln1",
02756 20320,-320,20000);
02757 hSigCorAllPln1->SetFillColor(0);
02758 hSigCorAllPln1->SetTitle("SigCor of All Planes in Slice");
02759 hSigCorAllPln1->GetXaxis()->SetTitle("SigCor");
02760 hSigCorAllPln1->GetXaxis()->CenterTitle();
02761
02762
02763
02764 hSigLinWinPln1=new TH1F("hSigLinWinPln1","hSigLinWinPln1",
02765 20320,-320,20000);
02766 hSigLinWinPln1->SetFillColor(0);
02767 hSigLinWinPln1->SetTitle("SigLin of Planes in Window (in Slice)");
02768 hSigLinWinPln1->GetXaxis()->SetTitle("SigLin");
02769 hSigLinWinPln1->GetXaxis()->CenterTitle();
02770
02771
02772 hSigLinAllPln1=new TH1F("hSigLinAllPln1","hSigLinAllPln1",
02773 20320,-320,20000);
02774 hSigLinAllPln1->SetFillColor(0);
02775 hSigLinAllPln1->SetTitle("SigLin of All Planes in Slice");
02776 hSigLinAllPln1->GetXaxis()->SetTitle("SigLin");
02777 hSigLinAllPln1->GetXaxis()->CenterTitle();
02778
02779
02780
02781 hAdcWinPln1=new TH1F("hAdcWinPln1","hAdcWinPln1",
02782 20320,-320,20000);
02783 hAdcWinPln1->SetFillColor(0);
02784 hAdcWinPln1->SetTitle("Adc of Planes in Window (in Slice)");
02785 hAdcWinPln1->GetXaxis()->SetTitle("Adc");
02786 hAdcWinPln1->GetXaxis()->CenterTitle();
02787
02788
02789 hAdcAllPln1=new TH1F("hAdcAllPln1","hAdcAllPln1",
02790 20320,-320,20000);
02791 hAdcAllPln1->SetFillColor(0);
02792 hAdcAllPln1->SetTitle("Adc of All Planes in Slice");
02793 hAdcAllPln1->GetXaxis()->SetTitle("Adc");
02794 hAdcAllPln1->GetXaxis()->CenterTitle();
02795
02796
02797
02798 hPeWinPln1=new TH1F("hPeWinPln1","hPeWinPln1",
02799 10160,-4,250);
02800 hPeWinPln1->SetFillColor(0);
02801 hPeWinPln1->SetTitle("Pe of Planes in Window (in Slice)");
02802 hPeWinPln1->GetXaxis()->SetTitle("Pe");
02803 hPeWinPln1->GetXaxis()->CenterTitle();
02804
02805
02806 hPeAllPln1=new TH1F("hPeAllPln1","hPeAllPln1",
02807 10160,-4,250);
02808 hPeAllPln1->SetFillColor(0);
02809 hPeAllPln1->SetTitle("Pe of All Planes in Slice");
02810 hPeAllPln1->GetXaxis()->SetTitle("Pe");
02811 hPeAllPln1->GetXaxis()->CenterTitle();
02812
02813
02814
02815 hDigitWinPln1=new TH1F("hDigitWinPln1","hDigitWinPln1",
02816 4020,-2,400);
02817 hDigitWinPln1->SetFillColor(0);
02818 hDigitWinPln1->SetTitle("Digits of Planes in Window (in Slice)");
02819 hDigitWinPln1->GetXaxis()->SetTitle("Digits");
02820 hDigitWinPln1->GetXaxis()->CenterTitle();
02821
02822
02823 hDigitAllPln1=new TH1F("hDigitAllPln1","hDigitAllPln1",
02824 4020,-2,400);
02825 hDigitAllPln1->SetFillColor(0);
02826 hDigitAllPln1->SetTitle("Digits of All Planes in Slice");
02827 hDigitAllPln1->GetXaxis()->SetTitle("Digits");
02828 hDigitAllPln1->GetXaxis()->CenterTitle();
02829
02830
02831 hStripWinPln1=new TH1F("hStripWinPln1","hStripWinPln1",
02832 4020,-2,400);
02833 hStripWinPln1->SetFillColor(0);
02834 hStripWinPln1->SetTitle("Stp/Pln of Plns in Window");
02835 hStripWinPln1->GetXaxis()->SetTitle("Strips");
02836 hStripWinPln1->GetXaxis()->CenterTitle();
02837
02838
02839 hStripAllPln1=new TH1F("hStripAllPln1","hStripAllPln1",
02840 4020,-2,400);
02841 hStripAllPln1->SetFillColor(0);
02842 hStripAllPln1->SetTitle("Stp/Pln of All Plns in Trk");
02843 hStripAllPln1->GetXaxis()->SetTitle("Strips");
02844 hStripAllPln1->GetXaxis()->CenterTitle();
02845
02846
02847
02848
02849 hSigMapWinPln2=new TH1F("hSigMapWinPln2","hSigMapWinPln2",
02850 20320,-320,20000);
02851 hSigMapWinPln2->SetFillColor(0);
02852 hSigMapWinPln2->SetTitle("SigMap of Planes in Window (in Slice)");
02853 hSigMapWinPln2->GetXaxis()->SetTitle("SigMap");
02854 hSigMapWinPln2->GetXaxis()->CenterTitle();
02855
02856
02857 hSigMapAllPln2=new TH1F("hSigMapAllPln2","hSigMapAllPln2",
02858 20320,-320,20000);
02859 hSigMapAllPln2->SetFillColor(0);
02860 hSigMapAllPln2->SetTitle("SigMap of All Planes in Slice");
02861 hSigMapAllPln2->GetXaxis()->SetTitle("SigMap");
02862 hSigMapAllPln2->GetXaxis()->CenterTitle();
02863
02864
02865
02866 hSigCorWinPln2=new TH1F("hSigCorWinPln2","hSigCorWinPln2",
02867 20320,-320,20000);
02868 hSigCorWinPln2->SetFillColor(0);
02869 hSigCorWinPln2->SetTitle("SigCor of Planes in Window (in Slice)");
02870 hSigCorWinPln2->GetXaxis()->SetTitle("SigCor");
02871 hSigCorWinPln2->GetXaxis()->CenterTitle();
02872
02873
02874 hSigCorAllPln2=new TH1F("hSigCorAllPln2","hSigCorAllPln2",
02875 20320,-320,20000);
02876 hSigCorAllPln2->SetFillColor(0);
02877 hSigCorAllPln2->SetTitle("SigCor of All Planes in Slice");
02878 hSigCorAllPln2->GetXaxis()->SetTitle("SigCor");
02879 hSigCorAllPln2->GetXaxis()->CenterTitle();
02880
02881
02882
02883 hSigLinWinPln2=new TH1F("hSigLinWinPln2","hSigLinWinPln2",
02884 20320,-320,20000);
02885 hSigLinWinPln2->SetFillColor(0);
02886 hSigLinWinPln2->SetTitle("SigLin of Planes in Window (in Slice)");
02887 hSigLinWinPln2->GetXaxis()->SetTitle("SigLin");
02888 hSigLinWinPln2->GetXaxis()->CenterTitle();
02889
02890
02891 hSigLinAllPln2=new TH1F("hSigLinAllPln2","hSigLinAllPln2",
02892 20320,-320,20000);
02893 hSigLinAllPln2->SetFillColor(0);
02894 hSigLinAllPln2->SetTitle("SigLin of All Planes in Slice");
02895 hSigLinAllPln2->GetXaxis()->SetTitle("SigLin");
02896 hSigLinAllPln2->GetXaxis()->CenterTitle();
02897
02898
02899
02900 hAdcWinPln2=new TH1F("hAdcWinPln2","hAdcWinPln2",
02901 20320,-320,20000);
02902 hAdcWinPln2->SetFillColor(0);
02903 hAdcWinPln2->SetTitle("Adc of Planes in Window (in Slice)");
02904 hAdcWinPln2->GetXaxis()->SetTitle("Adc");
02905 hAdcWinPln2->GetXaxis()->CenterTitle();
02906
02907
02908 hAdcAllPln2=new TH1F("hAdcAllPln2","hAdcAllPln2",
02909 20320,-320,20000);
02910 hAdcAllPln2->SetFillColor(0);
02911 hAdcAllPln2->SetTitle("Adc of All Planes in Slice");
02912 hAdcAllPln2->GetXaxis()->SetTitle("Adc");
02913 hAdcAllPln2->GetXaxis()->CenterTitle();
02914
02915
02916
02917 hPeWinPln2=new TH1F("hPeWinPln2","hPeWinPln2",
02918 10160,-4,250);
02919 hPeWinPln2->SetFillColor(0);
02920 hPeWinPln2->SetTitle("Pe of Planes in Window (in Slice)");
02921 hPeWinPln2->GetXaxis()->SetTitle("Pe");
02922 hPeWinPln2->GetXaxis()->CenterTitle();
02923
02924
02925 hPeAllPln2=new TH1F("hPeAllPln2","hPeAllPln2",
02926 10160,-4,250);
02927 hPeAllPln2->SetFillColor(0);
02928 hPeAllPln2->SetTitle("Pe of All Planes in Slice");
02929 hPeAllPln2->GetXaxis()->SetTitle("Pe");
02930 hPeAllPln2->GetXaxis()->CenterTitle();
02931
02932
02933
02934 hDigitWinPln2=new TH1F("hDigitWinPln2","hDigitWinPln2",
02935 4020,-2,400);
02936 hDigitWinPln2->SetFillColor(0);
02937 hDigitWinPln2->SetTitle("Digits of Planes in Window (in Slice)");
02938 hDigitWinPln2->GetXaxis()->SetTitle("Digits");
02939 hDigitWinPln2->GetXaxis()->CenterTitle();
02940
02941
02942 hDigitAllPln2=new TH1F("hDigitAllPln2","hDigitAllPln2",
02943 4020,-2,400);
02944 hDigitAllPln2->SetFillColor(0);
02945 hDigitAllPln2->SetTitle("Digits of All Planes in Slice");
02946 hDigitAllPln2->GetXaxis()->SetTitle("Digits");
02947 hDigitAllPln2->GetXaxis()->CenterTitle();
02948
02949
02950 hStripWinPln2=new TH1F("hStripWinPln2","hStripWinPln2",
02951 4020,-2,400);
02952 hStripWinPln2->SetFillColor(0);
02953 hStripWinPln2->SetTitle("Stp/Pln of Plns in Window");
02954 hStripWinPln2->GetXaxis()->SetTitle("Strips");
02955 hStripWinPln2->GetXaxis()->CenterTitle();
02956
02957
02958 hStripAllPln2=new TH1F("hStripAllPln2","hStripAllPln2",
02959 4020,-2,400);
02960 hStripAllPln2->SetFillColor(0);
02961 hStripAllPln2->SetTitle("Stp/Pln of All Plns in Trk");
02962 hStripAllPln2->GetXaxis()->SetTitle("Strips");
02963 hStripAllPln2->GetXaxis()->CenterTitle();
02964
02965 }
02966
02967
02968 if (s.EndPlane>=0){
02969
02970
02971
02972 Int_t pl=s.EndPlane;
02973 Int_t tmpVtxPlane=s.VtxPlane;
02974 if (posDir) tmpVtxPlane--;
02975 else tmpVtxPlane++;
02976
02977
02978 while (pl!=tmpVtxPlane){
02979
02980
02981 const MeuHitInfo& cp=plInfo.find(pl)->second;
02982 Float_t plMaterial=cp.PLCor*material01Pl;
02983
02984
02985 Float_t correctionFactor=0.0134*cp.Y+1;
02986 if (cp.Y<1) correctionFactor=0.005*cp.Y+1;
02987 Float_t meuSigMapCor=(cp.SigMap*correctionFactor)/cp.PLCor;
02988 Float_t meuSigMap=cp.SigMap/cp.PLCor;
02989 MAXMSG("MeuCuts",logLevel,100)
02990 <<"pl="<<pl<<", y="<<cp.Y<<", corFact="<<correctionFactor
02991 <<", meuSigMap="<<meuSigMap
02992 <<", meuSigMapCor="<<meuSigMapCor<<endl;
02993
02994
02995
02996 Float_t toPlotMaterial=plMaterial*0.5+materialFromTrkEnd;
02997
02998
02999
03000 if (s.WinSigMap>0 && s.WinSigMap<2000 &&
03001 fabs(s.WinAvCosThetaZ)>0.5 &&
03002 s.TotalMatTraversed>3 &&
03003 cp.SigMap/cp.PLCor<3000){
03004
03005 pEnVsDist->Fill((toPlotMaterial/Munits::cm),meuSigMap);
03006
03007
03008
03009 pEnVsDist0->Fill((materialFromTrkEnd/Munits::cm),meuSigMap);
03010 pEnVsDist0Cor->Fill((materialFromTrkEnd/Munits::cm),
03011 meuSigMapCor);
03012
03013 if (pl!=s.EndPlane) {
03014 pEnVsDistNoEnd->Fill((toPlotMaterial/Munits::cm),meuSigMap);
03015 }
03016 }
03017 materialFromTrkEnd+=plMaterial;
03018
03019 MAXMSG("MeuCuts",logLevel,1000)
03020 <<"plMat="<<plMaterial
03021 <<", matTrk="<<materialFromTrkEnd
03022 <<", toPlotMat="<<(toPlotMaterial/Munits::cm)
03023 <<"meuSigMap="<<meuSigMap<<", PLCor="<<cp.PLCor<<endl;
03024
03025
03027
03028 hSigMapAllPln->Fill(cp.SigMap);
03029 hSigCorAllPln->Fill(cp.SigCor);
03030 hSigLinAllPln->Fill(cp.SigLin);
03031 hAdcAllPln->Fill(cp.Adc);
03032 hPeAllPln->Fill(cp.Pe);
03033 hDigitAllPln->Fill(cp.NumDigits);
03034 hStripAllPln->Fill(cp.NumStrips);
03035
03036 hSigMapAllPln1->Fill(cp.SigMap1);
03037 hSigCorAllPln1->Fill(cp.SigCor1);
03038 hSigLinAllPln1->Fill(cp.SigLin1);
03039 hAdcAllPln1->Fill(cp.Adc1);
03040 hPeAllPln1->Fill(cp.Pe1);
03041
03042
03043
03044 hSigMapAllPln2->Fill(cp.SigMap2);
03045 hSigCorAllPln2->Fill(cp.SigCor2);
03046 hSigLinAllPln2->Fill(cp.SigLin2);
03047 hAdcAllPln2->Fill(cp.Adc2);
03048 hPeAllPln2->Fill(cp.Pe2);
03049
03050
03051
03052
03053
03054 if ((cp.Plane>=s.WinVtxSidePl && cp.Plane<=s.WinStopSidePl) ||
03055 (cp.Plane<=s.WinVtxSidePl && cp.Plane>=s.WinStopSidePl)) {
03056 hSigMapWinPln->Fill(cp.SigMap);
03057 hSigCorWinPln->Fill(cp.SigCor);
03058 hSigLinWinPln->Fill(cp.SigLin);
03059 hAdcWinPln->Fill(cp.Adc);
03060 hPeWinPln->Fill(cp.Pe);
03061 hDigitWinPln->Fill(cp.NumDigits);
03062 hStripWinPln->Fill(cp.NumStrips);
03063
03064 hSigMapWinPln1->Fill(cp.SigMap1);
03065 hSigCorWinPln1->Fill(cp.SigCor1);
03066 hSigLinWinPln1->Fill(cp.SigLin1);
03067 hAdcWinPln1->Fill(cp.Adc1);
03068 hPeWinPln1->Fill(cp.Pe1);
03069
03070
03071
03072 hSigMapWinPln2->Fill(cp.SigMap2);
03073 hSigCorWinPln2->Fill(cp.SigCor2);
03074 hSigLinWinPln2->Fill(cp.SigLin2);
03075 hAdcWinPln2->Fill(cp.Adc2);
03076 hPeWinPln2->Fill(cp.Pe2);
03077
03078
03079 }
03080
03081
03082 if (posDir) pl--;
03083 else pl++;
03084 }
03085 }
03086 }
03087
03088
03089
03090 void MeuCuts::FillEnVsDist2
03091 (const MeuSummary& s,const std::map<Int_t,MeuHitInfo>& plInfo) const
03092 {
03093
03094 Msg::LogLevel_t logLevel=Msg::kDebug;
03095
03096 Bool_t posDir=(s.EndPlane>s.VtxPlane);
03097 Float_t materialFromTrkEnd=0;
03098
03099 const Float_t material01Pl=(5.94)*Munits::cm;
03100
03101 if (s.TotalMatTraversed<(60*material01Pl) ||
03102 fabs(s.WinAvCosThetaZ)<0.3) {
03103 MAXMSG("MeuCuts",Msg::kDebug,10)
03104 <<"MatTrav="<<s.TotalMatTraversed/Munits::cm<<" cm, 60 planes="
03105 <<(60*material01Pl)/Munits::cm<<" cm"
03106 <<", cosThZ="<<s.WinAvCosThetaZ<<endl;
03107 return;
03108 }
03109
03111
03112 Float_t containDisttoEdge=0.1;
03113 Int_t pl2=s.EndPlane;
03114 Int_t tmpVtxPlane2=s.VtxPlane;
03115 if (posDir) tmpVtxPlane2--;
03116 else tmpVtxPlane2++;
03117
03118 Float_t matTrav=0;
03119 Bool_t allDeep=true;
03120
03121
03122 while (pl2!=tmpVtxPlane2){
03123
03124
03125 const MeuHitInfo& cp=plInfo.find(pl2)->second;
03126 Float_t plMaterial=cp.PLCor*material01Pl;
03127 matTrav+=plMaterial;
03128
03129 static MeuReco reco;
03130 Float_t distToEdge=reco.CheckDeepDistToEdge(cp);
03131
03132
03133 if (matTrav<(60*material01Pl)){
03134 if(distToEdge<containDisttoEdge) {
03135 allDeep=false;
03136 MAXMSG("MeuCuts",Msg::kDebug,500)
03137 <<"NOT DEEPLY CONTAINED: pl="<<cp.Plane<<", st="<<cp.Strip
03138 <<", d="<<distToEdge<<", bool="<<allDeep<<endl;
03139 return;
03140 }
03141 }
03142
03143 MAXMSG("MeuCuts",Msg::kDebug,500)
03144 <<"pl="<<pl2<<", matTrav="<<matTrav<<", plMat="<<plMaterial
03145 <<", d="<<distToEdge<<", bool="<<allDeep<<endl;
03146
03147
03148 if (posDir) pl2--;
03149 else pl2++;
03150 }
03151
03153 static TProfile* pEnVsDist2=0;
03154 static TProfile* pEnVsDist20=0;
03155 static TProfile* pEnVsDist20Cor=0;
03156 static TProfile* pEnVsDist2NoEnd=0;
03157 static TH1F* hEnVsDistMeu=0;
03158 if (pEnVsDist2==0) {
03159 pEnVsDist2=new TProfile("pEnVsDist2","pEnVsDist2",60,0,
03160 (60*material01Pl)/Munits::cm);
03161 pEnVsDist2NoEnd=new TProfile("pEnVsDist2NoEnd","pEnVsDist2NoEnd",60,0,
03162 (60*material01Pl)/Munits::cm);
03163 pEnVsDist20=new TProfile("pEnVsDist20","pEnVsDist20",60,0,
03164 (60*material01Pl)/Munits::cm);
03165 pEnVsDist20Cor=new TProfile("pEnVsDist20Cor","pEnVsDist20Cor",60,0,
03166 (60*material01Pl)/Munits::cm);
03167 hEnVsDistMeu=new TH1F("hEnVsDistMeu","hEnVsDistMeu",3001,-1,3000);
03168 hEnVsDistMeu->SetFillColor(0);
03169 hEnVsDistMeu->SetTitle("MEU in SigMap Units");
03170 hEnVsDistMeu->GetXaxis()->SetTitle("MEU");
03171 hEnVsDistMeu->GetXaxis()->CenterTitle();
03172
03173 MSG("MeuCuts",logLevel)
03174 <<"Creating TProfile, pEnVsDist2, max bin="
03175 <<(60*material01Pl)/Munits::cm<<endl;
03176 }
03177
03178
03179 if (s.EndPlane>=0){
03180
03181 hEnVsDistMeu->Fill(s.WinSigMap);
03182
03183
03184
03185 Int_t pl=s.EndPlane;
03186 Int_t tmpVtxPlane=s.VtxPlane;
03187 if (posDir) tmpVtxPlane--;
03188 else tmpVtxPlane++;
03189
03190
03191 while (pl!=tmpVtxPlane){
03192
03193
03194 const MeuHitInfo& cp=plInfo.find(pl)->second;
03195 Float_t plMaterial=cp.PLCor*material01Pl;
03196
03197
03198 Float_t correctionFactor=0.0134*cp.Y+1;
03199 if (cp.Y<1) correctionFactor=0.005*cp.Y+1;
03200 Float_t meuSigMapCor=(cp.SigMap*correctionFactor)/cp.PLCor;
03201 Float_t meuSigMap=cp.SigMap/cp.PLCor;
03202 MAXMSG("MeuCuts",logLevel,100)
03203 <<"pl="<<pl<<", y="<<cp.Y<<", corFact="<<correctionFactor
03204 <<", meuSigMap="<<meuSigMap
03205 <<", meuSigMapCor="<<meuSigMapCor<<endl;
03206
03207
03208
03209 Float_t toPlotMaterial=plMaterial*0.5+materialFromTrkEnd;
03210
03211
03212
03213 if (s.WinSigMap>0 && s.WinSigMap<2000 &&
03214 fabs(s.WinAvCosThetaZ)>0.5 &&
03215 s.TotalMatTraversed>3 &&
03216 cp.SigMap/cp.PLCor<3000){
03217
03218 pEnVsDist2->Fill((toPlotMaterial/Munits::cm),meuSigMap);
03219
03220
03221
03222 pEnVsDist20->Fill((materialFromTrkEnd/Munits::cm),meuSigMap);
03223 pEnVsDist20Cor->Fill((materialFromTrkEnd/Munits::cm),
03224 meuSigMapCor);
03225
03226 if (pl!=s.EndPlane) {
03227 pEnVsDist2NoEnd->Fill((toPlotMaterial/Munits::cm),meuSigMap);
03228 }
03229 }
03230 materialFromTrkEnd+=plMaterial;
03231
03232 MAXMSG("MeuCuts",logLevel,1000)
03233 <<"plMat="<<plMaterial
03234 <<", matTrk="<<materialFromTrkEnd
03235 <<", toPlotMat="<<(toPlotMaterial/Munits::cm)
03236 <<"meuSigMap="<<meuSigMap<<", PLCor="<<cp.PLCor<<endl;
03237
03238
03239 if (posDir) pl--;
03240 else pl++;
03241 }
03242 }
03243 }
03244
03245
03246
03247 Bool_t MeuCuts::IsGoodDataQuality(const NtpStRecord& ntp) const
03248 {
03249 Bool_t pass = true;
03250
03251 if (ntp.GetHeader().GetVldContext().GetDetector() != Detector::kFar)
03252 return pass;
03253
03254 NtpSRDataQuality dataqua = ntp.dataquality;
03255
03256
03257 pass = pass && dataqua.cratemask == 16;
03258
03259
03260 pass = pass && dataqua.trigsource < 130000;
03261
03262
03263 pass = pass && !( dataqua.litrigger == 1
03264 && dataqua.lisubtractedtime > -40000.0
03265 && dataqua.lisubtractedtime < +1000.0 );
03266
03267
03268 pass = pass && (dataqua.busychips + dataqua.coldchips) < 20;
03269
03270 return pass;
03271
03272 }
03273
03274
