loon: MadQuantities.cxx Source File

00001 #ifndef madquantities_cxx
00002 #define madquantities_cxx
00003 #include <iostream>
00004 #include <cmath>
00005 #include "TClonesArray.h"
00006 #include "TMath.h"
00007 #include "TFile.h"
00008 #include "TH1.h"
00009 #include "TH2.h"
00010 #include "Conventions/Munits.h"
00011 #include "Mad/MadQuantities.h"
00012 #include "Mad/SpillInfo.h"
00013 #include "Mad/SpillInfoBlock.h"
00014 #include "BeamDataUtil/BeamMonSpill.h"
00015 #include "BeamDataUtil/BDSpillAccessor.h"
00016 #include "SpillTiming/SpillTimeFinder.h"
00017 #include "DataUtil/infid_sr_interface.h"
00018 //#include "DataUtil/EnergyCorrections.h"
00019 
00020 
00021 using namespace EnergyCorrections;
00022 
00023 //********************************************************
00024 
00025 
00026 MadQuantities::MadQuantities(TChain *chainSR,TChain *chainMC,
00027                              TChain *chainTH,TChain *chainEM) 
00028    : MadBase(chainSR, chainMC, chainTH, chainEM)
00029 {
00030 
00031   if(!chainSR) {
00032     record = 0;
00033     strecord = 0;
00034     emrecord = 0;
00035     mcrecord = 0;
00036     threcord = 0;
00037     Clear();
00038     whichSource = -1;
00039     isMC = true;
00040     isTH = true;
00041     isEM = true;
00042     Nentries = -1;
00043     return;
00044   }
00045 
00046   fNumFSGeantino = 0;
00047   fGeantinoEnergy = 0;
00048   fNumFSProton = 0;
00049   fTotFSProtonEnergy = 0;
00050   fMaxFSProtonEnergy = 0;
00051   fNumFSNeutron = 0;
00052   fTotFSNeutronEnergy = 0;
00053   fMaxFSNeutronEnergy = 0;
00054   fNumFSPiPlus = 0;
00055   fTotFSPiPlusEnergy = 0;
00056   fMaxFSPiPlusEnergy = 0;
00057   fNumFSPiMinus = 0;
00058   fTotFSPiMinusEnergy = 0;
00059   fMaxFSPiMinusEnergy = 0;
00060   fNumFSPiZero = 0;
00061   fTotFSPiZeroEnergy = 0;
00062   fMaxFSPiZeroEnergy = 0;
00063   fNumFSKaon = 0;
00064   fTotFSKaonEnergy = 0;
00065   fMaxFSKaonEnergy = 0;
00066   fCurStdHepEntry = -1;
00067   fCurTruthIndex = -1;
00068 
00069   InitChain(chainSR,chainMC,chainTH,chainEM);
00070   whichSource = 0;
00071 
00072 }
00073 
00074 //********************************************************
00075 
00076 MadQuantities::MadQuantities(JobC *j,string path,int entries)
00077    : MadBase(j,path,entries)
00078 {
00079   record = 0;
00080   strecord = 0;
00081   emrecord = 0;
00082   mcrecord = 0;
00083   threcord = 0;
00084   Clear();
00085   isMC = true;
00086   isTH = true;
00087   isEM = true;
00088   Nentries = entries;
00089   whichSource = 1;
00090   jcPath = path;
00091   fJC = j;
00092   fChain = NULL;
00093 
00094   fNumFSGeantino = 0;
00095   fGeantinoEnergy = 0;
00096   fNumFSProton = 0;
00097   fTotFSProtonEnergy = 0;
00098   fMaxFSProtonEnergy = 0;
00099   fNumFSNeutron = 0;
00100   fTotFSNeutronEnergy = 0;
00101   fMaxFSNeutronEnergy = 0;
00102   fNumFSPiPlus = 0;
00103   fTotFSPiPlusEnergy = 0;
00104   fMaxFSPiPlusEnergy = 0;
00105   fNumFSPiMinus = 0;
00106   fTotFSPiMinusEnergy = 0;
00107   fMaxFSPiMinusEnergy = 0;
00108   fNumFSPiZero = 0;
00109   fTotFSPiZeroEnergy = 0;
00110   fMaxFSPiZeroEnergy = 0;
00111   fNumFSKaon = 0;
00112   fTotFSKaonEnergy = 0;
00113   fMaxFSKaonEnergy = 0;
00114 
00115 }
00116 
00117 //********************************************************
00118 
00119 MadQuantities::~MadQuantities()
00120 {
00121 
00122 }
00123 
00124 //********************************************************
00125 
00126 Bool_t MadQuantities::IsTrack(){
00127   if(eventSummary->ntrack>0) return true;
00128   return false;
00129 }
00130 
00131 //********************************************************
00132 
00133 Bool_t MadQuantities::IsSingleTrack(){
00134   if(eventSummary->ntrack==1) return true;
00135   return false;
00136 }
00137 
00138 //********************************************************
00139 
00140 Bool_t MadQuantities::PassTrackCuts(Int_t itrk){
00141   if(!LoadTrack(itrk)) return false; //load track
00142   if(ntpTrack->fit.pass==0) return false; //track fit passes
00143   if(ntpTrack->momentum.qp!=ntpTrack->momentum.qp) return false; //no nans
00144   if(fabs(ntpTrack->momentum.qp)<0.0005) return false; //no moms>2000
00145   return true;
00146 }
00147 
00148 //********************************************************
00149 
00150 Bool_t MadQuantities::PassCuts(Int_t itrk){
00151   //check basic fit track quality + check for contained vtx:
00152   if(!PassTrackCuts(itrk)) return false;
00153   if(!IsFidVtx(itrk)) return false;
00154   return true;
00155 }
00156 
00157 //********************************************************
00158 Bool_t MadQuantities::PassCuts(){  //specifically for event display
00159   //demand PassCuts for the first track
00160   if(!PassCuts(0)) return false;
00161   //make any cuts you like to define events to be displayed:
00162   if(TrueMuEnergy(0)==0) return false;
00163   return true;
00164 }
00165 
00166 //********************************************************
00167 Bool_t MadQuantities::IsFidVtx(Int_t itrk){
00168   if(!LoadTrack(itrk)) return false;
00169   
00170   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00171     
00172     if(ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>29.4 ||   //ends
00173        (ntpTrack->vtx.z<16.5&&ntpTrack->vtx.z>14.5) ||  //between SMs
00174        sqrt((ntpTrack->vtx.x*ntpTrack->vtx.x)           //radial cut
00175             +(ntpTrack->vtx.y*ntpTrack->vtx.y))>3.5 ||
00176        sqrt((ntpTrack->vtx.x*ntpTrack->vtx.x)           //radial cut
00177             +(ntpTrack->vtx.y*ntpTrack->vtx.y))<0.4) return false;
00178     
00179   }
00180   else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00181     if(ntpTrack->vtx.z<1 || ntpTrack->vtx.z>5 ||
00182        sqrt(((ntpTrack->vtx.x-1.4828)*(ntpTrack->vtx.x-1.4828)) +
00183             ((ntpTrack->vtx.y-0.2384)*(ntpTrack->vtx.y-0.2384)))>1) return false;
00184     //      sqrt(((ntpTrack->vtx.x-1.4885)*(ntpTrack->vtx.x-1.4885)) +
00185     //      ((ntpTrack->vtx.y-0.1397)*(ntpTrack->vtx.y-0.1397)))>1) return false;
00186   }
00187   
00188   return true;
00189 }
00190 
00191 //********************************************************
00192 Bool_t MadQuantities::IsFidVtxEvt(Int_t ievt){
00193   if(!LoadEvent(ievt)) return false;
00194   
00195   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00196     
00197     if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||   //ends
00198        (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||  //between SMs
00199        sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
00200             +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5 ||
00201        sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
00202             +(ntpEvent->vtx.y*ntpEvent->vtx.y))<0.4) return false;
00203     
00204   }
00205   else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00206     if(ntpEvent->vtx.z<1 || ntpEvent->vtx.z>5 ||
00207        sqrt(((ntpEvent->vtx.x-1.4828)*(ntpEvent->vtx.x-1.4828)) +
00208             ((ntpEvent->vtx.y-0.2384)*(ntpEvent->vtx.y-0.2384)))>1) return 
00209 false;
00210     //       sqrt(((ntpEvent->vtx.x-1.4885)*(ntpEvent->vtx.x-1.4885)) +
00211     //      ((ntpEvent->vtx.y-0.1397)*(ntpEvent->vtx.y-0.1397)))>1) return 
00212     //false;
00213   }
00214   
00215   return true;
00216 }
00217 
00218 
00219 //********************************************************
00220 Bool_t MadQuantities::IsFidVtxTrk(Int_t itrk){
00221   if(!LoadTrack(itrk)) return false;
00222   
00223   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00224     
00225     if(ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>29.4 ||   //ends
00226        (ntpTrack->vtx.z<16.5&&ntpTrack->vtx.z>14.5) ||  //between SMs
00227        sqrt((ntpTrack->vtx.x*ntpTrack->vtx.x)           //radial cut
00228             +(ntpTrack->vtx.y*ntpTrack->vtx.y))>3.5 ||
00229        sqrt((ntpTrack->vtx.x*ntpTrack->vtx.x)           //radial cut
00230             +(ntpTrack->vtx.y*ntpTrack->vtx.y))<0.4) return false;
00231     
00232   }
00233   else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00234     if(ntpTrack->vtx.z<1 || ntpTrack->vtx.z>5 ||
00235        sqrt(((ntpTrack->vtx.x-1.4828)*(ntpTrack->vtx.x-1.4828)) +
00236             ((ntpTrack->vtx.y-0.2384)*(ntpTrack->vtx.y-0.2384)))>1) return false;
00237     //       sqrt(((ntpTrack->vtx.x-1.4885)*(ntpTrack->vtx.x-1.4885)) +
00238     //      ((ntpTrack->vtx.y-0.1397)*(ntpTrack->vtx.y-0.1397)))>1) return false;
00239   }
00240   return true;
00241 }
00242 
00243 
00244 //********************************************************
00245 
00246 Bool_t MadQuantities::IsFidFD_AB(Int_t event)
00247 { 
00248 
00249   if(!LoadEvent(event)) return false;
00250 
00251   Int_t track = -1;
00252   LoadLargestTrackFromEvent(event,track);
00253 
00254   // only check for FD
00255 
00256   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00257  
00258 
00259     // for events with no track, use event vertex
00260     
00261     if (track==-1 && (ntpEvent->vtx.plane<=4 || 
00262                       ntpEvent->vtx.plane>=466 ||     // ends
00263                       (ntpEvent->vtx.plane<=253 && 
00264                        ntpEvent->vtx.plane>=241) ||   // between SMs
00265                       pow(ntpEvent->vtx.x,2)+         // radial cut
00266                       pow(ntpEvent->vtx.y,2)>14)) return false;
00267     
00268     
00269     // otherwise, use track vertex
00270     
00271     if (track!=-1 && (ntpTrack->vtx.plane<=4 || 
00272                       ntpTrack->vtx.plane>=466 ||   //ends
00273                       (ntpTrack->vtx.plane<=253 && 
00274                        ntpTrack->vtx.plane>=241) ||  //between SMs
00275                       pow(ntpTrack->vtx.x,2)+           //radial cut
00276                       pow(ntpTrack->vtx.y,2)>14)) return false;
00277 
00278   }
00279 
00280   return true; 
00281 
00282 }
00283 
00284 
00285 //********************************************************
00286 Bool_t MadQuantities::IsFid_2008(Int_t event)
00287 { 
00288 
00289 //first up initialise the cuts to be the cc2008 version
00290 //just need to do this once
00291   static Bool_t infid_initialised=false;
00292   if (!infid_initialised) {
00293      infid_initialised=true;
00294      choose_infid_set("cc2008");
00295   }
00296 
00297   if(!LoadEvent(event)) return false;
00298 
00299   Int_t track = -1;
00300   LoadLargestTrackFromEvent(event,track);
00301 
00302 
00303   bool IsInFid=false;
00304 
00305 //  Then find out if the track is in the fiducial volume:
00306   if(track!=-1){
00307      //if it has a track, use the track vertex
00308      IsInFid=infid(*strecord,*ntpTrack);
00309   }
00310   else{
00311      // for events with no track, use event vertex
00312      IsInFid=infid(*strecord,*ntpEvent);
00313   }
00314 
00315   return IsInFid;
00316 
00317 }
00318 
00319 
00320 //********************************************************
00321 
00322 
00323 
00324 Bool_t MadQuantities::IsFidAll(Int_t itrk){
00325   
00326   if(!LoadTrack(itrk)) return false;
00327 
00328   // 'empirical' definiton of CEV-like events
00329  
00330 
00331   // new definition - coil hole cuts removed for cedar
00332   bool contained=true;
00333   if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) {//near det
00334     
00335     if (!(ntpTrack->end.z<15 && 
00336           ntpTrack->end.x<2.7 && ntpTrack->end.x>-1.65 && 
00337           ntpTrack->end.y<1.65 && ntpTrack->end.y>-1.65 &&
00338           ntpTrack->end.y>(-ntpTrack->end.x)-1.65 && 
00339           ntpTrack->end.y<ntpTrack->end.x+1.65 &&
00340           ntpTrack->end.y<(-ntpTrack->end.x)+3.55 && 
00341           ntpTrack->end.y>ntpTrack->end.x-3.55)) {contained=false;}
00342     
00343     //updated for 2008 analysis
00344     if (ntpTrack->end.x>1.3) contained=ntpTrack->contained;
00345   }
00346   else if(pow(ntpTrack->end.x,2)+pow(ntpTrack->end.y,2)>14 || ntpTrack->end.plane>475 || ntpTrack->end.plane<5) contained=false;
00347 
00348   return contained;
00349 
00350 
00351 
00352 
00353 // old definition - cuts out ND & FD events close to coil in cedar...
00354 
00355 /*
00356  if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) {//near det
00357     
00358     if (!(ntpTrack->end.z<16 && sqrt(pow(ntpTrack->end.x,2)+
00359                                      pow(ntpTrack->end.y,2))>0.4 &&
00360           ntpTrack->end.x<2.7 && ntpTrack->end.x>-1.65 && 
00361           ntpTrack->end.y<1.65 && ntpTrack->end.y>-1.65 &&
00362           ntpTrack->end.y>(-ntpTrack->end.x)-1.65 && 
00363           ntpTrack->end.y<ntpTrack->end.x+1.65 &&
00364           ntpTrack->end.y<(-ntpTrack->end.x)+3.55 && 
00365           ntpTrack->end.y>ntpTrack->end.x-3.55)) {return false;}
00366      }      
00367      else if(ntpTrack->fidall.dr<0.5||ntpTrack->fidall.dz<0.5) return false;
00368  */
00369 
00370 
00371   return true;
00372 
00373 }
00374 
00375 //********************************************************
00376 Bool_t MadQuantities::IsFidAllEvt(Int_t ievt){
00377   Int_t track = -1;
00378   if(!LoadLargestTrackFromEvent(ievt,track)) return false;
00379   
00380   // 'empirical' definiton of CEV-like events
00381 
00382   // new definition - coil hole cuts removed for cedar
00383 
00384   if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) {//near det
00385     
00386     if (!(ntpTrack->end.z<15 && 
00387           ntpTrack->end.x<2.7 && ntpTrack->end.x>-1.65 && 
00388           ntpTrack->end.y<1.65 && ntpTrack->end.y>-1.65 &&
00389           ntpTrack->end.y>(-ntpTrack->end.x)-1.65 && 
00390           ntpTrack->end.y<ntpTrack->end.x+1.65 &&
00391           ntpTrack->end.y<(-ntpTrack->end.x)+3.55 && 
00392           ntpTrack->end.y>ntpTrack->end.x-3.55)) {return false;}
00393     
00394 
00395   }
00396   else if(pow(ntpTrack->end.x,2)+pow(ntpTrack->end.y,2)>14 || ntpTrack->end.plane>475 || ntpTrack->end.plane<5) return false;
00397 
00398 
00399 
00400 
00401 // old definition - cuts out ND & FD events close to coil in cedar...
00402 
00403 /*
00404  if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) {//near det
00405     
00406     if (!(ntpTrack->end.z<16 && sqrt(pow(ntpTrack->end.x,2)+
00407                                      pow(ntpTrack->end.y,2))>0.4 &&
00408           ntpTrack->end.x<2.7 && ntpTrack->end.x>-1.65 && 
00409           ntpTrack->end.y<1.65 && ntpTrack->end.y>-1.65 &&
00410           ntpTrack->end.y>(-ntpTrack->end.x)-1.65 && 
00411           ntpTrack->end.y<ntpTrack->end.x+1.65 &&
00412           ntpTrack->end.y<(-ntpTrack->end.x)+3.55 && 
00413           ntpTrack->end.y>ntpTrack->end.x-3.55)) {return false;}
00414      }      
00415      else if(ntpTrack->fidall.dr<0.5||ntpTrack->fidall.dz<0.5) return false;
00416  */
00417 
00418   return true;
00419 
00420 }
00421 
00422 //*******************************************************************
00423 
00424 
00425 Float_t* MadQuantities::CCNCSepVars(Int_t itrk){
00426   
00427   Float_t *sepVars = new Float_t[3];
00428   sepVars[0] = sepVars[1] = sepVars[2] = 0.0;
00429   
00430   if(!LoadTrack(itrk)) return sepVars;
00431 
00432   Int_t numtrkstp = ntpTrack->nstrip;
00433   Int_t *trkstrips = ntpTrack->stp;
00434   Float_t planeCharge[500] = {};
00435   Float_t totTrkCharge = 0;
00436   Int_t firstPlane = 500;
00437   Int_t lastPlane = 0;
00438 
00439   TClonesArray* pointStripArray = NULL;
00440   if(isST) pointStripArray = (strecord->stp);
00441   else pointStripArray = (record->stp);
00442   TClonesArray& stripArray = *pointStripArray;
00443 
00444   for(Int_t i=0;i<numtrkstp;i++){
00445     Int_t index = trkstrips[i];
00446     ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
00447     totTrkCharge += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
00448     Int_t tmpPln = ntpStrip->plane;
00449     planeCharge[tmpPln] += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor; 
00450     if(ntpStrip->plane<firstPlane) firstPlane = ntpStrip->plane;
00451     if(ntpStrip->plane>lastPlane) lastPlane = ntpStrip->plane;
00452   }
00453   
00454   Int_t numPlanes = (lastPlane-firstPlane+1);
00455   Float_t chargePerPlane = totTrkCharge/ntpTrack->ds;
00456   Float_t firstHalf = 0;
00457   Float_t lastHalf = 0.00001;
00458   
00459   for(Int_t i=firstPlane;i<=lastPlane;i++){
00460     if(numPlanes%2==0){
00461       if(i<firstPlane+Int_t(numPlanes/2)) firstHalf+=planeCharge[i];
00462       else lastHalf+=planeCharge[i];
00463     }
00464     else {
00465       if(i<firstPlane+Int_t(numPlanes/2)) firstHalf+=planeCharge[i];
00466       else if(i==firstPlane+Int_t(numPlanes/2)+1) {
00467         firstHalf+=planeCharge[i];
00468         lastHalf+=planeCharge[i];
00469       }
00470       else lastHalf+=planeCharge[i];
00471     }
00472   }
00473   Float_t halfRatio = firstHalf/lastHalf;
00474   
00475   sepVars[0] = ntpTrack->ds;
00476   sepVars[1] = chargePerPlane;
00477   sepVars[2] = halfRatio;
00478   
00479   return sepVars;
00480 
00481 }
00482 
00483 //********************************************************
00484 Bool_t MadQuantities::PassQuasiCuts(Int_t itrk,Float_t frac){
00485   if(!PassCuts(itrk)) return false;
00486   if(RecoShwEnergy()>frac*(RecoShwEnergy()+RecoMuEnergy(0,itrk))) return false;
00487   return true;
00488 }
00489 
00490 //********************************************************
00491 
00492 Int_t MadQuantities::INu(Int_t itr){
00493   if(!LoadTruth(itr)) return 0;
00494   return ntpTruth->inu;
00495 }
00496 
00497 //********************************************************
00498 
00499 Int_t MadQuantities::INuNoOsc(Int_t itr){
00500   if(!LoadTruth(itr)) return 0;
00501   return ntpTruth->inunoosc;
00502 }
00503 
00504 //********************************************************
00505 
00506 Int_t MadQuantities::Flavour(Int_t itr){
00507   Int_t f2=abs(INu(itr));
00508   Int_t flavcode=0;
00509   switch (f2) {
00510   case 12:
00511     flavcode=1; // nu_e
00512     break;
00513   case 14:
00514     flavcode=2; // nu_mu
00515     break;
00516   case 16:
00517     flavcode=3; // nu_tau
00518     break;
00519   default:
00520     flavcode=0;
00521     break;
00522   }
00523   return flavcode;
00524 }
00525 
00526 //********************************************************
00527 
00528 Int_t MadQuantities::IAction(Int_t itr){
00529   if(!LoadTruth(itr)) return 0;
00530   return ntpTruth->iaction;
00531 }
00532 
00533 //********************************************************
00534 
00535 Int_t MadQuantities::IResonance(Int_t itr){
00536   if(!LoadTruth(itr)) return 0;
00537   return ntpTruth->iresonance;
00538 }
00539 
00540 //******STDHEP********************************************
00541 
00542 void MadQuantities::SetStdHepVars(Int_t itr){
00543 
00544   fNumFSGeantino = 0;
00545   fGeantinoEnergy = 0;  
00546   fNumFSProton = 0;
00547   fTotFSProtonEnergy = 0;
00548   fMaxFSProtonEnergy = 0;
00549   fNumFSNeutron = 0;
00550   fTotFSNeutronEnergy = 0;
00551   fMaxFSNeutronEnergy = 0;
00552   fNumFSPiPlus = 0;
00553   fTotFSPiPlusEnergy = 0;
00554   fMaxFSPiPlusEnergy = 0;
00555   fNumFSPiMinus = 0;
00556   fTotFSPiMinusEnergy = 0;
00557   fMaxFSPiMinusEnergy = 0;
00558   fNumFSPiZero = 0;
00559   fTotFSPiZeroEnergy = 0;
00560   fMaxFSPiZeroEnergy = 0;
00561   fNumFSKaon = 0;
00562   fTotFSKaonEnergy = 0;
00563   fMaxFSKaonEnergy = 0;
00564 
00565   fCurStdHepEntry = lastEntry;
00566   fCurTruthIndex = itr;
00567   
00568   if(!LoadTruth(itr)) return;
00569 
00570   //otherwise start calculating things
00571 
00572   //First get indices to use:
00573   TClonesArray* pointStdhepArray = NULL;
00574   if(isST) pointStdhepArray = (strecord->stdhep);
00575   else pointStdhepArray = (mcrecord->stdhep);
00576   TClonesArray& stdhepArray = *pointStdhepArray;
00577   Int_t nStdHep = stdhepArray.GetEntries();
00578   Int_t *indicesToUse = new Int_t[nStdHep];
00579   Int_t cnt = 0;
00580   for(int i=0;i<nStdHep;i++){
00581     LoadStdHep(i);
00582     if(ntpStdHep->mc==itr){
00583       indicesToUse[cnt] = i; 
00584       cnt++;
00585     }
00586   }
00587 
00588   //now loop through these and fill variables
00589   for(int i=0;i<cnt;i++){
00590     LoadStdHep(indicesToUse[i]);
00591     if(ntpStdHep->IstHEP==1){ //final state particles:
00592       Float_t mom = TMath::Sqrt(TMath::Abs(ntpStdHep->p4[3]*ntpStdHep->p4[3] - 
00593                                            ntpStdHep->mass*ntpStdHep->mass));
00594       if(ntpStdHep->IdHEP==28) {
00595         fNumFSGeantino+=1;
00596         fGeantinoEnergy+=mom;
00597       }
00598       else if(ntpStdHep->IdHEP==2212) {
00599         fNumFSProton+=1;
00600         fTotFSProtonEnergy+=mom;
00601         if(mom>fMaxFSProtonEnergy) fMaxFSProtonEnergy = mom;
00602       }
00603       else if(ntpStdHep->IdHEP==2112) {
00604         fNumFSNeutron+=1;
00605         fTotFSNeutronEnergy+=mom;
00606         if(mom>fMaxFSNeutronEnergy) fMaxFSNeutronEnergy = mom;
00607       }
00608       else if(ntpStdHep->IdHEP==211) {
00609         fNumFSPiPlus+=1;
00610         fTotFSPiPlusEnergy+=mom;
00611         if(mom>fMaxFSPiPlusEnergy) fMaxFSPiPlusEnergy = mom;
00612       }
00613       else if(ntpStdHep->IdHEP==-211) {
00614         fNumFSPiMinus+=1;
00615         fTotFSPiMinusEnergy+=mom;
00616         if(mom>fMaxFSPiMinusEnergy) fMaxFSPiMinusEnergy = mom;
00617       }
00618       else if(ntpStdHep->IdHEP==111) {
00619         fNumFSPiZero+=1;
00620         fTotFSPiZeroEnergy+=mom;
00621         if(mom>fMaxFSPiZeroEnergy) fMaxFSPiZeroEnergy = mom;
00622       }
00623       else if(abs(ntpStdHep->IdHEP)==321) {
00624         fNumFSKaon+=1;
00625         fTotFSKaonEnergy+=mom;
00626         if(mom>fMaxFSKaonEnergy) fMaxFSPiZeroEnergy = mom;
00627       }
00628     }
00629   }
00630 
00631   delete[] indicesToUse;
00632 }
00633 
00634 Float_t MadQuantities::GeantinoEnergy(Int_t itr){
00635   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00636     SetStdHepVars(itr);
00637   return fGeantinoEnergy;
00638 }
00639 
00640 Int_t MadQuantities::NumFSGeantino(Int_t itr){
00641   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00642     SetStdHepVars(itr);
00643   return fNumFSGeantino;
00644 }
00645 
00646 Int_t MadQuantities::NumFSProton(Int_t itr){
00647   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00648     SetStdHepVars(itr);
00649   return fNumFSProton;
00650 }
00651 
00652 Float_t MadQuantities::TotFSProtonEnergy(Int_t itr){
00653   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00654     SetStdHepVars(itr);
00655   return fTotFSProtonEnergy;
00656 }
00657 
00658 Float_t MadQuantities::MaxFSProtonEnergy(Int_t itr){
00659   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00660     SetStdHepVars(itr);
00661   return fMaxFSProtonEnergy;
00662 }
00663 
00664 Int_t MadQuantities::NumFSNeutron(Int_t itr){
00665   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00666     SetStdHepVars(itr);
00667   return fNumFSNeutron;
00668 }
00669 
00670 Float_t MadQuantities::TotFSNeutronEnergy(Int_t itr){
00671   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00672     SetStdHepVars(itr);
00673   return fTotFSNeutronEnergy;
00674 }
00675 
00676 Float_t MadQuantities::MaxFSNeutronEnergy(Int_t itr){
00677   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00678     SetStdHepVars(itr);
00679   return fMaxFSNeutronEnergy;
00680 }
00681 
00682 Int_t MadQuantities::NumFSPion(Int_t itr,Int_t pm){
00683   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00684     SetStdHepVars(itr);
00685   if(pm==-1) return fNumFSPiMinus;
00686   else if(pm==1) return fNumFSPiPlus;
00687   return fNumFSPiPlus+fNumFSPiMinus;
00688 }
00689 
00690 Float_t MadQuantities::TotFSPionEnergy(Int_t itr,Int_t pm){
00691   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00692     SetStdHepVars(itr);
00693   if(pm==-1) return fTotFSPiMinusEnergy;
00694   else if(pm==1) return fTotFSPiPlusEnergy;
00695   return fTotFSPiPlusEnergy+fTotFSPiMinusEnergy;
00696 }
00697 
00698 Float_t MadQuantities::MaxFSPionEnergy(Int_t itr,Int_t pm){
00699   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00700     SetStdHepVars(itr);
00701   if(pm==-1) return fMaxFSPiMinusEnergy;
00702   else if(pm==1) return fMaxFSPiPlusEnergy;  
00703   if(fMaxFSPiPlusEnergy>fMaxFSPiMinusEnergy) return fMaxFSPiPlusEnergy;
00704   return fMaxFSPiMinusEnergy;
00705 }
00706 
00707 Int_t MadQuantities::NumFSKaon(Int_t itr){
00708   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00709     SetStdHepVars(itr);
00710   return fNumFSKaon;
00711 }
00712 
00713 Float_t MadQuantities::TotFSKaonEnergy(Int_t itr){
00714   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00715     SetStdHepVars(itr);
00716   return fTotFSKaonEnergy;
00717 }
00718 
00719 Float_t MadQuantities::MaxFSKaonEnergy(Int_t itr){
00720   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00721     SetStdHepVars(itr);
00722   return fMaxFSKaonEnergy;
00723 }
00724 
00725 Int_t MadQuantities::NumFSPiZero(Int_t itr){
00726   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00727     SetStdHepVars(itr);
00728   return fNumFSPiZero;
00729 }
00730 
00731 Float_t MadQuantities::TotFSPiZeroEnergy(Int_t itr){
00732   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00733     SetStdHepVars(itr);
00734   return fTotFSPiZeroEnergy;
00735 }
00736 
00737 Float_t MadQuantities::MaxFSPiZeroEnergy(Int_t itr){
00738   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00739     SetStdHepVars(itr);
00740   return fMaxFSPiZeroEnergy;
00741 }
00742 
00743 //********************************************************
00744 
00745 Float_t MadQuantities::Y(Int_t itr){
00746   if(!LoadTruth(itr)) return 0;
00747   return ntpTruth->y;
00748 }
00749 
00750 //********************************************************
00751 
00752 Float_t MadQuantities::X(Int_t itr){
00753   if(!LoadTruth(itr)) return 0;
00754   return ntpTruth->x;
00755 }
00756 
00757 //********************************************************
00758 
00759 Float_t MadQuantities::Q2(Int_t itr){
00760   if(!LoadTruth(itr)) return 0;
00761   return ntpTruth->q2;
00762 }
00763 
00764 //********************************************************
00765 
00766 Float_t MadQuantities::W2(Int_t itr){
00767   if(!LoadTruth(itr)) return 0;
00768   return ntpTruth->w2;
00769 }
00770 
00771 //********************************************************
00772 
00773 Int_t MadQuantities::Nucleus(Int_t itr){
00774  
00775   if(!LoadTruth(itr)) return 0;
00776 
00777   Int_t z=int(ntpTruth->z);
00778   Int_t a=int(ntpTruth->a);
00779   Int_t nucleus=1;
00780 
00781   switch (z) {
00782     //  case 1:
00783     //nucleus=0;   // free nucleon
00784     //break;
00785   case 1:
00786     switch (a) {
00787     case 1:
00788       nucleus=256;   // hydrogen1
00789       break;
00790     case 2:
00791       nucleus=257;   // hydrogen2
00792       break;
00793     case 3:
00794       nucleus=258;   // hydrogen2
00795       break;
00796     default:
00797       nucleus=256;   // hydrogen1
00798       break;
00799     }
00800     break;
00801   case 6:
00802     switch (a) {
00803     case 11:
00804       nucleus=274; // carbon11   
00805       break;
00806     case 12:
00807       nucleus=275; // carbon12
00808       break;
00809     case 13:
00810       nucleus=276; // carbon13
00811       break;
00812     case 14:
00813       nucleus=277; // carbon14
00814       break;
00815     default:
00816       nucleus=275; // carbon12
00817       break;
00818     }
00819     break;
00820   case 7:
00821     switch (a) {
00822     case 13:
00823       nucleus=278; // nitrogen13   
00824       break;
00825     case 14:
00826       nucleus=279; // nitrogen14
00827       break;
00828     case 15:
00829       nucleus=280; // nitrogen15
00830       break;
00831     case 16:
00832       nucleus=281; // nitrogen16
00833       break;
00834     case 17:
00835       nucleus=282; // nitrogen17
00836       break;
00837     default:
00838       nucleus=279; // nitrogen14
00839       break;
00840     }
00841     break;
00842   case 8:
00843     switch (a) {
00844     case 15:
00845       nucleus=283;   // oxygen15
00846       break;
00847     case 16:
00848       nucleus=284;   // oxygen16
00849       break;
00850     case 17:
00851       nucleus=285;   // oxygen17
00852       break;
00853     case 18:
00854       nucleus=286;   // oxygen18
00855       break;
00856     default:
00857       nucleus=284;   // oxygen16
00858       break;
00859     }
00860     break;
00861   case 13:
00862     switch (a) {
00863     case 26:
00864       nucleus=303;   // aluminium26
00865       break;
00866     case 27:
00867       nucleus=304;   // aluminium27
00868       break;
00869     case 28:
00870       nucleus=305;   // aluminium28
00871       break;
00872     case 29:
00873       nucleus=306;   // aluminium29
00874       break;
00875     default:
00876       nucleus=304;   // aluminium27
00877       break;
00878     }
00879     break;
00880   case 14:
00881     switch (a) {
00882     case 27:
00883       nucleus=307;   // silicon27
00884       break;
00885     case 28:
00886       nucleus=308;   // silicon28
00887       break;
00888     case 29:
00889       nucleus=309;   // silicon29
00890       break;
00891     case 30:
00892       nucleus=310;   // silicon30
00893       break;
00894     default:         
00895       nucleus=308;   // silicon28
00896       break;
00897     }
00898     break;
00899   case 15:
00900     switch (a) {
00901     case 30:
00902       nucleus=311;   //phosphorus30
00903       break;
00904     case 31:
00905       nucleus=312;   //phosphorus31
00906       break;
00907     case 32:
00908       nucleus=313;   //phosphorus32
00909       break;
00910     case 33:
00911       nucleus=314;   //phosphorus33
00912       break;
00913     default:
00914       nucleus=312;   //phosphorus31
00915       break;
00916     }
00917     break;
00918   case 16:
00919     switch (a) {
00920     case 31:
00921       nucleus=315;   //sulphur31
00922       break;
00923     case 32:
00924       nucleus=316;   //sulphur32
00925       break;
00926     case 33:
00927       nucleus=317;   //sulphur33
00928       break;
00929     case 34:
00930       nucleus=318;   //sulphur34
00931       break;
00932     case 35:
00933       nucleus=319;   //sulphur35
00934       break;
00935     case 36:
00936       nucleus=320;   //sulphur36
00937       break;
00938     default:
00939       nucleus=316;   //sulphur32
00940       break;
00941     }
00942     break;
00943   case 22:
00944     switch (a) {
00945     case 45:
00946       nucleus=347;   //titanium45
00947       break;
00948     case 46:
00949       nucleus=348;   //titanium46
00950       break;
00951     case 47:
00952       nucleus=349;   //titanium47
00953       break;
00954     case 48:
00955       nucleus=350;   //titanium48
00956       break;
00957     case 49:
00958       nucleus=351;   //titanium49
00959       break;
00960     case 50:
00961       nucleus=352;   //titanium50
00962       break;
00963     default:
00964       nucleus=350;   //titanium48
00965       break;
00966     }
00967     break;
00968   case 23:
00969     switch (a) {
00970     case 49:
00971       nucleus=353;   //vanadium49
00972       break;
00973     case 50:
00974       nucleus=354;   //vanadium50
00975       break;
00976     case 51:
00977       nucleus=355;   //vanadium51
00978       break;
00979     case 52:
00980       nucleus=356;   //vanadium52
00981       break;
00982     case 53:
00983       nucleus=357;   //vanadium53
00984       break;
00985     default:
00986       nucleus=355;   //vanadium51
00987       break;
00988     }
00989     break;
00990   case 24:
00991     switch (a) {
00992     case 49:
00993       nucleus=358;   //chromium49
00994       break;
00995     case 50:
00996       nucleus=359;   //chromium50
00997       break;
00998     case 51:
00999       nucleus=360;   //chromium51
01000       break;
01001     case 52:
01002       nucleus=361;   //chromium52
01003       break;
01004     case 53:
01005       nucleus=362;   //chromium53
01006       break;
01007     case 54:
01008       nucleus=363;   //chromium54
01009       break;
01010     default:
01011       nucleus=361;   //chromium52
01012       break;
01013     }
01014     break;
01015   case 25:
01016     switch (a) {
01017     case 53:
01018       nucleus=364;   //manganese53
01019       break;
01020     case 54:
01021       nucleus=365;   //manganese54
01022       break;    
01023     case 55:
01024       nucleus=366;   //manganese55
01025       break;    
01026     case 56:
01027       nucleus=367;   //manganese56
01028       break;    
01029     case 57:
01030       nucleus=368;   //manganese57
01031       break;    
01032     default:
01033       nucleus=366;   //manganese55
01034       break;
01035     }
01036     break;
01037   case 26:
01038     switch (a) {
01039     case 53:
01040       nucleus=369;   //iron53
01041       break;
01042     case 54:
01043       nucleus=370;   //iron54
01044       break;
01045     case 55:
01046       nucleus=371;   //iron55
01047       break;
01048     case 56:
01049       nucleus=372;   //iron56
01050       break;
01051     case 57:
01052       nucleus=373;   //iron57
01053       break;
01054     case 58:
01055       nucleus=374;   //iron58
01056       break;
01057     default:
01058       nucleus=372;   //iron56
01059       break;
01060     }
01061     break;
01062   case 28:
01063     switch (a) {
01064     case 57:
01065       nucleus=382;   //nickel57
01066       break;
01067     case 58:
01068       nucleus=383;   //nickel58
01069       break;
01070     case 59:
01071       nucleus=384;   //nickel59
01072       break;
01073     case 60:
01074       nucleus=385;   //nickel60
01075       break;
01076     case 61:
01077       nucleus=386;   //nickel61
01078       break;
01079     case 62:
01080       nucleus=387;   //nickel62
01081       break;
01082     case 63:
01083       nucleus=388;   //nickel63
01084       break;
01085     case 64:
01086       nucleus=389;   //nickel64
01087       break;
01088     default:
01089       nucleus=383;   //nickel58
01090       break;
01091     }
01092     break;
01093   case 29:
01094     switch (a) {
01095     case 62:
01096       nucleus=390;   //copper62
01097       break;
01098     case 63:
01099       nucleus=391;   //copper63
01100       break;
01101     case 64:
01102       nucleus=392;   //copper64
01103       break;
01104     case 65:
01105       nucleus=393;   //copper65
01106       break;
01107     case 66:
01108       nucleus=394;   //copper66
01109       break;
01110     case 67:
01111       nucleus=395;   //copper67
01112       break;
01113     default:
01114       nucleus=392;   //copper64
01115       break;
01116     }
01117     break;
01118   default:
01119     nucleus=1;   // unknown
01120     break;
01121   }
01122 
01123   return nucleus;
01124 }
01125 
01126 //********************************************************
01127 
01128 Int_t MadQuantities::Initial_state(Int_t itr){
01129   
01130   if(!LoadTruth(itr)) return 0;
01131   
01132   Int_t initial_state=0;  
01133   TClonesArray* pointStdhepArray = NULL;
01134   if(isST) pointStdhepArray = (strecord->stdhep);
01135   else pointStdhepArray = (mcrecord->stdhep);
01136   TClonesArray& stdhepArray = *pointStdhepArray;
01137   Int_t nStdHep = stdhepArray.GetEntries();
01138   
01139   int protneut = -1;  // 0 = neutron, 1 = proton, 2 = N, 3 = A
01140   int nubarnu = 0;    // +1 = neutrino, -1 = antineutrino
01141   
01142   for(int i=0;i<nStdHep;i++){    
01143     LoadStdHep(i);
01144     if(ntpStdHep->mc==itr){
01145       
01146       if(ntpStdHep->IstHEP==0){  //initial state particle
01147         if(abs(ntpStdHep->IdHEP)==12 || 
01148            abs(ntpStdHep->IdHEP)==14 || 
01149            abs(ntpStdHep->IdHEP)==16){   //(anti)neutrino         
01150           nubarnu = ntpStdHep->IdHEP/abs(ntpStdHep->IdHEP);  //get sign
01151         }
01152       }
01153       if(ntpStdHep->IstHEP==11){    //target nucleon
01154         if(ntpStdHep->IdHEP==2212) protneut = 1;  //proton
01155         else if(ntpStdHep->IdHEP==2112) protneut = 0;  //neutron
01156         else if(abs(ntpStdHep->IdHEP)>1000000000) protneut = 2;  //nucleus
01157         else protneut = 3; //atom - probably never get here since IdHEP A>N?
01158       }
01159     }
01160   }
01161 
01162   if(protneut==1 && nubarnu==1)  initial_state=1;  //p + v
01163   if(protneut==0 && nubarnu==1)  initial_state=2;  //n + v
01164   if(protneut==1 && nubarnu==-1) initial_state=3;  //p + vbar
01165   if(protneut==0 && nubarnu==-1) initial_state=4;  //n + vbar
01166   if(protneut==2 && nubarnu==1)  initial_state=5;  //N + v
01167   if(protneut==3 && nubarnu==1)  initial_state=6;  //A + v
01168   if(protneut==2 && nubarnu==-1) initial_state=7;  //N + vbar
01169   if(protneut==3 && nubarnu==-1) initial_state=8;  //A + vbar
01170 
01171   return initial_state;
01172 }
01173 
01174 //********************************************************
01175 Int_t MadQuantities::HadronicFinalState(Int_t itr){
01176   if(!LoadTruth(itr)) return 0;
01177   Int_t hfs=0;
01178   Int_t proc = IResonance(itr);
01179   TClonesArray* pointStdhepArray = NULL;
01180   if(isST) pointStdhepArray = (strecord->stdhep);
01181   else pointStdhepArray = (mcrecord->stdhep);
01182   TClonesArray& stdhepArray = *pointStdhepArray;
01183   Int_t nStdHep = stdhepArray.GetEntries();
01184   if(proc==1002){
01185     for(int i=0;i<nStdHep;i++){    
01186       LoadStdHep(i);
01187       if(ntpStdHep->mc==itr && ntpStdHep->IstHEP==3 &&
01188          !(abs(ntpStdHep->IdHEP)==15)){  //not a tau lepton
01189         hfs = ntpStdHep->IdHEP;
01190         break;
01191       }
01192     }
01193     hfs = hfs%1000;
01194   }
01195   else {
01196     for(int i=0;i<nStdHep;i++){    
01197       LoadStdHep(i);
01198       if(ntpStdHep->mc==itr && ntpStdHep->IstHEP==3){
01199         hfs = ntpStdHep->IdHEP;
01200         break;
01201       }
01202     }
01203     hfs = hfs%1000;
01204   }
01205   return hfs;
01206 }
01207 
01208 //********************************************************
01209 Float_t *MadQuantities::Target4Vector(Int_t itr){
01210   Float_t *vec = new Float_t[4];
01211   vec[0] = 0.; vec[1] = 0.; vec[2] = 0.; vec[3] = 0.;
01212   if(!LoadTruth(itr)) return vec;
01213   vec[0] = ntpTruth->p4tgt[0];
01214   vec[1] = ntpTruth->p4tgt[1];
01215   vec[2] = ntpTruth->p4tgt[2];
01216   vec[3] = ntpTruth->p4tgt[3];
01217   return vec;
01218 }
01219 
01220 //********************************************************
01221 
01222 Float_t MadQuantities::TrueNuEnergy(Int_t itr){
01223   if(!LoadTruth(itr)) return 0.;
01224   return ntpTruth->p4neu[3];
01225 }
01226 
01227 //********************************************************
01228 Float_t *MadQuantities::TrueNuMom(Int_t itr){
01229   Float_t *mom = new Float_t[3];  
01230   mom[0] = 0.; 
01231   mom[1] = 0.; 
01232   mom[2] = 0.;
01233   if(!LoadTruth(itr)) return mom;
01234   mom[0] = ntpTruth->p4neu[0];
01235   mom[1] = ntpTruth->p4neu[1];  
01236   mom[2] = ntpTruth->p4neu[2];
01237   return mom;
01238 }
01239 
01240 //********************************************************
01241 
01242 Float_t *MadQuantities::TrueVtx(Int_t itr){
01243   Float_t *vtx = new Float_t[3];  
01244   vtx[0] = 0.; 
01245   vtx[1] = 0.; 
01246   vtx[2] = 0.;
01247   if(!LoadTruth(itr)) return vtx;
01248   vtx[0] = ntpTruth->vtxx;
01249   vtx[1] = ntpTruth->vtxy;  
01250   vtx[2] = ntpTruth->vtxz;
01251   return vtx;
01252 }
01253 
01254 //********************************************************
01255 
01256 Float_t MadQuantities::TrueMuEnergy(Int_t itr){
01257   if(!LoadTruth(itr)) return 0;
01258   if(ntpTruth->p4mu1[3]!=0) return fabs(ntpTruth->p4mu1[3]);
01259   return 0;
01260 }
01261 
01262 //********************************************************
01263 Float_t MadQuantities::TrueLeptonEnergy(Int_t itr){
01264   if(!LoadTruth(itr)) return 0;
01265   if(ntpTruth->iaction==0) return 0;
01266   else if(ntpTruth->inu==12) return fabs(ntpTruth->p4el1[3]);
01267   else if(ntpTruth->inu==14) return fabs(ntpTruth->p4mu1[3]);
01268   else if(ntpTruth->inu==16) return fabs(ntpTruth->p4tau[3]);
01269   return 0;
01270 }
01271 
01272 //********************************************************
01273 
01274 Float_t MadQuantities::TrueMuQP(Int_t itr){
01275   if(!LoadTruth(itr)) return 0;
01276   if(ntpTruth->p4mu1[3]!=0) {
01277     Float_t p = 1./sqrt(ntpTruth->p4mu1[3]*ntpTruth->p4mu1[3]-0.10555*0.10555);
01278     Float_t q = ntpTruth->p4mu1[3]/fabs(ntpTruth->p4mu1[3]);
01279     return q*p;
01280   }
01281   return 0;
01282 }
01283 
01284 //********************************************************
01285 
01286 Float_t MadQuantities::TrueShwEnergy(Int_t itr){
01287   if(!LoadTruth(itr)) return 0.;
01288   return ntpTruth->p4shw[3];
01289 }
01290 
01291 //********************************************************
01292 
01293 Float_t MadQuantities::TrueMuDCosZ(Int_t itr){ //dz/ds
01294   if(!LoadTruth(itr)) return 0.;
01295   Float_t mom = sqrt((ntpTruth->p4mu1[0]*ntpTruth->p4mu1[0]) + 
01296                      (ntpTruth->p4mu1[1]*ntpTruth->p4mu1[1]) + 
01297                      (ntpTruth->p4mu1[2]*ntpTruth->p4mu1[2]));
01298   if(mom==0) return 0;
01299   return ntpTruth->p4mu1[2]/mom;
01300 }
01301 
01302 //********************************************************
01303 
01304 Float_t MadQuantities::TrueMuDCosNeu(Int_t itr){ //ds_mu/ds_nu
01305   if(!LoadTruth(itr)) return 0.;
01306   TVector3 *nuvec = new TVector3(ntpTruth->p4neu[0],
01307                                  ntpTruth->p4neu[1],
01308                                  ntpTruth->p4neu[2]);
01309   TVector3 *muvec = new TVector3(ntpTruth->p4mu1[0],
01310                                  ntpTruth->p4mu1[1],
01311                                  ntpTruth->p4mu1[2]);  
01312   Float_t MuAng = nuvec->Angle(*muvec); //angle in rads
01313   delete nuvec;
01314   delete muvec;
01315   return TMath::Cos(MuAng);
01316 }
01317 
01318 //********************************************************
01319 
01320 Float_t MadQuantities::GetTrueShwSC(Int_t itr){
01321 
01322   if(!LoadTruth(itr)) return 0.;
01323   //Requires truth digit information which is not present in ntuples  
01324   Float_t Summed_Shw_SC = 0;
01325   
01326   return Summed_Shw_SC;
01327 }
01328 
01329 //********************************************************
01330 Float_t MadQuantities::RecoMuEnergy(Int_t opt,Int_t itrk,Bool_t isdata,Detector::Detector_t det){
01331   const float mumass=0.10566;
01332 
01333   if(LoadTrack(itrk)){
01334     Float_t mr=ntpTrack->momentum.range;
01335     if(mr>0) mr=CorrectMomentumFromRange(mr,isdata,det);
01336     float mc=0.0; // signed momentum from curvature
01337     if(ntpTrack->momentum.qp!=0) mc = 1.0/(ntpTrack->momentum.qp);
01338     mc=CorrectSignedMomentumFromCurvature(mc,isdata,det);
01339     
01340     if(opt==0){  //return the most appropriate measure of momentum
01341       if(IsFidAll(itrk) || ntpTrack->momentum.qp==0) {
01342         return sqrt(mr*mr+ mumass*mumass);
01343       }
01344       else {
01345         if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
01346         else return sqrt(mc*mc+ mumass*mumass);
01347       }
01348     }
01349     else if(opt==1) {
01350       //return curvature measurement
01351       // in R1.9 the tracker will apparently return (q/p)=0.0
01352       // maybe it's when a track looks perfectly rigid?
01353       // if so, we have to do something
01354       // I don't want to use the range, that could be very wrong
01355       // so, we'll return an obviously ridiculous value of 10 TeV
01356       if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
01357       else return sqrt(mc*mc+ mumass*mumass);
01358     }
01359     else if(opt==2) //return range measurement
01360       return sqrt(mr*mr+ mumass*mumass);
01361     else return 0;
01362   }
01363   return 0.;
01364 }
01365 
01366 //********************************************************
01367 
01368 Float_t MadQuantities::RecoMuEnergyNew(Int_t opt, Int_t itrk, VldContext cx, ReleaseType::Release_t reltype, EnergyCorrections::WhichCorrection_t corrver) {
01369 
01370   // using the new version of energy corrections developed for Cedar/Daikon
01371 
01372   const float mumass=0.10566;
01373 
01374   if(LoadTrack(itrk)){
01375 
01376     Float_t mr=ntpTrack->momentum.range;
01377 
01378     if(mr>0) mr=FullyCorrectMomentumFromRange(mr,cx,reltype,corrver);
01379 
01380 
01381     float mc=0.0; // signed momentum from curvature
01382     if(ntpTrack->momentum.qp!=0) mc = 1.0/(ntpTrack->momentum.qp);
01383 
01384     mc=FullyCorrectSignedMomentumFromCurvature(mc,cx,reltype,corrver);
01385 
01386    
01387     if(opt==0){  //return the most appropriate measure of momentum
01388       if(IsFidAll(itrk) || ntpTrack->momentum.qp==0) {
01389         return sqrt(mr*mr+ mumass*mumass);
01390       }
01391       else {
01392         if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
01393         else return sqrt(mc*mc+ mumass*mumass);
01394       }
01395     }
01396     else if(opt==1) {
01397       //return curvature measurement
01398       // in R1.9 the tracker will apparently return (q/p)=0.0
01399       // maybe it's when a track looks perfectly rigid?
01400       // if so, we have to do something
01401       // I don't want to use the range, that could be very wrong
01402       // so, we'll return an obviously ridiculous value of 10 TeV
01403       if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
01404       else return sqrt(mc*mc+ mumass*mumass);
01405     }
01406     else if(opt==2) //return range measurement
01407       return sqrt(mr*mr+ mumass*mumass);
01408     else return 0;
01409   }
01410   return 0.;
01411 }
01412 
01413 //********************************************************
01414 
01415 Float_t MadQuantities::RecoMuDCosZ(Int_t itrk){ //dz/ds
01416   if(!LoadTrack(itrk)) return 0.;
01417   return ntpTrack->vtx.dcosz;
01418 }
01419 
01420 //********************************************************
01421 
01422 Float_t MadQuantities::RecoMuDCosNeu(Int_t itrk){ //ds_mu/ds_neu
01423   if(!LoadTrack(itrk)) return 0.;
01424   
01425   Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.); //3degree beam
01426   Float_t bl_y = sqrt(1. - bl_z*bl_z);
01427   Float_t costhbl = ntpTrack->vtx.dcosz*bl_z + ntpTrack->vtx.dcosy*bl_y;
01428    
01429   return  costhbl;
01430 }
01431 
01432 //********************************************************
01433 
01434 Float_t MadQuantities::RecoMuQP(Int_t itrk){
01435   if(LoadTrack(itrk)) return ntpTrack->momentum.qp;
01436   return 0.;
01437 }
01438 
01439 //********************************************************
01440 
01441 Float_t MadQuantities::GetNonTrkSC(Int_t itrk,Int_t npln){
01442   
01443   Float_t Summed_Stp_SC = 0;
01444   Float_t Summed_Trk_SC = 0;
01445   
01446   bool useAll = false;
01447   if(npln<=0) useAll=true;  //Most basic thing in the world ever: 
01448                             //if it's not in the track it's in the shower
01449   
01450   TClonesArray* pointStripArray = NULL;
01451   if(isST) pointStripArray = (strecord->stp);
01452   else pointStripArray = (record->stp);
01453   TClonesArray& stripArray = *pointStripArray;
01454   if(LoadTrack(itrk)) {
01455     Int_t numtrkstp = ntpTrack->nstrip;
01456     for(Int_t i=0;i<numtrkstp;i++){
01457       Int_t index = ntpTrack->stp[i];
01458       ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
01459       if((ntpStrip->plane <= ntpTrack->vtx.plane + npln
01460           && ntpStrip->plane >= ntpTrack->vtx.plane - npln)||useAll){
01461         Summed_Trk_SC += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
01462       }
01463     }
01464   }
01465   
01466   for(Int_t i=0;i<int(eventSummary->nstrip);i++) {
01467     if(!LoadStrip(i)) continue;
01468     if((ntpStrip->plane <= ntpTrack->vtx.plane + npln
01469         && ntpStrip->plane >= ntpTrack->vtx.plane - npln)||useAll){
01470       Summed_Stp_SC += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
01471     }
01472   }
01473   return Summed_Stp_SC - Summed_Trk_SC;
01474 }
01475 
01476 //********************************************************
01477 
01478 Float_t MadQuantities::GetShwSC(Int_t itrk){
01479 
01480   //Try something a bit smarter:
01481   //Look for a reconstructed shower near to the track vtx and sum the SigCors
01482   if(!LoadTrack(itrk)) return 0.;
01483   
01484   Float_t Summed_Shw_SC = 0;
01485   Int_t bestshower = -1;
01486   Float_t bestdistance = 1000.;
01487   if(eventSummary->nshower!=0){
01488     for(Int_t i=0;i<eventSummary->nshower;i++){
01489       if(!LoadShower(i)) continue;
01490       Float_t distance = sqrt((ntpTrack->vtx.x - ntpShower->vtx.x)
01491                               *(ntpTrack->vtx.x - ntpShower->vtx.x)
01492                               +(ntpTrack->vtx.y - ntpShower->vtx.y)
01493                               *(ntpTrack->vtx.y - ntpShower->vtx.y)
01494                               +(ntpTrack->vtx.z - ntpShower->vtx.z)
01495                               *(ntpTrack->vtx.z - ntpShower->vtx.z));
01496       
01497       if(distance<0.5){
01498         if(bestshower==-1) {
01499           bestshower = i;
01500           bestdistance = distance;
01501         }
01502         else if(distance<bestdistance) bestshower = i;
01503       }
01504     }
01505     
01506     TClonesArray* pointStripArray = NULL;
01507     if(isST) pointStripArray = (strecord->stp);
01508     else pointStripArray = (record->stp);
01509     TClonesArray& stripArray = *pointStripArray;
01510     if(bestshower!=-1){
01511       LoadShower(bestshower);
01512       Int_t *shwstrips = ntpShower->stp;
01513       for(Int_t j=0;j<ntpShower->nstrip;j++){
01514         Int_t index = shwstrips[j];
01515         ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
01516         Summed_Shw_SC += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
01517       }
01518       return Summed_Shw_SC;
01519     }
01520     return 0;
01521   }
01522   else return 0; //no reconstructed shower
01523 
01524 }
01525 
01526 //********************************************************
01527 
01528 Float_t MadQuantities::RecoShwEnergy(Int_t ishw, Int_t opt, Int_t det, bool isdata, int mode){
01529   // mode is used to deterimine if Niki or Andy's tuning should be used
01530   // they are more-or-less equivalent but we have choosen Niki's
01531   // mode=1 is set in the function declaration, so this choice will
01532   // be automatic for most people.
01533 
01534   Float_t shower_ene=0; 
01535 
01536   if(LoadShower(ishw)) { 
01537   if(det==1){
01538     if(opt==-1) shower_ene = ntpShower->ph.gev;
01539     if(opt==0)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.linCCgev,CandShowerHandle::kCC,mode,isdata);
01540     if(opt==1)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.wtCCgev,CandShowerHandle::kWtCC,mode,isdata);
01541     if(opt==2)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.linNCgev,CandShowerHandle::kNC,mode,isdata);
01542     if(opt==3)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.wtNCgev,CandShowerHandle::kWtNC,mode,isdata);
01543     if(opt==4)  shower_ene = ntpShower->shwph.EMgev;
01544     return shower_ene;
01545   }
01546    if(det==2){
01547     if(opt==-1) shower_ene = ntpShower->ph.gev;
01548     if(opt==0)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.linCCgev,CandShowerHandle::kCC,mode,isdata);
01549     if(opt==1)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.wtCCgev,CandShowerHandle::kWtCC,mode,isdata);
01550     if(opt==2)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.linNCgev,CandShowerHandle::kNC,mode,isdata);
01551     if(opt==3)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.wtNCgev,CandShowerHandle::kWtNC,mode,isdata);
01552     if(opt==4)  shower_ene = ntpShower->shwph.EMgev;
01553     return shower_ene;
01554 
01555   }
01556 } 
01557   return 0;
01558 }
01559 
01560 
01561 //********************************************************
01562 
01563 Float_t MadQuantities::RecoShwEnergyNew(Int_t ishw, Int_t opt, VldContext cx, ReleaseType::Release_t reltype, EnergyCorrections::WhichCorrection_t corrver) {
01564 
01565 
01566   Float_t shower_ene=0; 
01567 
01568   if(LoadShower(ishw)) { 
01569 
01570     if(opt==-1) shower_ene = ntpShower->ph.gev;
01571     if(opt==0)  shower_ene = FullyCorrectShowerEnergy(ntpShower->shwph.linCCgev,CandShowerHandle::kCC,cx,reltype,corrver);
01572     if(opt==1)  shower_ene = FullyCorrectShowerEnergy(ntpShower->shwph.wtCCgev,CandShowerHandle::kWtCC,cx,reltype,corrver);
01573     if(opt==2)  shower_ene = FullyCorrectShowerEnergy(ntpShower->shwph.linNCgev,CandShowerHandle::kNC,cx,reltype,corrver);
01574     if(opt==3)  shower_ene = FullyCorrectShowerEnergy(ntpShower->shwph.wtNCgev,CandShowerHandle::kWtNC,cx,reltype,corrver);
01575     if(opt==4)  shower_ene = ntpShower->shwph.EMgev;
01576     return shower_ene;
01577   }
01578  
01579   return 0;
01580 }
01581 
01582 //********************************************************
01583 //Andy Culling's method:
01584 Float_t MadQuantities::RecoShwEnergySqrt(Int_t ievt){
01585   Float_t sqrtSC = this->GetSqrtShwSC(ievt)+this->GetSqrtTrkSkimSC(ievt);
01586   if(sqrtSC==0) return 0;
01587   Float_t shwEnergy = (sqrtSC<3400) ? 0.108 + sqrtSC*0.0017 + 
01588     sqrtSC*sqrtSC*2.07e-7 : -3.76 + sqrtSC*0.00354;
01589   return shwEnergy;
01590 }
01591 
01592 //********************************************************
01593 //Needed for RecoShwEnergySqrt():
01594 Float_t MadQuantities::GetSqrtTrkSkimSC(Int_t ievt){
01595 
01596   Float_t skimmedSqrtSC = 0;
01597   Int_t track = -1;
01598   if(!LoadLargestTrackFromEvent(ievt,track)) return 0;
01599 
01600   if(ntpTrack->ph.gev>0){
01601     TClonesArray* pointStripArray = NULL;
01602     if(isST) pointStripArray = (strecord->stp);
01603     else pointStripArray = (record->stp);
01604     TClonesArray& stripArray = *pointStripArray;
01605     Int_t numtrkstp = ntpTrack->nstrip;
01606     for(Int_t i=0;i<numtrkstp;i++){
01607       Int_t index = ntpTrack->stp[i];
01608       ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
01609       if(ntpStrip->plane<ntpTrack->vtx.plane+11){
01610         skimmedSqrtSC += ((ntpStrip->ph0.sigcor + 
01611                            ntpStrip->ph1.sigcor) > 1600/ntpTrack->vtx.dcosz) ? 
01612           sqrt((ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor) - 
01613                (1600/ntpTrack->vtx.dcosz)) : 0;
01614       }
01615     }
01616   }
01617   return skimmedSqrtSC;  
01618 }
01619 
01620 //********************************************************
01621 //Needed for RecoShwEnergySqrt():
01622 Float_t MadQuantities::GetSqrtShwSC(Int_t ievt){
01623 
01624   Int_t shower = -1;
01625   if(!LoadLargestShowerFromEvent(ievt,shower)) return 0;
01626 
01627   Float_t summedSqrtSC = 0;
01628   TClonesArray* pointStripArray = NULL;
01629   if(isST) pointStripArray = (strecord->stp);
01630   else pointStripArray = (record->stp);
01631   TClonesArray& stripArray = *pointStripArray;
01632   Int_t *shwstrips = ntpShower->stp;
01633   for(Int_t j=0;j<ntpShower->nstrip;j++){
01634     Int_t index = shwstrips[j];
01635     ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
01636     summedSqrtSC += sqrt(ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor);
01637   }
01638   return summedSqrtSC;  
01639 }
01640 
01641 //********************************************************
01642 
01643 Int_t MadQuantities::GetChargeSign(Int_t itrk){  
01644   if(RecoMuQP(itrk)>0) return 1;
01645   return -1;
01646 }
01647 
01648 //********************************************************
01649 
01650 Float_t MadQuantities::RecoQENuEnergy(Int_t itrk){
01651   
01652   if(!LoadTrack(itrk)) return 0.;
01653   Float_t nucleonMass = 0.93956563; //mass of neutron by default
01654   if(GetChargeSign(itrk)==1) nucleonMass = 0.93827231; //proton mass for nubar
01655   Float_t muonEnergy = RecoMuEnergy(0,itrk);
01656   Float_t muonMass = 0.10555;
01657   Float_t muonMomentum = sqrt(muonEnergy*muonEnergy - muonMass*muonMass);
01658   Float_t costhbl = RecoMuDCosNeu(itrk);
01659 
01660   Float_t Eneu = (nucleonMass*muonEnergy - muonMass*muonMass/2.)
01661     /(nucleonMass - muonEnergy + muonMomentum*costhbl);  
01662 
01663   return Eneu;
01664 }
01665 
01666 //********************************************************
01667 
01668 Float_t MadQuantities::RecoQEQ2(Int_t itrk){
01669   
01670   if(!LoadTrack(itrk)) return 0.;
01671   Float_t Eneu = RecoQENuEnergy(itrk);
01672   Float_t muonEnergy = RecoMuEnergy(0,itrk);
01673   Float_t muonMass = 0.10555;
01674   Float_t muonMomentum = sqrt(muonEnergy*muonEnergy - muonMass*muonMass);
01675   Float_t costhbl = RecoMuDCosNeu(itrk);
01676   return -2.*Eneu*(muonEnergy-muonMomentum*costhbl)+muonMass*muonMass;
01677 
01678 }
01679 
01680 //*********************************************************
01681 Double_t MadQuantities::RecoEMFrac(Int_t shower,Int_t method,
01682                                    Double_t &EMFrac0,Double_t &EMFrac1)
01683 {
01684 
01685   if(!LoadShower(shower)) return 0;
01686   
01687   float EMFrac = 0;
01688   float EMFrac_end0 = 0;
01689   float EMFrac_end1 = 0;
01690   float EnTot = 0;
01691   float EnTot_end0 = 0;
01692   float EnTot_end1 = 0;
01693   for(int k=0; k<ntpShower->ncluster; k++){
01694     int index = ntpShower->clu[k];
01695     if(!LoadCluster(index)) continue;
01696     if(ntpCluster->id==0 && 
01697        (method==0 || ntpCluster->probem>0.2)){
01698       EMFrac+=ntpCluster->ph.gev;
01699       EnTot+=ntpCluster->ph.gev;
01700       if(ntpCluster->planeview==2){
01701         EMFrac_end0+=ntpCluster->ph.gev;
01702         EnTot_end0+=ntpCluster->ph.gev;
01703       }
01704       if(ntpCluster->planeview==3){
01705         EMFrac_end1+=ntpCluster->ph.gev;
01706         EnTot_end1+=ntpCluster->ph.gev;
01707       }
01708     }
01709     else {
01710       EnTot+=ntpCluster->ph.gev;
01711       if(ntpCluster->planeview==2){
01712         EnTot_end0+=ntpCluster->ph.gev;
01713       }
01714       if(ntpCluster->planeview==3){
01715         EnTot_end1+=ntpCluster->ph.gev;
01716       }
01717     }
01718   }
01719 
01720   if(EnTot==0) return 0;
01721   if(EnTot_end0==0) EMFrac = EMFrac_end1/EnTot_end1;
01722   else if(EnTot_end1==0) EMFrac = EMFrac_end0/EnTot_end0;
01723   else if(EMFrac_end0/EnTot_end0>EMFrac_end1/EnTot_end1) 
01724     EMFrac = EMFrac_end0/EnTot_end0;
01725   else EMFrac = EMFrac_end1/EnTot_end1;
01726   
01727   if(EnTot_end0>0) EMFrac0 = EMFrac_end0/EnTot_end0;
01728   if(EnTot_end1>1) EMFrac1 = EMFrac_end1/EnTot_end1;
01729   
01730   return EMFrac;
01731 
01732 }
01733 
01734 //-------------------------------------------------------------
01735 Int_t *MadQuantities::GetNShowerStrip(Int_t shower,Int_t opt)
01736 {
01737   Int_t *NShowerStrip = new Int_t[2];
01738   NShowerStrip[0] = 0; NShowerStrip[1] = 0;
01739   if(!LoadShower(shower)) {
01740     NShowerStrip[0] = -10; NShowerStrip[1] = -10;
01741     return NShowerStrip;
01742   }
01743   Int_t *clusters = ntpShower->clu;  
01744   for(int l=0;l<ntpShower->ncluster;l++){    
01745     Int_t index = clusters[l];
01746     if(!LoadCluster(index)) continue;
01747     if(opt==0||ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3) {
01748       if(ntpCluster->planeview==2) NShowerStrip[0] += ntpCluster->nstrip;
01749       else if(ntpCluster->planeview==3) NShowerStrip[1] += ntpCluster->nstrip;
01750     }
01751   }
01752   return NShowerStrip;
01753 }
01754 
01755 //******************************************************
01756 Float_t *MadQuantities::ClosestDistanceToEvent(Int_t evt,Float_t timeWeight)
01757 {
01758   Float_t *values = new Float_t[3];
01759   values[0] = values[1] = values[2] = -1;
01760   if(!LoadEvent(evt)) return values;
01761   Float_t smallestMetricValue = 0;
01762   Float_t smallestDeltaPos = 0;
01763   Float_t smallestDeltaTime = 0;
01764   for(int i=0;i<ntpEvent->nstrip;i++){
01765     if(!LoadStrip(ntpEvent->stp[i])) continue;
01766     Float_t z_val = ntpStrip->z;
01767     Float_t tpos_val = ntpStrip->tpos;
01768     Double_t time0_val = 1e9*ntpStrip->time0;
01769     Double_t time1_val = 1e9*ntpStrip->time1;    
01770     Double_t time_val = 0;
01771     if(time0_val>0 && time1_val>0) time_val = (time0_val + time1_val)/2.;
01772     else if(time0_val>0) time_val = time0_val;
01773     else if(time1_val>0) time_val = time1_val;
01774     Int_t planeview = ntpStrip->planeview;
01775     Float_t minDeltaStrip = 999999;
01776     Float_t minDeltaStripPos = 999999;
01777     Double_t minDeltaStripTime = 999999;
01778     for(UInt_t j=0;j<eventSummary->nstrip;j++){
01779       Bool_t useStrip = true;
01780       for(int k=0;k<ntpEvent->nstrip;k++){
01781         if(ntpEvent->stp[k]==Int_t(j)) {
01782           useStrip=false;
01783           break;
01784         }
01785       }
01786       if(!useStrip) continue;
01787       if(!LoadStrip(j)) continue;
01788       if(ntpStrip->planeview!=planeview) continue;
01789       Float_t deltaStripPos = TMath::Sqrt(TMath::Power(z_val - ntpStrip->z,2) +
01790                                           TMath::Power(tpos_val = ntpStrip->tpos,2));
01791       Double_t t0_val = 1e9*ntpStrip->time0; 
01792       Double_t t1_val = 1e9*ntpStrip->time1;
01793       Double_t t_val = 0;
01794       if(t0_val>0 && t1_val>0) t_val = (t0_val + t1_val)/2.;
01795       else if(t0_val>0) t_val = t0_val;
01796       else if(t1_val>0) t_val = t1_val;      
01797       Double_t deltaStripTime = TMath::Abs(time_val-t_val);
01798       Double_t metricValue = deltaStripPos + timeWeight*Float_t(deltaStripTime)*0.3;
01799       if(metricValue<minDeltaStrip) {
01800         minDeltaStrip = metricValue;
01801         minDeltaStripTime = deltaStripTime;
01802         minDeltaStripPos = deltaStripPos;
01803       }
01804     }
01805     smallestMetricValue += minDeltaStrip;
01806     smallestDeltaTime += Float_t(minDeltaStripTime);
01807     smallestDeltaPos += minDeltaStripPos;
01808   }
01809   values[0] = smallestMetricValue/Float_t(ntpEvent->nstrip);
01810   values[1] = smallestDeltaPos/Float_t(ntpEvent->nstrip);
01811   values[2] = smallestDeltaTime/Float_t(ntpEvent->nstrip);
01812   return values;
01813 }
01814 
01815 //*******************************************************
01816 void MadQuantities::MakeValidationFile(std::string tag){
01817 
01818   std::string savename = "RecoHistos_" + tag + ".root";
01819   TFile save(savename.c_str(),"RECREATE"); 
01820   save.cd();
01821   
01822   //#declare lots of histos:
01823   TH1F *TrueEnergy = new TH1F("TrueEnergy","True Neutrino Energy",200,0,50);
01824 
01825   TH1F *RecoEnergy = new TH1F("RecoEnergy","Reco Neutrino Energy",200,0,50);
01826 
01827   TH1F *TrueMuEn = new TH1F("TrueMuEn","True Muon Energy",200,0,50);
01828 
01829   TH1F *RecoMuEn = new TH1F("RecoMuEn","Reco Muon Energy",200,0,50);
01830   
01831   TH1F *TrueShwEn = new TH1F("TrueShwEn","True Shower Energy",200,0,50);
01832 
01833   TH1F *RecoShwEn = new TH1F("RecoShwEn","Reco Shower Energy",200,0,50);
01834   
01835   TH1F *TrueMuQP = new TH1F("TrueMuQP","True Muon q/p",300,-3,3);
01836 
01837   TH1F *RecoMuQP = new TH1F("RecoMuQP","Reconstructed Muon q/p",300,-3,3);
01838   
01839   TH2F *TrueVsReco = new TH2F("TrueVsReco",
01840                               "Reco Vs True Neutrino Energy",
01841                               200,0,50,200,0,50);
01842 
01843   TH2F *TrueVsRecoMu = new TH2F("TrueVsRecoMu",
01844                                 "Reco Vs True Muon Energy",200,0,50,200,0,50);
01845 
01846   TH2F *TrueVsRecoShw = new TH2F("TrueVsRecoShw",
01847                                  "Reco Vs True Shower Energy",
01848                                  200,0,50,200,0,50);
01849   
01850   TH2F *TrueRecoDiff 
01851     = new TH2F("TrueRecoDiff",
01852                "(True - Reco)/True Neutrino Energy vs True Neutrino Energy",
01853                200,0,50,1000,-9,1);
01854 
01855   TH2F *TrueRecoDiffMuCurv 
01856     = new TH2F("TrueRecoDiffMuCurv",
01857                "(True - Reco_{Curv})/True Muon Energy vs True Muon Energy",
01858                200,0,50,1000,-9,1);
01859   
01860   TH2F *TrueRecoDiffMuRange 
01861     = new TH2F("TrueRecoDiffMuRange",
01862                "(True - Reco_{Range})/True Muon Energy vs True Muon Energy",
01863                200,0,50,1000,-9,1);
01864 
01865   TH2F *TrueRecoDiffMuRangeHardCuts 
01866     = new TH2F("TrueRecoDiffMuRangeHardCuts",
01867                "(True - Reco_{Range})/True Muon Energy vs True Muon Energy with Hard Fiducial Cuts",200,0,50,1000,-9,1);
01868 
01869   TH2F *TrueRecoDiffShw 
01870     = new TH2F("TrueRecoDiffShw",
01871                "(True - Reco)/True Shower Energy vs True Shower Energy",
01872                200,0,50,1000,-9,1);
01873   
01874   TH2F *TrueRangeMu 
01875     = new TH2F("TrueRangeMu",
01876                "Reco Muon Energy from Range vs True Muon Energy",
01877                200,0,50,200,0,50); 
01878 
01879   TH2F *TrueCurvMu 
01880     = new TH2F("TrueCurvMu",
01881                "Reco Muon Energy from Curvature vs True Muon Energy",
01882                200,0,50,200,0,50); 
01883 
01884   TH2F *RecoRangeCurvMu 
01885     =new TH2F("RecoRangeCurvMu",
01886               "Reco Muon Energy from Range vs Reco Muon Energy from Curvature",
01887               200,0,50,200,0,50);
01888   
01889   TH2F *ShwSCVsTrueEn = new TH2F("ShwSCVsTrueEn",
01890                                  "Summed Shower SCs vs True Shower Energy",
01891                                  60,0,30,1000,0,100000);
01892 
01893   TH2F *ShwSCVsTrueEnZoom 
01894     = new TH2F("ShwSCVsTrueEnZoom",
01895                "Zoomed Summed Shower SCs vs True Shower Energy",
01896                50,0,5,1000,0,50000);
01897 
01898   TH1F *ShwSCperGeV = new TH1F("ShwSCperGeV",
01899                                "Summed Shower SCs per GeV Energy (from Truth)",
01900                                1000,0,10000);
01901 
01902   TH2F *RecoShwSCVsTrueShwEn 
01903     = new TH2F("RecoShwSCVsTrueShwEn",
01904                "Reconstructed Shower SCs Vs True Shower Energy",
01905                120,0,30,1000,0,100000);
01906     
01907   //Attempting to estimate neutrino energy in QE events using track angle:
01908   TH1F *TrueQENuEn = new TH1F("TrueQENuEn",
01909                               "True Neutrino Energy for QE-like Events",
01910                               200,0,50);
01911 
01912   TH1F *TrueNuEn_QEBackGround 
01913     = new TH1F("TrueNuEn_QEBackGround",
01914                "True Neutrino Energy for Background QE-like NC Events",
01915                200,0,50);
01916 
01917   TH1F *RecoQENuEn = new TH1F("RecoQENuEn","Reconstructed Neutrino Energy for QE-like Events using Muon Energy and Angle",200,0,50);
01918 
01919   TH1F *RecoQEMuEn = new TH1F("RecoQEMuEn",
01920                               "Reconstructed Muon Energy for QE-like Events",
01921                               200,0,50);
01922 
01923   TH2F *TrueRecoQEDiff 
01924     = new TH2F("TrueRecoQEDiff","(True - Reco_{QE})/True Neutrino Energy vs True Neutrino Energy for QE-like Events",200,0,50,1000,-9,1);
01925 
01926   TH2F *TrueRecoQEDiff_TrueQE 
01927     = new TH2F("TrueRecoQEDiff_TrueQE","(True - Reco_{QE})/True Neutrino Energy vs True Neutrino Energy for True QE Events",200,0,50,1000,-9,1);
01928 
01929   TH2F *TrueNuRecoMuDiff 
01930     = new TH2F("TrueNuRecoMuDiff","(TrueNu - RecoMu_{QE})/TrueNu Energy vs True Neutrino Energy for QE-like Events",200,0,50,1000,-9,1);
01931 
01932   TH2F *TrueNuRecoMuDiff_TrueQE 
01933     = new TH2F("TrueNuRecoMuDiff_TrueQE","(TrueNu - RecoMu_{QE})/TrueNu Energy vs True Neutrino Energy for True QE Events",200,0,50,1000,-9,1);
01934   
01935   //#Efficiency histos:  
01936   TH2F *MuEff = new TH2F("MuEff","MuEff",100,0,100,10,0,1);
01937   TH2F *MuEff_all = new TH2F("MuEff_all","MuEff_all",100,0,100,10,0,1);
01938   TH2F *MuEff_fid_all = new TH2F("MuEff_fid_all",
01939                                  "MuEff_fid_all",100,0,100,10,0,1);
01940 
01941   TH2F *NCMisEff = new TH2F("NCMisEff","NCMisEff",100,0,100,10,0,1);
01942   TH2F *NCMisEff_all = new TH2F("NCMisEff_all",
01943                                 "NCMisEff_all",100,0,100,10,0,1);
01944   TH2F *NCMisEff_fid_all = new TH2F("NCMisEff_fid_all",
01945                                     "NCMisEff_fid_all",100,0,100,10,0,1);
01946   
01947   //#Pi background histos:
01948   TH1F *TrkLenCC = new TH1F("TrkLenCC","Track length - NuMu CC",50,0,50);
01949   TH1F *TrkLenNC = new TH1F("TrkLenNC","Track length - NC",50,0,50);
01950   
01951   TH1F *dedxCC = new TH1F("dedxCC","Average dEdx - NuMu CC",1000,0,2000);
01952   TH1F *dedxNC = new TH1F("dedxNC","Average dEdx - NC",1000,0,2000);
01953 
01954   TH1F *HalfRatioCC 
01955     = new TH1F("HalfRatioCC",
01956                "Charge Ratio: First Half/Second Half of Track - NuMu CC",
01957                150,-1,14);
01958   TH1F *HalfRatioNC 
01959     = new TH1F("HalfRatioNC",
01960                "Charge Ratio: First Half/Second Half of Track - NuMu NC",
01961                150,-1,14);
01962   
01963   TH2F *MuRangeCurvDiffVsTrkLen 
01964       = new TH2F("MuRangeCurvDiffVsTrkLen",
01965                  "Diff in Momentum from Range and Curv Vs Track Length (CC)",
01966                  100,0,30,200,-3,17);
01967   TH2F *PiRangeCurvDiffVsTrkLen 
01968     = new TH2F("PiRangeCurvDiffVsTrkLen",
01969                "Diff in Momentum from Range and Curv Vs Track Length (NC)",
01970                100,0,30,200,-3,17);
01971   
01972   //#quasi plots
01973   TH2F *TrueShwSCVsQuasiNuEn
01974     = new TH2F("TrueShwSCVsQuasiNuEn",
01975                "True Shower SCs Vs Neutrino Energy for True QE events",
01976                60,0,30,500,0,10000);
01977   
01978   TH2F *TrueShwSCVsNonQENuEn
01979     = new TH2F("TrueShwSCVsNonQENuEn",
01980                "True Shower SCs Vs Neutrino Energy for True Non-QE events",
01981                60,0,30,500,0,10000);
01982   
01983   TH1F *TrueQENuEn_Pass 
01984     = new TH1F("TrueQENuEn_Pass",
01985                "True Neutrino Energy for QE Events that pass Reco Cuts",
01986                200,0,50);
01987   
01988   TH1F *SRMCVtxDistHist
01989     = new TH1F("SRMCVtxDistHist",
01990                "Distance between SR and MC Event Vertex",
01991                1000,0,5);
01992 
01993   TH1F *geantinoEnergy_CC[3];
01994   geantinoEnergy_CC[0] = new TH1F("geantinoEnergy_QECC",
01995                                   "Geantino Energy in Event",1000,0,5);
01996   geantinoEnergy_CC[1] = new TH1F("geantinoEnergy_RESCC",
01997                                   "Geantino Energy in Event",1000,0,5);
01998   geantinoEnergy_CC[2] = new TH1F("geantinoEnergy_DISCC",
01999                                   "Geantino Energy in Event",1000,0,5);
02000 
02001   TH1F *geantinoEnergy_NC[3];
02002   geantinoEnergy_NC[0] = new TH1F("geantinoEnergy_QENC",
02003                                   "Geantino Energy in Event",1000,0,5);
02004   geantinoEnergy_NC[1] = new TH1F("geantinoEnergy_RESNC",
02005                                   "Geantino Energy in Event",1000,0,5);
02006   geantinoEnergy_NC[2] = new TH1F("geantinoEnergy_DISNC",
02007                                   "Geantino Energy in Event",1000,0,5);
02008 
02009   TH1F *fracGeantinoEnergy_CC[3];
02010   fracGeantinoEnergy_CC[0] = new TH1F("fracGeantinoEnergy_QECC",
02011                                       "Fractional Geantino Energy in Event",
02012                                       1000,0,1);
02013   fracGeantinoEnergy_CC[1] = new TH1F("fracGeantinoEnergy_RESCC",
02014                                       "Fractional Geantino Energy in Event",
02015                                       1000,0,1);
02016   fracGeantinoEnergy_CC[2] = new TH1F("fracGeantinoEnergy_DISCC",
02017                                       "Fractional Geantino Energy in Event",
02018                                       1000,0,1);
02019 
02020   TH1F *fracGeantinoEnergy_NC[3];
02021   fracGeantinoEnergy_NC[0] = new TH1F("fracGeantinoEnergy_QENC",
02022                                       "Fractional Geantino Energy in Event",
02023                                       1000,0,1);
02024   fracGeantinoEnergy_NC[1] = new TH1F("fracGeantinoEnergy_RESNC",
02025                                       "Fractional Geantino Energy in Event",
02026                                       1000,0,1);
02027   fracGeantinoEnergy_NC[2] = new TH1F("fracGeantinoEnergy_DISNC",
02028                                       "Fractional Geantino Energy in Event",
02029                                       1000,0,1);
02030 
02031   TH1F *numGeantino_CC[3];
02032   numGeantino_CC[0] = new TH1F("numGeantino_QECC",
02033                                "Number of Geantino in Event",10,0,10);
02034   numGeantino_CC[1] = new TH1F("numGeantino_RESCC",
02035                                "Number of Geantino in Event",10,0,10);
02036   numGeantino_CC[2] = new TH1F("numGeantino_DISCC",
02037                                "Number of Geantino in Event",10,0,10);
02038 
02039   TH1F *numGeantino_NC[3];
02040   numGeantino_NC[0] = new TH1F("numGeantino_QENC",
02041                                "Number of Geantino in Event",10,0,10);
02042   numGeantino_NC[1] = new TH1F("numGeantino_RESNC",
02043                                "Number of Geantino in Event",10,0,10);
02044   numGeantino_NC[2] = new TH1F("numGeantino_DISNC",
02045                                "Number of Geantino in Event",10,0,10);
02046 
02047 
02048   for(int i=0;i<Nentries;i++){
02049     
02050     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
02051                             << "% done" << std::endl;
02052     
02053     if(!GetEntry(i)) continue;
02054 
02055     for(int j=0;j<eventSummary->nevent;j++){
02056       
02057       if(!LoadEvent(j)) continue;
02058       int mcevent = -1;
02059       if(!LoadTruthForRecoTH(j,mcevent)) continue;
02060 
02061       if(IAction(mcevent)==1) {
02062         Int_t index = IResonance(mcevent)-1001;
02063         if(index==3) index = 2;
02064         if(NumFSGeantino(mcevent)>0){
02065           geantinoEnergy_CC[index]->Fill(GeantinoEnergy(mcevent));
02066           fracGeantinoEnergy_CC[index]->Fill(GeantinoEnergy(mcevent) / 
02067                                              TrueNuEnergy(mcevent));
02068         }       
02069         numGeantino_CC[index]->Fill(NumFSGeantino(mcevent));
02070       }
02071       else {
02072         Int_t index = IResonance(mcevent)-1001;
02073         if(index==3) index = 2;
02074         if(NumFSGeantino(mcevent)>0){
02075           geantinoEnergy_NC[index]->Fill(GeantinoEnergy(mcevent));
02076           fracGeantinoEnergy_NC[index]->Fill(GeantinoEnergy(mcevent) / 
02077                                              TrueNuEnergy(mcevent));
02078         }
02079         numGeantino_NC[index]->Fill(NumFSGeantino(mcevent));
02080       }
02081       
02082       int track = -1;
02083       LoadLargestTrackFromEvent(j,track);
02084       int shower = -1;
02085       LoadLargestShowerFromEvent(j,shower);
02086       
02087       Float_t *vtx = this->TrueVtx(mcevent);
02088       float dist = sqrt(TMath::Power(vtx[0]-ntpEvent->vtx.x,2)+
02089                         TMath::Power(vtx[1]-ntpEvent->vtx.y,2)+
02090                         TMath::Power(vtx[2]-ntpEvent->vtx.z,2));
02091       SRMCVtxDistHist->Fill(dist);
02092       
02093       if(abs(this->INu(mcevent))==14&&this->IAction(mcevent)==1){
02094         
02095         MuEff_all->Fill(this->TrueNuEnergy(mcevent),this->Y(mcevent));
02096         
02097         if(sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1])<=3.5
02098            &&vtx[2]>=0.5){
02099           MuEff_fid_all->Fill(this->TrueNuEnergy(mcevent),this->Y(mcevent));
02100           if(this->PassCuts(track)) 
02101             MuEff->Fill(this->TrueNuEnergy(mcevent),this->Y(mcevent));
02102         }
02103         
02104         if(sqrt(this->W2(mcevent))<1.0) {
02105           TrueShwSCVsQuasiNuEn->Fill(this->TrueNuEnergy(mcevent),
02106                                      this->GetTrueShwSC(mcevent));
02107         }
02108         else TrueShwSCVsNonQENuEn->Fill(this->TrueNuEnergy(mcevent),
02109                                         this->GetTrueShwSC(mcevent));
02110         
02111         if(this->PassCuts(track)){
02112           //neutrino
02113           TrueEnergy->Fill(this->TrueNuEnergy(mcevent));
02114           RecoEnergy->Fill(this->RecoMuEnergy(1,track)
02115                            +this->RecoShwEnergy(shower));
02116           TrueVsReco->Fill(this->TrueNuEnergy(mcevent),
02117                            this->RecoMuEnergy(1,track) 
02118                            + this->RecoShwEnergy(shower));
02119           TrueRecoDiff->Fill(this->TrueNuEnergy(mcevent),
02120                              (this->TrueNuEnergy(mcevent)
02121                               -(this->RecoMuEnergy(1,track)
02122                                 +this->RecoShwEnergy(shower)))
02123                              /this->TrueNuEnergy(mcevent));
02124           
02125           //shower        
02126           TrueShwEn->Fill(this->TrueShwEnergy(mcevent));
02127           RecoShwEn->Fill(this->RecoShwEnergy(shower));
02128           TrueVsRecoShw->Fill(this->TrueShwEnergy(mcevent),
02129                               this->RecoShwEnergy(shower));
02130           TrueRecoDiffShw->Fill(this->TrueShwEnergy(mcevent),
02131                                 (this->TrueShwEnergy(mcevent)
02132                                  -this->RecoShwEnergy(shower))
02133                                 /this->TrueShwEnergy(mcevent));
02134           
02135           ShwSCVsTrueEn->Fill(this->TrueShwEnergy(mcevent),
02136                               this->GetTrueShwSC(mcevent));
02137           ShwSCVsTrueEnZoom->Fill(this->TrueShwEnergy(mcevent),
02138                                   this->GetTrueShwSC(mcevent));
02139           ShwSCperGeV->Fill(this->GetTrueShwSC(mcevent)
02140                             /this->TrueShwEnergy(mcevent));
02141           RecoShwSCVsTrueShwEn->Fill(this->TrueShwEnergy(mcevent),
02142                                      this->GetShwSC(track));
02143           
02144           //muon
02145           TrueMuEn->Fill(this->TrueMuEnergy(mcevent));
02146           RecoMuEn->Fill(this->RecoMuEnergy(0,track));
02147           TrueVsRecoMu->Fill(this->TrueMuEnergy(mcevent),
02148                              this->RecoMuEnergy(1,track));
02149             
02150           //curvature
02151           TrueRecoDiffMuCurv->Fill(this->TrueMuEnergy(mcevent),
02152                                    (this->TrueMuEnergy(mcevent)
02153                                     -this->RecoMuEnergy(1,track))
02154                                    /this->TrueMuEnergy(mcevent));
02155           TrueCurvMu->Fill(this->TrueMuEnergy(mcevent),
02156                            this->RecoMuEnergy(1,track));
02157           
02158           TrueMuQP->Fill(this->TrueMuQP(mcevent));
02159           RecoMuQP->Fill(this->RecoMuQP(track));
02160           
02161           //range
02162           if(this->IsFidAll(track)){
02163             TrueRecoDiffMuRange->Fill(this->TrueMuEnergy(mcevent),
02164                                       (this->TrueMuEnergy(mcevent)
02165                                        -this->RecoMuEnergy(2,track))
02166                                       /this->TrueMuEnergy(mcevent));      
02167               
02168             if(sqrt(TMath::Power(ntpTrack->vtx.x,2)
02169                     +TMath::Power(ntpTrack->vtx.y,2))>0.5
02170                &&sqrt(TMath::Power(ntpTrack->end.x,2)
02171                       +TMath::Power(ntpTrack->end.y,2))>0.5
02172                &&ntpTrack->end.z<25.) {
02173               TrueRecoDiffMuRangeHardCuts->Fill(this->TrueMuEnergy(mcevent),
02174                                                 (this->TrueMuEnergy(mcevent)
02175                                                  -this->RecoMuEnergy(2,track))
02176                                                 /this->TrueMuEnergy(mcevent));    
02177             }
02178             
02179             TrueRangeMu->Fill(this->TrueMuEnergy(mcevent),
02180                               this->RecoMuEnergy(2,track));
02181             
02182             RecoRangeCurvMu->Fill(this->RecoMuEnergy(1,track),
02183                                   this->RecoMuEnergy(2,track));
02184             MuRangeCurvDiffVsTrkLen->Fill(ntpTrack->ds,
02185                                           (this->RecoMuEnergy(1,track)
02186                                            -this->RecoMuEnergy(2,track))
02187                                           /this->RecoMuEnergy(1,track));
02188           }
02189           
02190           //QE-like
02191           if(this->PassQuasiCuts(track)){
02192             RecoQENuEn->Fill(this->RecoQENuEnergy(track));
02193             RecoQEMuEn->Fill(this->RecoMuEnergy(0,track));
02194             TrueQENuEn->Fill(this->TrueNuEnergy(mcevent));
02195             TrueRecoQEDiff->Fill(this->TrueNuEnergy(mcevent),
02196                                  (this->TrueNuEnergy(mcevent)
02197                                   -this->RecoQENuEnergy(track))
02198                                  /this->TrueNuEnergy(mcevent));
02199             TrueNuRecoMuDiff->Fill(this->TrueNuEnergy(mcevent),
02200                                    (this->TrueNuEnergy(mcevent)
02201                                     -this->RecoMuEnergy(0,track))
02202                                    /this->TrueNuEnergy(mcevent));
02203             if(this->W2(mcevent)<1.0) {
02204               TrueQENuEn_Pass->Fill(this->TrueNuEnergy(mcevent));
02205               TrueRecoQEDiff_TrueQE->Fill(this->TrueNuEnergy(mcevent),
02206                                           (this->TrueNuEnergy(mcevent)
02207                                            -this->RecoQENuEnergy(track))
02208                                           /this->TrueNuEnergy(mcevent));
02209               TrueNuRecoMuDiff_TrueQE->Fill(this->TrueNuEnergy(mcevent),
02210                                             (this->TrueNuEnergy(mcevent)
02211                                              -this->RecoMuEnergy(0,track))
02212                                             /this->TrueNuEnergy(mcevent));
02213             }
02214           }
02215           
02216           Float_t *CCNCVars = this->CCNCSepVars(track);
02217           TrkLenCC->Fill(CCNCVars[0]);
02218           dedxCC->Fill(CCNCVars[1]/100.);
02219           HalfRatioCC->Fill(CCNCVars[2]);
02220           delete CCNCVars;
02221         }
02222       }
02223       delete vtx;
02224       if(this->IAction(mcevent)==0){
02225         
02226         NCMisEff_all->Fill(this->TrueNuEnergy(mcevent),this->Y(mcevent));
02227         
02228         Float_t *vtx = this->TrueVtx(mcevent);
02229         if(sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1])<=3.5
02230            &&vtx[2]>=0.5){
02231           NCMisEff_fid_all->Fill(this->TrueNuEnergy(mcevent),this->Y(mcevent));
02232           if(this->PassCuts(0)) NCMisEff->Fill(this->TrueNuEnergy(mcevent),
02233                                                this->Y(mcevent));
02234         }
02235         delete vtx;
02236         
02237         if(this->PassCuts(track)) {
02238           Float_t *CCNCVars = this->CCNCSepVars(track);
02239           TrkLenNC->Fill(CCNCVars[0]);
02240           dedxNC->Fill(CCNCVars[1]/100.);
02241           HalfRatioNC->Fill(CCNCVars[2]);
02242           delete CCNCVars;
02243           if(this->IsFidAll(track)){
02244             PiRangeCurvDiffVsTrkLen->Fill(ntpTrack->ds,
02245                                           (this->RecoMuEnergy(1,track)
02246                                            -this->RecoMuEnergy(2,track))
02247                                           /this->RecoMuEnergy(1,track));
02248           }
02249           
02250           if(this->PassQuasiCuts(track)){
02251             TrueNuEn_QEBackGround->Fill(this->TrueNuEnergy(mcevent));
02252           }
02253         }
02254       }
02255     }
02256   }
02257 
02258   //#close file
02259   save.Write();
02260   save.Close();
02261  
02262 }
02263 
02264 void MadQuantities::ShowerValidation(char *fileName,Int_t startEnt)
02265 {
02266 
02267   this->GetEntry(0);
02268   const bool file_is_mc
02269     =(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC);
02270   Double_t pot = 0;
02271 
02272   char savename[256];
02273   sprintf(savename,"ShowerValidation_%s.root",fileName);
02274   TFile *save = new TFile(savename,"RECREATE");
02275   save->cd();
02276 
02277   //strip level:
02278   TH1F *stripPH = new TH1F("stripPH","stripPH",1000,0,100);
02279   TH1F *stripTime = new TH1F("stripTime","stripTime",1000,-1000,1000);
02280   TH1F *StripPHStripCut = new TH1F("StripPHStripCut","StripPHStripCut",1000,0,100);
02281   TH1F *StripTimeStripCut = new TH1F("StripTimeStripCut","StripTimeStripCut",1000,-1000,1000);
02282   TH1F *StripPHSSCut = new TH1F("StripPHSSCut","StripPHSSCut",1000,0,100);
02283   TH1F *StripTimeSSCut = new TH1F("StripTimeSSCut","StripTimeSSCut",1000,-1000,1000);
02284   TH1F *StripPHBothCut = new TH1F("StripPHBothCut","StripPHBothCut",1000,0,100);
02285   TH1F *StripTimeBothCut = new TH1F("StripTimeBothCut","StripTimeBothCut",1000,-1000,1000);
02286 
02287   TH1F *stripPH_Phys = new TH1F("stripPH_Phys","stripPH_Phys",1000,0,100);
02288   TH1F *stripTime_Phys = new TH1F("stripTime_Phys","stripTime_Phys",1000,-1000,1000);
02289   TH1F *StripPHStripCut_Phys = new TH1F("StripPHStripCut_Phys","StripPHStripCut_Phys",1000,0,100);
02290   TH1F *StripTimeStripCut_Phys = new TH1F("StripTimeStripCut_Phys","StripTimeStripCut_Phys",1000,-1000,1000);
02291   TH1F *StripPHSSCut_Phys = new TH1F("StripPHSSCut_Phys","StripPHSSCut_Phys",1000,0,100);
02292   TH1F *StripTimeSSCut_Phys = new TH1F("StripTimeSSCut_Phys","StripTimeSSCut_Phys",1000,-1000,1000);
02293   TH1F *StripPHBothCut_Phys = new TH1F("StripPHBothCut_Phys","StripPHBothCut_Phys",1000,0,100);
02294   TH1F *StripTimeBothCut_Phys = new TH1F("StripTimeBothCut_Phys","StripTimeBothCut_Phys",1000,-1000,1000);
02295 
02296   TH2F *stripTimeDiff = new TH2F("stripTimeDiff","stripTimeDiff",1000,0,1000,100,0,10);
02297 
02298   TH1F *multiplyHitStripFraction = new TH1F("multiplyHitStripFraction",
02299                                             "multiplyHitStripFraction",
02300                                             100,0,1.01);
02301   TH1F *multiplyHitStripFractionPhys = new TH1F("multiplyHitStripFractionPhys",
02302                                             "multiplyHitStripFractionPhys",
02303                                             100,0,1.01);
02304 
02305   //shower level:  
02306   TH1F *recoEshw = new TH1F("recoEshw","recoEshw",100,0,50);
02307   TH1F *recoEshwStripCut = new TH1F("recoEshwStripCut","recoEshwStripCut",100,0,50);
02308   TH1F *recoEshwSSCut = new TH1F("recoEshwSSCut","recoEshwSSCut",100,0,50);
02309   TH1F *recoEshwBothCut = new TH1F("recoEshwBothCut","recoEshwBothCut",100,0,50);
02310 
02311   TH1F *shwLengthTime = new TH1F("shwLengthTime","shwLengthTime",1000,0,1000);
02312   TH1F *shwLengthPlane = new TH1F("shwLengthPlane","shwLengthPlane",60,0,60);
02313   TH1F *shwStripN = new TH1F("shwStripN","shwStripN",200,0,200);
02314   TH1F *shwPhysStpN = new TH1F("shwPhysStpN","shwPhysStpN",200,0,200);
02315   TH1F *shwClusterN = new TH1F("shwClusterN","shwClusterN",60,0,60);
02316   TH1F *shwPhysCluN = new TH1F("shwPhysCluN","shwPhysCluN",60,0,60);
02317 
02318   //Truth Helper histos:
02319   TH1F *eventCompleteAll = new TH1F("eventCompleteAll",
02320                                     "Event Completeness",100,0,1);
02321   TH1F *eventCompleteSlc = new TH1F("eventCompleteSlc",
02322                                     "Event Completeness in Slice",100,0,1);
02323   TH1F *eventPurity = new TH1F("eventPurity","eventPurity",100,0,1);
02324   TH1F *eventEnergyRes = new TH1F("eventEnergyRes","eventEnergyRes",100,-1,9);
02325   TH1F *eventClosestDistance = new TH1F("eventClosestDistance","eventClosestDistance",
02326                                         1000,0,200);
02327   TH1F *eventLowResClosestDistance = new TH1F("eventLowResClosestDistance",
02328                                               "eventLowResClosestDistance",
02329                                               1000,0,200);
02330   TH2F *eventClosestDistTime = new TH2F("eventClosestDistTime","eventClosestDistTime",
02331                                         200,0,40,200,0,40);
02332   TH2F *eventLRClosestDistTime = new TH2F("eventLRClosestDistTime","eventLRClosestDistTime",
02333                                           200,0,40,200,0,40);
02334   TH2F *eventLRTrkShw = new TH2F("eventLRTrkShw","eventLRTrkShw",4,-0.5,3.5,4,-0.5,3.5);
02335   
02336 
02337   TH1F *trackCompleteAll = new TH1F("trackCompleteAll",
02338                                     "Track Completeness",100,0,1);
02339   TH1F *trackCompleteSlc = new TH1F("trackCompleteSlc",
02340                                     "Track Completeness in Slice",100,0,1);
02341   TH1F *trackPurity = new TH1F("trackPurity","trackPurity",100,0,1);
02342   TH1F *trackEnergyRes = new TH1F("trackEnergyRes","trackEnergyRes",100,-1,9);
02343 
02344   TH1F *showerCompleteAll = new TH1F("showerCompleteAll",
02345                                      "Shower Completeness",100,0,1);
02346   TH1F *showerCompleteSlc = new TH1F("showerCompleteSlc",
02347                                      "Shower Completeness in Slice",100,0,1);
02348   TH1F *showerPurity = new TH1F("showerPurity","showerPurity",100,0,1);
02349   TH1F *showerEnergyRes = new TH1F("showerEnergyRes","showerEnergyRes",100,-1,9);
02350   TH1F *showerVtxDist = new TH1F("showerVtxDist","showerVtxDist",100,0,6);
02351 
02352   //investigate low completeness showers:
02353   TH1F *showerLowCompleteProcess = new TH1F("showerLowCompleteProcess",
02354                                             "showerLowCompleteProcess",4,-0.5,3.5);
02355   TH1F *showerLowCompleteEnergy = new TH1F("showerLowCompleteEnergy",
02356                                            "showerLowCompleteEnergy",100,0,10);
02357   TH1F *showerLowCompleteCluID = new TH1F("showerLowCompleteCluID",
02358                                           "showerLowCompleteCluID",14,-0.5,13.5);
02359   TH1F *showerCompleteSlcSSOnly = new TH1F("showerCompleteSlcSSOnly",
02360                                           "showerCompleteSlcSSOnly",100,0,1);
02361   TH1F *showerLowCompleteVtxDist = new TH1F("showerLowCompleteVtxDist",
02362                                             "showerLowCompleteVtxDist",
02363                                             100,0,6);
02364   TH1F *showerLowCompleteProcessSSOnly = new TH1F("showerLowCompleteProcessSSOnly",
02365                                                   "showerLowCompleteProcessSSOnly",
02366                                                   4,-0.5,3.5);
02367 
02368   bool beamdb=false;
02369   bool checkeddb=false;
02370   for(int i=startEnt;i<Nentries;i++){
02371     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
02372                             << "% done" << std::endl;
02373     if(!this->GetEntry(i)) continue;
02374     if(!checkeddb){
02375       DbiResultPtr<BeamMonSpill> testbs(ntpHeader->GetVldContext(), 
02376                                         Dbi::kDefaultTask, 
02377                                         Dbi::kDisabled);
02378       checkeddb=true;
02379       if(testbs.GetNumRows()>0){
02380         beamdb=true;
02381         cout<<"Found the Beam Monitoring Database Tables"<<endl;
02382       }
02383       else{     
02384         cout<<"Didn't find the Beam Monitoring Database Tables"<<endl;
02385         cout<<"Resorting to old beamsummary ntuples"<<endl;
02386       }
02387     }
02388 
02389     if(!file_is_mc){ //file is mc
02390       Bool_t goodbeam = false;
02391       VldContext vc = ntpHeader->GetVldContext();
02392       VldTimeStamp vts = vc.GetTimeStamp();
02393       if(beamdb){ //beam db
02394         BDSpillAccessor &sa = BDSpillAccessor::Get();
02395         const BeamMonSpill *bs = sa.LoadSpill(vts);
02396         SpillTimeFinder &sfind= SpillTimeFinder::Instance();
02397         Double_t hpos2=0.;
02398         Double_t vpos2=0.;
02399         Double_t btint=0;
02400         for(int bt=0;bt<6;bt++){
02401           hpos2+=(bs->fTargBpmX[bt]*bs->fBpmInt[bt]);
02402           vpos2+=(bs->fTargBpmY[bt]*bs->fBpmInt[bt]);
02403           btint+=bs->fBpmInt[bt];
02404         }
02405         if(btint!=0){
02406           hpos2=hpos2*1000./btint;//make it in mm
02407           vpos2=vpos2*1000./btint;//make it in mm
02408         }
02409         else{
02410           hpos2=-999.;
02411           vpos2=-999.;
02412         }
02413         /*
02414           cout << bs->GetStatusBits().horn_on << " "
02415           << bs->GetStatusBits().target_in << " "
02416           << bs->fTortgt << " "
02417           << bs->fProfWidX*1000. << " "
02418           << bs->fProfWidY*1000. << " "
02419           << hpos2 << " " << vpos2 << " "
02420           << fabs(sfind.GetTimeToNearestSpill(vc)) << endl;
02421         */
02422         if(bs->GetStatusBits().horn_on && 
02423            bs->GetStatusBits().target_in &&
02424            bs->fTortgt > 0.1 &&
02425            bs->fProfWidX*1000.<1.5 &&
02426            bs->fProfWidY*1000.<2.0 &&
02427            hpos2>-2 && hpos2<0 &&    
02428            vpos2>0 && vpos2<2 &&
02429            fabs(sfind.GetTimeToNearestSpill(vc))<2 ) goodbeam=true;
02430         if(goodbeam) pot += bs->fTortgt;
02431       }
02432       if(!goodbeam) continue;
02433     }
02434 
02435     //Still have all good snarls here    
02436 
02437     if(eventSummary->nevent==0) continue;
02438 
02439     //find and store earliest time of hit strip
02440     Double_t stripArrayTime[500][192][2];
02441     for(int ii=0;ii<500;ii++) { 
02442       for(int jj=0;jj<192;jj++) {
02443         stripArrayTime[ii][jj][0] = stripArrayTime[ii][jj][1] = -1.0;
02444       }
02445     }
02446     for(int j=0;j<int(eventSummary->nstrip);j++){
02447       if(!LoadStrip(j)) continue;
02448       Int_t pl = ntpStrip->plane;
02449       Int_t st = ntpStrip->strip;
02450       if(ntpStrip->time1<stripArrayTime[pl][st][1] || 
02451          stripArrayTime[pl][st][1]<=0) {
02452         stripArrayTime[pl][st][1] = ntpStrip->time1;
02453       }
02454     }
02455 
02456     Int_t is_fid = 0; 
02457     for(int j=0;j<eventSummary->nevent;j++){ 
02458       if(!LoadEvent(j)) continue; //no event found
02459 
02460       int track_index   = -1;
02461       int shower_index  = -1;
02462       Bool_t vertex_shower = true;
02463       if(LoadLargestTrackFromEvent(j,track_index)) {
02464         if(!LoadShowerAtTrackVertex(j,track_index,shower_index)) {        
02465           LoadLargestShowerFromEvent(j,shower_index);
02466           vertex_shower = false;
02467         }
02468       }
02469       else LoadLargestShowerFromEvent(j,shower_index);
02470 
02471       is_fid = 1;
02472       if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
02473       if(ntpEvent->ntrack>0&&!IsFidVtxTrk(track_index)) is_fid = 0;
02474       if(is_fid==0) continue;
02475 
02476       //Still have all fiducial events from good snarls here
02477       Float_t phwCluID = 0;
02478       Float_t phClu = 0;
02479       if(ntpEvent->nshower>0) { //is shower?
02480 
02481         recoEshw->Fill(ntpShower->shwph.linCCgev);
02482         shwLengthPlane->Fill(ntpShower->plane.n);
02483         Double_t shwstptimemax = 0;
02484         Double_t shwstptimemin = 0;
02485         for(int k=0;k<ntpShower->nstrip;k++){
02486           if(!LoadStrip(ntpShower->stp[k])) continue;
02487           Double_t t0 = ntpStrip->time0;
02488           Double_t t1 = ntpStrip->time1;
02489           Double_t avet = 0;
02490           if(t0>0&&t1>0) avet = (t0+t1)/2.;
02491           else if(t0>0) avet = t0;
02492           else if(t1>0) avet = t1;
02493           if(k==0) shwstptimemin = shwstptimemax = avet;
02494           else {
02495             if(avet>shwstptimemax) shwstptimemax = avet;
02496             if(avet<shwstptimemin) shwstptimemin = avet;
02497           }
02498         }
02499         shwLengthTime->Fill(1e9*(shwstptimemax-shwstptimemin));
02500 
02501         Int_t nMultiplyHitStrips = 0;
02502         Int_t nPhysClusterU = 0;
02503         Int_t nPhysClusterV = 0;        
02504         Int_t *clusters = ntpShower->clu;
02505         for(int k=0; k<ntpShower->ncluster; k++){
02506           Int_t index = clusters[k];
02507           if(!LoadCluster(index)) continue;
02508           phwCluID += ntpCluster->id*ntpCluster->ph.gev;
02509           phClu += ntpCluster->ph.gev;
02510           for(int l=0;l<ntpCluster->nstrip;l++){
02511             if(!LoadStrip(ntpCluster->stp[l])) continue;
02512             Int_t pl = ntpStrip->plane;
02513             Int_t st = ntpStrip->strip;
02514             if(ntpStrip->time1>stripArrayTime[pl][st][1]) {       
02515               stripTimeDiff->Fill((ntpStrip->time1 - 
02516                                    stripArrayTime[pl][st][1])*1e9,
02517                                   ntpStrip->ph1.pe);
02518               nMultiplyHitStrips+=1;
02519             }
02520           }
02521           if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02522             if(ntpCluster->planeview==2) nPhysClusterU+=1;
02523             if(ntpCluster->planeview==3) nPhysClusterV+=1;
02524           }
02525         }
02526 
02527         if(phClu>0) phwCluID/=phClu;
02528 
02529         shwClusterN->Fill(ntpShower->nUcluster+ntpShower->nVcluster);
02530         shwPhysCluN->Fill(nPhysClusterU+nPhysClusterV);
02531         
02532         Int_t *nss = GetNShowerStrip(shower_index);
02533         Int_t *npss = GetNShowerStrip(shower_index,1);
02534     
02535         shwStripN->Fill(nss[0]+nss[1]);
02536         shwPhysStpN->Fill(npss[0]+npss[1]);
02537 
02538         multiplyHitStripFraction->Fill(Float_t(nMultiplyHitStrips)/
02539                                        Float_t(ntpShower->nstrip));
02540         
02541         if(Float_t(nMultiplyHitStrips)/Float_t(ntpShower->nstrip)<0.3) {
02542           recoEshwStripCut->Fill(ntpShower->shwph.linCCgev);
02543           if(1.-float(npss[0]+npss[1])/float(nss[0]+nss[1]) < 0.85) {
02544             recoEshwBothCut->Fill(ntpShower->shwph.linCCgev);
02545           }
02546         }
02547         if(1.-float(npss[0]+npss[1])/float(nss[0]+nss[1]) < 0.85) {
02548           recoEshwSSCut->Fill(ntpShower->shwph.linCCgev);
02549         }
02550         
02551         multiplyHitStripFractionPhys->Fill(Float_t(nMultiplyHitStrips)/
02552                                            Float_t(ntpShower->nstrip));
02553         
02554         
02555         for(int k=0; k<ntpShower->ncluster; k++){
02556           Int_t index = clusters[k];
02557           if(!LoadCluster(index)) continue;
02558           for(int l=0;l<ntpCluster->nstrip;l++){
02559             if(!LoadStrip(ntpCluster->stp[l])) continue;
02560             stripPH->Fill(ntpStrip->ph1.pe);
02561             stripTime->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));       
02562             if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02563               stripPH_Phys->Fill(ntpStrip->ph1.pe);
02564               stripTime_Phys->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));
02565             }
02566             
02567             if(Float_t(nMultiplyHitStrips)/Float_t(ntpShower->nstrip)<0.3) {
02568               StripPHStripCut->Fill(ntpStrip->ph1.pe);
02569               StripTimeStripCut->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));               
02570               if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02571                 StripPHStripCut_Phys->Fill(ntpStrip->ph1.pe);
02572                 StripTimeStripCut_Phys->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));
02573               }
02574             }
02575             
02576             if(1.-float(npss[0]+npss[1])/float(nss[0]+nss[1]) < 0.85) {
02577               StripPHSSCut->Fill(ntpStrip->ph1.pe);
02578               StripTimeSSCut->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));                  
02579               if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02580                 StripPHSSCut_Phys->Fill(ntpStrip->ph1.pe);
02581                 StripTimeSSCut_Phys->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));
02582               }
02583               
02584               if(Float_t(nMultiplyHitStrips)/Float_t(ntpShower->nstrip)<0.3) {
02585                 StripPHBothCut->Fill(ntpStrip->ph1.pe);
02586                 StripTimeBothCut->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));              
02587                 if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02588                   StripPHBothCut_Phys->Fill(ntpStrip->ph1.pe);
02589                   StripTimeBothCut_Phys->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));
02590                 }
02591               }
02592             }
02593           }
02594         }
02595         delete [] nss;
02596         delete [] npss;
02597       }
02598       //if this is MC:
02599       Int_t mcevent = -1;
02600       if(file_is_mc&&LoadTruthForRecoTH(j,mcevent)) {
02601         if(LoadTHEvent(j)) {
02602           eventCompleteAll->Fill(ntpTHEvent->completeall);
02603           eventCompleteSlc->Fill(ntpTHEvent->completeslc);
02604           eventPurity->Fill(ntpTHEvent->purity);
02605           Float_t true_en_to_use = ntpTruth->p4neu[3];
02606           if(ntpTruth->iaction==0) true_en_to_use = ntpTruth->p4shw[3];
02607           Float_t reco_en_to_use = RecoMuEnergy(track_index,0) + 
02608             RecoShwEnergy(shower_index,0);
02609           reco_en_to_use = ntpEvent->ph.gev;
02610           eventEnergyRes->Fill((reco_en_to_use - 
02611                                 true_en_to_use)/true_en_to_use);
02612           Float_t *closestDistanceToOtherEvent = 
02613             ClosestDistanceToEvent(ntpEvent->index,0.033); 
02614           eventClosestDistance->Fill(closestDistanceToOtherEvent[0]);
02615           eventClosestDistTime->Fill(closestDistanceToOtherEvent[2],
02616                                      closestDistanceToOtherEvent[1]);
02617           if((reco_en_to_use - true_en_to_use)/true_en_to_use<-0.9) {
02618             eventLowResClosestDistance->Fill(closestDistanceToOtherEvent[0]);
02619             eventLRClosestDistTime->Fill(closestDistanceToOtherEvent[2],
02620                                          closestDistanceToOtherEvent[1]);
02621             eventLRTrkShw->Fill(ntpEvent->ntrack,ntpEvent->nshower);
02622           }
02623           delete [] closestDistanceToOtherEvent;
02624         }
02625         if(LoadTHTrack(track_index)){
02626           trackCompleteAll->Fill(ntpTHTrack->completeall);
02627           trackCompleteSlc->Fill(ntpTHTrack->completeslc);
02628           trackPurity->Fill(ntpTHTrack->purity);
02629           Float_t true_en_to_use = ntpTruth->p4mu1[3];
02630           Float_t reco_en_to_use = RecoMuEnergy(track_index,0);
02631           if(true_en_to_use > 0)
02632             trackEnergyRes->Fill((reco_en_to_use - 
02633                                   true_en_to_use)/true_en_to_use);
02634           else trackEnergyRes->Fill(-1);
02635         }
02636         if(LoadTHShower(shower_index)&&vertex_shower){
02637           showerCompleteAll->Fill(ntpTHShower->completeall);
02638           showerCompleteSlc->Fill(ntpTHShower->completeslc);
02639           showerPurity->Fill(ntpTHShower->purity);
02640           Float_t true_en_to_use = ntpTruth->p4neu[3]*ntpTruth->y;
02641           if(true_en_to_use == 0) true_en_to_use = ntpTruth->p4shw[3];
02642           Float_t reco_en_to_use = RecoShwEnergy(shower_index,0);
02643           showerEnergyRes->Fill((reco_en_to_use - 
02644                                  true_en_to_use)/true_en_to_use);
02645           if(ntpTrack) showerVtxDist->Fill(TMath::Abs(ntpTrack->vtx.z - 
02646                                                       ntpShower->vtx.z));
02647           if(phwCluID>=5) {
02648             showerCompleteSlcSSOnly->Fill(ntpTHShower->completeslc);
02649             if(ntpTHShower->completeslc<0.1){
02650               showerLowCompleteProcessSSOnly->Fill((ntpTruth->iresonance-1000) * 
02651                                                    ntpTruth->iaction);
02652             }
02653           }
02654           if(ntpTHShower->completeslc<0.1){
02655             showerLowCompleteProcess->Fill((ntpTruth->iresonance-1000) * 
02656                                            ntpTruth->iaction);
02657             showerLowCompleteEnergy->Fill(ntpTruth->p4shw[3]);
02658             showerLowCompleteCluID->Fill(phwCluID);
02659             if(ntpTrack) showerLowCompleteVtxDist->Fill(TMath::Abs(ntpTrack->vtx.z - 
02660                                                                    ntpShower->vtx.z));
02661           }
02662         }
02663       }
02664     }
02665   }
02666 
02667   recoEshw->Write();
02668   recoEshwSSCut->Write();
02669   recoEshwStripCut->Write();
02670   recoEshwBothCut->Write();
02671 
02672   stripPH->Write();
02673   stripTime->Write();  
02674   stripPH_Phys->Write();
02675   stripTime_Phys->Write();  
02676 
02677   StripPHStripCut->Write();
02678   StripTimeStripCut->Write();
02679   StripPHStripCut_Phys->Write();
02680   StripTimeStripCut_Phys->Write();
02681 
02682   StripPHSSCut->Write();
02683   StripTimeSSCut->Write();
02684   StripPHSSCut_Phys->Write();
02685   StripTimeSSCut_Phys->Write();
02686 
02687   StripPHBothCut->Write();
02688   StripTimeBothCut->Write();
02689   StripPHBothCut_Phys->Write();
02690   StripTimeBothCut_Phys->Write();
02691 
02692   stripTimeDiff->Write();
02693   multiplyHitStripFraction->Write();
02694   multiplyHitStripFractionPhys->Write();
02695 
02696   shwLengthTime->Write();
02697   shwLengthPlane->Write();
02698   shwStripN->Write();
02699   shwPhysStpN->Write();
02700   shwClusterN->Write();
02701   shwPhysCluN->Write();
02702 
02703   eventCompleteAll->Write();
02704   eventCompleteSlc->Write();
02705   eventPurity->Write();
02706   eventEnergyRes->Write();
02707   eventClosestDistance->Write();
02708   eventLowResClosestDistance->Write();
02709   eventClosestDistTime->Write();
02710   eventLRClosestDistTime->Write();
02711   eventLRTrkShw->Write();
02712 
02713   trackCompleteAll->Write();
02714   trackCompleteSlc->Write();
02715   trackPurity->Write();
02716   trackEnergyRes->Write();
02717 
02718   showerCompleteAll->Write();
02719   showerCompleteSlc->Write();
02720   showerPurity->Write();
02721   showerEnergyRes->Write();
02722   showerVtxDist->Write();
02723 
02724   showerCompleteSlcSSOnly->Write();
02725   showerLowCompleteProcess->Write();
02726   showerLowCompleteEnergy->Write();
02727   showerLowCompleteCluID->Write();
02728   showerLowCompleteVtxDist->Write();
02729   showerLowCompleteProcessSSOnly->Write();
02730 
02731   save->Close();
02732   cout << "POT = " << pot << endl;
02733   return;
02734 }
02735 
02736 #endif // #ifdef madquantities_cxx