00001 #ifndef madcluanalysis_cxx
00002 #define madcluanalysis_cxx
00003 #include "MadCluAnalysis.h"
00004 #include "TH2.h"
00005 #include "TH3.h"
00006 #include "TProfile.h"
00007 #include "TStyle.h"
00008 #include "TCanvas.h"
00009 #include "TLatex.h"
00010 #include "TMath.h"
00011 #include "TTree.h"
00012 #include "TMatrixD.h"
00013 #include "TVectorD.h"
00014 #include "TROOT.h"
00015 #include "TGraphErrors.h"
00016 #include "TF1.h"
00017 #include "TF2.h"
00018 #include "TVector2.h"
00019 #include "TString.h"
00020 #include "TClonesArray.h"
00021 #include <sstream>
00022 #include <iostream>
00023 #include <cmath>
00024 #include "Mad/macros/OscReweight.C"
00025
00026 #include "Mad/SpillInfo.h"
00027 #include "Mad/SpillInfoBlock.h"
00028 #include "BeamDataUtil/BeamMonSpill.h"
00029 #include "BeamDataUtil/BDSpillAccessor.h"
00030 #include "SpillTiming/SpillTimeFinder.h"
00031
00032 using namespace std;
00033
00034
00035
00036 MadCluAnalysis::MadCluAnalysis(TChain *chainSR,TChain *chainMC,TChain *chainTH,TChain *chainEM)
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 }
00099
00100 MadCluAnalysis::MadCluAnalysis(JobC *j,string path,int entries)
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 }
00160
00161
00162 MadCluAnalysis::~MadCluAnalysis()
00163 {
00164
00165 }
00166
00167
00168
00169 void MadCluAnalysis::Plot(char *fileName,Int_t startEnt)
00170 {
00171
00172 char savename[256];
00173 sprintf(savename,"Clu_%s.root",fileName);
00174 TFile *savefile = new TFile(savename,"RECREATE");
00175
00176
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
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
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
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
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
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
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
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
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
00325 for(int j=0;j<eventSummary->nevent;j++){
00326 if(!LoadEvent(j)) continue;
00327 if(ntpEvent->nshower==0) continue;
00328
00329
00330 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00331
00332 is_fid = 1;
00333 if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
00334 }
00335 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00336
00337 is_fid = 1;
00338 if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
00339 }
00340
00341
00342 Int_t mcevent = -1;
00343
00344 if(LoadTruthForRecoTH(ntpEvent->index,mcevent)) {
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){
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
00366
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
00414
00415 if(ntpTruth->iaction==1){
00416
00417 yDist->Fill(ntpTruth->y);
00418
00419
00420
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
00428
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
00485
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
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
00523
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
00580
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
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
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
00638
00639
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
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 }
00679
00680 void MadCluAnalysis::Plotter(char *filename){
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 }
00865
00866 void MadCluAnalysis::DataDistributions(char *filename){
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
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
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
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
01095 if(eventSummary->nevent==0) continue;
01096
01097
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
01119 for(int j=0;j<eventSummary->nevent;j++){
01120 if(!LoadEvent(j)) continue;
01121 if(ntpEvent->nshower==0) continue;
01122
01123
01124 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01125
01126 is_fid = 1;
01127 if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
01128 }
01129 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01130
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
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
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 }
01279
01280 void MadCluAnalysis::DataMCComp(char *dataname,char *mcname){
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
01451
01452
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
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
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
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
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 }
01818
01819
01820 Bool_t MadCluAnalysis::PassDataCleaningCuts(Int_t shower)
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 }
01843
01844
01845 Bool_t MadCluAnalysis::IsTrackInSlice(Int_t slice){
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 }
01858
01859
01860 Double_t MadCluAnalysis::FOM()
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
01998 for(int j=0;j<eventSummary->nevent;j++){
01999 if(!LoadEvent(j)) continue;
02000 if(ntpEvent->nshower==0) continue;
02001
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
02009
02010
02011
02012
02013
02014
02015 }
02016
02017
02018
02019
02020 if(ntpTruth->emfrac>0.6){
02021
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 }
02232
02233
02234 Bool_t MadCluAnalysis::PassBasicCuts(Int_t event){
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 }
02282
02283 Bool_t MadCluAnalysis::PassAnalysisCuts(Int_t event){
02284 if(!PassBasicCuts(event)) return false;
02285 if(PID(event,1)>0. || PID(event,2)>0.) return true;
02286 return false;
02287 }
02288
02289 Bool_t MadCluAnalysis::PassTruthSignal(Int_t mcevent)
02290 {
02291 if(!LoadTruth(mcevent)) return false;
02292 if(ntpTruth->inu==12 && ntpTruth->iaction==1) return true;
02293 return false;
02294 }
02295
02296 Float_t MadCluAnalysis::PID(Int_t event,Int_t method)
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
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
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 }
02454
02455
02456 Float_t MadCluAnalysis::RecoAnalysisEnergy(Int_t event)
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 }
02469
02470
02471 Float_t MadCluAnalysis::GetWeight(Int_t )
02472 {
02473 return 1;
02474 }
02475
02476
02477 Double_t *MadCluAnalysis::PHWProbEM(Int_t event)
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 }
02512
02513 Bool_t MadCluAnalysis::FillPHWProbEM(Int_t event)
02514 {
02515 Double_t *phwProbEM = PHWProbEM(event);
02516 PAN_PHWProbEMU = phwProbEM[0];
02517 PAN_PHWProbEMV = phwProbEM[1];
02518 delete [] phwProbEM;
02519 return true;
02520 }
02521
02522 Double_t *MadCluAnalysis::PHWCluID(Int_t event)
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 }
02556
02557 Bool_t MadCluAnalysis::FillPHWCluID(Int_t event)
02558 {
02559 Double_t *phwCluID = PHWCluID(event);
02560 PAN_PHWCluIDU = phwCluID[0];
02561 PAN_PHWCluIDV = phwCluID[1];
02562 delete [] phwCluID;
02563 return true;
02564 }
02565
02566 Double_t *MadCluAnalysis::PHWAvgDev(Int_t event)
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 }
02600
02601 Bool_t MadCluAnalysis::FillPHWAvgDev(Int_t event)
02602 {
02603 Double_t *phwAvgDev = PHWAvgDev(event);
02604 PAN_PHWAvgDevU = phwAvgDev[0];
02605 PAN_PHWAvgDevV = phwAvgDev[1];
02606 delete [] phwAvgDev;
02607 return true;
02608 }
02609
02610 Double_t *MadCluAnalysis::ChargeFracRMS(Int_t event,Int_t method)
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 }
02654
02655 Bool_t MadCluAnalysis::FillChargeFracRMS(Int_t event,Int_t method)
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 }
02663
02664 Double_t *MadCluAnalysis::AveNumSSPlane(Int_t event)
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 }
02702
02703 Bool_t MadCluAnalysis::FillAveNumSSPlane(Int_t event)
02704 {
02705 Double_t *aveNumSSPlane = AveNumSSPlane(event);
02706 PAN_AveNumSSPlaneU = aveNumSSPlane[0];
02707 PAN_AveNumSSPlaneV = aveNumSSPlane[1];
02708 delete [] aveNumSSPlane;
02709 return true;
02710 }
02711
02712 Bool_t MadCluAnalysis::AddUserBranches(TTree *tree)
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 }
02756
02757 void MadCluAnalysis::CallUserFunctions(Int_t event)
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
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
02832
02833
02834
02835
02836
02837
02838
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 }
02874
02875 Bool_t MadCluAnalysis::ReadPIDFile(char *filename){
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 }
02896
02897 Bool_t MadCluAnalysis::ReadNNFile(char *filename){
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 }
02907
02908 Bool_t MadCluAnalysis::IsNueFidVtxEvt(Int_t ievt){
02909 if(!LoadEvent(ievt)) return false;
02910 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
02911 if(ntpEvent->vtx.z<1 || ntpEvent->vtx.z>29 ||
02912 (ntpEvent->vtx.z<17&&ntpEvent->vtx.z>14) ||
02913 sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)
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 }
02923
02924
02925 Bool_t MadCluAnalysis::FillNCluster(Int_t event){
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 }
02949
02950 Bool_t MadCluAnalysis::FillBestProbEM(Int_t event){
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 }
02977
02978 #endif //madcluanalysis_cxx
