loon: MeuCuts Class Reference

00052 {
00053   MSG("MeuCuts",Msg::kDebug)
00054     <<"Running MeuCuts Constructor..."<<endl;
00055 
00056 
00057   MSG("MeuCuts",Msg::kDebug)
00058     <<"Finished MeuCuts Constructor"<<endl;
00059 }

00064 {
00065   MSG("MeuCuts",Msg::kDebug)
00066     <<"Running MeuCuts Destructor..."<<endl;
00067   
00068 
00069   MSG("MeuCuts",Msg::kDebug)
00070     <<"Finished MeuCuts Destructor"<<endl;
00071 }

01497 {
01498   Bool_t awayFromCoil=false;  
01499   std::string sLogLevel="kDebug";
01500   //set up the log level (can't use msg in .h file???)
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   //loop over the strips in the snarl
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     //cut out the non-tracked strips
01526     if (!st.Trk) continue;
01527     
01528     //calc the x,y,z positions
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       //this->PrintSnarlInfo(ev,sLogLevel);
01539       MSG("MeuCuts",logLevel)
01540         <<"Already filled en dep for this plane="<<st.Plane<<endl;
01541     }
01542 
01543     //add the positions to the map as long as it hasn't already 
01544     //been done
01545     if (plPosition.count(st.Plane)) {
01546       //this->PrintSnarlInfo(ev,sLogLevel);
01547       MSG("MeuCuts",logLevel)<<"already filled plane"<<endl;
01548     }
01549     else  plPosition[st.Plane]=TVector3(x,y,z);
01550     
01551     //check if radius is smaller than current smallest
01552     Float_t radius=sqrt(pow(x,2)+pow(y,2));
01553     if (radius<smallestRadius) smallestRadius=radius;
01554 
01555   }//end of loop over strips
01556 
01557   //analyse coil proximity
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   //set the output variable
01567   if (r) *r=smallestRadius;
01568 
01569   //print out info below for these cases
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     //calc direction of muon
01577     Bool_t tmpPosDir=true;
01578     if (s.VtxPlane-s.EndPlane>0) tmpPosDir=false;
01579     const Bool_t posDir=tmpPosDir;
01580 
01582     //loop over the track, starting at the vtx (not the end)
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       //look at position if it exists
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         //analyse coil proximity
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       //increment the plane
01612       if (posDir) pl++;
01613       else pl--;
01614     }
01615   }
01616   
01617   return awayFromCoil;
01618 }

01468 {
01469   VldTimeStamp vldts;
01470   //VldTimeStamp vldts(ev.UnixTime,0);
01471   //always kFar here because AtmosEvent
01472   VldContext vc(Detector::kFar,SimFlag::kData,vldts);
01473   //get the ugh
01474   UgliGeomHandle ugh(vc);
01475   TVector3 uvz;
01476   uvz.SetZ(st.Z);
01477   if ((st.View+2)==PlaneView::kU){
01478     uvz.SetX(st.T);//tpos measures u
01479     uvz.SetY(st.L);//lpos measures v
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 }

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     //a gap of 1 is expected
00119     //a gap of 3 means that 1 plane is missing at the end of 1 view
00120     //a gap of 5 means that 2 planes are missing at the end of 1 view
00121     //a gap of 7 means that 3 planes are missing at the end of 1 view
00122     //a gap of 9 means that 4 planes are missing at the end of 1 view
00123       
00124     //allow for 2 planes to be missing at the 
00125     //end of one view
00126       
00127     //in the near detector the full planes mean that 4 planes
00128     //are missing between views even with 100% efficiency
00129 
00130     if (s.Detector==Detector::kFar){
00131       if (diffAtLow>=7 || diffAtHigh>=7){//equivalent to >5
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     //write out the info if requested
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 }

01732 {
01733   //variable to turn on all the useful messages if required
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   //get variables for loops
01743   Int_t pl=s.VtxPlane;
01744   Int_t tmpEndPlane=s.EndPlane;
01745   if (posDir) tmpEndPlane++;//make sure last plane is used
01746   else tmpEndPlane--;
01747   
01748   VldTimeStamp vldts;
01749   //VldTimeStamp vldts(ev.UnixTime,0);
01750   VldContext vc(static_cast<Detector::Detector_t>(s.Detector),
01751                 SimFlag::kData,vldts);
01752   //get the ugh
01753   UgliGeomHandle ugh(vc);
01754   TVector3 uvz;
01755 
01756   //loop and calc xyz using tpos, lpos and z in each plane
01757   //while (pl!=tmpEndPlane){
01758   typedef map<Int_t,MeuHitInfo>::iterator plIt;
01759   for (plIt it=plInfo.begin();it!=plInfo.end();++it){
01760     
01761     //cache the current cp
01762     MeuHitInfo& cp=it->second;
01763 
01764     if (cp.LPos==-999) continue;
01765     
01766     //set z first
01767     uvz.SetZ(cp.Z);
01768     
01769     //now decide what tpos and lpos measure depending on the view
01770     if (cp.View==PlaneView::kU){
01771       uvz.SetX(cp.TPos);//in u-view tpos measures u
01772       uvz.SetY(cp.LPos);//in u-view lpos measures v
01773     }
01774     else if (cp.View==PlaneView::kV){
01775       uvz.SetY(cp.TPos);//in v-view tpos measures v
01776       uvz.SetX(cp.LPos);//in v-view lpos measures u
01777     }
01778     else cout<<"ahhhhhhh no plane view"<<endl;
01779 
01780     //now calc xyz
01781     TVector3 xyz=ugh.uvz2xyz(uvz);
01782     
01783     //now set the cp x and y
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     //now zero the LPos value since later algorithms use it
01792     //to look for gaps. This is bad coding on my part!
01793     cp.LPos=initValue;
01794     MAXMSG("MeuCuts",Msg::kInfo,1)
01795       <<"Doing nasty hack for AtmosEvent X and Y pos"<<endl;
01796   }
01797 }

01803 {
01804   //variable to turn on all the useful messages if required
01805   //Msg::LogLevel_t logLevel=Msg::kDebug;
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   //loop over the strips in the snarl
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];//get a reference to save looking up later
01822 
01823     //set the position variables that are plane based
01824     cp.Plane=pl;
01825     cp.View=(st.View+2);
01826     cp.Z=st.Z;
01827 
01828     //also set the lpos so you can get at x and y
01829     cp.LPos=st.L;
01830     
01831     //sum the energy deposition
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]);//not right
01839     cp.SigLin1+=st.Qadc[0];//not right
01840     cp.SigLin2+=st.Qadc[1];//not right
01841     cp.SigCor+=(st.QPEcorr[0]+st.QPEcorr[1]);
01842     cp.SigCor1+=st.QPEcorr[0];
01843     cp.SigCor2+=st.QPEcorr[1];
01844 
01845     //sum the number of digits/buckets
01846     cp.NumDigits+=st.Ndigits;
01847     cp.NumStrips++;
01848 
01849     //cut out the non-tracked strips
01850     if (!st.Trk) continue;
01851 
01852     //get the min and max planes for each view
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     //record the best strip found in each plane
01864     MeuHitInfo& tmpcp=plSigCorLast[pl];
01865     //check if already found a tracked strip on this plane
01866     if (cp.Strip<0){//no strip found yet
01867       //set the variables
01868       cp.Strip=st.Strip;
01869       cp.TPos=st.T;
01870 
01871       //store the sigcor of both ends of the tracked strip
01872       //this means that fSigCor!=fSigCor1+fSigCor2 always
01873       cp.SigCorTrk1=st.QPEcorr[0];
01874       cp.SigCorTrk2=st.QPEcorr[1];
01875 
01876       //store as current best
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{//already have a strip
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       //decide whether to use current strip or best previous strip
01896       if ((st.QPEcorr[0]+st.QPEcorr[1])>tmpcp.SigCor){
01897         cp.Strip=st.Strip;
01898         cp.TPos=st.T;
01899 
01900         //store the sigcor of both ends of the tracked strip
01901         //this means that fSigCor!=fSigCor1+fSigCor2 always
01902         cp.SigCorTrk1=st.QPEcorr[0];
01903         cp.SigCorTrk2=st.QPEcorr[1];
01904 
01905         //store this strip as the best current strip
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{ //already using best previous strip
01915         MAXMSG("MeuCuts",Msg::kDebug,1000)
01916           <<"...using previous strip"<<endl;
01917       }
01918     }
01919   }//end of loop over strips
01920 
01921   //now calc the xy info
01922   this->ConvertToXY(s,plInfo);
01923 }

01930 {
01931   //variable to turn on all the useful messages if required
01932   Msg::LogLevel_t logLevel=Msg::kDebug;
01933 
01934   //calibration test histos
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     //stripend 1
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     //hAdcPerSigLinNtp1->SetBit(TH1::kCanRebin);
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     //hSigLinPerSigCorNtp1->SetBit(TH1::kCanRebin);
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     //hSigCorPerSigMapNtp1->SetBit(TH1::kCanRebin);
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     //hSigMapPerMipNtp1->SetBit(TH1::kCanRebin);
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     //hMipPerGeVNtp1->SetBit(TH1::kCanRebin);
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     //hAdcPerPeNtp1->SetBit(TH1::kCanRebin);
02001 
02002 
02003     //stripend 2
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     //hAdcPerSigLinNtp2->SetBit(TH1::kCanRebin);
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     //hSigLinPerSigCorNtp2->SetBit(TH1::kCanRebin);
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     //hSigCorPerSigMapNtp2->SetBit(TH1::kCanRebin);
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     //hSigMapPerMipNtp2->SetBit(TH1::kCanRebin);
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     //hMipPerGeVNtp2->SetBit(TH1::kCanRebin);
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     //hAdcPerPeNtp2->SetBit(TH1::kCanRebin);
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   //Int_t numStps=stpTca.GetEntries();
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   //get the slice associated with this event
02068   const NtpSRSlice* pslc=
02069     dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]);
02070   const NtpSRSlice& slc=(*pslc);
02071 
02072   //loop over strips in slc
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     //get the time0, if between 0 and 1 switch to time1
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;//messes up beg/end planes (if not filtered bad trks)!
02091     }
02092     
02093     Int_t pl=stp.plane;
02094     MeuHitInfo& cp=plInfo[pl];//get a reference
02095       
02096     //set the position variables that are plane based
02097     cp.Plane=pl;
02098     cp.View=stp.planeview;
02099     //cp.Z=stp.Z;
02100       
02101     //sum the energy deposition
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     //fill calibration test histos
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     //fill calibration test histos
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     //sum the number of digits/buckets
02131     cp.NumDigits+=stp.ndigit;
02132     //sum the number of strips
02133     cp.NumStrips++;
02134   }
02135   
02136   TClonesArray& trkTca=(*ntp.trk);
02137   //Int_t numTrks=tcaTk.GetEntries();
02138   Int_t numTrks=evt.ntrack;
02139   
02140   //Loop over tracks
02141   //for (Int_t itrk=0;itrk<numTrks;itrk++){
02142   //const NtpSRTrack* ptrk=
02143   //  dynamic_cast<NtpSRTrack*>(tcaTk[itrk]);
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     //this loop is just over the tracked strips
02153     for (Int_t i=0;i<trk.nstrip;i++){
02154 
02155       //check for bug where strip index is -1
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       //MAXMSG("NDMeuCuts",Msg::kInfo,200)
02166       //  <<"Strip="<<stp.strip<<", pl="<<stp.plane
02167       //  <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
02168 
02169 
02170       //fill calibration test histos
02171       //note these 3 are only filled for tracked strips
02172       //stripend 1
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       //stripend 2
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       //check times
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       //sanity check that no trk times are outside window
02197       //trks with large dts should have been filtered
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;//messes up beg/end planes (if not already cut)
02205       }
02206       
02207       Int_t pl=stp.plane;
02208         
02209       //get the min and max planes for each view
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       //fill the calibpos
02224       MeuHitInfo& cp=plInfo[pl];
02225 
02226       //record the best strip found in each plane
02227       MeuHitInfo& tmpcp=plSigCorLast[pl];
02228       //check if already found a tracked strip on this plane
02229       if (cp.Strip<0){//no strip found yet
02230         //set the variables
02231         cp.Strip=stp.strip;
02232         cp.TPos=stp.tpos;
02233             
02234         //store the sigcor of both ends of the tracked strip
02235         //this means that fSigCor!=fSigCor1+fSigCor2 always
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         //store as current best
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         //trk.stpph0sigmap[i];
02251         //trk.stpph0mip[i];
02252         //trk.stpph0gev[i];
02253         //trk.stpattn0c0[i];
02254         time0[pl]=stp.time0;
02255         time1[pl]=stp.time1;
02256       }
02257       else{//already have a strip
02258         MAXMSG("MeuCuts",logLevel,100)
02259           <<"Choosing between strips... run="//<<ev.Run
02260           //<<", sn="<<ev.Snarl<<endl
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         //decide whether to use current strip or best previous strip
02271         if (stp.ph0.sigcor+stp.ph1.sigcor>tmpcp.SigCor){
02272           cp.Strip=stp.strip;
02273           cp.TPos=stp.tpos;
02274                 
02275           //store the sigcor of both ends of the tracked strip
02276           //this means that fSigCor!=fSigCor1+fSigCor2 always
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           //store this strip as the best current strip
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{ //already using best previous strip
02297           MAXMSG("MeuCuts",logLevel,100)
02298             <<"...using previous strip="<<cp.Strip<<endl;
02299         }
02300       }
02301     }
02302   }
02303 }

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   //initialise
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;// (m/GeV)
02376   
02377   //have a vector for efficiency (expensive to search plInfo for dshs)
02378   vector<Float_t> vEn(1000);
02379   vector<Float_t> vEnSupp(1000);
02380 
02381   //loop over scint hits
02382   TClonesArray& scintHits=(*ev.ScintHitList);
02383   Int_t numDshs=scintHits.GetEntries();
02384   //loop over the truth hits
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     //get the highest and lowest muon energy and associated planes
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;//complete saturation
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     //reset to -1 if not found to be consistent with highEn
02417     if (lowEn==999999) lowEn=-1;
02418     
02419     //set all the variables
02420     s.MCHighEn=highEn;
02421     s.MCLowEn=lowEn;
02422     s.MCParticleId=mainParticle;
02423     s.MCVtxPlane=highEnPl;
02424     s.MCEndPlane=lowEnPl;
02425     
02426     //now transfer the cache of variables to the map
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 }

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   //initialise
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;// (m/GeV)
02461   
02462   //have a vector for efficiency (expensive to search plInfo for dshs)
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   //loop over dshs
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     //get the highest and lowest muon energy and associated planes
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;//complete saturation
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     //reset to -1 if not found to be consistent with highEn
02507     if (lowEn==999999) lowEn=-1;
02508     
02509     //set all the variables
02510     s.MCHighEn=highEn;
02511     s.MCLowEn=lowEn;
02512     s.MCParticleId=mainParticle;
02513     s.MCVtxPlane=highEnPl;
02514     s.MCEndPlane=lowEnPl;
02515     
02516     //now transfer the cache of variables to the map
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 }

00575 {
00576   s.BDtortgt=bd.tortgt;
00577   s.BDtor101=bd.tor101;
00578 }

02535 {
02536   //variable to turn on all the useful messages if required
02537   Msg::LogLevel_t logLevel=Msg::kDebug;
02538 
02539   Bool_t posDir=(s.EndPlane>s.VtxPlane);
02540   Float_t materialFromTrkEnd=0;
02541   //const Float_t material01Pl=(2.5+1)*Munits::cm;
02542   const Float_t material01Pl=(5.94)*Munits::cm;//including air gap
02543 
02544   static TProfile* pEnVsDist=0;
02545   static TProfile* pEnVsDist0=0;
02546   static TProfile* pEnVsDist0Cor=0;
02547   static TProfile* pEnVsDistNoEnd=0;
02548 
02549   //these are nothing to do with "time" but don't want to have a 
02550   //separate loop
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     //hSigMapWinPln->SetBit(TH1::kCanRebin);
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     //hSigMapAllPln->SetBit(TH1::kCanRebin);
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     //hSigCorWinPln->SetBit(TH1::kCanRebin);
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     //hSigCorAllPln->SetBit(TH1::kCanRebin);
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     //hSigLinWinPln->SetBit(TH1::kCanRebin);
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     //hSigLinAllPln->SetBit(TH1::kCanRebin);
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     //hAdcWinPln->SetBit(TH1::kCanRebin);
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     //hAdcAllPln->SetBit(TH1::kCanRebin);
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     //hPeWinPln->SetBit(TH1::kCanRebin);
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     //hPeAllPln->SetBit(TH1::kCanRebin);
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     //hDigitAllPln->SetBit(TH1::kCanRebin);
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     //hStripAllPln->SetBit(TH1::kCanRebin);
02727 
02728 
02729     //stripend 1
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     //hSigMapWinPln1->SetBit(TH1::kCanRebin);
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     //hSigMapAllPln1->SetBit(TH1::kCanRebin);
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     //hSigCorWinPln1->SetBit(TH1::kCanRebin);
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     //hSigCorAllPln1->SetBit(TH1::kCanRebin);
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     //hSigLinWinPln1->SetBit(TH1::kCanRebin);
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     //hSigLinAllPln1->SetBit(TH1::kCanRebin);
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     //hAdcWinPln1->SetBit(TH1::kCanRebin);
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     //hAdcAllPln1->SetBit(TH1::kCanRebin);
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     //hPeWinPln1->SetBit(TH1::kCanRebin);
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     //hPeAllPln1->SetBit(TH1::kCanRebin);
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     //hDigitWinPln1->SetBit(TH1::kCanRebin);
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     //hDigitAllPln1->SetBit(TH1::kCanRebin);
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     //hStripWinPln1->SetBit(TH1::kCanRebin);
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     //hStripAllPln1->SetBit(TH1::kCanRebin);
02846 
02847 
02848     //stripend 2
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     //hSigMapWinPln2->SetBit(TH1::kCanRebin);
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     //hSigMapAllPln2->SetBit(TH1::kCanRebin);
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     //hSigCorWinPln2->SetBit(TH1::kCanRebin);
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     //hSigCorAllPln2->SetBit(TH1::kCanRebin);
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     //hSigLinWinPln2->SetBit(TH1::kCanRebin);
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     //hSigLinAllPln2->SetBit(TH1::kCanRebin);
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     //hAdcWinPln2->SetBit(TH1::kCanRebin);
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     //hAdcAllPln2->SetBit(TH1::kCanRebin);
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     //hPeWinPln2->SetBit(TH1::kCanRebin);
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     //hPeAllPln2->SetBit(TH1::kCanRebin);
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     //hDigitWinPln2->SetBit(TH1::kCanRebin);
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     //hDigitAllPln2->SetBit(TH1::kCanRebin);
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     //hStripWinPln2->SetBit(TH1::kCanRebin);
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     //hStripAllPln2->SetBit(TH1::kCanRebin);
02965   }
02966 
02967   //calculate energy deposition in the window
02968   if (s.EndPlane>=0){
02969     
02970     //get variables for loops
02971     //NOTE: start at the end of the track NOT the vtx
02972     Int_t pl=s.EndPlane;
02973     Int_t tmpVtxPlane=s.VtxPlane;
02974     if (posDir) tmpVtxPlane--;//make sure last plane is used
02975     else tmpVtxPlane++;
02976     
02977     //loop over the track starting at the end
02978     while (pl!=tmpVtxPlane){
02979       
02980       //cache the current cp
02981       const MeuHitInfo& cp=plInfo.find(pl)->second;
02982       Float_t plMaterial=cp.PLCor*material01Pl;
02983       //Float_t correctionFactor=0.0052*cp.Y+1;//fit to all y
02984       //the gradient changes with y, so do a simplistic split
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       //plot the value half way through the range of dE/dx it effectively
02995       //integrates over
02996       Float_t toPlotMaterial=plMaterial*0.5+materialFromTrkEnd;
02997 
02998       //put a check on the crazy brems and cut out the really steep ones
02999       //can also get high landau-tail events in individual planes
03000       if (s.WinSigMap>0 && s.WinSigMap<2000 &&
03001           fabs(s.WinAvCosThetaZ)>0.5 &&
03002           s.TotalMatTraversed>3 &&//material traversed more than 3 m
03003           cp.SigMap/cp.PLCor<3000){//3000 is 8 sigma from the peak
03004         
03005         pEnVsDist->Fill((toPlotMaterial/Munits::cm),meuSigMap);
03006         
03007         //fill using 0 for the first pl
03008         //don't try and get an average position
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       //section to make plots of all hit strip distributions
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       //hDigitAllPln1->Fill(cp.NumDigits);
03042       //hStripAllPln1->Fill(cp.NumStrips);
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       //hDigitAllPln2->Fill(cp.NumDigits);
03050       //hStripAllPln2->Fill(cp.NumStrips);
03051 
03052 
03053       //check if plane is in window
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         //hDigitWinPln1->Fill(cp.NumDigits);
03070         //hStripWinPln1->Fill(cp.NumStrips);
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         //hDigitWinPln2->Fill(cp.NumDigits);
03078         //hStripWinPln2->Fill(cp.NumStrips);
03079       }
03080 
03081       //increment the plane
03082       if (posDir) pl--;//working backwards from end to vtx
03083       else pl++;
03084     }
03085   }
03086 }

03092 {
03093   //variable to turn on all the useful messages if required
03094   Msg::LogLevel_t logLevel=Msg::kDebug;
03095 
03096   Bool_t posDir=(s.EndPlane>s.VtxPlane);
03097   Float_t materialFromTrkEnd=0;
03098   //const Float_t material01Pl=(2.5+1)*Munits::cm;
03099   const Float_t material01Pl=(5.94)*Munits::cm;//including air gap
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   //check that all hits on track up to 60 planes are deeply contained
03112   Float_t containDisttoEdge=0.1;//fiducial distance from edge for win
03113   Int_t pl2=s.EndPlane;
03114   Int_t tmpVtxPlane2=s.VtxPlane;
03115   if (posDir) tmpVtxPlane2--;//make sure last plane is used
03116   else tmpVtxPlane2++;
03117 
03118   Float_t matTrav=0;
03119   Bool_t allDeep=true;
03120   
03121   //loop over the track starting at the end
03122   while (pl2!=tmpVtxPlane2){
03123 
03124     //cache the current cp
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     //check that the last 60 planes of track are deeply contained
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     //increment the plane
03148     if (posDir) pl2--;//working backwards from end to vtx
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   //calculate energy deposition in the window
03179   if (s.EndPlane>=0){
03180     
03181     hEnVsDistMeu->Fill(s.WinSigMap);//use to test if meu is good
03182     
03183     //get variables for loops
03184     //NOTE: start at the end of the track NOT the vtx
03185     Int_t pl=s.EndPlane;
03186     Int_t tmpVtxPlane=s.VtxPlane;
03187     if (posDir) tmpVtxPlane--;//make sure last plane is used
03188     else tmpVtxPlane++;
03189     
03190     //loop over the track starting at the end
03191     while (pl!=tmpVtxPlane){
03192       
03193       //cache the current cp
03194       const MeuHitInfo& cp=plInfo.find(pl)->second;
03195       Float_t plMaterial=cp.PLCor*material01Pl;
03196       //Float_t correctionFactor=0.0052*cp.Y+1;//fit to all y
03197       //the gradient changes with y, so do a simplistic split
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       //plot the value half way through the range of dE/dx it effectively
03208       //integrates over
03209       Float_t toPlotMaterial=plMaterial*0.5+materialFromTrkEnd;
03210 
03211       //put a check on the crazy brems and cut out the really steep ones
03212       //can also get high landau-tail events in individual planes
03213       if (s.WinSigMap>0 && s.WinSigMap<2000 &&
03214           fabs(s.WinAvCosThetaZ)>0.5 &&
03215           s.TotalMatTraversed>3 &&//material traversed more than 3 m
03216           cp.SigMap/cp.PLCor<3000){//3000 is 8 sigma from the peak
03217         
03218         pEnVsDist2->Fill((toPlotMaterial/Munits::cm),meuSigMap);
03219         
03220         //fill using 0 for the first pl
03221         //don't try and get an average position
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       //increment the plane
03239       if (posDir) pl--;//working backwards from end to vtx
03240       else pl++;
03241     }
03242   }
03243 }

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   //s.NStrip=tcaStp.GetEntries();
00559   s.NStrip=(*ev.StripList).GetEntries();
00560   s.Detector=Detector::kFar;//always kFar for AtmosEvent
00561   //have to work out the simflag
00562   if (ev.NScintHits<=0) s.SimFlag=SimFlag::kData;
00563   else s.SimFlag=SimFlag::kMC;
00564   s.TrigSrc=ev.TrigSrc;
00565 
00566   //these don't matter for the FD
00567   s.Evt=-1;
00568   s.Slc=-1;
00569 }

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   //s.NStrip=tcaStp.GetEntries();
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   //a hack to correct the trigger bit of the cosmic MC in the ND
00523   if (s.SimFlag==SimFlag::kMC && s.Detector==Detector::kNear) {
00524     //the cosmic MC doesn't have any entries
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;//a plane trigger
00535     }
00536   }
00537 
00538   //apply data quality cuts to FD data, Tingjun 03/27/07
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 }

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 }

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   //these are nothing to do with "time" but don't want to have a 
00651   //separate loop
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     //hSigMapWinStp->SetBit(TH1::kCanRebin);
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     //hSigMapAllStp->SetBit(TH1::kCanRebin);
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     //hSigCorWinStp->SetBit(TH1::kCanRebin);
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     //hSigCorAllStp->SetBit(TH1::kCanRebin);
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     //hSigLinWinStp->SetBit(TH1::kCanRebin);
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     //hSigLinAllStp->SetBit(TH1::kCanRebin);
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     //hAdcWinStp->SetBit(TH1::kCanRebin);
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     //hAdcAllStp->SetBit(TH1::kCanRebin);
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     //hPeWinStp->SetBit(TH1::kCanRebin);
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     //hPeAllStp->SetBit(TH1::kCanRebin);
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     //hDigitAllStp->SetBit(TH1::kCanRebin);
00881 
00882 
00883     //stripend 1
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     //hSigMapWinStp1->SetBit(TH1::kCanRebin);
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     //hSigMapAllStp1->SetBit(TH1::kCanRebin);
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     //hSigCorWinStp1->SetBit(TH1::kCanRebin);
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     //hSigCorAllStp1->SetBit(TH1::kCanRebin);
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     //hSigLinWinStp1->SetBit(TH1::kCanRebin);
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     //hSigLinAllStp1->SetBit(TH1::kCanRebin);
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     //hAdcWinStp1->SetBit(TH1::kCanRebin);
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     //hAdcAllStp1->SetBit(TH1::kCanRebin);
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     //hPeWinStp1->SetBit(TH1::kCanRebin);
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     //hPeAllStp1->SetBit(TH1::kCanRebin);
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     //hDigitWinStp1->SetBit(TH1::kCanRebin);
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     //hDigitAllStp1->SetBit(TH1::kCanRebin);
00984 
00985 
00986     //stripend 2
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     //hSigMapWinStp2->SetBit(TH1::kCanRebin);
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     //hSigMapAllStp2->SetBit(TH1::kCanRebin);
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     //hSigCorWinStp2->SetBit(TH1::kCanRebin);
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     //hSigCorAllStp2->SetBit(TH1::kCanRebin);
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     //hSigLinWinStp2->SetBit(TH1::kCanRebin);
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     //hSigLinAllStp2->SetBit(TH1::kCanRebin);
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     //hAdcWinStp2->SetBit(TH1::kCanRebin);
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     //hAdcAllStp2->SetBit(TH1::kCanRebin);
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     //hPeWinStp2->SetBit(TH1::kCanRebin);
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     //hPeAllStp2->SetBit(TH1::kCanRebin);
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     //hDigitWinStp2->SetBit(TH1::kCanRebin);
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     //hDigitAllStp2->SetBit(TH1::kCanRebin);
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   //this loop is over ALL the strips in the whole SNARL
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     //Int_t pl=stp.plane;
01103 
01104     //get earliest and latest times in snarl
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   //get the median time from the map
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   //these get filled for each selected event in the snarl
01133   //so there is potentially double counting if more than one good 
01134   //stopping muon per snarl
01135   hTimeAfter->Fill((latestTime-medianTime)/Munits::ns);
01136   hTimeBefore->Fill((medianTime-earliestTime)/Munits::ns);
01137   hTimeLength->Fill((latestTime-earliestTime)/Munits::ns);
01138 
01140   //now look at slice times
01142 
01143   //reset variables
01144   earliestTime=999;
01145   latestTime=-999;
01146   times.clear();
01147   
01148   TClonesArray& slcTca=(*ntp.slc);
01149   
01150   //get the slice associated with this event
01151   const NtpSRSlice* pslc=
01152     dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]);
01153   const NtpSRSlice& slc=(*pslc);
01154   
01155   //loop over the strips in the slc that the event was found in
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     //get earliest and latest times in snarl
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     //make a plot of how the energy deposition changes with time cut
01175     //I don't think this actually works...
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     //section to make plots of all hit strip distributions
01205     //hSigMapWinStp->Fill();//no sigmap data
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     //hSigMapAllStp1->Fill();//no sigmap data
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     //hDigitAllStp1->Fill(stp.ndigit);//no data for individual stripend 
01218 
01219     //hSigMapAllStp2->Fill();//no sigmap data
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     //hDigitAllStp2->Fill(stp.ndigit);//no data for individual stripend 
01225 
01226     
01227     //check if plane is in window
01228     if ((stp.plane>=s.WinVtxSidePl && stp.plane<=s.WinStopSidePl) ||
01229         (stp.plane<=s.WinVtxSidePl && stp.plane>=s.WinStopSidePl)) {
01230       //hSigMapWinStp->Fill();//no sigmap data
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       //hSigMapWinStp1->Fill();//no sigmap data
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       //hDigitWinStp1->Fill(stp.ndigit);//no data 4 individual stripend 
01243       
01244       //hSigMapWinStp2->Fill();//no sigmap data
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       //hDigitWinStp2->Fill(stp.ndigit);//no data 4 individual stripend 
01250     }
01251   }
01252 
01253   //get the median time from the map
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   //now look at track times
01265 
01266   //reset variables
01267   earliestTime=999;
01268   latestTime=-999;
01269   times.clear();
01270   
01271   TClonesArray& trkTca=(*ntp.trk);
01272   Int_t numTrks=evt.ntrack;
01273   
01274   //loop over the tracks in the event
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     //this loop is just over the tracked strips
01283     for (Int_t i=0;i<trk.nstrip;i++){
01284 
01285       //check for bug where strip index is -1
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       //get earliest and latest times in snarl
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   //get the median time from the map
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 }

00413 {
00414   Bool_t filter=false;
00415 
00416   //only run for ND at present
00417   //if (s.Detector!=Detector::kNear) return filter;
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   //check that planes exist
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   //loop over planes in window
00437   for (map<Int_t,MeuHitInfo>::const_iterator it=plInfo.find(lowPl);
00438        it!=++plInfo.find(highPl);++it){//stop at 1 past highPl
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   //write out minDist if requested
00458   if (minDistEndStrip) *minDistEndStrip=minDist;
00459   
00460   return filter;
00461 }

00342 {
00343   Bool_t posDir=(s.EndPlane>s.VtxPlane);
00344   Bool_t badGap=false;
00345   Float_t sumAdcInGap=0;
00346 
00347   //static TH1F* hSumTrkEndGapAdc=0;
00348   //if (hSumTrkEndGapAdc==0) {
00349   //  hSumTrkEndGapAdc=new TH1F("hSumTrkEndGapAdc",
00350   //                  "hSumTrkEndGapAdc",1000,-1,5000);
00351   //}
00352 
00353   //look in the last 7 planes
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     //check if there was a tracked strip in this plane
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   //output value
00403   if (adcInGap) *adcInGap=sumAdcInGap;
00404 
00405   return badGap;
00406 }

00235 {
00236   static Int_t lastEntry=-1;
00237   static map<Int_t,Int_t> evtPerSlc;
00238 
00239   //build the map if new entry
00240   if (entry!=lastEntry){
00241     MAXMSG("MeuCuts",Msg::kDebug,30)
00242       <<"Building map..."<<endl;
00243 
00244     evtPerSlc.clear();
00245     
00246     //cache the entry the map was built with
00247     lastEntry=entry;
00248     
00249     TClonesArray& evtTca=(*ntp.evt);
00250     const Int_t numEvts=evtTca.GetEntriesFast();
00251     
00252     //loop over the events
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       //count up the number of events per slice
00259       evtPerSlc[evt.slc]++;      
00260     }
00261 
00262     Bool_t print=false;
00263     if (print){
00264       //make this into a function
00265       //Bool_t FilterMultiEvtPerSlc(ntp,slc,entry)
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;//true if evts != exactly 1 in this slc
00285 }

01406 {
01407   //only run this for the ND at present
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   //this algorithm assumes forward going showers so
01416   //only run filter for spills (don't worry about shws for cosmics)
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   //loop over showers in evt
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     //define a vtx shower to be one that starts within 5 pl
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   //if a vtx shower was found then can check if it was close to window
01451   if (vtxShw){
01452     winPlFromShwEnd=s.WinVtxSidePl-shwEndPl;
01453     //require that shw stops 10 planes before window starts
01454     if (winPlFromShwEnd<10) badShwDist=true;
01455   }
01456 
01457   //write out shw dist if requested
01458   if (shwDist) *shwDist=winPlFromShwEnd;
01459   
01460   return badShwDist;
01461 }

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     //work out if the plane is outside the track range on the stop
00305     //end side
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     //if there is more than 15 pe beyond the end of a stopping muon
00323     //track it is almost certainly misreconstructed.
00324     //this is dangerous for containment and it also affect the position
00325     //of the track window
00326     //15 pe = 15*70 = 1050
00327     //but changed to be 500, since the distribution falls off strongly
00328     if (sumSigCorBeyondEnd>500 || sumAdcBeyondEnd>500) badTracking=true;
00329   }
00330   
00331   //write out signal quantity if requested
00332   if (sigBeyondEnd!=0) *sigBeyondEnd=sumAdcBeyondEnd;
00333 
00334   return badTracking;
00335 }

01332 {
01333   Bool_t badTrkTimes=false;
01334 
01335   //reset variables
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     //this loop is just over the tracked strips
01351     for (Int_t i=0;i<trk.nstrip;i++){
01352 
01353       //check for bug where strip index is -1
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       //get earliest and latest times in snarl
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   //get the median time from the map
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   //set the median time in the summary
01396   s.MedianTime=medianTimeTrk;
01397   
01398   return badTrkTimes;
01399 }

00181 {
00182   const AtmosTrack* track=dynamic_cast<const AtmosTrack*>
00183     (ev.TrackList->At(0));
00184   //cut out events with no track
00185   if (!track) {
00186     //fAwayFromCoil=false;//set to false by default
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     //fAwayFromCoil=false;//set to false by default
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     //fAwayFromCoil=false;//set to false by default
00203     return true;
00204   } 
00205 
00207   //loop over the strips in the snarl
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     //cut out the non-tracked strips
00214     if (!st.Trk) continue;
00215     
00216     //calc the x,y,z positions
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;//true=filter this one
00225     }    
00226   }//end of loop over strips
00227 
00228   return false;//no need to filter
00229 }

00466 {
00467   const Float_t material01Pl=(5.94)*Munits::cm;//including air gap
00468   const Float_t material16Pl=16*material01Pl;
00469   const Float_t material14Pl=14*material01Pl;
00470 
00471   Bool_t lowMatTrav=true;
00472   
00473   //if a window was completed
00474   if (s.TotalMatTraversed>material16Pl+material14Pl){
00475     lowMatTrav=false;
00476   }
00477   
00478   return lowMatTrav;
00479 }

01692 {
01693   static BeamMonSpill spill;
01694   static BMSpillAna bmsa;
01695 
01696   // change default values
01697   static Registry r;
01698   r.UnLockValues();
01699   r.Set("TargetIn",1);
01700   r.Set("BeamType",-1);//don't care if ME,LE,etc
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     //only count new spills, as to events in good spills
01723     if (entry!=lastEntry) goodSpillCounter++;
01724   }
01725   lastEntry=entry;
01726 }

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   //loop over the planes
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     //check if radius is smaller than current smallest
01667     Float_t radius=sqrt(pow(cp.X,2)+pow(cp.Y,2));
01668     if (radius<smallestRadius) smallestRadius=radius;
01669 
01670   }//end of loop over strips
01671 
01672   //analyse coil proximity
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   //set the output variable
01682   if (r) *r=smallestRadius;
01683 
01684   return awayFromCoil;
01685 }

01623 {
01624   //cut out FD spills
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 }

03248 {
03249   Bool_t pass = true;
03250   //Only check FD
03251   if (ntp.GetHeader().GetVldContext().GetDetector() != Detector::kFar)
03252     return pass;
03253 
03254   NtpSRDataQuality dataqua = ntp.dataquality;
03255   // Cuts applied to the far detector data in the atmospheric analysis:
03256   // require full readout
03257   pass = pass && dataqua.cratemask == 16;
03258 
03259   // remove spill triggers
03260   pass = pass && dataqua.trigsource < 130000;
03261 
03262   // remove light injection
03263   pass = pass && !( dataqua.litrigger == 1 
03264                     && dataqua.lisubtractedtime > -40000.0 
03265                     && dataqua.lisubtractedtime < +1000.0 );
03266 
03267   // remove events with lots of busy or cold chips
03268   pass = pass && (dataqua.busychips + dataqua.coldchips) < 20;  
03269 
03270   return pass;
03271 
03272 }

00077 {
00078   // Is the event vertex, (x,y,u,v,z) inside the fiducial volume
00079   // used by the pittsburgh group?
00080   // code from Mike
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 }

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   //Loop over tracks
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     //this loop is just over the tracked strips
02327     for (Int_t i=0;i<trk.nstrip;i++){
02328 
02329       //check for bug where strip index is -1
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   //this loop is over ALL the strips
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 }

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 }


Public Member Functions
	MeuCuts ()
	~MeuCuts ()
Bool_t	AnalyseCoilProximity (AtmosEvent &ev, const MeuSummary &s, Float_t *r=0) const
void	CalcXYZ (const AtmosEvent &ev, const AtmosStrip &st, Float_t &x, Float_t &y, Float_t &z) const
Bool_t	CheckTrkViews (const MeuSummary &s, Int_t diffAtVtx, Int_t diffAtEnd) const
void	ConvertToXY (const MeuSummary &s, std::map< Int_t, MeuHitInfo > &plInfo) const
void	ExtractPlInfo (const NtpStRecord &ntp, MeuSummary &s, std::map< Int_t, MeuHitInfo > &plInfo, const NtpSREvent &evt) const
void	ExtractPlInfo (const AtmosEvent &ntp, MeuSummary &s, std::map< Int_t, MeuHitInfo > &plInfo) const
void	ExtractTruthInfo (const NtpStRecord &ntp, MeuSummary &s, std::map< Int_t, MeuHitInfo > &plInfo) const
void	ExtractTruthInfo (const AtmosEvent &ev, MeuSummary &s, std::map< Int_t, MeuHitInfo > &plInfo) const
void	FillBeamMonDetails (const NtpBDLiteRecord &bd, MeuSummary &s) const
void	FillEnVsDist (const MeuSummary &s, const std::map< Int_t, MeuHitInfo > &plInfo) const
void	FillEnVsDist2 (const MeuSummary &s, const std::map< Int_t, MeuHitInfo > &plInfo) const
void	FillSTSumDetails (const NtpStRecord &ntp, MeuSummary &s, Int_t evt, Int_t slc) const
void	FillSTSumDetails (const AtmosEvent &ev, MeuSummary &s) const
void	FillSTSumRecoDetails (MeuSummary &s, const std::map< Int_t, MeuHitInfo > &plInfo) const
void	FillTimeHistos (const NtpStRecord &ntp, const NtpSREvent &evt, const MeuSummary &s) const
Bool_t	FilterBadEndGaps (const MeuSummary &s, const std::map< Int_t, MeuHitInfo > &plInfo, Float_t *adcInGap=0) const
Bool_t	FilterBadEvtPerSlc (const NtpStRecord &ntp, Int_t slc, Int_t entry, Int_t *evts=0) const
Bool_t	FilterBadTrackEnd (const MeuSummary &s, const std::map< Int_t, MeuHitInfo > &plInfo, Float_t *sigBeyondEnd=0) const
Bool_t	FilterBadDistEndStrip (const MeuSummary &s, const std::map< Int_t, MeuHitInfo > &plInfo, Float_t *minDistEndStrip=0) const
Bool_t	FilterBadShwDist (const NtpStRecord &ntp, MeuSummary &s, const NtpSREvent &evt, Int_t *shwDist=0) const
Bool_t	FilterBadTrkTimes (const NtpStRecord &ntp, MeuSummary &s, const NtpSREvent &evt) const
Bool_t	FilterBadXY (const AtmosEvent &ev, const MeuSummary &s) const
Bool_t	FilterLowMatTraversed (const MeuSummary &s) const
void	GetBDSelectSpillInfo (const NtpBDLiteRecord &ntpBD, MeuSummary &s, Int_t entry, Int_t &goodSpillCounter) const
Bool_t	IsAwayFromCoil (const MeuSummary &s, const std::map< Int_t, MeuHitInfo > &plInfo, Float_t *r=0) const
Bool_t	IsCorrectTrigSrc (const MeuSummary &s) const
Bool_t	IsInPittFidVol (Float_t x, Float_t y, Float_t z, Float_t u, Float_t v) const
Bool_t	IsGoodDataQuality (const NtpStRecord &ntp) const
void	PrintNtpSt (const NtpStRecord &ntp) const
void	SetSpecificCuts (MeuSummary &s) const

MeuCuts Class Reference

Public Member Functions

Constructor & Destructor Documentation

Member Function Documentation