MadCluAnalysis Class Reference

#include <MadCluAnalysis.h>

Inheritance diagram for MadCluAnalysis:
MadAnalysis MadQuantities MadBase

List of all members.

Public Member Functions

 MadCluAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
 MadCluAnalysis (JobC *, string, int)
 ~MadCluAnalysis ()
void Plot (char *fileName, Int_t startEnt=0)
void Plotter (char *fileName)
Double_t FOM ()
void DataDistributions (char *filename)
void DataMCComp (char *, char *)
Bool_t PassDataCleaningCuts (Int_t shower)
Bool_t IsTrackInSlice (Int_t slice)
Bool_t IsNueFidVtxEvt (Int_t event)
Double_t * PHWProbEM (Int_t event)
Double_t * PHWCluID (Int_t event)
Double_t * PHWAvgDev (Int_t event)
Double_t * ChargeFracRMS (Int_t event, Int_t method=0)
Double_t * AveNumSSPlane (Int_t event)
Bool_t FillPHWProbEM (Int_t event)
Bool_t FillPHWCluID (Int_t event)
Bool_t FillPHWAvgDev (Int_t event)
Bool_t FillChargeFracRMS (Int_t event, Int_t method=0)
Bool_t FillAveNumSSPlane (Int_t event)
Bool_t FillNCluster (Int_t event)
Bool_t FillBestProbEM (Int_t event)
Bool_t ReadPIDFile (char *filename)
Bool_t ReadNNFile (char *filename)

Public Attributes

Double_t PAN_PHWCluIDU
Double_t PAN_PHWCluIDV
Double_t PAN_PHWCluIDBest
Double_t PAN_PHWProbEMU
Double_t PAN_PHWProbEMV
Double_t PAN_PHWProbEMBest
Double_t PAN_PHWAvgDevU
Double_t PAN_PHWAvgDevV
Double_t PAN_PHWAvgDevBest
Double_t PAN_ChargeFracRMSU
Double_t PAN_ChargeFracRMSV
Double_t PAN_ChargeFracRMSBest
Double_t PAN_AveNumSSPlaneU
Double_t PAN_AveNumSSPlaneV
Double_t PAN_AveNumSSPlaneBest
Int_t PAN_NShowerStripU
Int_t PAN_NShowerStripV
Int_t PAN_NShowerStripBest
Int_t PAN_NPhysShowerStripU
Int_t PAN_NPhysShowerStripV
Double_t PAN_BestProbEM
Double_t PAN_RatioProbEM
Double_t PAN_Weight
Float_t fPID3
Float_t fPID4
Float_t fPID5
Int_t fis_nue_pid
Int_t PAN_NClusterU
Int_t PAN_NClusterV
Int_t PAN_NPhysClusterU
Int_t PAN_NPhysClusterV
TH1F * flikelihoodDists [12]
TMultiLayerPerceptron * fneuralNet

Protected Member Functions

Bool_t PassTruthSignal (Int_t mcevent=0)
Bool_t PassAnalysisCuts (Int_t event=0)
Bool_t PassBasicCuts (Int_t event=0)
Float_t PID (Int_t event=0, Int_t method=0)
Float_t RecoAnalysisEnergy (Int_t event=0)
Float_t GetWeight (Int_t event=0)
Bool_t AddUserBranches (TTree *)
void CallUserFunctions (Int_t)

Detailed Description

Definition at line 8 of file MadCluAnalysis.h.


Constructor & Destructor Documentation

MadCluAnalysis::MadCluAnalysis ( TChain *  chainSR = 0,
TChain *  chainMC = 0,
TChain *  chainTH = 0,
TChain *  chainEM = 0 
)

Definition at line 36 of file MadCluAnalysis.cxx.

References MadBase::Clear(), MadBase::emrecord, fis_nue_pid, fneuralNet, fPID3, fPID4, fPID5, MadBase::InitChain(), MadBase::isEM, MadBase::isMC, MadBase::isTH, MadBase::mcrecord, MadBase::Nentries, PAN_AveNumSSPlaneBest, PAN_AveNumSSPlaneU, PAN_AveNumSSPlaneV, PAN_BestProbEM, PAN_ChargeFracRMSBest, PAN_ChargeFracRMSU, PAN_ChargeFracRMSV, PAN_NClusterU, PAN_NClusterV, PAN_NPhysClusterU, PAN_NPhysClusterV, PAN_NPhysShowerStripU, PAN_NPhysShowerStripV, PAN_NShowerStripBest, PAN_NShowerStripU, PAN_NShowerStripV, PAN_PHWAvgDevBest, PAN_PHWAvgDevU, PAN_PHWAvgDevV, PAN_PHWCluIDBest, PAN_PHWCluIDU, PAN_PHWCluIDV, PAN_PHWProbEMBest, PAN_PHWProbEMU, PAN_PHWProbEMV, PAN_RatioProbEM, PAN_Weight, MadBase::record, MadBase::threcord, and MadBase::whichSource.

00037 {
00038 
00039   if(!chainSR) {
00040     record = 0;
00041     mcrecord = 0;
00042     threcord = 0;
00043     emrecord = 0;
00044     Clear();
00045     whichSource = -1;
00046     isMC = true;
00047     isTH = true;
00048     isEM = true;
00049     Nentries = -1;
00050     return;
00051   }
00052   
00053   InitChain(chainSR,chainMC,chainTH,chainEM);
00054   whichSource = 0;
00055 
00056   PAN_PHWCluIDU = 0;
00057   PAN_PHWCluIDV = 0;
00058   PAN_PHWCluIDBest = 0;
00059 
00060   PAN_PHWProbEMU = 0;
00061   PAN_PHWProbEMV = 0;
00062   PAN_PHWProbEMBest = 0;
00063 
00064   PAN_PHWAvgDevU = 0;
00065   PAN_PHWAvgDevV = 0;
00066   PAN_PHWAvgDevBest = 0;
00067 
00068   PAN_ChargeFracRMSU = 0;
00069   PAN_ChargeFracRMSV = 0;
00070   PAN_ChargeFracRMSBest = 0;
00071 
00072   PAN_AveNumSSPlaneU = 0;
00073   PAN_AveNumSSPlaneV = 0;
00074   PAN_AveNumSSPlaneBest = 0;
00075 
00076   PAN_NShowerStripU = 0;
00077   PAN_NShowerStripV = 0;
00078   PAN_NShowerStripBest = 0;
00079 
00080   PAN_NPhysShowerStripU = 0;
00081   PAN_NPhysShowerStripV = 0;
00082 
00083   PAN_BestProbEM = 0;
00084   PAN_RatioProbEM = 0;
00085   PAN_Weight = 0;
00086   fPID3 = 0;
00087   fPID4 = 0;
00088   fPID5 = 0;
00089   fis_nue_pid = 0;
00090   
00091   PAN_NClusterU = 0;
00092   PAN_NClusterV = 0;
00093   PAN_NPhysClusterU = 0;
00094   PAN_NPhysClusterV = 0;
00095 
00096   fneuralNet = NULL;
00097 
00098 }

MadCluAnalysis::MadCluAnalysis ( JobC j,
string  path,
int  entries 
)

Definition at line 100 of file MadCluAnalysis.cxx.

References MadBase::Clear(), MadBase::emrecord, MadBase::fChain, fis_nue_pid, MadBase::fJC, fneuralNet, fPID3, fPID4, fPID5, MadBase::isEM, MadBase::isMC, MadBase::isTH, MadBase::jcPath, MadBase::mcrecord, MadBase::Nentries, PAN_AveNumSSPlaneBest, PAN_AveNumSSPlaneU, PAN_AveNumSSPlaneV, PAN_BestProbEM, PAN_ChargeFracRMSBest, PAN_ChargeFracRMSU, PAN_ChargeFracRMSV, PAN_NClusterU, PAN_NClusterV, PAN_NPhysClusterU, PAN_NPhysClusterV, PAN_NPhysShowerStripU, PAN_NPhysShowerStripV, PAN_NShowerStripBest, PAN_NShowerStripU, PAN_NShowerStripV, PAN_PHWAvgDevBest, PAN_PHWAvgDevU, PAN_PHWAvgDevV, PAN_PHWCluIDBest, PAN_PHWCluIDU, PAN_PHWCluIDV, PAN_PHWProbEMBest, PAN_PHWProbEMU, PAN_PHWProbEMV, PAN_RatioProbEM, PAN_Weight, MadBase::record, MadBase::threcord, and MadBase::whichSource.

00101 {
00102 
00103   record = 0;
00104   mcrecord = 0;
00105   threcord = 0;
00106   emrecord = 0;
00107   Clear();
00108   isMC = true;
00109   isTH = true;
00110   isEM = true;
00111   Nentries = entries;
00112   whichSource = 1;
00113   jcPath = path;
00114   fJC = j;
00115   fChain = NULL;
00116 
00117   PAN_PHWCluIDU = 0;
00118   PAN_PHWCluIDV = 0;
00119   PAN_PHWCluIDBest = 0;
00120 
00121   PAN_PHWProbEMU = 0;
00122   PAN_PHWProbEMV = 0;
00123   PAN_PHWProbEMBest = 0;
00124 
00125   PAN_PHWAvgDevU = 0;
00126   PAN_PHWAvgDevV = 0;
00127   PAN_PHWAvgDevBest = 0;
00128 
00129   PAN_ChargeFracRMSU = 0;
00130   PAN_ChargeFracRMSV = 0;
00131   PAN_ChargeFracRMSBest = 0;
00132 
00133   PAN_AveNumSSPlaneU = 0;
00134   PAN_AveNumSSPlaneV = 0;
00135   PAN_AveNumSSPlaneBest = 0;
00136 
00137   PAN_NShowerStripU = 0;
00138   PAN_NShowerStripV = 0;
00139   PAN_NShowerStripBest = 0;
00140 
00141   PAN_NPhysShowerStripU = 0;
00142   PAN_NPhysShowerStripV = 0;
00143 
00144   PAN_BestProbEM = 0;
00145   PAN_RatioProbEM = 0;
00146   PAN_Weight = 0;
00147   fPID3 = 0;
00148   fPID4 = 0;
00149   fPID5 = 0;
00150   fis_nue_pid = 0;
00151 
00152   PAN_NClusterU = 0;
00153   PAN_NClusterV = 0;
00154   PAN_NPhysClusterU = 0;
00155   PAN_NPhysClusterV = 0;
00156 
00157   fneuralNet = NULL;
00158 
00159 }

MadCluAnalysis::~MadCluAnalysis (  ) 

Definition at line 162 of file MadCluAnalysis.cxx.

00163 {
00164   
00165 }


Member Function Documentation

Bool_t MadCluAnalysis::AddUserBranches ( TTree *  tree  )  [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 2712 of file MadCluAnalysis.cxx.

References fis_nue_pid, fPID3, fPID4, fPID5, PAN_AveNumSSPlaneBest, PAN_AveNumSSPlaneU, PAN_AveNumSSPlaneV, PAN_BestProbEM, PAN_ChargeFracRMSBest, PAN_ChargeFracRMSU, PAN_ChargeFracRMSV, PAN_NClusterU, PAN_NClusterV, PAN_NPhysClusterU, PAN_NPhysClusterV, PAN_NPhysShowerStripU, PAN_NPhysShowerStripV, PAN_NShowerStripBest, PAN_NShowerStripU, PAN_NShowerStripV, PAN_PHWAvgDevBest, PAN_PHWAvgDevU, PAN_PHWAvgDevV, PAN_PHWCluIDBest, PAN_PHWCluIDU, PAN_PHWCluIDV, PAN_PHWProbEMBest, PAN_PHWProbEMU, PAN_PHWProbEMV, PAN_RatioProbEM, and PAN_Weight.

02713 {
02714   if(!tree) return false;
02715   tree->Branch("ss_PhwIDU",&PAN_PHWCluIDU,"ss_PhwIDU/D",32000);
02716   tree->Branch("ss_PhwIDV",&PAN_PHWCluIDV,"ss_PhwIDV/D",32000);
02717   tree->Branch("ss_PhwIDBest",&PAN_PHWCluIDBest,"ss_PhwIDBest/D",32000);
02718 
02719   tree->Branch("ss_PhwProbEMU",&PAN_PHWProbEMU,"ss_PhwProbEMU/D",32000);
02720   tree->Branch("ss_PhwProbEMV",&PAN_PHWProbEMV,"ss_PhwProbEMV/D",32000);
02721   tree->Branch("ss_PhwProbEMBest",&PAN_PHWProbEMBest,"ss_PhwProbEMBest/D",32000);
02722 
02723   tree->Branch("ss_PHWAvgDevU",&PAN_PHWAvgDevU,"ss_PHWAvgDevU/D",32000);
02724   tree->Branch("ss_PHWAvgDevV",&PAN_PHWAvgDevV,"ss_PHWAvgDevV/D",32000);
02725   tree->Branch("ss_PHWAvgDevBest",&PAN_PHWAvgDevBest,"ss_PHWAvgDevBest/D",32000);
02726 
02727   tree->Branch("ss_ChargeFracRMSU",&PAN_ChargeFracRMSU,"ss_ChargeFracRMSU/D",32000);
02728   tree->Branch("ss_ChargeFracRMSV",&PAN_ChargeFracRMSV,"ss_ChargeFracRMSV/D",32000);
02729   tree->Branch("ss_ChargeFracRMSBest",&PAN_ChargeFracRMSBest,"ss_ChargeFracRMSBest/D",32000);
02730 
02731   tree->Branch("ss_AveNumSSPlaneU",&PAN_AveNumSSPlaneU,"ss_AveNumSSPlaneU/D",32000);
02732   tree->Branch("ss_AveNumSSPlaneV",&PAN_AveNumSSPlaneV,"ss_AveNumSSPlaneV/D",32000);
02733   tree->Branch("ss_AveNumSSPlaneBest",&PAN_AveNumSSPlaneBest,"ss_AveNumSSPlaneBest/D",32000);
02734 
02735   tree->Branch("ss_NShowerStripU",&PAN_NShowerStripU,"ss_NShowerStripU/I",32000);
02736   tree->Branch("ss_NShowerStripV",&PAN_NShowerStripV,"ss_NShowerStripV/I",32000);
02737   tree->Branch("ss_NShowerStripBest",&PAN_NShowerStripBest,"ss_NShowerStripBest/I",32000);
02738 
02739   tree->Branch("ss_NPhysShowerStripU",&PAN_NPhysShowerStripU,"ss_NPhysShowerStripU/I",32000);
02740   tree->Branch("ss_NPhysShowerStripV",&PAN_NPhysShowerStripV,"ss_NPhysShowerStripV/I",32000);
02741 
02742   tree->Branch("ss_BestProbEM",&PAN_BestProbEM,"ss_BestProbEM/D",32000);
02743   tree->Branch("ss_RatioProbEM",&PAN_RatioProbEM,"ss_RatioProbEM/D",32000);
02744   tree->Branch("ss_Weight",&PAN_Weight,"ss_Weight/D",32000);
02745 
02746   tree->Branch("pid3",&fPID3,"pid3/F",32000);
02747   tree->Branch("pid4",&fPID4,"pid4/F",32000);
02748   tree->Branch("pid5",&fPID5,"pid5/F",32000);
02749   tree->Branch("is_nue_fid",&fis_nue_pid,"is_nue_fid/I",32000);
02750   tree->Branch("NClusterU",&PAN_NClusterU,"NClusterU/I",32000);
02751   tree->Branch("NClusterV",&PAN_NClusterV,"NClusterV/I",32000);
02752   tree->Branch("NPhysClusterU",&PAN_NPhysClusterU,"NPhysClusterU/I",32000);
02753   tree->Branch("NPhysClusterV",&PAN_NPhysClusterV,"NPhysClusterV/I",32000);
02754   return true;
02755 }

Double_t * MadCluAnalysis::AveNumSSPlane ( Int_t  event  ) 

Definition at line 2664 of file MadCluAnalysis.cxx.

References NtpSRPlane::beg, NtpSRShower::clu, NtpSRCluster::id, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadStrip(), Munits::m, NtpSRShower::ncluster, NtpSRCluster::nstrip, MadBase::ntpCluster, MadBase::ntpShower, MadBase::ntpStrip, NtpSRPlane::nu, NtpSRPlane::nv, NtpSRShower::plane, NtpSRStrip::plane, NtpSRCluster::planeview, and NtpSRCluster::stp.

Referenced by FillAveNumSSPlane().

02665 {
02666   Double_t *aveNumSSPlane = new Double_t[2];
02667   aveNumSSPlane[0] = 0; aveNumSSPlane[1] = 0;
02668   if(!LoadEvent(event)) {
02669     aveNumSSPlane[0] = -10.; aveNumSSPlane[1] = -10.;
02670     return aveNumSSPlane;
02671   }
02672   Int_t shower = -1;
02673   if(!LoadLargestShowerFromEvent(event,shower)) {
02674     aveNumSSPlane[0] = -10.; aveNumSSPlane[1] = -10.;
02675     return aveNumSSPlane;
02676   }  
02677   Int_t *clusters = ntpShower->clu;  
02678   for(int k=ntpShower->plane.beg; k<=ntpShower->plane.end; k++){
02679     Double_t planeAveU = 0;
02680     Double_t planeAveV = 0;
02681     for(int l=0;l<ntpShower->ncluster;l++){    
02682       Int_t index = clusters[l];
02683       if(!LoadCluster(index)) continue;
02684       if(!(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3)) continue;
02685       Int_t *strip = ntpCluster->stp;
02686       for(int m=0;m<ntpCluster->nstrip;m++){
02687         if(!LoadStrip(strip[m])) continue;
02688         if(ntpStrip->plane==k){
02689           if(ntpCluster->planeview==2) planeAveU += 1;
02690           else if(ntpCluster->planeview==3) planeAveV += 1;
02691           break;
02692         }
02693       }
02694     }
02695     aveNumSSPlane[0] += planeAveU;
02696     aveNumSSPlane[1] += planeAveV;
02697   }
02698   aveNumSSPlane[0]/=ntpShower->plane.nu;
02699   aveNumSSPlane[1]/=ntpShower->plane.nv;
02700   return aveNumSSPlane;
02701 }

void MadCluAnalysis::CallUserFunctions ( Int_t  event  )  [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 2757 of file MadCluAnalysis.cxx.

References FillAveNumSSPlane(), FillBestProbEM(), FillChargeFracRMS(), FillNCluster(), FillPHWAvgDev(), FillPHWCluID(), FillPHWProbEM(), fis_nue_pid, fPID3, fPID4, fPID5, MadQuantities::GetNShowerStrip(), NtpMCTruth::iaction, NtpSREvent::index, NtpMCTruth::inu, NtpMCTruth::inunoosc, IsNueFidVtxEvt(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadTruthForRecoTH(), max, min, MadBase::ntpEvent, MadBase::ntpTruth, PAN_AveNumSSPlaneBest, PAN_AveNumSSPlaneU, PAN_AveNumSSPlaneV, PAN_BestProbEM, PAN_ChargeFracRMSBest, PAN_ChargeFracRMSU, PAN_ChargeFracRMSV, PAN_NClusterU, PAN_NClusterV, PAN_NPhysClusterU, PAN_NPhysClusterV, PAN_NPhysShowerStripU, PAN_NPhysShowerStripV, PAN_NShowerStripBest, PAN_NShowerStripU, PAN_NShowerStripV, PAN_PHWAvgDevBest, PAN_PHWAvgDevU, PAN_PHWAvgDevV, PAN_PHWCluIDBest, PAN_PHWCluIDU, PAN_PHWCluIDV, PAN_PHWProbEMBest, PAN_PHWProbEMU, PAN_PHWProbEMV, PAN_RatioProbEM, PAN_Weight, and PID().

02758 {
02759 
02760   PAN_PHWCluIDU = 0;
02761   PAN_PHWCluIDV = 0;
02762   PAN_PHWCluIDBest = 0;
02763   PAN_PHWProbEMU = 0;
02764   PAN_PHWProbEMV = 0;
02765   PAN_PHWProbEMBest = 0;
02766   PAN_PHWAvgDevU = 0;
02767   PAN_PHWAvgDevV = 0;
02768   PAN_PHWAvgDevBest = 0;
02769   PAN_ChargeFracRMSU = 0;
02770   PAN_ChargeFracRMSV = 0;
02771   PAN_ChargeFracRMSBest = 0;
02772   PAN_AveNumSSPlaneU = 0;
02773   PAN_AveNumSSPlaneV = 0;
02774   PAN_AveNumSSPlaneBest = 0;
02775   PAN_NShowerStripU = 0;
02776   PAN_NShowerStripV = 0;
02777   PAN_NShowerStripBest = 0;
02778   PAN_NPhysShowerStripU = 0;
02779   PAN_NPhysShowerStripV = 0;
02780   PAN_BestProbEM = 0;
02781   PAN_RatioProbEM = 0;
02782   PAN_Weight = 0;
02783   fPID3 = 0;
02784   fPID4 = 0;
02785   fPID5 = 0;
02786   fis_nue_pid = 0;
02787   PAN_NClusterU = 0;
02788   PAN_NClusterV = 0;
02789   PAN_NPhysClusterU = 0;
02790   PAN_NPhysClusterV = 0;
02791 
02792   if(!LoadEvent(event)) return;
02793 
02794   fis_nue_pid = this->IsNueFidVtxEvt(event);
02795 
02796   this->FillNCluster(event);
02797   this->FillBestProbEM(event);
02798   this->FillPHWCluID(event);
02799   this->FillPHWProbEM(event);
02800   this->FillPHWAvgDev(event);
02801   this->FillChargeFracRMS(event,1);
02802   this->FillAveNumSSPlane(event);
02803   
02804   //set ratio
02805   if(PAN_PHWProbEMU>0 && PAN_PHWProbEMV>0) {
02806     if(PAN_PHWProbEMU>PAN_PHWProbEMV) 
02807       PAN_RatioProbEM = PAN_PHWProbEMV/PAN_PHWProbEMU;
02808     else PAN_RatioProbEM = PAN_PHWProbEMU/PAN_PHWProbEMV;
02809   }
02810 
02811   Int_t shower = -1;
02812   LoadLargestShowerFromEvent(event,shower);
02813 
02814   Int_t *nshowerstrip = this->GetNShowerStrip(shower);
02815   PAN_NShowerStripU = nshowerstrip[0];
02816   PAN_NShowerStripV = nshowerstrip[1];
02817   delete [] nshowerstrip;
02818  
02819   Int_t *nphysshowerstrip = this->GetNShowerStrip(shower,1);
02820   PAN_NPhysShowerStripU = nphysshowerstrip[0];
02821   PAN_NPhysShowerStripV = nphysshowerstrip[1];
02822   delete [] nphysshowerstrip;
02823 
02824   double prob = 1.;
02825   Int_t mcevent = -1;
02826   if(LoadTruthForRecoTH(ntpEvent->index,mcevent)) {
02827     Double_t nNuMuFiles = 40;
02828     Double_t nNuEFiles = 39;
02829     Double_t nNuTauFiles = 39;
02830     Double_t POT = 7.4;
02831     //prob = Oscillate(ntpTruth->inu,ntpTruth->inunoosc,
02832     //       ntpTruth->p4neu[3],735.,0.002175,
02833     //       TMath::ASin(TMath::Sqrt(0.925))/2.,
02834     //       0.01);
02835     //prob = Oscillate(ntpTruth->inu,ntpTruth->inunoosc,
02836     //       ntpTruth->p4neu[3],735.,0.0025,
02837     //       TMath::ASin(TMath::Sqrt(1.0))/2.,
02838     //       0.01);      
02839 
02840     if(ntpTruth->iaction==1){
02841       if(abs(ntpTruth->inunoosc)==12 && abs(ntpTruth->inu)==12)
02842         PAN_Weight = prob*(POT/((nNuEFiles+nNuMuFiles)*6.5));
02843       else if(abs(ntpTruth->inu)==12) 
02844         PAN_Weight = prob*(POT/(nNuEFiles*6.5));
02845       else if(abs(ntpTruth->inu)==14)
02846         PAN_Weight = prob*(POT/(nNuMuFiles*6.5));
02847       else if(abs(ntpTruth->inu)==16)
02848         PAN_Weight = prob*(POT/(nNuTauFiles*6.5));
02849     }
02850     else if(ntpTruth->iaction==0) {
02851       if(abs(ntpTruth->inu)==12)
02852         PAN_Weight = prob*(POT/((nNuEFiles)*6.5));      
02853       else if(abs(ntpTruth->inu)==14)
02854         PAN_Weight = prob*(POT/((nNuMuFiles)*6.5));
02855       else if(abs(ntpTruth->inu)==16)
02856         PAN_Weight = prob*(POT/((nNuTauFiles)*6.5));
02857     }
02858   }
02859   else PAN_Weight = 1.;
02860 
02861   fPID3 = PID(event,3);
02862   fPID4 = PID(event,4);
02863   fPID5 = PID(event,5);
02864 
02865   PAN_ChargeFracRMSBest = max(PAN_ChargeFracRMSU,PAN_ChargeFracRMSV);
02866   PAN_PHWProbEMBest     = max(PAN_PHWProbEMU,PAN_PHWProbEMV);
02867   PAN_PHWCluIDBest      = min(PAN_PHWCluIDU,PAN_PHWCluIDV);
02868   PAN_NShowerStripBest  = min(PAN_NShowerStripU,PAN_NShowerStripV);
02869   PAN_AveNumSSPlaneBest = min(PAN_AveNumSSPlaneU,PAN_AveNumSSPlaneV);
02870   PAN_PHWAvgDevBest     = min(fabs(PAN_PHWAvgDevU-0.02),
02871                               fabs(PAN_PHWAvgDevV-0.02));
02872 
02873 }

Double_t * MadCluAnalysis::ChargeFracRMS ( Int_t  event,
Int_t  method = 0 
)

Definition at line 2610 of file MadCluAnalysis.cxx.

References NtpSRShower::clu, NtpSRStripPulseHeight::gev, NtpSRCluster::id, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), NtpSRShower::ncluster, MadBase::ntpCluster, MadBase::ntpShower, NtpSRShower::nUcluster, NtpSRShower::nVcluster, NtpSRCluster::ph, and NtpSRCluster::planeview.

Referenced by FillChargeFracRMS().

02611 {
02612   Double_t *chargeFracRMS = new Double_t[2];
02613   chargeFracRMS[0] = 0; chargeFracRMS[1] = 0;
02614   
02615   if(!LoadEvent(event)) {
02616     chargeFracRMS[0] = -10.; chargeFracRMS[1] = -10.;
02617     return chargeFracRMS;
02618   }
02619   Int_t shower = -1;
02620   if(!LoadLargestShowerFromEvent(event,shower)) {
02621     chargeFracRMS[0] = -10.; chargeFracRMS[1] = -10.;
02622     return chargeFracRMS;
02623   }
02624   Double_t totE_U = 0;
02625   Double_t totE_V = 0;
02626   Int_t *clusters = ntpShower->clu;
02627   for(int k=0; k<ntpShower->ncluster; k++){
02628     Int_t index = clusters[k];
02629     if(!LoadCluster(index)) continue;
02630     if(method==1 && !(ntpCluster->id==0 || 
02631                       ntpCluster->id==1 || 
02632                       ntpCluster->id==3)) continue;
02633     if(ntpCluster->planeview==2){
02634       chargeFracRMS[0] += ntpCluster->ph.gev*ntpCluster->ph.gev;
02635       totE_U += ntpCluster->ph.gev;
02636     }
02637     else if(ntpCluster->planeview==3){
02638       chargeFracRMS[1] += ntpCluster->ph.gev*ntpCluster->ph.gev;
02639       totE_V += ntpCluster->ph.gev;
02640     }
02641   }
02642   if(totE_U!=0) {
02643     chargeFracRMS[0] /= (totE_U*totE_U*ntpShower->nUcluster);
02644     chargeFracRMS[0] -= TMath::Power(1./ntpShower->nUcluster,2);
02645     chargeFracRMS[0]  = TMath::Sqrt(chargeFracRMS[0]);
02646   }
02647   if(totE_V!=0) {
02648     chargeFracRMS[1] /= (totE_V*totE_V*ntpShower->nVcluster);
02649     chargeFracRMS[1] -= TMath::Power(1./ntpShower->nVcluster,2);
02650     chargeFracRMS[1]  = TMath::Sqrt(chargeFracRMS[1]);
02651   }
02652   return chargeFracRMS;
02653 }

void MadCluAnalysis::DataDistributions ( char *  filename  ) 

Definition at line 866 of file MadCluAnalysis.cxx.

References NtpSRShower::clu, NtpSREventSummary::date, MadBase::eventSummary, VldContext::GetDetector(), MadBase::GetEntry(), SpillInfoBlock::GetHorn(), SpillInfoBlock::GetHpos(), SpillInfoBlock::GetMagnet(), VldTimeStamp::GetNanoSec(), SpillInfoBlock::GetPot(), VldTimeStamp::GetSec(), VldContext::GetSimFlag(), SpillInfo::GetSpillInfo(), SpillInfoBlock::GetTgtpos(), VldContext::GetTimeStamp(), RecHeader::GetVldContext(), SpillInfoBlock::GetVpos(), NtpSRStripPulseHeight::gev, NtpMCTruth::iaction, NtpSRCluster::id, NtpSREvent::index, NtpMCTruth::inu, NtpMCTruth::inunoosc, MadQuantities::IsFidVtxEvt(), Detector::kFar, Detector::kNear, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadTruthForRecoTH(), month, NtpSRDate::month, NtpSRShower::ncluster, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, MadBase::ntpCluster, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpTruth, PassDataCleaningCuts(), MadQuantities::PassTrackCuts(), NtpSRShower::ph, NtpSRCluster::ph, NtpSREvent::ph, NtpSRCluster::planeview, NtpSRCluster::probem, MadQuantities::RecoEMFrac(), MadQuantities::RecoMuEnergy(), NtpSRDate::year, and Munits::year.

00866                                                     {
00867 
00868   SpillInfo pppinfo;
00869   Double_t datapot = 0;
00870 
00871   char savename[256];
00872   sprintf(savename,"DataClu_%s.root",filename);
00873   TFile *save = new TFile(savename,"RECREATE");
00874   save->cd();
00875 
00876   string SplitName[6] = {"beamnue","nue","numu","nutau","nc","unknown"};
00877 
00878   TH1F *IDHistTrack = new TH1F("IDHistTrack","Cluster ID for Track Events",6,0,5);
00879   IDHistTrack->SetXTitle("Cluster ID");
00880   TH1F *IDHistShower = new TH1F("IDHistShower","Cluster ID for Shower Events",6,0,5);
00881   IDHistShower->SetXTitle("Cluster ID");
00882 
00883   TH1F *IDHistTrackSplit[6] = {0};
00884   for(int i=0;i<6;i++) {
00885     char name[256];
00886     sprintf(name,"%s_%s",IDHistTrack->GetName(),SplitName[i].c_str());
00887     IDHistTrackSplit[i] = new TH1F(name,"Cluster ID for Track Events",6,0,5);
00888   }
00889   TH1F *IDHistShowerSplit[6] = {0};
00890   for(int i=0;i<6;i++) {
00891     char name[256];
00892     sprintf(name,"%s_%s",IDHistShower->GetName(),SplitName[i].c_str());
00893     IDHistShowerSplit[i] = new TH1F(name,"Cluster ID for Shower Events",6,0,5);
00894   }
00895 
00896   TH1F *RecoyTrack = new TH1F("RecoyTrack","Reco y Distribution for Track Events",
00897                               101,-0.005,1.005);
00898   RecoyTrack->SetXTitle("Reco y");
00899   TH1F *RecoyShower = new TH1F("RecoyShower","Reco y Distribution for Shower Events",
00900                                101,-0.005,1.005);
00901   RecoyShower->SetXTitle("Reco y");
00902 
00903   TH1F *RecoyTrackSplit[6] = {0};
00904   for(int i=0;i<6;i++) {
00905     char name[256];
00906     sprintf(name,"%s_%s",RecoyTrack->GetName(),SplitName[i].c_str());
00907     RecoyTrackSplit[i] = new TH1F(name,"Reco y Distribution for Track Events",
00908                                101,-0.005,1.005);
00909   }
00910   TH1F *RecoyShowerSplit[6] = {0};
00911   for(int i=0;i<6;i++) {
00912     char name[256];
00913     sprintf(name,"%s_%s",RecoyShower->GetName(),SplitName[i].c_str());
00914     RecoyShowerSplit[i] = new TH1F(name,"Reco y Distribution for Shower Events",
00915                                 101,-0.005,1.005);
00916   }
00917 
00918   TH1F *RecoEMFracTrack = new TH1F("RecoEMFracTrack","Reco EMFrac for Track Events",
00919                                    101,-0.005,1.005);
00920   RecoEMFracTrack->SetXTitle("Reco EMFrac");
00921   TH1F *RecoEMFracShower = new TH1F("RecoEMFracShower","Reco EMFrac for Shower Events",
00922                                     101,-0.005,1.005);
00923   RecoEMFracShower->SetXTitle("Reco EMFrac");
00924 
00925   TH1F *RecoEMFracTrackSplit[6] = {0};
00926   for(int i=0;i<6;i++) {
00927     char name[256];
00928     sprintf(name,"%s_%s",RecoEMFracTrack->GetName(),SplitName[i].c_str());
00929     RecoEMFracTrackSplit[i] = new TH1F(name,"Reco EMFrac for Track Events",
00930                                     101,-0.005,1.005);
00931   }
00932   TH1F *RecoEMFracShowerSplit[6] = {0};
00933   for(int i=0;i<6;i++) {
00934     char name[256];
00935     sprintf(name,"%s_%s",RecoEMFracShower->GetName(),SplitName[i].c_str());
00936     RecoEMFracShowerSplit[i] = new TH1F(name,"Reco EMFrac for Shower Events",
00937                                     101,-0.005,1.005);
00938   }
00939 
00940   TH1F *RecoEMFracProbTrack = new TH1F("RecoEMFracProbTrack","Reco EMFrac for Track Events",
00941                                        101,-0.005,1.005);
00942   RecoEMFracProbTrack->SetXTitle("Reco EMFrac");
00943   TH1F *RecoEMFracProbShower = new TH1F("RecoEMFracProbShower","Reco EMFrac for Shower Events",
00944                                         101,-0.005,1.005);
00945   RecoEMFracProbShower->SetXTitle("Reco EMFrac");
00946 
00947   TH1F *RecoEMFracProbTrackSplit[6] = {0};
00948   for(int i=0;i<6;i++) {
00949     char name[256];
00950     sprintf(name,"%s_%s",RecoEMFracProbTrack->GetName(),SplitName[i].c_str());
00951     RecoEMFracProbTrackSplit[i] = new TH1F(name,"Reco EMFrac for Track Events",
00952                                         101,-0.005,1.005);
00953   }
00954   TH1F *RecoEMFracProbShowerSplit[6] = {0};
00955   for(int i=0;i<6;i++) {
00956     char name[256];
00957     sprintf(name,"%s_%s",RecoEMFracProbShower->GetName(),SplitName[i].c_str());
00958     RecoEMFracProbShowerSplit[i] = new TH1F(name,"Reco EMFrac for Shower Events",
00959                                          101,-0.005,1.005);
00960   }
00961 
00962   //primary cluster IDs for NC and CC
00963   TH1F *UClusterIDTrack = new TH1F("UClusterIDTrack","U View Primary Cluster ID for Track Events",6,-0.5,5.5);
00964   UClusterIDTrack->SetXTitle("Cluster ID");
00965   TH1F *UClusterIDShower = new TH1F("UClusterIDShower","U View Primary Cluster ID for Shower Events",6,-0.5,5.5);
00966   UClusterIDShower->SetXTitle("Cluster ID");
00967   TH1F *VClusterIDTrack = new TH1F("VClusterIDTrack","V View Primary Cluster ID for Track Events",6,-0.5,5.5);
00968   VClusterIDTrack->SetXTitle("Cluster ID");
00969   TH1F *VClusterIDShower = new TH1F("VClusterIDShower","V View Primary Cluster ID for Shower Events",6,-0.5,5.5);
00970   VClusterIDShower->SetXTitle("Cluster ID");
00971 
00972   TH1F *UClusterIDTrackSplit[6] = {0};
00973   TH1F *VClusterIDTrackSplit[6] = {0};
00974   for(int i=0;i<6;i++) {
00975     char name[256];
00976     sprintf(name,"%s_%s",UClusterIDTrack->GetName(),SplitName[i].c_str());
00977     UClusterIDTrackSplit[i] = new TH1F(name,"U View Primary Cluster ID for Track Events",6,-0.5,5.5);
00978     sprintf(name,"%s_%s",VClusterIDTrack->GetName(),SplitName[i].c_str());
00979     VClusterIDTrackSplit[i] = new TH1F(name,"V View Primary Cluster ID for Track Events",6,-0.5,5.5);
00980   }
00981   TH1F *UClusterIDShowerSplit[6] = {0};
00982   TH1F *VClusterIDShowerSplit[6] = {0};
00983   for(int i=0;i<6;i++) {
00984     char name[256];
00985     sprintf(name,"%s_%s",UClusterIDShower->GetName(),SplitName[i].c_str());
00986     UClusterIDShowerSplit[i] = new TH1F(name,"U View Primary Cluster ID for Shower Events",6,-0.5,5.5);
00987     sprintf(name,"%s_%s",VClusterIDShower->GetName(),SplitName[i].c_str());
00988     VClusterIDShowerSplit[i] = new TH1F(name,"V View Primary Cluster ID for Shower Events",6,-0.5,5.5);
00989   }
00990 
00991     //primary cluster ProbEMs
00992   TH1F *UClusterProbEMTrack = new TH1F("UClusterProbEMTrack","U View Primary Cluster ProbEM for Track Events",100,0,1);
00993   UClusterProbEMTrack->SetXTitle("Cluster ProbEM");
00994   TH1F *UClusterProbEMShower = new TH1F("UClusterProbEMShower","U View Primary Cluster ProbEM for Shower Events",100,0,1);
00995   UClusterProbEMShower->SetXTitle("Cluster ProbEM");
00996   TH1F *VClusterProbEMTrack = new TH1F("VClusterProbEMTrack","V View Primary Cluster ProbEM for Track Events",100,0,1);
00997   VClusterProbEMTrack->SetXTitle("Cluster ProbEM");
00998   TH1F *VClusterProbEMShower = new TH1F("VClusterProbEMShower","V View Primary Cluster ProbEM for Shower Events",100,0,1);
00999   VClusterProbEMShower->SetXTitle("Cluster ProbEM");
01000 
01001   TH1F *UClusterProbEMTrackSplit[6] = {0};
01002   TH1F *VClusterProbEMTrackSplit[6] = {0};
01003   for(int i=0;i<6;i++) {
01004     char name[256];
01005     sprintf(name,"%s_%s",UClusterProbEMTrack->GetName(),SplitName[i].c_str());
01006     UClusterProbEMTrackSplit[i] = new TH1F(name,"U View Primary Cluster ProbEM for Track Events",100,0,1);
01007     sprintf(name,"%s_%s",VClusterProbEMTrack->GetName(),SplitName[i].c_str());
01008     VClusterProbEMTrackSplit[i] = new TH1F(name,"V View Primary Cluster ProbEM for Track Events",100,0,1);
01009   }
01010   TH1F *UClusterProbEMShowerSplit[6] = {0};
01011   TH1F *VClusterProbEMShowerSplit[6] = {0};
01012   for(int i=0;i<6;i++) {
01013     char name[256];
01014     sprintf(name,"%s_%s",UClusterProbEMShower->GetName(),SplitName[i].c_str());
01015     UClusterProbEMShowerSplit[i] = new TH1F(name,"U View Primary Cluster ProbEM for Shower Events",100,0,1);
01016     sprintf(name,"%s_%s",VClusterProbEMShower->GetName(),SplitName[i].c_str());
01017     VClusterProbEMShowerSplit[i] = new TH1F(name,"V View Primary Cluster ProbEM for Shower Events",100,0,1);
01018   }
01019 
01020 
01021   //Average #cluster with energy
01022   TH2F *AvgCluEnTrack = new TH2F("AvgCluEnTrack",
01023                                  "Average Number of Clusters with Reco Shower Energy",
01024                                  100,0,25,20,0,20);
01025   AvgCluEnTrack->SetXTitle("Reco Shower Energy (GeV)");
01026   AvgCluEnTrack->SetYTitle("# of Clusters");
01027   TH2F *AvgCluEnShower = new TH2F("AvgCluEnShower",
01028                                   "Average Number of Clusters with Reco Shower Energy",
01029                                   100,0,25,20,0,20);
01030   AvgCluEnShower->SetXTitle("Reco Shower Energy (GeV)");
01031   AvgCluEnShower->SetYTitle("# of Clusters");
01032   TH2F *AvgCluEnTracknox = 
01033     new TH2F("AvgCluEnTracknox",
01034              "Average Number of Clusters with Reco Shower Energy (excl xtalk clu) - CC",
01035              100,0,25,20,0,20);
01036   AvgCluEnTracknox->SetXTitle("Reco Shower Energy (GeV)");
01037   AvgCluEnTracknox->SetYTitle("# of Clusters");
01038   TH2F *AvgCluEnShowernox = 
01039     new TH2F("AvgCluEnShowernox",
01040              "Average Number of Clusters with Reco Shower Energy (excl xtalk clu) - NC",
01041              100,0,25,20,0,20);
01042   AvgCluEnShowernox->SetXTitle("Reco Shower Energy (GeV)");
01043   AvgCluEnShowernox->SetYTitle("# of Clusters");
01044 
01045   TH2F *AvgCluEnTrackSplit[6] = {0};
01046   TH2F *AvgCluEnnoxTrackSplit[6] = {0};
01047   for(int i=0;i<6;i++) {
01048     char name[256];
01049     sprintf(name,"%s_%s",AvgCluEnTrack->GetName(),SplitName[i].c_str());
01050     AvgCluEnTrackSplit[i] = new TH2F(name,"Average Number of Clusters with Reco Shower Energy",
01051                                   100,0,25,20,0,20);
01052     sprintf(name,"%s_%s",AvgCluEnTracknox->GetName(),SplitName[i].c_str());
01053     AvgCluEnnoxTrackSplit[i] = new TH2F(name,
01054                     "Average Number of Clusters with Reco Shower Energy (excl xtalk clu)",
01055                                      100,0,25,20,0,20);
01056   }
01057   TH2F *AvgCluEnShowerSplit[6] = {0};
01058   TH2F *AvgCluEnnoxShowerSplit[6] = {0};
01059   for(int i=0;i<6;i++) {
01060     char name[256];
01061     sprintf(name,"%s_%s",AvgCluEnShower->GetName(),SplitName[i].c_str());
01062     AvgCluEnShowerSplit[i] = new TH2F(name,"Average Number of Clusters with Reco Shower Energy",
01063                                    100,0,25,20,0,20);
01064     sprintf(name,"%s_%s",AvgCluEnShowernox->GetName(),SplitName[i].c_str());
01065     AvgCluEnnoxShowerSplit[i] = new TH2F(name,
01066              "Average Number of Clusters with Reco Shower Energy (excl xtalk clu)",
01067                                    100,0,25,20,0,20);
01068   }
01069 
01070 
01071   TH1F *EMFracUVDiffTrack = new TH1F("EMFracUVDiffTrack","U/V Diff in Reco EMFrac",200,-2,2);
01072   EMFracUVDiffTrack->SetXTitle("U-V EMFrac");
01073   TH1F *EMFracUVDiffShower = new TH1F("EMFracUVDiffShower","U/V Diff in Reco EMFrac",200,-2,2);
01074   EMFracUVDiffShower->SetXTitle("U-V EMFrac");
01075   
01076   TH1F *EMFracUVDiffTrackSplit[6] = {0};
01077   for(int i=0;i<6;i++) {
01078     char name[256];
01079     sprintf(name,"%s_%s",EMFracUVDiffTrack->GetName(),SplitName[i].c_str());
01080     EMFracUVDiffTrackSplit[i] = new TH1F(name,"U/V Diff in Reco EMFrac",200,-2,2);
01081   }
01082   TH1F *EMFracUVDiffShowerSplit[6] = {0};
01083   for(int i=0;i<6;i++) {
01084     char name[256];
01085     sprintf(name,"%s_%s",EMFracUVDiffShower->GetName(),SplitName[i].c_str());
01086     EMFracUVDiffShowerSplit[i] = new TH1F(name,"U/V Diff in Reco EMFrac",200,-2,2);
01087   }
01088 
01089   for(int i=0;i<Nentries;i++){
01090     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
01091                             << "% done" << std::endl;
01092 
01093     if(!this->GetEntry(i)) continue;
01094     //cout << ntpHeader->GetRun() << endl;
01095     if(eventSummary->nevent==0) continue;
01096 
01097     //if not MC:
01098     if(ntpHeader->GetVldContext().GetSimFlag()!=4) {
01099       Double_t snarltime = (ntpHeader->GetVldContext().GetTimeStamp().GetSec() + 
01100                             (ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9));
01101       Int_t    month     = eventSummary->date.month; 
01102       Int_t    year      = eventSummary->date.year;
01103 
01104       SpillInfoBlock spillinfo;
01105 
01106       pppinfo.GetSpillInfo(month,year,snarltime,spillinfo);
01107       if(spillinfo.GetTgtpos()>1000) continue;
01108       if(spillinfo.GetHorn()<16000)  continue;
01109       if(spillinfo.GetHpos()<-2 || 
01110          spillinfo.GetHpos()>-1.6)   continue;
01111       if(spillinfo.GetVpos()<0.8 ||
01112          spillinfo.GetVpos()>1.2) continue;
01113       if(spillinfo.GetMagnet()>-4800.) continue;
01114       datapot += spillinfo.GetPot();      
01115     }
01116 
01117     Int_t is_fid = 0;
01118     //fill tree once for each reconstructed event
01119     for(int j=0;j<eventSummary->nevent;j++){ 
01120       if(!LoadEvent(j)) continue; //no event found
01121       if(ntpEvent->nshower==0) continue; //no shower found
01122       
01123       //get detector type:
01124       if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01125         //fiducial criteria on vtx for far detector
01126         is_fid = 1;
01127         if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
01128       }
01129       else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01130         //fiducial criteria on vtx for near detector
01131         is_fid = 1;
01132         if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
01133       }
01134       if(is_fid==0) continue;
01135       
01136       Int_t shower = -1;
01137       if(!LoadLargestShowerFromEvent(ntpEvent->index,shower)) continue;
01138       
01139       if(!PassDataCleaningCuts(shower)) continue;
01140 
01141       Int_t noxClu = 0;
01142       double emfracU = 0;
01143       double emfracV = 0;
01144       double emfrac = RecoEMFrac(shower,0,emfracU,emfracV);
01145       double emfracprob = RecoEMFrac(shower,1,emfracU,emfracV);
01146       Int_t primaryUID = -1;
01147       Float_t primaryUprobem = 0;
01148       Double_t maxEnU = 0;
01149       Int_t primaryVID = -1;
01150       Float_t primaryVprobem = 0;
01151       Double_t maxEnV = 0;
01152 
01153       for(int k=0; k<ntpShower->ncluster; k++){
01154         int index = ntpShower->clu[k];
01155         if(!LoadCluster(index)) continue;
01156         if(ntpCluster->ph.gev>0.5 && 
01157            (ntpCluster->id==0 || ntpCluster->id==1 || 
01158             ntpCluster->id==3)) noxClu++; 
01159         if(ntpCluster->planeview==2 && 
01160            ntpCluster->ph.gev>maxEnU) {
01161           maxEnU = ntpCluster->ph.gev;
01162           primaryUID = ntpCluster->id;
01163           primaryUprobem = ntpCluster->probem;
01164         }
01165         if(ntpCluster->planeview==3 && 
01166            ntpCluster->ph.gev>maxEnV) {
01167           maxEnV = ntpCluster->ph.gev;
01168           primaryVID = ntpCluster->id;
01169           primaryVprobem = ntpCluster->probem;
01170         }
01171       }
01172       
01173       Int_t track = -1;      
01174       LoadLargestTrackFromEvent(ntpEvent->index,track);
01175       
01176       if(track>=0 && PassTrackCuts(track)) {
01177         //Fill track events:
01178         double y = 1. - RecoMuEnergy(0,track)/ntpEvent->ph.gev;
01179         for(int k=0; k<ntpShower->ncluster; k++){
01180           int index = ntpShower->clu[k];
01181           if(!LoadCluster(index)) continue;
01182           IDHistTrack->Fill(ntpCluster->id);
01183         }
01184         RecoyTrack->Fill(y);
01185         RecoEMFracTrack->Fill(emfrac);
01186         RecoEMFracProbTrack->Fill(emfracprob);
01187         UClusterIDTrack->Fill(primaryUID);
01188         VClusterIDTrack->Fill(primaryVID);
01189         if(primaryUID==0) UClusterProbEMTrack->Fill(primaryUprobem);
01190         if(primaryVID==0) VClusterProbEMTrack->Fill(primaryVprobem);
01191         AvgCluEnTrack->Fill(ntpShower->ph.gev,ntpShower->ncluster);
01192         AvgCluEnTracknox->Fill(ntpShower->ph.gev,noxClu);
01193         EMFracUVDiffTrack->Fill(emfracU-emfracV);
01194 
01195         Int_t mcevent = -1;
01196         Int_t hist_index = 5;
01197         if(LoadTruthForRecoTH(ntpEvent->index,mcevent)){          
01198           if(ntpTruth->iaction==1){
01199             if(ntpTruth->inu==12 && ntpTruth->inunoosc==12) hist_index = 0;
01200             if(ntpTruth->inu==12 && ntpTruth->inunoosc==14) hist_index = 1;
01201             if(ntpTruth->inu==14 && ntpTruth->inunoosc==14) hist_index = 2;
01202             if(ntpTruth->inu==16 && ntpTruth->inunoosc==14) hist_index = 3;
01203           }
01204           else hist_index = 4;
01205         }
01206         for(int k=0; k<ntpShower->ncluster; k++){
01207           int index = ntpShower->clu[k];
01208           if(!LoadCluster(index)) continue;
01209           IDHistTrackSplit[hist_index]->Fill(ntpCluster->id);
01210         }
01211         RecoyTrackSplit[hist_index]->Fill(y);
01212         RecoEMFracTrackSplit[hist_index]->Fill(emfrac);
01213         RecoEMFracProbTrackSplit[hist_index]->Fill(emfracprob);
01214         UClusterIDTrackSplit[hist_index]->Fill(primaryUID);
01215         VClusterIDTrackSplit[hist_index]->Fill(primaryVID);
01216         if(primaryUID==0) UClusterProbEMTrackSplit[hist_index]->Fill(primaryUprobem);
01217         if(primaryVID==0) VClusterProbEMTrackSplit[hist_index]->Fill(primaryVprobem);
01218         AvgCluEnTrackSplit[hist_index]->Fill(ntpShower->ph.gev,ntpShower->ncluster);
01219         AvgCluEnnoxTrackSplit[hist_index]->Fill(ntpShower->ph.gev,noxClu);
01220         EMFracUVDiffTrackSplit[hist_index]->Fill(emfracU-emfracV);
01221       }
01222       else {
01223         //Fill shower events:
01224         for(int k=0; k<ntpShower->ncluster; k++){
01225           int index = ntpShower->clu[k];
01226           if(!LoadCluster(index)) continue;
01227           IDHistShower->Fill(ntpCluster->id);
01228         }
01229         if(primaryVprobem>0.2 && primaryUprobem>0.2 &&
01230            primaryUID==0  && primaryVID==0) RecoyShower->Fill(emfracprob);
01231         RecoEMFracShower->Fill(emfrac);
01232         RecoEMFracProbShower->Fill(emfracprob);
01233         UClusterIDShower->Fill(primaryUID);
01234         VClusterIDShower->Fill(primaryVID);
01235         if(primaryUID==0) UClusterProbEMShower->Fill(primaryUprobem);
01236         if(primaryVID==0) VClusterProbEMShower->Fill(primaryVprobem);
01237         AvgCluEnShower->Fill(ntpShower->ph.gev,ntpShower->ncluster);
01238         AvgCluEnShowernox->Fill(ntpShower->ph.gev,noxClu);
01239         EMFracUVDiffShower->Fill(emfracU-emfracV);      
01240 
01241         Int_t mcevent = -1;
01242         Int_t hist_index = 5;
01243         if(LoadTruthForRecoTH(ntpEvent->index,mcevent)){          
01244           if(ntpTruth->iaction==1){
01245             if(ntpTruth->inu==12 && ntpTruth->inunoosc==12) hist_index = 0;
01246             if(ntpTruth->inu==12 && ntpTruth->inunoosc==14) hist_index = 1;
01247             if(ntpTruth->inu==14 && ntpTruth->inunoosc==14) hist_index = 2;
01248             if(ntpTruth->inu==16 && ntpTruth->inunoosc==14) hist_index = 3;
01249           }
01250           else hist_index = 4;
01251         }
01252         for(int k=0; k<ntpShower->ncluster; k++){
01253           int index = ntpShower->clu[k];
01254           if(!LoadCluster(index)) continue;
01255           IDHistShowerSplit[hist_index]->Fill(ntpCluster->id);
01256         }
01257         if(primaryVprobem>0.2 && primaryUprobem>0.2 &&
01258            primaryUID==0  && primaryVID==0) {
01259           RecoyShowerSplit[hist_index]->Fill(emfrac);
01260         }
01261         RecoEMFracShowerSplit[hist_index]->Fill(emfrac);
01262         RecoEMFracProbShowerSplit[hist_index]->Fill(emfracprob);
01263         UClusterIDShowerSplit[hist_index]->Fill(primaryUID);
01264         VClusterIDShowerSplit[hist_index]->Fill(primaryVID);
01265         if(primaryUID==0) UClusterProbEMShowerSplit[hist_index]->Fill(primaryUprobem);
01266         if(primaryVID==0) VClusterProbEMShowerSplit[hist_index]->Fill(primaryVprobem);
01267         AvgCluEnShowerSplit[hist_index]->Fill(ntpShower->ph.gev,ntpShower->ncluster);
01268         AvgCluEnnoxShowerSplit[hist_index]->Fill(ntpShower->ph.gev,noxClu);
01269         EMFracUVDiffShowerSplit[hist_index]->Fill(emfracU-emfracV);
01270       }
01271     }
01272   }
01273   cout << "DATAPOT = " << datapot << endl;
01274   save->Write();
01275   save->Close();
01276   return;
01277 
01278 }

void MadCluAnalysis::DataMCComp ( char *  dataname,
char *  mcname 
)

Definition at line 1280 of file MadCluAnalysis.cxx.

01280                                                           {
01281 
01282   TFile *datafile = new TFile(dataname,"READ");
01283   TH1F *DataIDHistTrack = (TH1F*) datafile->Get("IDHistTrack");
01284   TH1F *DataRecoyTrack = (TH1F*) datafile->Get("RecoyTrack");
01285   TH1F *DataRecoEMFracTrack = (TH1F*) datafile->Get("RecoEMFracTrack");
01286   TH1F *DataRecoEMFracProbTrack = (TH1F*) datafile->Get("RecoEMFracProbTrack");
01287   TH1F *DataUClusterIDTrack = (TH1F*) datafile->Get("UClusterIDTrack");
01288   TH1F *DataVClusterIDTrack = (TH1F*) datafile->Get("VClusterIDTrack");
01289   TH1F *DataUClusterProbEMTrack = (TH1F*) datafile->Get("UClusterProbEMTrack");
01290   TH1F *DataVClusterProbEMTrack = (TH1F*) datafile->Get("VClusterProbEMTrack");
01291   TH2F *DataAvgCluEnTrack = (TH2F*) datafile->Get("AvgCluEnTrack");
01292   TH2F *DataAvgCluEnTracknox = (TH2F*) datafile->Get("AvgCluEnTracknox");
01293   TH1F *DataEMFracUVDiffTrack = (TH1F*) datafile->Get("EMFracUVDiffTrack");
01294 
01295   TH1F *DataIDHistShower = (TH1F*) datafile->Get("IDHistShower");
01296   TH1F *DataRecoyShower = (TH1F*) datafile->Get("RecoyShower");
01297   TH1F *DataRecoEMFracShower = (TH1F*) datafile->Get("RecoEMFracShower");
01298   TH1F *DataRecoEMFracProbShower = (TH1F*) datafile->Get("RecoEMFracProbShower");
01299   TH1F *DataUClusterIDShower = (TH1F*) datafile->Get("UClusterIDShower");
01300   TH1F *DataVClusterIDShower = (TH1F*) datafile->Get("VClusterIDShower");
01301   TH1F *DataUClusterProbEMShower = (TH1F*) datafile->Get("UClusterProbEMShower");
01302   TH1F *DataVClusterProbEMShower = (TH1F*) datafile->Get("VClusterProbEMShower");
01303   TH2F *DataAvgCluEnShower = (TH2F*) datafile->Get("AvgCluEnShower");
01304   TH2F *DataAvgCluEnShowernox = (TH2F*) datafile->Get("AvgCluEnShowernox");
01305   TH1F *DataEMFracUVDiffShower = (TH1F*) datafile->Get("EMFracUVDiffShower");
01306 
01307   TFile *mcfile = new TFile(mcname,"READ");  
01308   TH1F *IDHistTrack = (TH1F*) mcfile->Get("IDHistTrack");
01309   TH1F *RecoyTrack = (TH1F*) mcfile->Get("RecoyTrack");
01310   TH1F *RecoEMFracTrack = (TH1F*) mcfile->Get("RecoEMFracTrack");
01311   TH1F *RecoEMFracProbTrack = (TH1F*) mcfile->Get("RecoEMFracProbTrack");
01312   TH1F *UClusterIDTrack = (TH1F*) mcfile->Get("UClusterIDTrack");
01313   TH1F *VClusterIDTrack = (TH1F*) mcfile->Get("VClusterIDTrack");
01314   TH1F *UClusterProbEMTrack = (TH1F*) mcfile->Get("UClusterProbEMTrack");
01315   TH1F *VClusterProbEMTrack = (TH1F*) mcfile->Get("VClusterProbEMTrack");
01316   TH2F *AvgCluEnTrack = (TH2F*) mcfile->Get("AvgCluEnTrack");
01317   TH2F *AvgCluEnTracknox = (TH2F*) mcfile->Get("AvgCluEnTracknox");
01318   TH1F *EMFracUVDiffTrack = (TH1F*) mcfile->Get("EMFracUVDiffTrack");
01319 
01320   TH1F *IDHistShower = (TH1F*) mcfile->Get("IDHistShower");
01321   TH1F *RecoyShower = (TH1F*) mcfile->Get("RecoyShower");
01322   TH1F *RecoEMFracShower = (TH1F*) mcfile->Get("RecoEMFracShower");
01323   TH1F *RecoEMFracProbShower = (TH1F*) mcfile->Get("RecoEMFracProbShower");
01324   TH1F *UClusterIDShower = (TH1F*) mcfile->Get("UClusterIDShower");
01325   TH1F *VClusterIDShower = (TH1F*) mcfile->Get("VClusterIDShower");
01326   TH1F *UClusterProbEMShower = (TH1F*) mcfile->Get("UClusterProbEMShower");
01327   TH1F *VClusterProbEMShower = (TH1F*) mcfile->Get("VClusterProbEMShower");
01328   TH2F *AvgCluEnShower = (TH2F*) mcfile->Get("AvgCluEnShower");
01329   TH2F *AvgCluEnShowernox = (TH2F*) mcfile->Get("AvgCluEnShowernox");
01330   TH1F *EMFracUVDiffShower = (TH1F*) mcfile->Get("EMFracUVDiffShower");
01331 
01332   string SplitName[6] = {"beamnue","nue","numu","nutau","nc","unknown"};
01333 
01334   TH1F *IDHistTrackSplit[6] = {0};
01335   for(int i=0;i<6;i++) {
01336     char name[256];
01337     sprintf(name,"%s_%s",IDHistTrack->GetName(),SplitName[i].c_str());
01338     IDHistTrackSplit[i] = (TH1F*) mcfile->Get(name);
01339   }
01340   TH1F *IDHistShowerSplit[6] = {0};
01341   for(int i=0;i<6;i++) {
01342     char name[256];
01343     sprintf(name,"%s_%s",IDHistShower->GetName(),SplitName[i].c_str());
01344     IDHistShowerSplit[i] = (TH1F*) mcfile->Get(name);
01345   }
01346 
01347   TH1F *RecoyTrackSplit[6] = {0};
01348   for(int i=0;i<6;i++) {
01349     char name[256];
01350     sprintf(name,"%s_%s",RecoyTrack->GetName(),SplitName[i].c_str());
01351     RecoyTrackSplit[i] = (TH1F*) mcfile->Get(name);
01352   }
01353   TH1F *RecoyShowerSplit[6] = {0};
01354   for(int i=0;i<6;i++) {
01355     char name[256];
01356     sprintf(name,"%s_%s",RecoyShower->GetName(),SplitName[i].c_str());
01357     RecoyShowerSplit[i] = (TH1F*) mcfile->Get(name);
01358   }
01359 
01360   TH1F *RecoEMFracTrackSplit[6] = {0};
01361   for(int i=0;i<6;i++) {
01362     char name[256];
01363     sprintf(name,"%s_%s",RecoEMFracTrack->GetName(),SplitName[i].c_str());
01364     RecoEMFracTrackSplit[i] = (TH1F*) mcfile->Get(name);
01365   }
01366   TH1F *RecoEMFracShowerSplit[6] = {0};
01367   for(int i=0;i<6;i++) {
01368     char name[256];
01369     sprintf(name,"%s_%s",RecoEMFracShower->GetName(),SplitName[i].c_str());
01370     RecoEMFracShowerSplit[i] = (TH1F*) mcfile->Get(name);
01371   }
01372 
01373   TH1F *RecoEMFracProbTrackSplit[6] = {0};
01374   for(int i=0;i<6;i++) {
01375     char name[256];
01376     sprintf(name,"%s_%s",RecoEMFracProbTrack->GetName(),SplitName[i].c_str());
01377     RecoEMFracProbTrackSplit[i] = (TH1F*) mcfile->Get(name);
01378   }
01379   TH1F *RecoEMFracProbShowerSplit[6] = {0};
01380   for(int i=0;i<6;i++) {
01381     char name[256];
01382     sprintf(name,"%s_%s",RecoEMFracProbShower->GetName(),SplitName[i].c_str());
01383     RecoEMFracProbShowerSplit[i] = (TH1F*) mcfile->Get(name);
01384   }
01385 
01386   TH1F *UClusterIDTrackSplit[6] = {0};
01387   TH1F *VClusterIDTrackSplit[6] = {0};
01388   TH1F *UClusterProbEMTrackSplit[6] = {0};
01389   TH1F *VClusterProbEMTrackSplit[6] = {0};
01390   for(int i=0;i<6;i++) {
01391     char name[256];
01392     sprintf(name,"%s_%s",UClusterIDTrack->GetName(),SplitName[i].c_str());
01393     UClusterIDTrackSplit[i] = (TH1F*) mcfile->Get(name);
01394     sprintf(name,"%s_%s",VClusterIDTrack->GetName(),SplitName[i].c_str());
01395     VClusterIDTrackSplit[i] = (TH1F*) mcfile->Get(name);
01396     sprintf(name,"%s_%s",UClusterProbEMTrack->GetName(),SplitName[i].c_str());
01397     UClusterProbEMTrackSplit[i] = (TH1F*) mcfile->Get(name);
01398     sprintf(name,"%s_%s",VClusterProbEMTrack->GetName(),SplitName[i].c_str());
01399     VClusterProbEMTrackSplit[i] = (TH1F*) mcfile->Get(name);
01400   }
01401   TH1F *UClusterIDShowerSplit[6] = {0};
01402   TH1F *VClusterIDShowerSplit[6] = {0};
01403   TH1F *UClusterProbEMShowerSplit[6] = {0};
01404   TH1F *VClusterProbEMShowerSplit[6] = {0};
01405   for(int i=0;i<6;i++) {
01406     char name[256];
01407     sprintf(name,"%s_%s",UClusterIDShower->GetName(),SplitName[i].c_str());
01408     UClusterIDShowerSplit[i] = (TH1F*) mcfile->Get(name);
01409     sprintf(name,"%s_%s",VClusterIDShower->GetName(),SplitName[i].c_str());
01410     VClusterIDShowerSplit[i] = (TH1F*) mcfile->Get(name);
01411     sprintf(name,"%s_%s",UClusterProbEMShower->GetName(),SplitName[i].c_str());
01412     UClusterProbEMShowerSplit[i] = (TH1F*) mcfile->Get(name);
01413     sprintf(name,"%s_%s",VClusterProbEMShower->GetName(),SplitName[i].c_str());
01414     VClusterProbEMShowerSplit[i] = (TH1F*) mcfile->Get(name);
01415   }
01416 
01417   TH2F *AvgCluEnTrackSplit[6] = {0};
01418   TH2F *AvgCluEnnoxTrackSplit[6] = {0};
01419   for(int i=0;i<6;i++) {
01420     char name[256];
01421     sprintf(name,"%s_%s",AvgCluEnTrack->GetName(),SplitName[i].c_str());
01422     AvgCluEnTrackSplit[i] = (TH2F*) mcfile->Get(name);
01423     sprintf(name,"%s_%s",AvgCluEnTracknox->GetName(),SplitName[i].c_str());
01424     AvgCluEnnoxTrackSplit[i] = (TH2F*) mcfile->Get(name);
01425   }
01426   TH2F *AvgCluEnShowerSplit[6] = {0};
01427   TH2F *AvgCluEnnoxShowerSplit[6] = {0};
01428   for(int i=0;i<6;i++) {
01429     char name[256];
01430     sprintf(name,"%s_%s",AvgCluEnShower->GetName(),SplitName[i].c_str());
01431     AvgCluEnShowerSplit[i] = (TH2F*) mcfile->Get(name);
01432     sprintf(name,"%s_%s",AvgCluEnShowernox->GetName(),SplitName[i].c_str());
01433     AvgCluEnnoxShowerSplit[i] = (TH2F*) mcfile->Get(name);
01434   }
01435 
01436   TH1F *EMFracUVDiffTrackSplit[6] = {0};
01437   for(int i=0;i<6;i++) {
01438     char name[256];
01439     sprintf(name,"%s_%s",EMFracUVDiffTrack->GetName(),SplitName[i].c_str());
01440     EMFracUVDiffTrackSplit[i] = (TH1F*) mcfile->Get(name);
01441   }
01442   TH1F *EMFracUVDiffShowerSplit[6] = {0};
01443   for(int i=0;i<6;i++) {
01444     char name[256];
01445     sprintf(name,"%s_%s",EMFracUVDiffShower->GetName(),SplitName[i].c_str());
01446     EMFracUVDiffShowerSplit[i] = (TH1F*) mcfile->Get(name);
01447   }
01448 
01449 
01450   //Double_t mc_scale = DataUClusterIDTrack->GetEntries() + DataUClusterIDShower->GetEntries();
01451   //mc_scale /= (UClusterIDTrack->GetEntries() + UClusterIDShower->GetEntries());
01452   //cout << "mc_scale = " << mc_scale << endl;
01453   cout << "MD: " << 550*2.4e13*40 << " POT" << endl;
01454   cout << "DATA: " << 354657e12 << " POT" << endl;
01455   Double_t mc_scale = 3.54657e17/(550.*2.4e13*40.);
01456 
01457   IDHistTrack->Scale(mc_scale);
01458   RecoyTrack->Scale(mc_scale);
01459   RecoEMFracTrack->Scale(mc_scale);
01460   RecoEMFracProbTrack->Scale(mc_scale);
01461   UClusterIDTrack->Scale(mc_scale);
01462   VClusterIDTrack->Scale(mc_scale);
01463   UClusterProbEMTrack->Scale(mc_scale);
01464   VClusterProbEMTrack->Scale(mc_scale);
01465   EMFracUVDiffTrack->Scale(mc_scale);
01466 
01467   IDHistTrack->SetLineColor(2);
01468   RecoyTrack->SetLineColor(2);
01469   RecoEMFracTrack->SetLineColor(2);
01470   RecoEMFracProbTrack->SetLineColor(2);
01471   UClusterIDTrack->SetLineColor(2);
01472   VClusterIDTrack->SetLineColor(2);
01473   UClusterProbEMTrack->SetLineColor(2);
01474   VClusterProbEMTrack->SetLineColor(2);
01475   EMFracUVDiffTrack->SetLineColor(2);
01476 
01477   for(int i=0;i<6;i++){
01478     IDHistTrackSplit[i]->Scale(mc_scale);
01479     RecoyTrackSplit[i]->Scale(mc_scale);
01480     RecoEMFracTrackSplit[i]->Scale(mc_scale);
01481     RecoEMFracProbTrackSplit[i]->Scale(mc_scale);
01482     UClusterIDTrackSplit[i]->Scale(mc_scale);
01483     VClusterIDTrackSplit[i]->Scale(mc_scale);
01484     UClusterProbEMTrackSplit[i]->Scale(mc_scale);
01485     VClusterProbEMTrackSplit[i]->Scale(mc_scale);
01486     EMFracUVDiffTrackSplit[i]->Scale(mc_scale);
01487     Int_t col = i+3;
01488     if(col>=5) col+=1;
01489     IDHistTrackSplit[i]->SetLineColor(col);
01490     RecoyTrackSplit[i]->SetLineColor(col);
01491     RecoEMFracTrackSplit[i]->SetLineColor(col);
01492     RecoEMFracProbTrackSplit[i]->SetLineColor(col);
01493     UClusterIDTrackSplit[i]->SetLineColor(col);
01494     VClusterIDTrackSplit[i]->SetLineColor(col);
01495     UClusterProbEMTrackSplit[i]->SetLineColor(col);
01496     VClusterProbEMTrackSplit[i]->SetLineColor(col);
01497     EMFracUVDiffTrackSplit[i]->SetLineColor(col);
01498   }
01499 
01500   IDHistShower->Scale(mc_scale);
01501   RecoyShower->Scale(mc_scale);
01502   RecoEMFracShower->Scale(mc_scale);
01503   RecoEMFracProbShower->Scale(mc_scale);
01504   UClusterIDShower->Scale(mc_scale);
01505   VClusterIDShower->Scale(mc_scale);
01506   UClusterProbEMShower->Scale(mc_scale);
01507   VClusterProbEMShower->Scale(mc_scale);
01508   EMFracUVDiffShower->Scale(mc_scale);
01509 
01510   IDHistShower->SetLineColor(2);
01511   RecoyShower->SetLineColor(2);
01512   RecoEMFracShower->SetLineColor(2);
01513   RecoEMFracProbShower->SetLineColor(2);
01514   UClusterIDShower->SetLineColor(2);
01515   VClusterIDShower->SetLineColor(2);
01516   UClusterProbEMShower->SetLineColor(2);
01517   VClusterProbEMShower->SetLineColor(2);
01518   EMFracUVDiffShower->SetLineColor(2);
01519 
01520   for(int i=0;i<6;i++){
01521     IDHistShowerSplit[i]->Scale(mc_scale);
01522     RecoyShowerSplit[i]->Scale(mc_scale);
01523     RecoEMFracShowerSplit[i]->Scale(mc_scale);
01524     RecoEMFracProbShowerSplit[i]->Scale(mc_scale);
01525     UClusterIDShowerSplit[i]->Scale(mc_scale);
01526     VClusterIDShowerSplit[i]->Scale(mc_scale);
01527     UClusterProbEMShowerSplit[i]->Scale(mc_scale);
01528     VClusterProbEMShowerSplit[i]->Scale(mc_scale);
01529     EMFracUVDiffShowerSplit[i]->Scale(mc_scale);
01530     Int_t col = i+3;
01531     if(col>=5) col+=1;
01532     IDHistShowerSplit[i]->SetLineColor(col);
01533     RecoyShowerSplit[i]->SetLineColor(col);
01534     RecoEMFracShowerSplit[i]->SetLineColor(col);
01535     RecoEMFracProbShowerSplit[i]->SetLineColor(col);
01536     UClusterIDShowerSplit[i]->SetLineColor(col);
01537     VClusterIDShowerSplit[i]->SetLineColor(col);
01538     UClusterProbEMShowerSplit[i]->SetLineColor(col);
01539     VClusterProbEMShowerSplit[i]->SetLineColor(col);
01540     EMFracUVDiffShowerSplit[i]->SetLineColor(col);
01541   }
01542 
01543   TCanvas *can1 = new TCanvas();
01544   can1->Divide(1,2);
01545   can1->cd(1);
01546   DataIDHistTrack->Draw();
01547   IDHistTrack->Draw("sames");
01548   if(IDHistTrack->GetMaximum()>DataIDHistTrack->GetMaximum()){
01549     DataIDHistTrack->SetMaximum(IDHistTrack->GetMaximum()*1.1);
01550   }
01551   for(int i=0;i<6;i++){
01552     IDHistTrackSplit[i]->Draw("same");
01553   }
01554   can1->cd(2);
01555   DataIDHistShower->Draw();
01556   IDHistShower->Draw("sames");
01557   if(IDHistShower->GetMaximum()>DataIDHistShower->GetMaximum()){
01558     DataIDHistShower->SetMaximum(IDHistShower->GetMaximum()*1.1);
01559   }
01560   for(int i=0;i<6;i++){
01561     IDHistShowerSplit[i]->Draw("same");
01562   }
01563 
01564   TCanvas *can2 = new TCanvas();
01565   can2->Divide(1,2);
01566   can2->cd(1);
01567   DataRecoyTrack->Draw();
01568   RecoyTrack->Draw("sames");
01569   if(RecoyTrack->GetMaximum()>DataRecoyTrack->GetMaximum()){
01570     DataRecoyTrack->SetMaximum(RecoyTrack->GetMaximum()*1.1);
01571   }
01572   for(int i=0;i<6;i++){
01573     RecoyTrackSplit[i]->Draw("same");
01574   }
01575   can2->cd(2);
01576   DataRecoyShower->Draw();
01577   RecoyShower->Draw("sames");
01578   if(RecoyShower->GetMaximum()>DataRecoyShower->GetMaximum()){
01579     DataRecoyShower->SetMaximum(RecoyShower->GetMaximum()*1.1);
01580   }
01581   for(int i=0;i<6;i++){
01582     RecoyShowerSplit[i]->Draw("same");
01583   }
01584 
01585   TCanvas *can3 = new TCanvas();
01586   can3->Divide(1,2);
01587   can3->cd(1);
01588   DataRecoEMFracTrack->Draw();
01589   RecoEMFracTrack->Draw("sames");
01590   if(RecoEMFracTrack->GetMaximum()>DataRecoEMFracTrack->GetMaximum()){
01591     DataRecoEMFracTrack->SetMaximum(RecoEMFracTrack->GetMaximum()*1.1);
01592   }
01593   for(int i=0;i<6;i++){
01594     RecoEMFracTrackSplit[i]->Draw("same");
01595   }
01596   can3->cd(2);
01597   DataRecoEMFracShower->Draw();
01598   RecoEMFracShower->Draw("sames");
01599   if(RecoEMFracShower->GetMaximum()>DataRecoEMFracShower->GetMaximum()){
01600     DataRecoEMFracShower->SetMaximum(RecoEMFracShower->GetMaximum()*1.1);
01601   }
01602   for(int i=0;i<6;i++){
01603     RecoEMFracShowerSplit[i]->Draw("same");
01604   }
01605 
01606   TCanvas *can4 = new TCanvas();
01607   can4->Divide(1,2);
01608   can4->cd(1);
01609   DataRecoEMFracProbTrack->Draw();
01610   RecoEMFracProbTrack->Draw("sames");
01611   if(RecoEMFracProbTrack->GetMaximum()>DataRecoEMFracProbTrack->GetMaximum()){
01612     DataRecoEMFracProbTrack->SetMaximum(RecoEMFracProbTrack->GetMaximum()*1.1);
01613   }
01614   for(int i=0;i<6;i++){
01615     RecoEMFracProbTrackSplit[i]->Draw("same");
01616   }
01617 
01618   can4->cd(2);
01619   DataRecoEMFracProbShower->Draw();
01620   RecoEMFracProbShower->Draw("sames");
01621   if(RecoEMFracProbShower->GetMaximum()>DataRecoEMFracProbShower->GetMaximum()){
01622     DataRecoEMFracProbShower->SetMaximum(RecoEMFracProbShower->GetMaximum()*1.1);
01623   }
01624   for(int i=0;i<6;i++){
01625     RecoEMFracProbShowerSplit[i]->Draw("same");
01626   }
01627 
01628   TCanvas *can5 = new TCanvas();
01629   can5->Divide(1,2);
01630   can5->cd(1);
01631   DataUClusterIDTrack->Draw();
01632   UClusterIDTrack->Draw("sames");
01633   if(UClusterIDTrack->GetMaximum()>DataUClusterIDTrack->GetMaximum()){
01634     DataUClusterIDTrack->SetMaximum(UClusterIDTrack->GetMaximum()*1.1);
01635   }
01636   for(int i=0;i<6;i++){
01637     UClusterIDTrackSplit[i]->Draw("same");
01638   }
01639   can5->cd(2);
01640   DataUClusterIDShower->Draw();
01641   UClusterIDShower->Draw("sames");
01642   if(UClusterIDShower->GetMaximum()>DataUClusterIDShower->GetMaximum()){
01643     DataUClusterIDShower->SetMaximum(UClusterIDShower->GetMaximum()*1.1);
01644   }
01645   for(int i=0;i<6;i++){
01646     UClusterIDShowerSplit[i]->Draw("same");
01647   }
01648 
01649   TCanvas *can5b = new TCanvas();
01650   can5b->Divide(1,2);
01651   can5b->cd(1);
01652   DataUClusterProbEMTrack->Draw();
01653   UClusterProbEMTrack->Draw("sames");
01654   if(UClusterProbEMTrack->GetMaximum()>DataUClusterProbEMTrack->GetMaximum()){
01655     DataUClusterProbEMTrack->SetMaximum(UClusterProbEMTrack->GetMaximum()*1.1);
01656   }
01657   for(int i=0;i<6;i++){
01658     UClusterProbEMTrackSplit[i]->Draw("same");
01659   }
01660   can5b->cd(2);
01661   DataUClusterProbEMShower->Draw();
01662   UClusterProbEMShower->Draw("sames");
01663   if(UClusterProbEMShower->GetMaximum()>DataUClusterProbEMShower->GetMaximum()){
01664     DataUClusterProbEMShower->SetMaximum(UClusterProbEMShower->GetMaximum()*1.1);
01665   }
01666   for(int i=0;i<6;i++){
01667     UClusterProbEMShowerSplit[i]->Draw("same");
01668   }
01669 
01670 
01671   TCanvas *can6 = new TCanvas();
01672   can6->Divide(1,2);
01673   can6->cd(1);
01674   DataVClusterIDTrack->Draw();
01675   VClusterIDTrack->Draw("sames");
01676   if(VClusterIDTrack->GetMaximum()>DataVClusterIDTrack->GetMaximum()){
01677     DataVClusterIDTrack->SetMaximum(VClusterIDTrack->GetMaximum()*1.1);
01678   }
01679   for(int i=0;i<6;i++){
01680     VClusterIDTrackSplit[i]->Draw("same");
01681   }
01682   can6->cd(2);
01683   DataVClusterIDShower->Draw();
01684   VClusterIDShower->Draw("sames");
01685   if(VClusterIDShower->GetMaximum()>DataVClusterIDShower->GetMaximum()){
01686     DataVClusterIDShower->SetMaximum(VClusterIDShower->GetMaximum()*1.1);
01687   }
01688   for(int i=0;i<6;i++){
01689     VClusterIDShowerSplit[i]->Draw("same");
01690   }
01691 
01692   TCanvas *can6b = new TCanvas();
01693   can6b->Divide(1,2);
01694   can6b->cd(1);
01695   DataVClusterProbEMTrack->Draw();
01696   VClusterProbEMTrack->Draw("sames");
01697   if(VClusterProbEMTrack->GetMaximum()>DataVClusterProbEMTrack->GetMaximum()){
01698     DataVClusterProbEMTrack->SetMaximum(VClusterProbEMTrack->GetMaximum()*1.1);
01699   }
01700   for(int i=0;i<6;i++){
01701     VClusterProbEMTrackSplit[i]->Draw("same");
01702   }
01703   can6b->cd(2);
01704   DataVClusterProbEMShower->Draw();
01705   VClusterProbEMShower->Draw("sames");
01706   if(VClusterProbEMShower->GetMaximum()>DataVClusterProbEMShower->GetMaximum()){
01707     DataVClusterProbEMShower->SetMaximum(VClusterProbEMShower->GetMaximum()*1.1);
01708   }
01709   for(int i=0;i<6;i++){
01710     VClusterProbEMShowerSplit[i]->Draw("same");
01711   }
01712 
01713   TCanvas *can7 = new TCanvas();
01714   can7->Divide(3,2);
01715   can7->cd(1);
01716   DataAvgCluEnTrack->Draw("COLZ");
01717   can7->cd(2);
01718   AvgCluEnTrack->Draw("COLZ");
01719   can7->cd(3);
01720   TProfile *pDT = DataAvgCluEnTrack->ProfileX("pDT");
01721   TProfile *pMT = AvgCluEnTrack->ProfileX("pMT");
01722   pMT->SetLineColor(2);
01723   pMT->SetMarkerColor(2);
01724   pDT->Draw();
01725   pMT->Draw("sames");
01726   for(int i=0;i<6;i++){
01727     TProfile *temp = AvgCluEnTrackSplit[i]->ProfileX();
01728     Int_t col = i+3;
01729     if(col>=5) col+=1;
01730     temp->SetLineColor(col);
01731     temp->SetMarkerColor(col);
01732     //temp->Draw("same");
01733   }
01734 
01735   can7->cd(4);
01736   DataAvgCluEnShower->Draw("COLZ");
01737   can7->cd(5);
01738   AvgCluEnShower->Draw("COLZ");
01739   can7->cd(6);
01740   TProfile *pDS = DataAvgCluEnShower->ProfileX("pDS");
01741   TProfile *pMS = AvgCluEnShower->ProfileX("pMS");
01742   pMS->SetLineColor(2);
01743   pMS->SetMarkerColor(2);
01744   pDS->Draw();
01745   pMS->Draw("sames");
01746   for(int i=0;i<6;i++){
01747     TProfile *temp = AvgCluEnShowerSplit[i]->ProfileX();
01748     Int_t col = i+3;
01749     if(col>=5) col+=1;
01750     temp->SetLineColor(col);
01751     temp->SetMarkerColor(col);
01752     //temp->Draw("same");
01753   }
01754 
01755   TCanvas *can8 = new TCanvas();
01756   can8->Divide(3,2);
01757   can8->cd(1);
01758   DataAvgCluEnTracknox->Draw("COLZ");
01759   can8->cd(2);
01760   AvgCluEnTracknox->Draw("COLZ");
01761   can8->cd(3);
01762   TProfile *pDTx = DataAvgCluEnTracknox->ProfileX("pDTx");
01763   TProfile *pMTx = AvgCluEnTracknox->ProfileX("pMTx");
01764   pMTx->SetLineColor(2);
01765   pMTx->SetMarkerColor(2);
01766   pDTx->Draw();
01767   pMTx->Draw("sames");
01768   for(int i=0;i<6;i++){
01769     TProfile *temp = AvgCluEnnoxTrackSplit[i]->ProfileX();
01770     Int_t col = i+3;
01771     if(col>=5) col+=1;
01772     temp->SetLineColor(col);
01773     temp->SetMarkerColor(col);
01774     //temp->Draw("same");
01775   }
01776 
01777   can8->cd(4);
01778   DataAvgCluEnShowernox->Draw("COLZ");
01779   can8->cd(5);
01780   AvgCluEnShowernox->Draw("COLZ");
01781   can8->cd(6);
01782   TProfile *pDSx = DataAvgCluEnShowernox->ProfileX("pDSx");
01783   TProfile *pMSx = AvgCluEnShowernox->ProfileX("pMSx");
01784   pMSx->SetLineColor(2);
01785   pMSx->SetMarkerColor(2);
01786   pDSx->Draw();
01787   pMSx->Draw("sames");
01788   for(int i=0;i<6;i++){
01789     TProfile *temp = AvgCluEnnoxShowerSplit[i]->ProfileX();
01790     Int_t col = i+3;
01791     if(col>=5) col+=1;
01792     temp->SetLineColor(col);
01793     temp->SetMarkerColor(col);
01794     //temp->Draw("same");
01795   }
01796 
01797   TCanvas *can9 = new TCanvas();
01798   can9->Divide(1,2);
01799   can9->cd(1);
01800   DataEMFracUVDiffTrack->Draw();
01801   EMFracUVDiffTrack->Draw("sames");
01802   if(EMFracUVDiffTrack->GetMaximum()>DataEMFracUVDiffTrack->GetMaximum()){
01803     DataEMFracUVDiffTrack->SetMaximum(EMFracUVDiffTrack->GetMaximum()*1.1);
01804   }
01805   for(int i=0;i<6;i++){
01806     EMFracUVDiffTrackSplit[i]->Draw("same"); 
01807   }
01808   can9->cd(2);
01809   DataEMFracUVDiffShower->Draw();
01810   EMFracUVDiffShower->Draw("sames");
01811   if(EMFracUVDiffShower->GetMaximum()>DataEMFracUVDiffShower->GetMaximum()){
01812     DataEMFracUVDiffShower->SetMaximum(EMFracUVDiffShower->GetMaximum()*1.1);
01813   }
01814   for(int i=0;i<6;i++){
01815     EMFracUVDiffShowerSplit[i]->Draw("same"); 
01816   }
01817 }

Bool_t MadCluAnalysis::FillAveNumSSPlane ( Int_t  event  ) 

Definition at line 2703 of file MadCluAnalysis.cxx.

References AveNumSSPlane(), PAN_AveNumSSPlaneU, and PAN_AveNumSSPlaneV.

Referenced by CallUserFunctions().

02704 {
02705   Double_t *aveNumSSPlane = AveNumSSPlane(event);
02706   PAN_AveNumSSPlaneU = aveNumSSPlane[0];
02707   PAN_AveNumSSPlaneV = aveNumSSPlane[1];
02708   delete [] aveNumSSPlane;
02709   return true;
02710 }

Bool_t MadCluAnalysis::FillBestProbEM ( Int_t  event  ) 

Definition at line 2950 of file MadCluAnalysis.cxx.

References NtpSRShower::clu, NtpSRStripPulseHeight::gev, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), NtpSRShower::ncluster, MadBase::ntpCluster, MadBase::ntpShower, PAN_BestProbEM, NtpSRCluster::ph, NtpSRCluster::planeview, and NtpSRCluster::probem.

Referenced by CallUserFunctions().

02950                                                 {
02951   
02952   PAN_BestProbEM = 0;
02953   
02954   if(!LoadEvent(event)) return false;
02955   Int_t shower = -1;
02956   if(!LoadLargestShowerFromEvent(event,shower)) return false;
02957 
02958   Double_t bestEnU = 0.;
02959   Double_t bestEnV = 0.;
02960   Double_t bestProbU = 0;
02961   Double_t bestProbV = 0;
02962   Int_t *clusters = ntpShower->clu;
02963   for(int i=0;i<ntpShower->ncluster;i++){
02964     if(!LoadCluster(clusters[i])) continue;    
02965     if(ntpCluster->planeview==2&&ntpCluster->ph.gev>bestEnU) {
02966       bestProbU = ntpCluster->probem;
02967       bestEnU = ntpCluster->ph.gev;
02968     }
02969     else if(ntpCluster->planeview==3&&ntpCluster->ph.gev>bestEnV) {
02970       bestProbV += ntpCluster->probem;
02971       bestEnV = ntpCluster->ph.gev;
02972     }
02973   }
02974   PAN_BestProbEM = (bestProbU + bestProbV)/2.;
02975   return true;
02976 }

Bool_t MadCluAnalysis::FillChargeFracRMS ( Int_t  event,
Int_t  method = 0 
)

Definition at line 2655 of file MadCluAnalysis.cxx.

References ChargeFracRMS(), PAN_ChargeFracRMSU, and PAN_ChargeFracRMSV.

Referenced by CallUserFunctions().

02656 {
02657   Double_t *chargeFracRMS = ChargeFracRMS(event,method);
02658   PAN_ChargeFracRMSU = chargeFracRMS[0];
02659   PAN_ChargeFracRMSV = chargeFracRMS[1];
02660   delete [] chargeFracRMS;
02661   return true;
02662 }

Bool_t MadCluAnalysis::FillNCluster ( Int_t  event  ) 

Definition at line 2925 of file MadCluAnalysis.cxx.

References NtpSRShower::clu, NtpSRCluster::id, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), NtpSRShower::ncluster, MadBase::ntpCluster, MadBase::ntpShower, NtpSRShower::nUcluster, NtpSRShower::nVcluster, PAN_NClusterU, PAN_NClusterV, PAN_NPhysClusterU, PAN_NPhysClusterV, and NtpSRCluster::planeview.

Referenced by CallUserFunctions().

02925                                               {
02926 
02927   PAN_NClusterU = 0;
02928   PAN_NClusterV = 0;
02929   PAN_NPhysClusterU = 0;
02930   PAN_NPhysClusterV = 0;
02931 
02932   if(!LoadEvent(event)) return false;  
02933   Int_t shower = -1;
02934   if(!LoadLargestShowerFromEvent(event,shower)) return false;
02935 
02936   PAN_NClusterU = ntpShower->nUcluster;
02937   PAN_NClusterV = ntpShower->nVcluster;
02938 
02939   Int_t *clusters = ntpShower->clu;
02940   for(int i=0;i<ntpShower->ncluster;i++){
02941     if(!LoadCluster(clusters[i])) continue;
02942     if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02943       if(ntpCluster->planeview==2) PAN_NPhysClusterU+=1;
02944       else if(ntpCluster->planeview==3) PAN_NPhysClusterV+=1;
02945     }
02946   }
02947   return true;
02948 }

Bool_t MadCluAnalysis::FillPHWAvgDev ( Int_t  event  ) 

Definition at line 2601 of file MadCluAnalysis.cxx.

References PAN_PHWAvgDevU, PAN_PHWAvgDevV, and PHWAvgDev().

Referenced by CallUserFunctions().

02602 {
02603   Double_t *phwAvgDev = PHWAvgDev(event);
02604   PAN_PHWAvgDevU = phwAvgDev[0];
02605   PAN_PHWAvgDevV = phwAvgDev[1];
02606   delete [] phwAvgDev;
02607   return true;
02608 }

Bool_t MadCluAnalysis::FillPHWCluID ( Int_t  event  ) 

Definition at line 2557 of file MadCluAnalysis.cxx.

References PAN_PHWCluIDU, PAN_PHWCluIDV, and PHWCluID().

Referenced by CallUserFunctions().

02558 {
02559   Double_t *phwCluID = PHWCluID(event);
02560   PAN_PHWCluIDU = phwCluID[0];
02561   PAN_PHWCluIDV = phwCluID[1];
02562   delete [] phwCluID;
02563   return true;
02564 }

Bool_t MadCluAnalysis::FillPHWProbEM ( Int_t  event  ) 

Definition at line 2513 of file MadCluAnalysis.cxx.

References PAN_PHWProbEMU, PAN_PHWProbEMV, and PHWProbEM().

Referenced by CallUserFunctions().

02514 {
02515   Double_t *phwProbEM = PHWProbEM(event);
02516   PAN_PHWProbEMU = phwProbEM[0];
02517   PAN_PHWProbEMV = phwProbEM[1];
02518   delete [] phwProbEM;
02519   return true;
02520 }

Double_t MadCluAnalysis::FOM (  ) 

Definition at line 1860 of file MadCluAnalysis.cxx.

References NtpSRVertex::dcosz, NtpMCTruth::emfrac, MadBase::eventSummary, MadBase::GetEntry(), NtpSRStripPulseHeight::gev, NtpMCTruth::iaction, NtpSREvent::index, NtpMCTruth::inu, NtpMCTruth::inunoosc, NtpMCTruth::iresonance, MadQuantities::IsFidVtxEvt(), MadBase::LoadEvent(), MadBase::LoadTruthForRecoTH(), MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, NtpSRShower::nstrip, MadBase::ntpEvent, MadBase::ntpShower, MadBase::ntpTruth, PassBasicCuts(), NtpSRShower::ph, PHWAvgDev(), PHWCluID(), PHWProbEM(), NtpSRShower::vtx, and NtpMCTruth::y.

01861 {
01862 
01863   double NueCCn = 0.;
01864   double NueNCn = 0.;
01865   double NueBeamn = 0.;
01866   double NumuCCn = 0.;
01867   double NumuNCn = 0.;
01868   double NutauCCn = 0.;
01869   double NutauNCn = 0.;
01870   double NueQEn = 0.;
01871   double Unknown = 0;
01872 
01873   TH1F *totHist = new TH1F("totHist","totHist",10,0,10);
01874   totHist->SetLineColor(1);
01875   TH1F *beamnueHist = new TH1F("beamnueHist","beamnueHist",10,0,10);
01876   beamnueHist->SetLineColor(8);
01877   TH1F *nueCCHist = new TH1F("nueCCHist","nueCCHist",10,0,10);
01878   nueCCHist->SetLineColor(3);
01879   TH1F *numuCCHist = new TH1F("numuCCHist","numuCCHist",10,0,10);
01880   numuCCHist->SetLineColor(4);
01881   TH1F *nutauCCHist = new TH1F("nutauCCHist","nutauCCHist",10,0,10);
01882   nutauCCHist->SetLineColor(6);
01883   TH1F *NCHist = new TH1F("NCHist","NCHist",10,0,10);
01884   NCHist->SetLineColor(2);
01885 
01886   TH2F *totHist_nstp = new TH2F("totHist_nstp","totHist_nstp",10,0,10,300,0,300);
01887   totHist_nstp->SetLineColor(1);
01888   TH2F *beamnueHist_nstp = new TH2F("beamnueHist_nstp","beamnueHist_nstp",10,0,10,300,0,300);
01889   beamnueHist_nstp->SetLineColor(8);
01890   TH2F *nueCCHist_nstp = new TH2F("nueCCHist_nstp","nueCCHist_nstp",10,0,10,300,0,300);
01891   nueCCHist_nstp->SetLineColor(3);
01892   TH2F *numuCCHist_nstp = new TH2F("numuCCHist_nstp","numuCCHist_nstp",10,0,10,300,0,300);
01893   numuCCHist_nstp->SetLineColor(4);
01894   TH2F *nutauCCHist_nstp = new TH2F("nutauCCHist_nstp","nutauCCHist_nstp",10,0,10,300,0,300);
01895   nutauCCHist_nstp->SetLineColor(6);
01896   TH2F *NCHist_nstp = new TH2F("NCHist_nstp","NCHist_nstp",10,0,10,300,0,300);
01897   NCHist_nstp->SetLineColor(2);
01898 
01899   TH1F *totHist_dcz = new TH1F("totHist_dcz","totHist_dcz",100,0,1);
01900   totHist_dcz->SetLineColor(1);
01901   TH1F *beamnueHist_dcz = new TH1F("beamnueHist_dcz","beamnueHist_dcz",100,0,1);
01902   beamnueHist_dcz->SetLineColor(8);
01903   TH1F *nueCCHist_dcz = new TH1F("nueCCHist_dcz","nueCCHist_dcz",100,0,1);
01904   nueCCHist_dcz->SetLineColor(3);
01905   TH1F *numuCCHist_dcz = new TH1F("numuCCHist_dcz","numuCCHist_dcz",100,0,1);
01906   numuCCHist_dcz->SetLineColor(4);
01907   TH1F *nutauCCHist_dcz = new TH1F("nutauCCHist_dcz","nutauCCHist_dcz",100,0,1);
01908   nutauCCHist_dcz->SetLineColor(6);
01909   TH1F *NCHist_dcz = new TH1F("NCHist_dcz","NCHist_dcz",100,0,1);
01910   NCHist_dcz->SetLineColor(2);
01911 
01912 
01913   TH2F *totHist_phw = new TH2F("totHist_phw","totHist_phw",30,-0.5,5.5,20,0,1);
01914   totHist_phw->SetLineColor(1);
01915   TH2F *beamnueHist_phw = new TH2F("beamnueHist_phw","beamnueHist_phw",30,-0.5,5.5,20,0,1);
01916   beamnueHist_phw->SetLineColor(8);
01917   TH2F *nueCCHist_phw = new TH2F("nueCCHist_phw","nueCCHist_phw",30,-0.5,5.5,20,0,1);
01918   nueCCHist_phw->SetLineColor(3);
01919   TH2F *numuCCHist_phw = new TH2F("numuCCHist_phw","numuCCHist_phw",30,-0.5,5.5,20,0,1);
01920   numuCCHist_phw->SetLineColor(4);
01921   TH2F *nutauCCHist_phw = new TH2F("nutauCCHist_phw","nutauCCHist_phw",30,-0.5,5.5,20,0,1);
01922   nutauCCHist_phw->SetLineColor(6);
01923   TH2F *NCHist_phw = new TH2F("NCHist_phw","NCHist_phw",30,-0.5,5.5,20,0,1);
01924   NCHist_phw->SetLineColor(2);
01925   
01926   TH2F *highyCCHist_phw = new TH2F("highyHist_phw","highyHist_phw",
01927                                    30,-0.5,5.5,20,0,1);
01928   TH2F *midyCCHist_phw = new TH2F("midyHist_phw","midyHist_phw",
01929                                   30,-0.5,5.5,20,0,1);
01930   TH2F *lowyCCHist_phw = new TH2F("lowyHist_phw","lowyHist_phw",
01931                                   30,-0.5,5.5,20,0,1);
01932   TH2F *qenueCCHist_phw = new TH2F("qenueCCHist_phw","qenueCCHist_phw",
01933                                    30,-0.5,5.5,20,0,1);
01934 
01935   TH1F *totHist_avgdev = new TH1F("totHist_avgdev","totHist_avgdev",100,0,0.2);
01936   totHist_avgdev->SetLineColor(1);
01937   TH1F *beamnueHist_avgdev = new TH1F("beamnueHist_avgdev","beamnueHist_avgdev",100,0,0.2);
01938   beamnueHist_avgdev->SetLineColor(8);
01939   TH1F *nueCCHist_avgdev = new TH1F("nueCCHist_avgdev","nueCCHist_avgdev",100,0,0.2);
01940   nueCCHist_avgdev->SetLineColor(3);
01941   TH1F *numuCCHist_avgdev = new TH1F("numuCCHist_avgdev","numuCCHist_avgdev",100,0,0.2);
01942   numuCCHist_avgdev->SetLineColor(4);
01943   TH1F *nutauCCHist_avgdev = new TH1F("nutauCCHist_avgdev","nutauCCHist_avgdev",100,0,0.2);
01944   nutauCCHist_avgdev->SetLineColor(6);
01945   TH1F *NCHist_avgdev = new TH1F("NCHist_avgdev","NCHist_avgdev",100,0,0.2);
01946   NCHist_avgdev->SetLineColor(2);
01947   
01948   TH1F *highyCCHist_avgdev = new TH1F("highyHist_avgdev","highyHist_avgdev",
01949                                    100,0,0.2);
01950   highyCCHist_avgdev->SetLineColor(5);
01951   TH1F *midyCCHist_avgdev = new TH1F("midyHist_avgdev","midyHist_avgdev",
01952                                   100,0,0.2);
01953   midyCCHist_avgdev->SetLineColor(9);
01954   TH1F *lowyCCHist_avgdev = new TH1F("lowyHist_avgdev","lowyHist_avgdev",
01955                                   100,0,0.2);
01956   lowyCCHist_avgdev->SetLineColor(3);
01957   TH1F *qenueCCHist_avgdev = new TH1F("qenueCCHist_avgdev",
01958                                       "qenueCCHist_avgdev",
01959                                       100,0,0.2);
01960   qenueCCHist_avgdev->SetLineColor(4);
01961 
01962   
01963   TH2F *totHist_ratio = new TH2F("totHist_ratio","totHist_ratio",20,0,5,20,0,1);
01964   totHist_ratio->SetLineColor(1);
01965   TH2F *beamnueHist_ratio = new TH2F("beamnueHist_ratio","beamnueHist_ratio",20,0,5,20,0,1);
01966   beamnueHist_ratio->SetLineColor(8);
01967   TH2F *nueCCHist_ratio = new TH2F("nueCCHist_ratio","nueCCHist_ratio",20,0,5,20,0,1);
01968   nueCCHist_ratio->SetLineColor(3);
01969   TH2F *numuCCHist_ratio = new TH2F("numuCCHist_ratio","numuCCHist_ratio",20,0,5,20,0,1);
01970   numuCCHist_ratio->SetLineColor(4);
01971   TH2F *nutauCCHist_ratio = new TH2F("nutauCCHist_ratio","nutauCCHist_ratio",20,0,5,20,0,1);
01972   nutauCCHist_ratio->SetLineColor(6);
01973   TH2F *NCHist_ratio = new TH2F("NCHist_ratio","NCHist_ratio",20,0,5,20,0,1);
01974   NCHist_ratio->SetLineColor(2);
01975   
01976   TH2F *highyCCHist_ratio = new TH2F("highyHist_ratio","highyHist_ratio",
01977                                      20,0,5,20,0,1);
01978   TH2F *midyCCHist_ratio = new TH2F("midyHist_ratio","midyHist_ratio",
01979                                     20,0,5,20,0,1);
01980   TH2F *lowyCCHist_ratio = new TH2F("lowyHist_ratio","lowyHist_ratio",
01981                                     20,0,5,20,0,1);
01982   TH2F *qenueCCHist_ratio = new TH2F("qenueCCHist_ratio","qenueCCHist_ratio",
01983                                      20,0,5,20,0,1);
01984 
01985   Double_t nNuMuFiles = 10;
01986   Double_t nNuEFiles = 20;
01987   Double_t nNuTauFiles = 10;
01988   Double_t POT = 9.25;
01989 
01990   for(int i=0;i<Nentries;i++){
01991     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
01992                             << "% done" << std::endl;
01993     
01994     if(!this->GetEntry(i)) continue;
01995     if(eventSummary->nevent==0) continue;
01996     
01997     //fill tree once for each reconstructed event
01998     for(int j=0;j<eventSummary->nevent;j++){ 
01999       if(!LoadEvent(j)) continue; //no event found
02000       if(ntpEvent->nshower==0) continue; //no shower found
02001       //if(ntpEvent->ntrack!=0) continue; //no shower found
02002       if(!IsFidVtxEvt(ntpEvent->index)) continue;
02003       
02004       if(PassBasicCuts(ntpEvent->index)){
02005         double prob = 1.;
02006         Int_t mcevent = -1;
02007         if(LoadTruthForRecoTH(ntpEvent->index,mcevent)) {
02008           //prob = Oscillate(ntpTruth->inu,ntpTruth->inunoosc,
02009           //       ntpTruth->p4neu[3],735.,0.0025,
02010           //       TMath::ASin(1.0)/2.,0.01);
02011           //prob = Oscillate(ntpTruth->inu,ntpTruth->inunoosc,
02012           //   ntpTruth->p4neu[3],735.,0.002175,
02013           //   TMath::ASin(TMath::Sqrt(0.925))/2.,
02014           //   0.01);
02015         }
02016         //double emfracU = 0;
02017         //double emfracV = 0;     
02018         //double emfracprob = RecoEMFrac(ntpShower->index,0,emfracU,emfracV);
02019         //if(emfracprob>0.75){
02020         if(ntpTruth->emfrac>0.6){
02021         //if(PassAnalysisCuts(ntpEvent->index)){
02022           if( ntpShower->ph.gev>8. || ntpShower->ph.gev<2.) continue;
02023           if(mcevent>=0) {
02024             Double_t *phwcluid = this->PHWCluID(ntpEvent->index);
02025             Double_t *phwprobem = this->PHWProbEM(ntpEvent->index);
02026             Double_t rat_cluid = phwcluid[0]+phwcluid[1];
02027             Double_t rat_probem = 0;
02028             if(phwprobem[0]>0 && phwprobem[1]>0) {
02029               if(phwprobem[0]>phwprobem[1]) 
02030                 rat_probem = phwprobem[1]/phwprobem[0];
02031               else rat_probem = phwprobem[0]/phwprobem[1];
02032             }
02033             Double_t *avgdev = this->PHWAvgDev(ntpEvent->index);
02034             double w1 = 0;
02035             if(ntpTruth->iaction==1){
02036               if(abs(ntpTruth->inunoosc)==12 && abs(ntpTruth->inu)==12){
02037                 w1 = prob*(POT/((nNuEFiles+nNuMuFiles)*6.5));
02038                 NueBeamn+=w1;
02039                 beamnueHist->Fill(ntpShower->ph.gev,w1);
02040                 beamnueHist_nstp->Fill(ntpShower->ph.gev,ntpShower->nstrip,w1);
02041                 beamnueHist_dcz->Fill(ntpShower->vtx.dcosz,w1);
02042                 beamnueHist_phw->Fill(phwcluid[0],phwprobem[0],w1);
02043                 if(phwprobem[0]<0.05) beamnueHist_avgdev->Fill(avgdev[0],w1);
02044                 beamnueHist_ratio->Fill(rat_cluid,rat_probem,w1);
02045               }
02046               else if(abs(ntpTruth->inu)==12) {
02047                 w1 = prob*(POT/(nNuEFiles*6.5));
02048                 NueCCn+=w1;
02049                 nueCCHist->Fill(ntpShower->ph.gev,w1);
02050                 nueCCHist_nstp->Fill(ntpShower->ph.gev,ntpShower->nstrip,w1);
02051                 nueCCHist_dcz->Fill(ntpShower->vtx.dcosz,w1);
02052                 nueCCHist_phw->Fill(phwcluid[0],phwprobem[0],w1);
02053                 if(phwprobem[0]<0.05) nueCCHist_avgdev->Fill(avgdev[0],w1);
02054                 nueCCHist_ratio->Fill(rat_cluid,rat_probem,w1);
02055                 if(ntpTruth->iresonance==1001) {
02056                   NueQEn+=w1;
02057                   qenueCCHist_phw->Fill(phwcluid[0],phwprobem[0],w1);
02058                   if(phwprobem[0]<0.05) qenueCCHist_avgdev->Fill(avgdev[0],w1);
02059                   qenueCCHist_ratio->Fill(rat_cluid,rat_probem,w1);
02060                 }
02061               }
02062               else if(abs(ntpTruth->inu)==14) {
02063                 w1 = prob*(POT/(nNuMuFiles*6.5));
02064                 NumuCCn+=w1;
02065                 numuCCHist->Fill(ntpShower->ph.gev,w1);
02066                 numuCCHist_nstp->Fill(ntpShower->ph.gev,ntpShower->nstrip,w1);
02067                 numuCCHist_dcz->Fill(ntpShower->vtx.dcosz,w1);
02068                 numuCCHist_phw->Fill(phwcluid[0],phwprobem[0],w1);
02069                 if(phwprobem[0]<0.05) numuCCHist_avgdev->Fill(avgdev[0],w1);
02070                 numuCCHist_ratio->Fill(rat_cluid,rat_probem,w1);
02071               }
02072               else if(abs(ntpTruth->inu)==16) {
02073                 w1 = prob*(POT/(nNuTauFiles*6.5));
02074                 NutauCCn+=w1;
02075                 nutauCCHist->Fill(ntpShower->ph.gev,w1);
02076                 nutauCCHist_nstp->Fill(ntpShower->ph.gev,ntpShower->nstrip,w1);
02077                 nutauCCHist_dcz->Fill(ntpShower->vtx.dcosz,w1);
02078                 nutauCCHist_phw->Fill(phwcluid[0],phwprobem[0],w1);
02079                 if(phwprobem[0]<0.05) nutauCCHist_avgdev->Fill(avgdev[0],w1);
02080                 nutauCCHist_ratio->Fill(rat_cluid,rat_probem,w1);
02081               }
02082               if(ntpTruth->y<0.3) {
02083                 lowyCCHist_phw->Fill(phwcluid[0],phwprobem[0],w1);
02084                 if(phwprobem[0]<0.05) lowyCCHist_avgdev->Fill(avgdev[0],w1);
02085                 lowyCCHist_ratio->Fill(rat_cluid,rat_probem,w1);
02086               }
02087               else if(ntpTruth->y<0.6) {
02088                 midyCCHist_phw->Fill(phwcluid[0],phwprobem[0],w1);
02089                 if(phwprobem[0]<0.05) midyCCHist_avgdev->Fill(avgdev[0],w1);
02090                 midyCCHist_ratio->Fill(rat_cluid,rat_probem,w1);
02091               }
02092               else {
02093                 highyCCHist_phw->Fill(phwcluid[0],phwprobem[0],w1);
02094                 if(phwprobem[0]<0.05) highyCCHist_avgdev->Fill(avgdev[0],w1);
02095                 highyCCHist_ratio->Fill(rat_cluid,rat_probem,w1);
02096               }
02097             }
02098             else if(ntpTruth->iaction==0){
02099               if(abs(ntpTruth->inu)==12) {
02100                 w1 = prob*(POT/((nNuEFiles)*6.5));
02101                 NueNCn+=w1;
02102               }
02103               else if(abs(ntpTruth->inu)==14) {
02104                 w1 = prob*(POT/((nNuMuFiles)*6.5));
02105                 NumuNCn+=w1;
02106               }
02107               else if(abs(ntpTruth->inu)==16) {
02108                 w1 = prob*(POT/((nNuTauFiles)*6.5));
02109                 NutauNCn+=w1;
02110               }
02111               NCHist->Fill(ntpShower->ph.gev,w1);
02112               NCHist_nstp->Fill(ntpShower->ph.gev,ntpShower->nstrip,w1);
02113               NCHist_dcz->Fill(ntpShower->vtx.dcosz,w1);
02114               NCHist_phw->Fill(phwcluid[0],phwprobem[0],w1);
02115               if(phwprobem[0]<0.05) NCHist_avgdev->Fill(avgdev[0],w1);
02116               NCHist_ratio->Fill(rat_cluid,rat_probem,w1);
02117             }
02118             totHist->Fill(ntpShower->ph.gev,w1);
02119             totHist_nstp->Fill(ntpShower->ph.gev,ntpShower->nstrip,w1);
02120             totHist_dcz->Fill(ntpShower->vtx.dcosz,w1);
02121             totHist_phw->Fill(phwcluid[0],phwprobem[0],w1);
02122             if(phwprobem[0]<0.05) totHist_avgdev->Fill(avgdev[0],w1);
02123             totHist_ratio->Fill(rat_cluid,rat_probem,w1);
02124             
02125             delete [] phwprobem;
02126             delete [] phwcluid;
02127             delete [] avgdev;
02128           }
02129           else Unknown+=prob*(POT/((nNuEFiles+nNuMuFiles+nNuTauFiles)*6.5));
02130         }
02131       }
02132     }
02133   }
02134   
02135   Double_t fom = NueCCn/sqrt(NueBeamn+NumuCCn+NutauCCn+NueNCn+NumuNCn+NutauNCn+Unknown);
02136   cout<<"NueBeamn "<<NueBeamn<<" NueCCn "<<NueCCn<<" NumuCCn "<<NumuCCn<<" NutauCCn "<<NutauCCn
02137       <<" NueNCn "<<NueNCn<<" NumuNCn "<<NumuNCn<<" NutauNCn "<<NutauNCn<<" Unknown "<<Unknown 
02138       << " fom " << fom<<" NueQEn "<<NueQEn<<endl;
02139 
02140   TCanvas *can = new TCanvas();
02141   can->Divide(2,2);
02142   can->cd(1);
02143   
02144   totHist->Draw();
02145   NCHist->Draw("sames");
02146   numuCCHist->Draw("sames");
02147   nutauCCHist->Draw("sames");
02148   nueCCHist->Draw("sames");
02149   beamnueHist->Draw("sames");
02150   can->cd(2);
02151   totHist_dcz->Draw();
02152   NCHist_dcz->Draw("sames");
02153   numuCCHist_dcz->Draw("sames");
02154   nutauCCHist_dcz->Draw("sames");
02155   nueCCHist_dcz->Draw("sames");
02156   beamnueHist_dcz->Draw("sames");
02157   can->cd(3);
02158   NCHist_nstp->Draw("COLZ");
02159   can->cd(4);
02160   nueCCHist_nstp->Draw("COLZ");
02161 
02162   TCanvas *can2 = new TCanvas();
02163   can2->Divide(2,3);
02164   can2->cd(1);
02165   totHist_phw->Draw("COLZ");
02166   can2->cd(2);
02167   numuCCHist_phw->Draw("COLZ");
02168   can2->cd(3);
02169   nutauCCHist_phw->Draw("COLZ");
02170   can2->cd(4);
02171   nueCCHist_phw->Draw("COLZ");
02172   can2->cd(5);
02173   beamnueHist_phw->Draw("COLZ");
02174   can2->cd(6);
02175   NCHist_phw->Draw("COLZ");
02176 
02177   TCanvas *can3 = new TCanvas();
02178   can3->Divide(2,2);
02179   can3->cd(1);
02180   lowyCCHist_phw->Draw("LEGO2");
02181   can3->cd(2);
02182   midyCCHist_phw->Draw("LEGO2");
02183   can3->cd(3);
02184   highyCCHist_phw->Draw("LEGO2");
02185   can3->cd(4);
02186   qenueCCHist_phw->Draw("LEGO2");
02187 
02188   TCanvas *can4 = new TCanvas();
02189   can4->Divide(1,2);
02190   can4->cd(1);
02191   totHist_avgdev->Draw();
02192   numuCCHist_avgdev->Draw("sames");
02193   nutauCCHist_avgdev->Draw("sames");
02194   nueCCHist_avgdev->Draw("sames");
02195   beamnueHist_avgdev->Draw("sames");
02196   NCHist_avgdev->Draw("sames");
02197   can4->cd(2);
02198   lowyCCHist_avgdev->Draw();
02199   qenueCCHist_avgdev->Draw("sames");
02200   midyCCHist_avgdev->Draw("sames");
02201   highyCCHist_avgdev->Draw("sames");
02202   NCHist_avgdev->Draw("sames");
02203 
02204   TCanvas *can5 = new TCanvas();
02205   can5->Divide(2,3);
02206   can5->cd(1);
02207   totHist_ratio->Draw("COLZ");
02208   can5->cd(2);
02209   numuCCHist_ratio->Draw("COLZ");
02210   can5->cd(3);
02211   nutauCCHist_ratio->Draw("COLZ");
02212   can5->cd(4);
02213   nueCCHist_ratio->Draw("COLZ");
02214   can5->cd(5);
02215   beamnueHist_ratio->Draw("COLZ");
02216   can5->cd(6);
02217   NCHist_ratio->Draw("COLZ");
02218 
02219   TCanvas *can6 = new TCanvas();
02220   can6->Divide(2,2);
02221   can6->cd(1);
02222   lowyCCHist_ratio->Draw("LEGO2");
02223   can6->cd(2);
02224   midyCCHist_ratio->Draw("LEGO2");
02225   can6->cd(3);
02226   highyCCHist_ratio->Draw("LEGO2");
02227   can6->cd(4);
02228   qenueCCHist_ratio->Draw("LEGO2");
02229 
02230   return fom;
02231 }

Float_t MadCluAnalysis::GetWeight ( Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 2471 of file MadCluAnalysis.cxx.

02472 {
02473   return 1;
02474 }

Bool_t MadCluAnalysis::IsNueFidVtxEvt ( Int_t  event  ) 

Definition at line 2908 of file MadCluAnalysis.cxx.

References VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kFar, Detector::kNear, MadBase::LoadEvent(), MadBase::ntpEvent, MadBase::ntpHeader, NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by CallUserFunctions().

02908                                                {
02909   if(!LoadEvent(ievt)) return false;
02910   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
02911     if(ntpEvent->vtx.z<1 || ntpEvent->vtx.z>29 ||   //ends
02912        (ntpEvent->vtx.z<17&&ntpEvent->vtx.z>14) ||  //between SMs
02913        sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
02914             +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.87) return false;
02915   }
02916   else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
02917     if(ntpEvent->vtx.z<1 || ntpEvent->vtx.z>5 ||
02918        sqrt(((ntpEvent->vtx.x-1.4885)*(ntpEvent->vtx.x-1.4885)) +
02919             ((ntpEvent->vtx.y-0.1397)*(ntpEvent->vtx.y-0.1397)))>1) return false;
02920   }  
02921   return true;
02922 }

Bool_t MadCluAnalysis::IsTrackInSlice ( Int_t  slice  ) 

Definition at line 1845 of file MadCluAnalysis.cxx.

References MadBase::eventSummary, NtpSRTrack::index, MadBase::LoadSlice(), MadBase::LoadTrack(), MadBase::ntpTrack, NtpSREventSummary::ntrack, and NtpSRTrack::slc.

Referenced by PassDataCleaningCuts().

01845                                                 {
01846   Int_t cur_track = -1;
01847   if(ntpTrack) cur_track = ntpTrack->index;
01848   if(!LoadSlice(slice)) return false;
01849   for(int i=0;i<eventSummary->ntrack;i++){
01850     if(!LoadTrack(i)) continue;
01851     if(ntpTrack->slc==slice) {
01852       LoadTrack(cur_track);
01853       return true;
01854     }
01855   }
01856   return false;
01857 }

Bool_t MadCluAnalysis::PassAnalysisCuts ( Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 2283 of file MadCluAnalysis.cxx.

References PassBasicCuts(), and PID().

02283                                                   {
02284   if(!PassBasicCuts(event)) return false;
02285   if(PID(event,1)>0. || PID(event,2)>0.) return true;
02286   return false;
02287 }

Bool_t MadCluAnalysis::PassBasicCuts ( Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 2234 of file MadCluAnalysis.cxx.

References NtpSRPlane::beg, NtpSRPlane::end, NtpSRMomentum::eqp, NtpSRTrack::fit, MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), NtpSRTrack::momentum, NtpSRPlane::n, NtpSREvent::nshower, MadBase::ntpEvent, MadBase::ntpShower, MadBase::ntpTrack, NtpSREvent::ntrack, NtpSRFitTrack::pass, NtpSRShower::plane, NtpSRTrack::plane, and NtpSRMomentum::qp.

Referenced by FOM(), and PassAnalysisCuts().

02234                                                {
02235 
02236   if(!LoadEvent(event)) return false;
02237   if(ntpEvent->nshower!=1) return false;
02238   Int_t shower = -1;
02239   if(!LoadLargestShowerFromEvent(event,shower)) return false;
02240 
02241   if(ntpEvent->ntrack==0) return true; 
02242 
02243   Int_t ntrkplane = 0; 
02244   Int_t trkpass = 0;
02245   Double_t trkerror = 0;
02246   Int_t track = -1;
02247   if(LoadLargestTrackFromEvent(event,track)) {
02248     ntrkplane = ntpTrack->plane.end-ntpTrack->plane.beg+1;
02249     if (ntrkplane<ntpTrack->plane.n) ntrkplane = ntpTrack->plane.n;
02250     if(ntpTrack->momentum.qp!=0) 
02251       trkerror = ntpTrack->momentum.eqp/ntpTrack->momentum.qp;
02252     else trkerror = 1;
02253     trkpass = ntpTrack->fit.pass;
02254   }
02255   
02256   double trkshw_overlap = 0.;
02257   if (ntrkplane>0){
02258     if(ntpShower->plane.beg<=ntpTrack->plane.beg &&
02259        ntpShower->plane.end>=ntpTrack->plane.end)
02260       trkshw_overlap = 1.;
02261     else if(ntpShower->plane.beg>ntpTrack->plane.beg &&
02262             ntpShower->plane.end<ntpTrack->plane.end)
02263       trkshw_overlap = double(ntpShower->plane.end - 
02264                               ntpShower->plane.beg+1)/double(ntrkplane);
02265     else if(ntpShower->plane.beg<=ntpTrack->plane.beg &&
02266             ntpShower->plane.end>=ntpTrack->plane.beg)
02267       trkshw_overlap = double(ntpShower->plane.end - 
02268                               ntpTrack->plane.beg+1)/double(ntrkplane);
02269     else if(ntpShower->plane.beg<=ntpTrack->plane.end &&
02270             ntpShower->plane.end>=ntpTrack->plane.end)
02271       trkshw_overlap = double(ntpTrack->plane.end - 
02272                               ntpShower->plane.beg+1)/double(ntrkplane);        
02273   }
02274   
02275   if(trkshw_overlap>0.6 && 
02276      (trkerror>0.4 || trkpass==0 || 
02277       ntrkplane<20)) return true;
02278   
02279   return false;
02280 
02281 }

Bool_t MadCluAnalysis::PassDataCleaningCuts ( Int_t  shower  ) 

Definition at line 1820 of file MadCluAnalysis.cxx.

References NtpSRShower::clu, NtpSRStripPulseHeight::gev, NtpSRCluster::id, IsTrackInSlice(), MadBase::LoadCluster(), MadBase::LoadShower(), NtpSRShower::ncluster, MadBase::ntpCluster, MadBase::ntpShower, NtpSRShower::ph, NtpSRCluster::planeview, and NtpSRShower::slc.

Referenced by DataDistributions().

01821 {
01822   if(!LoadShower(shower)) return false;
01823   if(ntpShower->ph.gev<0.2) return false;
01824   Int_t nphyclusU = 0;
01825   Int_t nphyclusV = 0;
01826   for(int k=0; k<ntpShower->ncluster; k++){
01827     int index = ntpShower->clu[k];
01828     LoadCluster(index);
01829     if(ntpCluster->id==0 || 
01830        ntpCluster->id==1 || 
01831        ntpCluster->id==3) {
01832       if(ntpCluster->planeview==2) nphyclusU += 1;
01833       if(ntpCluster->planeview==3) nphyclusV += 1;
01834     }
01835   }
01836   if(nphyclusU==0&&nphyclusV==0) return 0;
01837   if(nphyclusU==0||nphyclusV==0) {
01838     if(IsTrackInSlice(ntpShower->slc) &&
01839        ntpShower->ph.gev<0.5) return false;
01840   }
01841   return true;
01842 }

Bool_t MadCluAnalysis::PassTruthSignal ( Int_t  mcevent = 0  )  [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 2289 of file MadCluAnalysis.cxx.

References NtpMCTruth::iaction, NtpMCTruth::inu, MadBase::LoadTruth(), and MadBase::ntpTruth.

02290 {
02291   if(!LoadTruth(mcevent)) return false;
02292   if(ntpTruth->inu==12 && ntpTruth->iaction==1) return true;
02293   return false;
02294 }

Double_t * MadCluAnalysis::PHWAvgDev ( Int_t  event  ) 

Definition at line 2566 of file MadCluAnalysis.cxx.

References NtpSRCluster::avgdev, NtpSRShower::clu, NtpSRStripPulseHeight::gev, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), NtpSRShower::ncluster, MadBase::ntpCluster, MadBase::ntpShower, NtpSRCluster::ph, and NtpSRCluster::planeview.

Referenced by FillPHWAvgDev(), and FOM().

02567 {
02568   Double_t *phwAvgDev = new Double_t[2];
02569   phwAvgDev[0] = 0; phwAvgDev[1] = 0;
02570   
02571   if(!LoadEvent(event)) {
02572     phwAvgDev[0] = -10.; phwAvgDev[1] = -10.;
02573     return phwAvgDev;
02574   }
02575   Int_t shower = -1;
02576   if(!LoadLargestShowerFromEvent(event,shower)) {
02577     phwAvgDev[0] = -10.; phwAvgDev[1] = -10.;
02578     return phwAvgDev;
02579   }
02580 
02581   Double_t totPH[2] = {0};
02582   Int_t *clusters = ntpShower->clu;
02583   for(int k=0; k<ntpShower->ncluster; k++){
02584     Int_t index = clusters[k];
02585     if(!LoadCluster(index)) continue;
02586     if(ntpCluster->planeview==2){
02587       phwAvgDev[0] += ntpCluster->ph.gev*ntpCluster->avgdev;
02588       totPH[0] += ntpCluster->ph.gev;
02589     }
02590     else if(ntpCluster->planeview==3){
02591       phwAvgDev[1] += ntpCluster->ph.gev*ntpCluster->avgdev;
02592       totPH[1] += ntpCluster->ph.gev;
02593     }
02594   }
02595   phwAvgDev[0]/=totPH[0];
02596   phwAvgDev[1]/=totPH[1];
02597   return phwAvgDev;
02598 
02599 }

Double_t * MadCluAnalysis::PHWCluID ( Int_t  event  ) 

Definition at line 2522 of file MadCluAnalysis.cxx.

References NtpSRShower::clu, NtpSRStripPulseHeight::gev, NtpSRCluster::id, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), NtpSRShower::ncluster, MadBase::ntpCluster, MadBase::ntpShower, NtpSRCluster::ph, and NtpSRCluster::planeview.

Referenced by FillPHWCluID(), and FOM().

02523 {
02524 
02525   Double_t *phwCluID = new Double_t[2];
02526   phwCluID[0] = 0; phwCluID[1] = 0;
02527   
02528   if(!LoadEvent(event)) {
02529     phwCluID[0] = -10.; phwCluID[1] = -10.;
02530     return phwCluID;
02531   }
02532   Int_t shower = -1;
02533   if(!LoadLargestShowerFromEvent(event,shower)) {
02534     phwCluID[0] = -10.; phwCluID[1] = -10.;
02535     return phwCluID;
02536   }
02537 
02538   Double_t totPH[2] = {0};
02539   Int_t *clusters = ntpShower->clu;
02540   for(int k=0; k<ntpShower->ncluster; k++){
02541     Int_t index = clusters[k];
02542     if(!LoadCluster(index)) continue;
02543     if(ntpCluster->planeview==2){
02544       phwCluID[0] += ntpCluster->ph.gev*ntpCluster->id;
02545       totPH[0] += ntpCluster->ph.gev;
02546     }
02547     else if(ntpCluster->planeview==3){
02548       phwCluID[1] += ntpCluster->ph.gev*ntpCluster->id;
02549       totPH[1] += ntpCluster->ph.gev;
02550     }
02551   }
02552   phwCluID[0]/=totPH[0];
02553   phwCluID[1]/=totPH[1];
02554   return phwCluID;
02555 }

Double_t * MadCluAnalysis::PHWProbEM ( Int_t  event  ) 

Definition at line 2477 of file MadCluAnalysis.cxx.

References NtpSRShower::clu, NtpSRStripPulseHeight::gev, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), NtpSRShower::ncluster, MadBase::ntpCluster, MadBase::ntpShower, NtpSRCluster::ph, NtpSRCluster::planeview, and NtpSRCluster::probem.

Referenced by FillPHWProbEM(), and FOM().

02478 {
02479   Double_t *phwProbEM = new Double_t[2];
02480   phwProbEM[0] = 0; phwProbEM[1] = 0;
02481   
02482   if(!LoadEvent(event)) {
02483     phwProbEM[0] = -10.; phwProbEM[1] = -10.;
02484     return phwProbEM;
02485   }
02486   Int_t shower = -1;
02487   if(!LoadLargestShowerFromEvent(event,shower)) {
02488     phwProbEM[0] = -10.; phwProbEM[1] = -10.;
02489     return phwProbEM;
02490   }
02491 
02492   Double_t totPH[2] = {0};
02493   Int_t *clusters = ntpShower->clu;
02494   for(int k=0; k<ntpShower->ncluster; k++){
02495     Int_t index = clusters[k];
02496     if(!LoadCluster(index)) continue;
02497     if(ntpCluster->planeview==2){
02498       if(ntpCluster->ph.gev>0.7) 
02499         phwProbEM[0] += ntpCluster->ph.gev*ntpCluster->probem;
02500       totPH[0] += ntpCluster->ph.gev;
02501     }
02502     else if(ntpCluster->planeview==3){
02503       if(ntpCluster->ph.gev>0.7) 
02504         phwProbEM[1] += ntpCluster->ph.gev*ntpCluster->probem;
02505       totPH[1] += ntpCluster->ph.gev;
02506     }
02507   }
02508   phwProbEM[0]/=totPH[0];
02509   phwProbEM[1]/=totPH[1];
02510   return phwProbEM;
02511 }

Float_t MadCluAnalysis::PID ( Int_t  event = 0,
Int_t  method = 0 
) [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 2296 of file MadCluAnalysis.cxx.

References NtpSRCluster::avgdev, NtpSRCluster::begplane, NtpSRPlane::begu, NtpSRPlane::begv, NtpSRShower::clu, flikelihoodDists, fneuralNet, NtpSRStripPulseHeight::gev, NtpSRCluster::id, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), NtpSRShower::ncluster, MadBase::ntpCluster, MadBase::ntpShower, PAN_ChargeFracRMSU, PAN_ChargeFracRMSV, PAN_PHWAvgDevU, PAN_PHWAvgDevV, PAN_PHWCluIDU, PAN_PHWCluIDV, PAN_PHWProbEMU, PAN_PHWProbEMV, NtpSRCluster::ph, NtpSRShower::plane, NtpSRCluster::planeview, and NtpSRCluster::probem.

Referenced by CallUserFunctions(), and PassAnalysisCuts().

02297 {
02298   if(!LoadEvent(event)) return -10.;
02299   Int_t shower = -1;
02300   Int_t track = -1;
02301   LoadLargestTrackFromEvent(event,track);
02302   if(!LoadLargestShowerFromEvent(event,shower)) return -10.;
02303 
02304   if(method==0){
02305     int goodU = 0;
02306     int goodV = 0;
02307     double maxU = 0;
02308     double maxV = 0;
02309     double avgdevU = 1.;
02310     double avgdevV = 1.;
02311     Int_t nemclu = 0;
02312     Int_t nhadclu = 0;
02313     Int_t nphysclu = 0;
02314     
02315     Int_t *clusters = ntpShower->clu;
02316     for(int k=0; k<ntpShower->ncluster; k++){
02317       Int_t index = clusters[k];
02318       if(!LoadCluster(index)) continue;
02319       if(ntpCluster->id==0) nemclu++;
02320       else if(ntpCluster->id==1) nhadclu++;
02321       else if(ntpCluster->id!=2 && ntpCluster->ph.gev>0.5) nphysclu++;
02322       if(ntpCluster->planeview==2){
02323         if(ntpCluster->ph.gev>maxU){
02324           if(ntpCluster->id==0 && 
02325              TMath::Abs(ntpCluster->begplane-ntpShower->plane.begu)<=4){
02326             goodU = 1;
02327             if(ntpCluster->probem>0.2) goodU=3;
02328             avgdevU = ntpCluster->avgdev;
02329             maxU = ntpCluster->ph.gev;
02330           }
02331           else if(ntpCluster->id==1 && ntpCluster->probem>0.5 &&
02332                   TMath::Abs(ntpCluster->begplane-ntpShower->plane.begu)<=4){
02333             goodU = 2;
02334             avgdevU = ntpCluster->avgdev;
02335             maxU = ntpCluster->ph.gev;
02336           }
02337           else {
02338             goodU = 0;
02339             avgdevU = ntpCluster->avgdev;
02340             maxU = ntpCluster->ph.gev;
02341           }
02342         }
02343       }
02344       else if(ntpCluster->planeview==3){
02345         if(ntpCluster->ph.gev>maxV){
02346           if(ntpCluster->id==0 &&
02347              TMath::Abs(ntpCluster->begplane-ntpShower->plane.begv)<=4){
02348             goodV = 1;
02349             if(ntpCluster->probem>0.2) goodV=3;
02350             avgdevV = ntpCluster->avgdev;
02351             maxV = ntpCluster->ph.gev;
02352           }
02353           else if(ntpCluster->id==1 && ntpCluster->probem>0.5 &&
02354                   TMath::Abs(ntpCluster->begplane-ntpShower->plane.begv)<=4){
02355             goodV = 2;
02356             avgdevV = ntpCluster->avgdev;
02357             maxV = ntpCluster->ph.gev;
02358           }
02359           else {
02360             goodV = 0;
02361             avgdevV = ntpCluster->avgdev;
02362             maxV = ntpCluster->ph.gev;
02363           }
02364         }
02365       }
02366     }
02367     
02368     if( ( (goodU==1 && avgdevU<0.05 && goodV==1 && avgdevV<0.05) || 
02369           goodU==3 || goodV==3 ) && nemclu < 3){
02370       return 1;
02371     }
02372   
02373     return 0.;
02374   }
02375   else if(method==1 || method==2){ 
02376 
02377     Int_t index = 0;
02378     Float_t id = PAN_PHWCluIDU;
02379     Float_t em = PAN_PHWProbEMU;
02380     Float_t cf = PAN_ChargeFracRMSU;
02381     if(method==2) {
02382       index = 6;
02383       id = PAN_PHWCluIDV;
02384       em = PAN_PHWProbEMV;
02385       cf = PAN_ChargeFracRMSV;
02386     }
02387 
02388     Float_t ProbNC = 1.;
02389     Float_t ProbQE = 1.;
02390     Float_t PidParQENC = 0;
02391 
02392     //NC:
02393     int bin1 = flikelihoodDists[index+1]->FindBin(id);
02394     ProbNC *= flikelihoodDists[index+1]->GetBinContent(bin1);
02395     
02396     bin1 = flikelihoodDists[index+3]->FindBin(em);
02397     ProbNC *= flikelihoodDists[index+3]->GetBinContent(bin1);
02398     
02399     bin1 = flikelihoodDists[index+5]->FindBin(cf);
02400     ProbNC *= flikelihoodDists[index+5]->GetBinContent(bin1);
02401       
02402     //QE:
02403     bin1 = flikelihoodDists[index]->FindBin(id);
02404     ProbQE *= flikelihoodDists[index]->GetBinContent(bin1);
02405     
02406     bin1 = flikelihoodDists[index+2]->FindBin(em);
02407     ProbQE *= flikelihoodDists[index+2]->GetBinContent(bin1);
02408     
02409     bin1 = flikelihoodDists[index+4]->FindBin(cf);
02410     ProbQE *= flikelihoodDists[index+4]->GetBinContent(bin1);
02411     
02412     if(ProbQE!=0&&ProbNC!=0)
02413       PidParQENC = TMath::Sqrt(-TMath::Log(ProbNC)) -
02414         TMath::Sqrt(-TMath::Log(ProbQE));
02415     else if(ProbQE==0 && ProbNC==0) PidParQENC=-10;
02416     else if(ProbQE==0) PidParQENC = TMath::Sqrt(-TMath::Log(ProbNC)) - 10;
02417     else if(ProbNC==0) PidParQENC = 10 - TMath::Sqrt( -TMath::Log(ProbQE));
02418     
02419     return PidParQENC;
02420   }
02421   else if(method==3){
02422     
02423     Double_t params[4];
02424     params[0] = PAN_ChargeFracRMSU*2.;
02425     params[1] = PAN_PHWProbEMU*PAN_PHWProbEMU;
02426     params[2] = PAN_PHWCluIDU;
02427     params[3] = PAN_PHWAvgDevU*2.;
02428     
02429     return fneuralNet->Evaluate(0,params);
02430   }
02431   else if(method==4){
02432 
02433     Double_t params[4];
02434     params[0] = PAN_ChargeFracRMSV*2.;
02435     params[1] = PAN_PHWProbEMV*PAN_PHWProbEMV;
02436     params[2] = PAN_PHWCluIDV;
02437     params[3] = PAN_PHWAvgDevV*2.;
02438     
02439     return fneuralNet->Evaluate(0,params);
02440   }
02441   else if(method==5){
02442 
02443     Double_t params[4];
02444     params[0] = PAN_ChargeFracRMSU + PAN_ChargeFracRMSV;
02445     params[1] = (PAN_PHWProbEMU*PAN_PHWProbEMU + 
02446                  PAN_PHWProbEMV*PAN_PHWProbEMV)/2.;
02447     params[2] = (PAN_PHWCluIDU + PAN_PHWCluIDV)/2.;
02448     params[3] = PAN_PHWAvgDevU + PAN_PHWAvgDevV;
02449     
02450     return fneuralNet->Evaluate(0,params);
02451   }
02452   return 0;
02453 }

void MadCluAnalysis::Plot ( char *  fileName,
Int_t  startEnt = 0 
)

Definition at line 169 of file MadCluAnalysis.cxx.

References NtpSRShower::clu, NtpMCTruth::emfrac, MadBase::eventSummary, VldContext::GetDetector(), MadBase::GetEntry(), RecHeader::GetVldContext(), NtpSRStripPulseHeight::gev, NtpMCTruth::iaction, NtpSRCluster::id, NtpMCStdHep::IdHEP, NtpSREvent::index, MadQuantities::IsFidVtxEvt(), MadBase::isST, NtpMCStdHep::IstHEP, Detector::kFar, Detector::kNear, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadShower(), MadBase::LoadStdHep(), MadBase::LoadTruthForRecoTH(), NtpMCStdHep::mc, MadBase::mcrecord, NtpSRShower::ncluster, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, MadBase::ntpCluster, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpStdHep, MadBase::ntpTruth, NtpSREvent::ntrack, NtpMCStdHep::p4, NtpMCTruth::p4neu, NtpMCTruth::p4shw, NtpSRCluster::ph, NtpSRCluster::planeview, NtpSRCluster::probem, NtpSREvent::shw, NtpStRecord::stdhep, NtpMCRecord::stdhep, MadBase::strecord, and NtpMCTruth::y.

00170 {
00171 
00172   char savename[256];
00173   sprintf(savename,"Clu_%s.root",fileName);
00174   TFile *savefile = new TFile(savename,"RECREATE");
00175 
00176   //all clusters ID
00177   TH1F *IDHist = new TH1F("IDHist","Cluster ID - CC",6,0,5);
00178   IDHist->SetXTitle("Cluster ID");
00179   TH1F *IDHistNC = new TH1F("IDHistNC","Cluster ID - NC",6,0,5);
00180   IDHistNC->SetXTitle("Cluster ID");
00181 
00182   //CC plots with y
00183   TH1F *yDist = new TH1F("yDist","y Distribution for CC Events",101,-0.005,1.005);  
00184   TH1F *EMyDist = new TH1F("EMyDist","Fraction of CC Events with y",101,-0.005,1.005);
00185   EMyDist->SetXTitle("y");
00186   TH1F *HADyDist = new TH1F("HADyDist","Fraction of CC Events with y",101,-0.005,1.005);
00187   HADyDist->SetXTitle("y");
00188   TH1F *EMProbyDist = new TH1F("EMProbyDist","Fraction of CC Events with y",101,-0.005,1.005);
00189   EMProbyDist->SetXTitle("y");
00190   TH1F *UVEMyDist = new TH1F("UVEMyDist","Fraction of CC Events with y",101,-0.005,1.005);
00191   UVEMyDist->SetXTitle("y");
00192   TH1F *UVHADyDist = new TH1F("UVHADyDist","Fraction of CC Events with y",101,-0.005,1.005);
00193   UVHADyDist->SetXTitle("y");
00194   TH1F *UVEMProbyDist = new TH1F("UVEMProbyDist","Fraction of CC Events with y",101,-0.005,1.005);
00195   UVEMProbyDist->SetXTitle("y");
00196 
00197   //NC plots with y
00198   TH1F *yDistNC = new TH1F("yDistNC","y Distribution for NC Events",101,-0.005,1.005); 
00199   TH1F *EMyDistNC = new TH1F("EMyDistNC","Fraction of NC Events with y",101,-0.005,1.005);
00200   EMyDistNC->SetXTitle("y");
00201   TH1F *HADyDistNC = new TH1F("HADyDistNC","Fraction of NC Events with y",
00202                               101,-0.005,1.005);
00203   HADyDistNC->SetXTitle("y");
00204   TH1F *EMProbyDistNC = new TH1F("EMProbyDistNC","Fraction of NC Events with y",101,-0.005,1.005);
00205   EMProbyDistNC->SetXTitle("y");
00206   TH1F *UVEMyDistNC = new TH1F("UVEMyDistNC","Fraction of NC Events with y",
00207                                101,-0.005,1.005);
00208   UVEMyDistNC->SetXTitle("y");
00209   TH1F *UVHADyDistNC = new TH1F("UVHADyDistNC","Fraction of NC Events with y",
00210                                 101,-0.005,1.005);
00211   UVHADyDistNC->SetXTitle("y");
00212   TH1F *UVEMProbyDistNC = new TH1F("UVEMProbyDistNC","Fraction of NC Events with y",
00213                                    101,-0.005,1.005);
00214   UVEMProbyDistNC->SetXTitle("y");
00215 
00216   //primary cluster IDs for NC and CC
00217   TH1F *UClusterIDCC = new TH1F("UClusterIDCC","U View Primary Cluster ID",6,-0.5,5.5);
00218   UClusterIDCC->SetXTitle("Cluster ID");
00219   TH1F *UClusterIDNC = new TH1F("UClusterIDNC","U View Primary Cluster ID",6,-0.5,5.5);
00220   UClusterIDNC->SetXTitle("Cluster ID");
00221   TH1F *VClusterIDCC = new TH1F("VClusterIDCC","V View Primary Cluster ID",6,-0.5,5.5);
00222   VClusterIDCC->SetXTitle("Cluster ID");
00223   TH1F *VClusterIDNC = new TH1F("VClusterIDNC","V View Primary Cluster ID",6,-0.5,5.5);
00224   VClusterIDNC->SetXTitle("Cluster ID");
00225 
00226   //Average #cluster with energy
00227   TH2F *AvgCluEnCC = new TH2F("AvgCluEnCC","Average Number of Clusters with Visible Energy - CC",
00228                       100,0,25,20,0,20);
00229   AvgCluEnCC->SetXTitle("Energy (GeV)");
00230   AvgCluEnCC->SetYTitle("# of Clusters");
00231   TH2F *AvgCluEnNC = new TH2F("AvgCluEnNC","Average Number of Clusters with Visible Energy - NC",
00232                       100,0,25,20,0,20);
00233   AvgCluEnNC->SetXTitle("Energy (GeV)");
00234   AvgCluEnNC->SetYTitle("# of Clusters");
00235   TH2F *AvgCluEnCCnox = new TH2F("AvgCluEnCCnox","Average Number of Clusters with Visible Energy (excl xtalk clu) - CC",
00236                        100,0,25,20,0,20);
00237   AvgCluEnCCnox->SetXTitle("Energy (GeV)");
00238   AvgCluEnCCnox->SetYTitle("# of Clusters");
00239   TH2F *AvgCluEnNCnox = new TH2F("AvgCluEnNCnox","Average Number of Clusters with Visible Energy (excl xtalk clu) - NC",
00240                        100,0,25,20,0,20);
00241   AvgCluEnNCnox->SetXTitle("Energy (GeV)");
00242   AvgCluEnNCnox->SetYTitle("# of Clusters");
00243 
00244   //Average #cluster with #stdhep final states
00245   TH2F *AvgCluSHCC = new TH2F("AvgCluSHCC","Average Number of Clusters Vs #StdHep Final States >0.1GeV - CC",
00246                       25,0,25,20,0,20);
00247   AvgCluSHCC->SetXTitle("#StdHep FS");
00248   AvgCluSHCC->SetYTitle("# of Clusters");
00249   TH2F *AvgCluSHNC = new TH2F("AvgCluSHNC","Average Number of Clusters Vs #StdHep Final States >0.1GeV - NC",
00250                       25,0,25,20,0,20);
00251   AvgCluSHNC->SetXTitle("#StdHep FS");
00252   AvgCluSHNC->SetYTitle("# of Clusters");
00253   TH2F *AvgCluSHCCnox = new TH2F("AvgCluSHCCnox","Average Number of Clusters Vs #StdHep Final States >0.5GeV (excl xtalk clu) - CC",
00254                        25,0,25,20,0,20);
00255   AvgCluSHCCnox->SetXTitle("#StdHep FS");
00256   AvgCluSHCCnox->SetYTitle("# of Clusters");
00257   TH2F *AvgCluSHNCnox = new TH2F("AvgCluSHNCnox","Average Number of Clusters Vs #StdHep Final States >0.5GeV (excl xtalk clu) - NC",
00258                        25,0,25,20,0,20);
00259   AvgCluSHNCnox->SetXTitle("#StdHep FS");
00260   AvgCluSHNCnox->SetYTitle("# of Clusters");
00261 
00262   //Average #cluster with #stdhep final states & energy
00263   TH3F *AvgCluSHEnergyCC = new TH3F("AvgCluSHEnergyCC","Average Number of Clusters Vs #StdHep Final States >0.1GeV vs Energy - CC",
00264                                     100,0,25,25,0,25,50,0,50);
00265   AvgCluSHEnergyCC->SetXTitle("Energy (GeV)");
00266   AvgCluSHEnergyCC->SetYTitle("#StdHep FS");
00267   AvgCluSHEnergyCC->SetZTitle("# of Clusters");
00268   TH3F *AvgCluSHEnergyNC = new TH3F("AvgCluSHEnergyNC","Average Number of Clusters Vs #StdHep Final States >0.1GeV vs Energy - NC",
00269                                     100,0,25,25,0,25,50,0,50);
00270   AvgCluSHEnergyNC->SetXTitle("Energy (GeV)");
00271   AvgCluSHEnergyNC->SetYTitle("#StdHep FS");
00272   AvgCluSHEnergyNC->SetZTitle("# of Clusters");
00273   TH3F *AvgCluSHEnergyCCnox = new TH3F("AvgCluSHEnergyCCnox","Average Number of Clusters Vs #StdHep Final States >0.5GeV (excl xtalk clu) vs Energy - CC",
00274                                        100,0,25,25,0,25,50,0,50);
00275   AvgCluSHEnergyCCnox->SetXTitle("Energy (GeV)");
00276   AvgCluSHEnergyCCnox->SetYTitle("#StdHep FS");
00277   AvgCluSHEnergyCCnox->SetZTitle("# of Clusters");
00278   TH3F *AvgCluSHEnergyNCnox = new TH3F("AvgCluSHEnergyNCnox","Average Number of Clusters Vs #StdHep Final States >0.5GeV (excl xtalk clu) vs Energy - NC",
00279                                        100,0,25,25,0,25,50,0,50);
00280   AvgCluSHEnergyNCnox->SetXTitle("Energy (GeV)");
00281   AvgCluSHEnergyNCnox->SetYTitle("#StdHep FS");
00282   AvgCluSHEnergyNCnox->SetZTitle("# of Clusters");
00283 
00284   //Reco'd EMFrac
00285   TH1F *EMFracCC = new TH1F("EMFracCC","Estimated EMFrac - CC Events",100,0,1);
00286   TH1F *EMFracNC = new TH1F("EMFracNC","Estimated EMFrac - NC Events",100,0,1);
00287 
00288   TH1F *EMFracCCDiff = new TH1F("EMFracCCDiff","(Truth - Estimated)/Truth EMFrac - CC Events",200,-2,2);
00289   TH1F *EMFracNCDiff = new TH1F("EMFracNCDiff","(Truth - Estimated)/Truth EMFrac - NC Events",200,-2,2);
00290 
00291   TH2F *EMFracCCY = new TH2F("EMFracCCY","Estimated EMFrac Vs Y - CC Events",100,0,1,100,0,1);
00292   TH2F *EMFracNCY = new TH2F("EMFracNCY","Estimated EMFrac Vs Y - NC Events",100,0,1,100,0,1);
00293 
00294   //2D reco vs true:
00295   TH2F *EMFrac_Truth_Reco = new TH2F("EMFracTrueReco","Reco vs Truth EMFrac",
00296                                      100,0,1,100,0,1);
00297   EMFrac_Truth_Reco->SetXTitle("True EMFrac");
00298   EMFrac_Truth_Reco->SetYTitle("Reco EMFrac");
00299   TH2F *EMFrac_Truth_Reco_end0 = new TH2F("EMFracTrueReco_end0","Reco vs Truth EMFrac",
00300                                           100,0,1,100,0,1);
00301   EMFrac_Truth_Reco_end0->SetXTitle("True EMFrac");
00302   EMFrac_Truth_Reco_end0->SetYTitle("Reco EMFrac");
00303   TH2F *EMFrac_Truth_Reco_end1 = new TH2F("EMFracTrueReco_end1","Reco vs Truth EMFrac",
00304                                           100,0,1,100,0,1);
00305   EMFrac_Truth_Reco_end1->SetXTitle("True EMFrac");
00306   EMFrac_Truth_Reco_end1->SetYTitle("Reco EMFrac");  
00307 
00308   TH2F *asymHistNC = new TH2F("asymHistNC","asymHistNC",1000,0,10,100,0,1);
00309   TH2F *asymHistCC = new TH2F("asymHistCC","asymHistCC",1000,0,10,100,0,1);
00310 
00311   TH2F *asymHistNC_NotPrime = new TH2F("asymHistNC_NotPrime","asymHistNC_NotPrime",1000,0,10,100,0,1);
00312   TH2F *asymHistCC_NotPrime = new TH2F("asymHistCC_NotPrime","asymHistCC_NotPrime",1000,0,10,100,0,1);
00313 
00314   TH2F *asymHist = new TH2F("asymHist","asymHist",1000,0,10,100,0,1);
00315 
00316   for(int i=startEnt;i<Nentries;i++){
00317     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
00318                             << "% done" << std::endl;
00319 
00320     if(!this->GetEntry(i)) continue;
00321     if(eventSummary->nevent==0) continue;
00322     
00323     Int_t is_fid = 0;
00324     //fill tree once for each reconstructed event
00325     for(int j=0;j<eventSummary->nevent;j++){ 
00326       if(!LoadEvent(j)) continue; //no event found
00327       if(ntpEvent->nshower==0) continue; //no shower found
00328       
00329       //get detector type:
00330       if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00331         //fiducial criteria on vtx for far detector
00332         is_fid = 1;
00333         if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
00334       }
00335       else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00336         //fiducial criteria on vtx for near detector
00337         is_fid = 1;
00338         if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
00339       }
00340       //if(is_fid==0) continue;
00341 
00342       Int_t mcevent = -1;
00343       //check for a corresponding mc event      
00344       if(LoadTruthForRecoTH(ntpEvent->index,mcevent)) {  //loadtruth
00345 
00346         float nBigStdHep = 0;
00347         TClonesArray* pointStdhepArray = NULL;
00348         if(isST) pointStdhepArray = (strecord->stdhep);
00349         else pointStdhepArray = (mcrecord->stdhep);
00350         TClonesArray& stdhepArray = *pointStdhepArray;
00351         Int_t nStdHep = stdhepArray.GetEntries();       
00352         for(int k=0;k<nStdHep;k++){
00353           LoadStdHep(k);
00354           if(ntpStdHep->mc==mcevent) {
00355             if(ntpStdHep->IstHEP==1){ //final state particles
00356               if(ntpStdHep->p4[3]>0.5&&!(abs(ntpStdHep->IdHEP)==12 || 
00357                                          abs(ntpStdHep->IdHEP)==14 || 
00358                                          abs(ntpStdHep->IdHEP)==16)) {
00359                 nBigStdHep+=1;
00360               }
00361             }
00362           }
00363         }
00364       
00365         //Int_t shower = -1;
00366         //if(!LoadLargestShowerFromEvent(ntpEvent->index,shower)) continue;
00367         if(!LoadShower(ntpEvent->shw[0])) continue;
00368 
00369         int goodU = 0;
00370         int goodV = 0;
00371         int bestU = 0;
00372         int bestV = 0;
00373         double enU = 0;
00374         double enV = 0;
00375         double totU = 0;
00376         double totV = 0;
00377         
00378         for(int k=0; k<ntpShower->ncluster; k++){
00379           int index = ntpShower->clu[k];
00380           LoadCluster(index);
00381           if(ntpCluster->planeview==2){
00382             totU+=ntpCluster->ph.gev;
00383             if(ntpCluster->ph.gev>enU) {
00384               enU = ntpCluster->ph.gev;
00385               bestU = index;
00386               if(ntpCluster->id==0) {
00387                 if(ntpCluster->probem>0.2) goodU = 2;
00388                 else goodU = 1;
00389               }
00390               else if(ntpCluster->id==1) goodU = -1;
00391               else goodU = 0;
00392             }
00393           }
00394           else if(ntpCluster->planeview==3){
00395             totV+=ntpCluster->ph.gev;
00396             if(ntpCluster->ph.gev>enV) {
00397               enV = ntpCluster->ph.gev;
00398               bestV = index;
00399               if(ntpCluster->id==0){
00400                 if(ntpCluster->probem>0.2) goodV = 2;
00401                 else goodV = 1;
00402               }
00403               else if(ntpCluster->id==1) goodV = -1;
00404               else goodV = 0;
00405             }
00406           }
00407         }
00408 
00409         double asymmetry = TMath::Abs(totU-totV);
00410         if((totU+totV)>0.) asymmetry /= (totU+totV);
00411         else asymmetry = 1;
00412 
00413         //cout << i << " " << j << " " << totU << " " << totV << " " << asymmetry << endl;
00414 
00415         if(ntpTruth->iaction==1){
00416 
00417           yDist->Fill(ntpTruth->y);
00418           
00419           //if(totU<totV) asymHistCC->Fill(totU,asymmetry);
00420           //else asymHistCC->Fill(totV,asymmetry);
00421           if(asymmetry>0.8){
00422             if(totU<totV) asymHistCC->Fill(totU+totV,totU);
00423             else asymHistCC->Fill(totU+totV,totV);
00424           }
00425 
00426           AvgCluEnCC->Fill(ntpTruth->p4neu[3],ntpShower->ncluster);
00427           //if(ntpTruth->p4neu[3]>2.5 && 
00428           // ntpTruth->p4neu[3]<3.5) 
00429           AvgCluSHCC->Fill(nBigStdHep,ntpShower->ncluster);
00430           AvgCluSHEnergyCC->Fill(ntpTruth->p4neu[3],nBigStdHep,
00431                                  ntpShower->ncluster);
00432 
00433           int cntU = 0;
00434           int cntV = 0;
00435           float EMFrac = 0;
00436           float EMFrac_end0 = 0;
00437           float EMFrac_end1 = 0;
00438           float EnTot = 0;
00439           float EnTot_end0 = 0;
00440           float EnTot_end1 = 0;
00441           
00442           for(int k=0; k<ntpShower->ncluster; k++){
00443             int index = ntpShower->clu[k];
00444             LoadCluster(index);
00445             IDHist->Fill(ntpCluster->id);
00446             if(ntpCluster->id==0 && ntpCluster->probem>0.2) {
00447               EMFrac+=ntpCluster->ph.gev;
00448               EnTot+=ntpCluster->ph.gev;
00449               if(ntpCluster->planeview==2){
00450                 EMFrac_end0+=ntpCluster->ph.gev;
00451                 EnTot_end0+=ntpCluster->ph.gev;
00452               }
00453               if(ntpCluster->planeview==3){
00454                 EMFrac_end1+=ntpCluster->ph.gev;
00455                 EnTot_end1+=ntpCluster->ph.gev;
00456               }
00457             }
00458             else {
00459               EnTot+=ntpCluster->ph.gev;
00460               if(ntpCluster->planeview==2){
00461                 EnTot_end0+=ntpCluster->ph.gev;
00462               }
00463               if(ntpCluster->planeview==3){
00464                 EnTot_end1+=ntpCluster->ph.gev;
00465               }
00466             }
00467             if(ntpCluster->ph.gev>0.25 && 
00468                (ntpCluster->id==0 || ntpCluster->id==1 || ntpCluster->id==3)) {
00469               if(ntpCluster->planeview==2) cntU++; 
00470               if(ntpCluster->planeview==3) cntV++;
00471             }
00472             if (index==bestU) UClusterIDCC->Fill(ntpCluster->id);
00473             if (index==bestV) VClusterIDCC->Fill(ntpCluster->id); 
00474           }
00475         
00476           if(goodU>=1 || goodV>=1) EMyDist->Fill(ntpTruth->y);
00477           if(goodU==-1 || goodV==-1) HADyDist->Fill(ntpTruth->y);
00478           if(goodU==-1 && goodV==-1) UVHADyDist->Fill(ntpTruth->y);
00479           if(goodU>=1 && goodV>=1) UVEMyDist->Fill(ntpTruth->y);
00480           if(goodU==2 || goodV==2) EMProbyDist->Fill(ntpTruth->y);
00481           if(goodU==2 && goodV==2) UVEMProbyDist->Fill(ntpTruth->y);
00482 
00483           AvgCluEnCCnox->Fill(ntpTruth->p4neu[3],(cntU+cntV)/2.);
00484           //if(ntpTruth->p4neu[3]>2.5 && 
00485           // ntpTruth->p4neu[3]<3.5) 
00486           AvgCluSHCCnox->Fill(nBigStdHep,(cntU+cntV)/2.);
00487           AvgCluSHEnergyCCnox->Fill(ntpTruth->p4neu[3],nBigStdHep,
00488                                     (cntU+cntV)/2.);
00489           
00490           if(EnTot_end0!=0 || EnTot_end1!=0){
00491             if(EnTot_end0==0) EMFrac = EMFrac_end1/EnTot_end1;
00492             else if(EnTot_end1==0) EMFrac = EMFrac_end0/EnTot_end0;
00493             else if(EMFrac_end0/EnTot_end0>EMFrac_end1/EnTot_end1) 
00494               EMFrac = EMFrac_end0/EnTot_end0;
00495             else EMFrac = EMFrac_end1/EnTot_end1;
00496             //else EMFrac = (EMFrac_end0/EnTot_end0 + EMFrac_end1/EnTot_end1)/2.;
00497             EnTot = 1;
00498           }
00499           if(EnTot!=0) {
00500             EMFracCC->Fill(EMFrac/EnTot);
00501             EMFracCCY->Fill(ntpTruth->y,EMFrac/EnTot);
00502             if(ntpTruth->emfrac!=0) EMFracCCDiff->Fill((ntpTruth->emfrac - 
00503                                                         (EMFrac/EnTot))/ntpTruth->emfrac);
00504             else if(EMFrac==0) EMFracCCDiff->Fill(0);
00505             else EMFracCCDiff->Fill(2-EMFrac/EnTot);
00506             EMFrac_Truth_Reco->Fill(ntpTruth->emfrac,EMFrac/EnTot);
00507             if(EnTot_end0>0) 
00508               EMFrac_Truth_Reco_end0->Fill(ntpTruth->emfrac,EMFrac_end0/EnTot_end0);
00509             if(EnTot_end1>0) 
00510               EMFrac_Truth_Reco_end1->Fill(ntpTruth->emfrac,EMFrac_end1/EnTot_end1);
00511           }
00512         }
00513         
00514         else if(ntpTruth->iaction==0){
00515           
00516           yDistNC->Fill(ntpTruth->y);
00517 
00518           if(totU<totV) asymHistNC->Fill(totU,asymmetry);
00519           else asymHistNC->Fill(totV,asymmetry);
00520           
00521           AvgCluEnNC->Fill(ntpTruth->p4shw[3],ntpShower->ncluster);
00522           //if(ntpTruth->p4neu[3]>3.0 && 
00523           // ntpTruth->p4neu[3]<3.1) 
00524           AvgCluSHNC->Fill(nBigStdHep,ntpShower->ncluster);
00525           AvgCluSHEnergyNC->Fill(ntpTruth->p4neu[3],nBigStdHep,
00526                                  ntpShower->ncluster);
00527 
00528           int cntU = 0;
00529           int cntV = 0;
00530           float EMFrac = 0;
00531           float EMFrac_end0 = 0;
00532           float EMFrac_end1 = 0;
00533           float EnTot = 0;
00534           float EnTot_end0 = 0;
00535           float EnTot_end1 = 0;
00536           for(int k=0; k<ntpShower->ncluster; k++){
00537             int index = ntpShower->clu[k];
00538             LoadCluster(index);
00539             IDHistNC->Fill(ntpCluster->id);
00540             if(ntpCluster->id==0 && ntpCluster->probem>0.2) {
00541               EMFrac+=ntpCluster->ph.gev;
00542               EnTot+=ntpCluster->ph.gev;
00543               if(ntpCluster->planeview==2){
00544                 EMFrac_end0+=ntpCluster->ph.gev;
00545                 EnTot_end0+=ntpCluster->ph.gev;
00546               }
00547               if(ntpCluster->planeview==3){
00548                 EMFrac_end1+=ntpCluster->ph.gev;
00549                 EnTot_end1+=ntpCluster->ph.gev;
00550               }
00551             }
00552             else {
00553               EnTot+=ntpCluster->ph.gev;
00554               if(ntpCluster->planeview==2){
00555                 EnTot_end0+=ntpCluster->ph.gev;
00556               }
00557               if(ntpCluster->planeview==3){
00558                 EnTot_end1+=ntpCluster->ph.gev;
00559               }
00560             }
00561             
00562             if(ntpCluster->ph.gev>0.5 && 
00563                (ntpCluster->id==0 || ntpCluster->id==1 || ntpCluster->id==3)) {
00564               if(ntpCluster->planeview==2) cntU++; 
00565               if(ntpCluster->planeview==3) cntV++; 
00566             }
00567             if (index==bestU) UClusterIDNC->Fill(ntpCluster->id);
00568             if (index==bestV) VClusterIDNC->Fill(ntpCluster->id);         
00569           }
00570         
00571           if(goodU>=1 || goodV>=1) EMyDistNC->Fill(ntpTruth->y);
00572           if(goodU==-1 || goodV==-1) HADyDistNC->Fill(ntpTruth->y);
00573           if(goodU==-1 && goodV==-1) UVHADyDistNC->Fill(ntpTruth->y);
00574           if(goodU>=1 && goodV>=1) UVEMyDistNC->Fill(ntpTruth->y);
00575           if(goodU==2 || goodV==2) EMProbyDistNC->Fill(ntpTruth->y);
00576           if(goodU==2 && goodV==2) UVEMProbyDistNC->Fill(ntpTruth->y);
00577 
00578           AvgCluEnNCnox->Fill(ntpTruth->p4shw[3],(cntU+cntV)/2.);
00579           //if(ntpTruth->p4neu[3]>3.0 && 
00580           // ntpTruth->p4neu[3]<3.1) 
00581           AvgCluSHNCnox->Fill(nBigStdHep,(cntU+cntV)/2.);
00582           AvgCluSHEnergyNCnox->Fill(ntpTruth->p4neu[3],nBigStdHep,
00583                                     (cntU+cntV)/2.);
00584           
00585           if(EnTot_end0!=0 || EnTot_end1!=0){
00586             if(EnTot_end0==0) EMFrac = EMFrac_end1/EnTot_end1;
00587             else if(EnTot_end1==0) EMFrac = EMFrac_end0/EnTot_end0;
00588             else if(EMFrac_end0/EnTot_end0>EMFrac_end1/EnTot_end1) 
00589               EMFrac = EMFrac_end0/EnTot_end0;
00590             else EMFrac = EMFrac_end1/EnTot_end1;
00591             //else EMFrac = (EMFrac_end0/EnTot_end0 + EMFrac_end1/EnTot_end1)/2.;
00592             EnTot = 1;
00593           }
00594           if(EnTot!=0) {
00595             EMFracNC->Fill(EMFrac/EnTot);
00596             EMFracNCY->Fill(ntpTruth->y,EMFrac/EnTot);
00597             if(ntpTruth->emfrac!=0) EMFracNCDiff->Fill((ntpTruth->emfrac - 
00598                                                         (EMFrac/EnTot))/ntpTruth->emfrac);
00599             else if(EMFrac==0) EMFracNCDiff->Fill(0);
00600             else EMFracNCDiff->Fill(2-EMFrac/EnTot);
00601             EMFrac_Truth_Reco->Fill(ntpTruth->emfrac,EMFrac/EnTot);
00602             if(EnTot_end0>0 && EnTot_end1==0) 
00603               EMFrac_Truth_Reco_end0->Fill(ntpTruth->emfrac,EMFrac_end0/EnTot_end0);
00604             if(EnTot_end1>0 && EnTot_end0==0) 
00605               EMFrac_Truth_Reco_end1->Fill(ntpTruth->emfrac,EMFrac_end1/EnTot_end1);
00606           }
00607         }
00608 
00609         //look at non primary showers
00610         for(int k=1;k<ntpEvent->nshower;k++){
00611           totU = 0;
00612           totV = 0;
00613           Bool_t UHalo = true;
00614           Bool_t VHalo = true;
00615           if(LoadShower(ntpEvent->shw[k])){
00616             for(int l=0; l<ntpShower->ncluster; l++){
00617               if(LoadCluster(ntpShower->clu[l])){
00618                 if(ntpCluster->planeview==2){
00619                   totU+=ntpCluster->ph.gev;
00620                   if(ntpCluster->id!=4) UHalo = false;
00621                 }
00622                 if(ntpCluster->planeview==3){
00623                   totV+=ntpCluster->ph.gev;
00624                   if(ntpCluster->id!=4) VHalo = false;
00625                 }
00626               }
00627             }
00628           }
00629           asymmetry = TMath::Abs(totU-totV);
00630           if((totU+totV)>0.) asymmetry /= (totU+totV);
00631           else asymmetry = 1;
00632           if(ntpTruth->iaction==1&&ntpEvent->ntrack==0){
00633             if(UHalo || VHalo) {
00634               if(totU<totV) asymHistCC_NotPrime->Fill(totU,asymmetry);
00635               else asymHistCC_NotPrime->Fill(totV,asymmetry);
00636             }
00637             //if(asymmetry>0.8){
00638             //if(totU<totV) asymHistCC_NotPrime->Fill(totU+totV,totU);
00639             //else asymHistCC_NotPrime->Fill(totU+totV,totV);
00640             //}
00641           }
00642           else if(ntpTruth->iaction==0){
00643             if(totU<totV) asymHistNC_NotPrime->Fill(totU,asymmetry);
00644             else asymHistNC_NotPrime->Fill(totV,asymmetry);
00645           }
00646         }
00647 
00648       }
00649       else {
00650 
00651         //look at asymmetry for showers not associated with MC event
00652         Double_t totU = 0;
00653         Double_t totV = 0;
00654         if(LoadShower(ntpEvent->shw[0])){
00655           for(int l=0; l<ntpShower->ncluster; l++){
00656             if(LoadCluster(ntpShower->clu[l])){
00657               if(ntpCluster->planeview==2){
00658                 totU+=ntpCluster->ph.gev;
00659               }
00660               if(ntpCluster->planeview==3){
00661                 totV+=ntpCluster->ph.gev;
00662               }
00663             }
00664           }
00665         }
00666         Double_t asymmetry = TMath::Abs(totU-totV);
00667         if((totU+totV)>0.) asymmetry /= (totU+totV);
00668         else asymmetry = 1;
00669         if(totU<totV) asymHist->Fill(totU,asymmetry);
00670         else asymHist->Fill(totV,asymmetry);
00671       }
00672     }
00673   }
00674 
00675   savefile->Write();
00676   savefile->Close();
00677   return;
00678 }

void MadCluAnalysis::Plotter ( char *  fileName  ) 

Definition at line 680 of file MadCluAnalysis.cxx.

00680                                           {
00681 
00682   TFile *file = new TFile(filename,"READ");
00683   file->cd();
00684 
00685   TH1F *IDHist = (TH1F*) file->Get("IDHist");
00686   TH1F *IDHistNC = (TH1F*) file->Get("IDHistNC");
00687   TH1F *yDist = (TH1F*) file->Get("yDist");
00688   TH1F *EMyDist = (TH1F*) file->Get("EMyDist");
00689   TH1F *HADyDist = (TH1F*) file->Get("HADyDist");
00690   TH1F *EMProbyDist = (TH1F*) file->Get("EMProbyDist");
00691   TH1F *UVHADyDist = (TH1F*) file->Get("UVHADyDist");
00692   TH1F *UVEMyDist = (TH1F*) file->Get("UVEMyDist");
00693   TH1F *UVEMProbyDist = (TH1F*) file->Get("UVEMProbyDist");
00694   TH1F *yDistNC = (TH1F*) file->Get("yDistNC");
00695   TH1F *EMyDistNC = (TH1F*) file->Get("EMyDistNC");
00696   TH1F *HADyDistNC = (TH1F*) file->Get("HADyDistNC");
00697   TH1F *EMProbyDistNC = (TH1F*) file->Get("EMProbyDistNC");
00698   TH1F *UVHADyDistNC = (TH1F*) file->Get("UVHADyDistNC");
00699   TH1F *UVEMyDistNC = (TH1F*) file->Get("UVEMyDistNC");
00700   TH1F *UVEMProbyDistNC = (TH1F*) file->Get("UVEMProbyDistNC");
00701   TH1F *UClusterIDCC = (TH1F*) file->Get("UClusterIDCC");
00702   TH1F *UClusterIDNC = (TH1F*) file->Get("UClusterIDNC");
00703   TH1F *VClusterIDCC = (TH1F*) file->Get("VClusterIDCC");
00704   TH1F *VClusterIDNC = (TH1F*) file->Get("VClusterIDNC");
00705   TH2F *AvgCluEnCC = (TH2F*) file->Get("AvgCluEnCC");
00706   TH2F *AvgCluEnNC = (TH2F*) file->Get("AvgCluEnNC");
00707   TH2F *AvgCluEnCCnox = (TH2F*) file->Get("AvgCluEnCCnox");
00708   TH2F *AvgCluEnNCnox = (TH2F*) file->Get("AvgCluEnNCnox");
00709   TH2F *AvgCluSHCC = (TH2F*) file->Get("AvgCluSHCC");
00710   TH2F *AvgCluSHNC = (TH2F*) file->Get("AvgCluSHNC");
00711   TH2F *AvgCluSHCCnox = (TH2F*) file->Get("AvgCluSHCCnox");
00712   TH2F *AvgCluSHNCnox = (TH2F*) file->Get("AvgCluSHNCnox");
00713   TH1F *EMFracCC = (TH1F*) file->Get("EMFracCC");
00714   TH1F *EMFracNC = (TH1F*) file->Get("EMFracNC");
00715   TH1F *EMFracCCDiff = (TH1F*) file->Get("EMFracCCDiff");
00716   TH1F *EMFracNCDiff = (TH1F*) file->Get("EMFracNCDiff");
00717   TH2F *EMFracCCY = (TH2F*) file->Get("EMFracCCY");
00718   TH2F *EMFracNCY = (TH2F*) file->Get("EMFracNCY");
00719   TH2F *EMFracTrueReco = (TH2F*) file->Get("EMFracTrueReco");
00720   TH2F *EMFracTrueReco_end0 = (TH2F*) file->Get("EMFracTrueReco_end0");
00721   TH2F *EMFracTrueReco_end1 = (TH2F*) file->Get("EMFracTrueReco_end1");
00722 
00723   TCanvas *can = new TCanvas();
00724   can->Divide(1,2);
00725   can->cd(1);
00726   EMyDist->Divide(yDist);
00727   HADyDist->Divide(yDist);
00728   EMProbyDist->Divide(yDist);
00729   UVEMyDist->Divide(yDist);
00730   UVHADyDist->Divide(yDist);
00731   UVEMProbyDist->Divide(yDist);
00732   EMyDist->SetLineColor(4);
00733   EMProbyDist->SetLineColor(6);
00734   HADyDist->SetLineColor(2);
00735   UVEMyDist->SetLineColor(3);  
00736   UVHADyDist->SetLineColor(1);
00737   UVEMProbyDist->SetLineColor(7);
00738   EMyDist->Draw();
00739   HADyDist->Draw("sames");
00740   EMProbyDist->Draw("sames");
00741   UVEMyDist->Draw("sames");
00742   UVHADyDist->Draw("sames");
00743   UVEMProbyDist->Draw("sames");
00744   can->cd(2);
00745   EMyDistNC->Divide(yDistNC);
00746   HADyDistNC->Divide(yDistNC);
00747   EMProbyDistNC->Divide(yDistNC);
00748   UVEMyDistNC->Divide(yDistNC);
00749   UVHADyDistNC->Divide(yDistNC);
00750   UVEMProbyDistNC->Divide(yDistNC);
00751   EMyDistNC->SetLineColor(4);
00752   HADyDistNC->SetLineColor(2);
00753   EMProbyDistNC->SetLineColor(6);
00754   UVEMyDistNC->SetLineColor(3);
00755   UVHADyDistNC->SetLineColor(1);
00756   UVEMProbyDistNC->SetLineColor(7);
00757   EMyDistNC->Draw();
00758   HADyDistNC->Draw("sames");
00759   EMProbyDistNC->Draw("sames");
00760   UVEMyDistNC->Draw("sames");
00761   UVHADyDistNC->Draw("sames");
00762   UVEMProbyDistNC->Draw("sames");
00763 
00764   TCanvas *can1 = new TCanvas();
00765   can1->Divide(1,2);
00766   can1->cd(1);
00767   UClusterIDCC->SetLineColor(4);
00768   float integ = UClusterIDCC->Integral();
00769   UClusterIDCC->Scale(1./integ);
00770   UClusterIDCC->Draw();
00771   UClusterIDNC->SetLineColor(2);
00772   integ = UClusterIDNC->Integral();
00773   UClusterIDNC->Scale(1./integ);
00774   UClusterIDNC->Draw("sames");
00775   can1->cd(2);
00776   VClusterIDCC->SetLineColor(4);
00777   integ = VClusterIDCC->Integral();
00778   VClusterIDCC->Scale(1./integ);
00779   VClusterIDCC->Draw();
00780   VClusterIDNC->SetLineColor(2);
00781   integ = VClusterIDNC->Integral();
00782   VClusterIDNC->Scale(1./integ);
00783   VClusterIDNC->Draw("sames");
00784 
00785   TCanvas *can2 = new TCanvas();
00786   can2->Divide(2,2);
00787   can2->cd(1);
00788   AvgCluEnCC->Draw("COLZ");
00789   TProfile *AvgCluEnCC_pfx = AvgCluEnCC->ProfileX();
00790   AvgCluEnCC_pfx->Draw("same");
00791   can2->cd(2);
00792   AvgCluEnNC->Draw("COLZ");
00793   TProfile *AvgCluEnNC_pfx = AvgCluEnNC->ProfileX();
00794   AvgCluEnNC_pfx->Draw("same");
00795   can2->cd(3);
00796   AvgCluEnCCnox->Draw("COLZ");
00797   TProfile *AvgCluEnCCnox_pfx = AvgCluEnCCnox->ProfileX();
00798   AvgCluEnCCnox_pfx->Draw("same");
00799   can2->cd(4);
00800   AvgCluEnNCnox->Draw("COLZ");
00801   TProfile *AvgCluEnNCnox_pfx = AvgCluEnNCnox->ProfileX();
00802   AvgCluEnNCnox_pfx->Draw("same");
00803 
00804   TCanvas *can3 = new TCanvas();
00805   can3->Divide(2,2);
00806   can3->cd(1);
00807   AvgCluSHCC->Draw("COLZ");
00808   TProfile *AvgCluSHCC_pfx = AvgCluSHCC->ProfileX();
00809   AvgCluSHCC_pfx->Draw("same");
00810   can3->cd(2);
00811   AvgCluSHNC->Draw("COLZ");
00812   TProfile *AvgCluSHNC_pfx = AvgCluSHNC->ProfileX();
00813   AvgCluSHNC_pfx->Draw("same");
00814   can3->cd(3);
00815   AvgCluSHCCnox->Draw("COLZ");
00816   TProfile *AvgCluSHCCnox_pfx = AvgCluSHCCnox->ProfileX();
00817   AvgCluSHCCnox_pfx->Draw("same");
00818   can3->cd(4);
00819   AvgCluSHNCnox->Draw("COLZ");
00820   TProfile *AvgCluSHNCnox_pfx = AvgCluSHNCnox->ProfileX();
00821   AvgCluSHNCnox_pfx->Draw("same");
00822 
00823   TCanvas *can4 = new TCanvas();
00824   can4->cd();
00825   IDHist->SetLineColor(4);
00826   IDHistNC->SetLineColor(2);
00827   IDHist->Draw();
00828   IDHistNC->Draw("sames");
00829 
00830   TCanvas *can5 = new TCanvas();
00831   can5->Divide(2,2);
00832   can5->cd(1);
00833   EMFracCC->Draw();
00834   can5->cd(2);
00835   EMFracNC->Draw();
00836   can5->cd(3);
00837   EMFracCCDiff->Draw();
00838   can5->cd(4);
00839   EMFracNCDiff->Draw();
00840 
00841   TCanvas *can6 = new TCanvas();
00842   can6->Divide(1,2);
00843   can6->cd(1);
00844   EMFracCCY->Draw("colz");
00845   can6->cd(2);
00846   EMFracNCY->Draw("colz");
00847 
00848   TCanvas *can7 = new TCanvas();
00849   can7->Divide(1,3);
00850   can7->cd(1);
00851   EMFracTrueReco->Draw("COLZ");
00852   TProfile *prof = EMFracTrueReco->ProfileX();
00853   prof->Draw("sames");
00854   can7->cd(2);
00855   EMFracTrueReco_end0->Draw("COLZ");
00856   prof = EMFracTrueReco_end0->ProfileX();
00857   prof->Draw("sames");
00858   can7->cd(3);
00859   EMFracTrueReco_end1->Draw("COLZ");
00860   prof = EMFracTrueReco_end1->ProfileX();
00861   prof->Draw("sames");
00862   
00863   return;
00864 }

Bool_t MadCluAnalysis::ReadNNFile ( char *  filename  ) 

Definition at line 2897 of file MadCluAnalysis.cxx.

References fneuralNet.

02897                                                {
02898 
02899   TFile *f = new TFile(filename,"READ");
02900   if(f->IsOpen()&&!f->IsZombie()){
02901     fneuralNet = (TMultiLayerPerceptron*) f->Get("TMultiLayerPerceptron");
02902     if(fneuralNet) cout << "Read NN file" << endl;
02903   }
02904   else return false;
02905   return true;
02906 }

Bool_t MadCluAnalysis::ReadPIDFile ( char *  filename  ) 

Definition at line 2875 of file MadCluAnalysis.cxx.

References flikelihoodDists.

02875                                                 {
02876 
02877   TFile *f = new TFile(filename,"READ");
02878   if(f->IsOpen()&&!f->IsZombie()){
02879     flikelihoodDists[0] = (TH1F*) f->Get("PhwIDU8");
02880     flikelihoodDists[1] = (TH1F*) f->Get("PhwIDU10");
02881     flikelihoodDists[2] = (TH1F*) f->Get("PhwProbEMU8");
02882     flikelihoodDists[3] = (TH1F*) f->Get("PhwProbEMU10");
02883     flikelihoodDists[4] = (TH1F*) f->Get("ChargeFracRMSU8");
02884     flikelihoodDists[5] = (TH1F*) f->Get("ChargeFracRMSU10");
02885     flikelihoodDists[6] = (TH1F*) f->Get("PhwIDV8");
02886     flikelihoodDists[7] = (TH1F*) f->Get("PhwIDV10");
02887     flikelihoodDists[8] = (TH1F*) f->Get("PhwProbEMV8");
02888     flikelihoodDists[9] = (TH1F*) f->Get("PhwProbEMV10");
02889     flikelihoodDists[10] = (TH1F*) f->Get("ChargeFracRMSV8");
02890     flikelihoodDists[11] = (TH1F*) f->Get("ChargeFracRMSV10");
02891     cout << "Read PID file" << endl;
02892   }
02893   else return false;
02894   return true;
02895 }

Float_t MadCluAnalysis::RecoAnalysisEnergy ( Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis.

Definition at line 2456 of file MadCluAnalysis.cxx.

References NtpSRStripPulseHeight::gev, MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::ntpShower, NtpSRShower::ph, and MadQuantities::RecoMuEnergy().

02457 {
02458   if(!LoadEvent(event)) return 0;
02459   Int_t shower = -1;
02460   LoadLargestShowerFromEvent(event,shower);
02461   Int_t track = -1;
02462   LoadLargestTrackFromEvent(event,track);
02463   Double_t shw_en = 0;
02464   Double_t trk_en = 0;
02465   if(shower>=0) shw_en = ntpShower->ph.gev;
02466   if(track>=0) trk_en = RecoMuEnergy(0,track);
02467   return shw_en + trk_en;
02468 }


Member Data Documentation

Definition at line 91 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 98 of file MadCluAnalysis.h.

Referenced by PID(), and ReadPIDFile().

TMultiLayerPerceptron* MadCluAnalysis::fneuralNet

Definition at line 99 of file MadCluAnalysis.h.

Referenced by MadCluAnalysis(), PID(), and ReadNNFile().

Definition at line 87 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 88 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 89 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 74 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 70 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 93 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), FillNCluster(), and MadCluAnalysis().

Definition at line 94 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), FillNCluster(), and MadCluAnalysis().

Definition at line 95 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), FillNCluster(), and MadCluAnalysis().

Definition at line 96 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), FillNCluster(), and MadCluAnalysis().

Definition at line 80 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 81 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 78 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 76 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 77 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 66 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 58 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 62 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 84 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().

Definition at line 85 of file MadCluAnalysis.h.

Referenced by AddUserBranches(), CallUserFunctions(), and MadCluAnalysis().


The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1