00001 #ifndef maddpanalysis_cxx
00002 #define maddpanalysis_cxx
00003 #include <iostream>
00004 #include <cmath>
00005 #include "Mad/MadDpAnalysis.h"
00006 #include "TClonesArray.h"
00007 #include "TMath.h"
00008 #include "Validity/VldContext.h"
00009 #include "Conventions/Detector.h"
00010 #include "Registry/Registry.h"
00011 #include "Mad/MadAnalysis.h"
00012 #include "MCReweight/MCReweight.h"
00013 #include "MCReweight/GnumiInterface.h"
00014 #include "MCReweight/NuParent.h"
00015 #include "Mad/SpillInfo.h"
00016 #include "Mad/SpillInfoBlock.h"
00017 #include "BeamDataUtil/BDSpillAccessor.h"
00018 #include "BeamDataUtil/BeamMonSpill.h"
00019 #include "SpillTiming/SpillTimeFinder.h"
00020 #include "Riostream.h"
00021 #include "Mad/ScanList.h"
00022 #include "Mad/fddataquality.h"
00023
00024 using namespace std;
00025
00026
00027 MadDpAnalysis::MadDpAnalysis(TChain *chainSR,TChain *chainMC,
00028 TChain *chainTH,TChain *chainEM)
00029 {
00030
00031 if(!chainSR) {
00032 record = 0;
00033 strecord = 0;
00034 emrecord = 0;
00035 mcrecord = 0;
00036 threcord = 0;
00037 Clear();
00038 whichSource = -1;
00039 isMC = true;
00040 isTH = true;
00041 isEM = true;
00042 Nentries = -1;
00043 return;
00044 }
00045
00046 InitChain(chainSR,chainMC,chainTH,chainEM);
00047 whichSource = 0;
00048 fLikeliFile = NULL;
00049 for(int i=0;i<6;i++) fLikeliHist[i] = NULL;
00050
00051 }
00052
00053
00054 MadDpAnalysis::MadDpAnalysis(JobC *j,string path,int entries)
00055 {
00056 record = 0;
00057 strecord = 0;
00058 emrecord = 0;
00059 mcrecord = 0;
00060 threcord = 0;
00061 Clear();
00062 isMC = true;
00063 isTH = true;
00064 isEM = true;
00065 Nentries = entries;
00066 whichSource = 1;
00067 jcPath = path;
00068 fJC = j;
00069 fLikeliFile = NULL;
00070 for(int i=0;i<6;i++) fLikeliHist[i] = NULL;
00071 }
00072
00073
00074 MadDpAnalysis::~MadDpAnalysis()
00075 {
00076 if(fLikeliFile) fLikeliFile->Close();
00077 }
00078
00079
00080 Bool_t MadDpAnalysis::PassTruthSignal(Int_t mcevent)
00081 {
00082 if(!LoadTruth(mcevent)) return false;
00083 if(abs(ntpTruth->inu)==14&&ntpTruth->iaction==1) return true;
00084 return false;
00085 }
00086
00087
00088 Bool_t MadDpAnalysis::PassAnalysisCuts(Int_t event)
00089 {
00090 if(!LoadEvent(event)) return false;
00091
00092 Float_t pidcut=-0.4;
00093 if(ntpHeader->GetVldContext().GetDetector()==1) {pidcut=-0.1;}
00094
00095 if(!(PID(event,0)>pidcut)) return false;
00096 return true;
00097 }
00098 Bool_t MadDpAnalysis::PassFidNew(Int_t event)
00099 {
00100
00101 if(!LoadEvent(event)) return false;
00102
00103 Int_t track = -1;
00104 LoadLargestTrackFromEvent(event,track);
00105
00106
00107
00108 if (track==-1 && (ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>28.0 ||
00109 (ntpEvent->vtx.z<16.2 && ntpEvent->vtx.z>14.3) ||
00110 (pow(ntpEvent->vtx.x,2)+
00111 pow(ntpEvent->vtx.y,2))>14)) return false;
00112
00113
00114
00115 if (track!=-1 && (ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>28.0 ||
00116 (ntpTrack->vtx.z<16.2 && ntpTrack->vtx.z>14.3) ||
00117 (pow(ntpTrack->vtx.x,2)+
00118 pow(ntpTrack->vtx.y,2))>14)) return false;
00119
00120 return true;
00121
00122 }
00123
00124
00125
00126 Bool_t MadDpAnalysis::PassFidNewN(Int_t event)
00127 {
00128
00129 if(!LoadEvent(event)) return false;
00130
00131 Int_t track = -1;
00132 LoadLargestTrackFromEvent(event,track);
00133
00134
00135
00136 if (track==-1 && (ntpEvent->vtx.z<1.0 || ntpEvent->vtx.z>27.0 ||
00137 (ntpEvent->vtx.z<17.0 && ntpEvent->vtx.z>14.0) ||
00138 (pow(ntpEvent->vtx.x,2)+
00139 pow(ntpEvent->vtx.y,2))>(3.5*3.5))) return false;
00140
00141
00142
00143 if (track!=-1 && (ntpTrack->vtx.z<1.0 || ntpTrack->vtx.z>27.0 ||
00144 (ntpTrack->vtx.z<17. && ntpTrack->vtx.z>14.0) ||
00145 (pow(ntpTrack->vtx.x,2)+
00146 pow(ntpTrack->vtx.y,2))>(3.5*3.5))) return false;
00147
00148 return true;
00149
00150 }
00151
00152 Bool_t MadDpAnalysis::PassBasicCuts(Int_t event)
00153 {
00154 if(!LoadEvent(event)) return false;
00155 if(ntpEvent->ntrack==0) return false;
00156
00157
00158
00159
00160
00161 return true;
00162 }
00163 Bool_t MadDpAnalysis::PassFid(Int_t event)
00164 {
00165 if(!LoadEvent(event)) return false;
00166 if(!IsFidVtxEvt(ntpEvent->index)) return false;
00167 return true;
00168 }
00169
00170
00171
00172 AnnInputBlock MadDpAnalysis::AnnVar(Int_t event)
00173 {
00174 AnnInputBlock anninput;
00175
00176 anninput.aTotrk =0.;
00177 anninput.aTotstp =0.;
00178 anninput.aTotph =0.;
00179 anninput.aTotlen =0.;
00180 anninput.aPhperpl =0.;
00181 anninput.aPhperstp =0.;
00182
00183
00184 anninput.aTrkpass =0.;
00185 anninput.aTrkph =0.;
00186 anninput.aTrklen =0.;
00187 anninput.aTrkphperpl =0.;
00188 anninput.aTrkphper =0.;
00189 anninput.aTrkplu =0.;
00190 anninput.aTrkplv =0.;
00191 anninput.aTrkstp =0.;
00192
00193 anninput.aShwph =0.;
00194 anninput.aShwstp =0.;
00195 anninput.aShwdig =0.;
00196 anninput.aShwpl =0.;
00197 anninput.aShwphper =0.;
00198 anninput.aShwphperpl =0.;
00199 anninput.aShwphperdig =0.;
00200 anninput.aShwphperstp =0.;
00201 anninput.aShwplu =0.;
00202 anninput.aShwplv =0.;
00203 anninput.aPh3 =0.;
00204 anninput.aPh6 =0.;
00205 anninput.aPhlast =0.;
00206 anninput.aPhcommon =0.;
00207 anninput.aTimemax =0.;
00208 anninput.aTimemin =0.;
00209
00210
00211
00212 Int_t detector = 0;
00213 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) detector=1;
00214 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar) detector=2;
00215
00216 LoadEvent(event) ;
00217 int track_index = -1;
00218 LoadLargestTrackFromEvent(event,track_index);
00219 int shower_index = -1;
00220
00221 LoadShower_Jim(event,shower_index,detector);
00222
00223 anninput.aTotrk =ntpEvent->ntrack;
00224 anninput.aTotstp =ntpEvent->nstrip;
00225 anninput.aTotph =ntpEvent->ph.sigcor;
00226 anninput.aTotlen =ntpEvent->plane.end-ntpEvent->plane.beg+1;
00227
00228 if (ntpEvent->plane.n>0) anninput.aPhperpl =ntpEvent->ph.sigcor/ntpEvent->plane.n;
00229 if (ntpEvent->nstrip>0) anninput.aPhperstp =ntpEvent->ph.sigcor/ntpEvent->nstrip;
00230
00231 if(track_index!=-1) {
00232 anninput.aTrkpass =ntpTrack->fit.pass;
00233 anninput.aTrkph =ntpTrack->ph.sigcor;
00234 anninput.aTrklen =ntpTrack->plane.ntrklike;
00235 if(ntpTrack->plane.ntrklike>0) anninput.aTrkphperpl =ntpTrack->ph.sigcor/ntpTrack->plane.ntrklike;
00236 if(ntpEvent->ph.sigcor>0) anninput.aTrkphper =ntpTrack->ph.sigcor/ntpEvent->ph.sigcor;
00237 anninput.aTrkplu =ntpTrack->plane.nu;
00238 anninput.aTrkplv =ntpTrack->plane.nv;
00239 anninput.aTrkstp =ntpTrack->nstrip;
00240 }
00241 if(shower_index!=-1) {
00242 anninput.aShwph =ntpShower->ph.sigcor;
00243 anninput.aShwstp =ntpShower->nstrip;
00244 anninput.aShwdig =ntpShower->ndigit;
00245 anninput.aShwpl =ntpShower->plane.end-ntpShower->plane.beg+1;
00246 if (ntpEvent->ph.sigcor>0) anninput.aShwphper =ntpShower->ph.sigcor/ntpEvent->ph.sigcor;
00247 anninput.aShwphperpl =ntpShower->ph.sigcor/(ntpShower->plane.end-ntpShower->plane.beg+1);
00248
00249 if (ntpShower->ndigit>0) anninput.aShwphperdig =ntpShower->ph.sigcor/ntpShower->ndigit;
00250 if (ntpShower->nstrip>0) anninput.aShwphperstp =ntpShower->ph.sigcor/ntpShower->nstrip;
00251 anninput.aShwplu =ntpShower->plane.nu;
00252 anninput.aShwplv =ntpShower->plane.nv;
00253 }
00254
00255 Int_t ssind,ssplind;
00256 Double_t ssphtot;
00257 Bool_t foundshw,foundtrk;
00258 Int_t planes=ntpEvent->plane.beg;
00259
00260 Double_t ph3,ph6,phlast,phcommon,phnowere,hitcommon,hitnowere;
00261 ph3=0.;ph6=0.;phlast=0.;phcommon=0.;phnowere=0.;hitcommon=0.;hitnowere=0.;
00262
00263 Double_t timemin=9.*10e30;
00264 Double_t timemax=-9999999;
00265
00266 for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
00267 Int_t stp_index=((ntpEvent->stp)[evss]);
00268 if(stp_index!=-1){
00269 LoadStrip(stp_index);
00270
00271
00272 if(detector==1){
00273 Double_t striptime=ntpStrip->time1;
00274 if(striptime<=timemin) timemin=striptime;
00275 if(striptime>=timemax) timemax=striptime;
00276 }
00277
00278 if(detector==2){
00279 Double_t striptime1=ntpStrip->time1;
00280 Double_t striptime0=ntpStrip->time0;
00281 Double_t striptime=0;
00282 if(striptime1>0 && striptime0<0) striptime=striptime1;
00283 if(striptime0>0 && striptime1<0) striptime=striptime0;
00284 if(striptime0>0 && striptime1>0) striptime=(striptime0+striptime1)/2.;
00285 if(striptime<=timemin) timemin=striptime;
00286 if(striptime>=timemax) timemax=striptime;
00287
00288 }
00289
00290 ssind=ntpStrip->strip;
00291 ssplind=ntpStrip->plane;
00292 ssphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
00293
00294 foundshw=false;
00295 foundtrk=false;
00296
00297 Double_t phstrips=0;
00298 Double_t phstript=0;
00299
00300 Int_t sshwind,sshwplind;
00301 Double_t sshwphtot;
00302
00303 if(shower_index!=-1) {
00304 for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
00305 Int_t stp_indexs=((ntpShower->stp)[jj]);
00306 if(stp_indexs!=-1){
00307 LoadStrip(stp_indexs);
00308
00309 if(!foundshw){
00310 sshwind=ntpStrip->strip;
00311 sshwplind=ntpStrip->plane;
00312 sshwphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
00313 if(sshwind==ssind && sshwplind==ssplind) {
00314 foundshw=true;
00315 phstrips=sshwphtot;
00316 }
00317 }
00318 }
00319 }
00320
00321 }
00322
00323 Int_t strkind,strkplind;
00324 Double_t strkphtot;
00325 if(track_index!=-1) {
00326 for(Int_t jj=0;jj<ntpTrack->nstrip;jj++){
00327 Int_t stp_indext=((ntpTrack->stp)[jj]);
00328
00329 if(stp_indext!=-1){
00330 LoadStrip(stp_indext);
00331 if(!foundtrk){
00332 strkind=ntpStrip->strip;
00333 strkplind=ntpStrip->plane;
00334 strkphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
00335 if(strkind==ssind && strkplind==ssplind) {
00336 foundtrk=true;
00337 phstript=strkphtot;
00338 }
00339 }
00340 }
00341
00342 }
00343 }
00344
00345 if(foundshw && foundtrk) {
00346 hitcommon=hitcommon+1;
00347 phcommon=phcommon+phstrips+phstript;
00348 }
00349 if(!foundshw && ! foundtrk) {
00350 hitnowere=hitnowere+1;
00351 phnowere=phnowere+ssphtot;
00352 }
00353 if(ssplind>=planes && ssplind<=(planes+3)){
00354 ph3=ph3+ssphtot;
00355 }
00356 else if(ssplind>(planes+3) && ssplind<=(planes+6)){
00357 ph6=ph6+ssphtot;
00358 }
00359 else {
00360 phlast=phlast+ssphtot;
00361 }
00362 }
00363 }
00364
00365 anninput.aTimemin = timemin;
00366 anninput.aTimemax = timemax;
00367
00368
00369 anninput.aPh3 =ph3;
00370 anninput.aPh6 =ph6;
00371 anninput.aPhlast =phlast;
00372 anninput.aPhcommon =phcommon;
00373 return anninput;
00374 }
00375
00376 Double_t MadDpAnalysis::Sigmoid(Double_t x)
00377 {
00378 double sig;
00379 if(x>37.){
00380 sig = 1.;
00381 }
00382 else if(x<-37.){
00383 sig = 0.;
00384 }
00385 else {
00386 sig = 1./(1.+exp(-x));
00387 }
00388 return sig;
00389 }
00390
00391
00392 Double_t MadDpAnalysis::GetAnnPid(AnnInputBlock anninput,Int_t det, Int_t tar,Int_t fid, Int_t prior, Bool_t first_ann)
00393 {
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404 int inneuron = 7;
00405 int hidneuron = 15;
00406 double var;
00407
00408 double out[15];
00409 double rin[15];
00410
00411
00412 double weight1[7][15];
00413 double constant1[15];
00414
00415 double weighto[15];
00416 double constanto[1];
00417 double prob;
00418
00419
00420
00421 for(Int_t i=0;i<hidneuron;i++) {
00422 out[i]=0.;
00423 rin[i]=0.;
00424 constant1[i]=0.;
00425 }
00426
00427 for(Int_t i=0;i<inneuron;i++) {
00428 for(Int_t j=0;j<hidneuron;j++){
00429 weight1[i][j]=0.;
00430 }
00431 }
00432
00433 constanto[0]=0.;
00434 for(Int_t i=0;i<hidneuron;i++){
00435 weighto[i]=0.;
00436 }
00437 Int_t all = inneuron*hidneuron+hidneuron+hidneuron+1;
00438 Int_t all1 = inneuron*hidneuron+hidneuron;
00439 Int_t n1 =0;
00440 Int_t nw1 =0;
00441 Int_t no =0;
00442 Int_t nwo =0;
00443 ifstream weightfile;
00444
00445
00446 if(first_ann){
00447
00448 if(det==2){
00449 if(prior==1){
00450 if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00451 if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00452 }
00453 if(prior==2){
00454 if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00455 if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00456 }
00457 if(prior==3){
00458 if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00459 if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00460 }
00461 }
00462
00463 if(det==1){
00464 if(tar==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearle_cedar_daikon7v2.dat");
00465 if(tar==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearpme_cedar_daikon7v2.dat");
00466 if(tar==3) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearphe_cedar_daikon7v2.dat");
00467 }
00468 for(Int_t k=0;k<all;k++){
00469 weightfile >> var ;
00470 Ann_weight_vector[k]=var;
00471 }
00472 weightfile.close();
00473 }
00474 if(det==2){
00475
00476 out[0] = (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/2000.;
00477 out[1] = (anninput.aPhperpl/10000.);
00478
00479 out[2] = (anninput.aTotlen/40);
00480 out[3] = (anninput.aPh3/(anninput.aTotph+1.));
00481 out[4] = (anninput.aPh6/(anninput.aTotph+1.));
00482 out[5] = (anninput.aPhlast/(anninput.aTotph+1.));
00483 out[6] = ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.;
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495 }
00496
00497
00498 if(det==1){
00499
00500 out[0] = (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/5000.;
00501 out[1] = (anninput.aPhperpl/10000.);
00502
00503 out[2] = (anninput.aTotlen/40);
00504 out[3] = (anninput.aPh3/(anninput.aTotph+1.));
00505 out[4] = (anninput.aPh6/(anninput.aTotph+1.));
00506 out[5] = (anninput.aPhlast/(anninput.aTotph+1.));
00507 out[6] = ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.;
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519 }
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529 for(Int_t k=0;k<all;k++){
00530 if(k<all1){
00531 if(k%(inneuron+1)==0){
00532 nw1=0;
00533 constant1[n1]=Ann_weight_vector[k];
00534 n1=n1+1;
00535 }
00536 else {
00537 weight1[nw1][n1-1]=Ann_weight_vector[k];
00538 nw1=nw1+1;
00539 }
00540 }
00541 else if(k>=all1){
00542 if((k-all1)%(hidneuron+1)==0){
00543 nwo=0;
00544 constanto[no]=Ann_weight_vector[k];
00545 no=no+1;
00546 }
00547 else {
00548 weighto[nwo]=Ann_weight_vector[k];
00549 nwo=nwo+1;
00550 }
00551 }
00552 }
00553
00554
00555
00556
00557
00558 for(Int_t i=0;i<hidneuron;i++){
00559 rin[i]=constant1[i];
00560 for(Int_t j=0;j<inneuron;j++){
00561 rin[i]=rin[i]+weight1[j][i]*out[j];
00562 }
00563 }
00564
00565 for(Int_t i=0;i<hidneuron;i++){
00566 out[i]=Sigmoid(rin[i]);
00567 }
00568
00569 prob=constanto[0];
00570 for(Int_t i=0;i<hidneuron;i++) prob=prob+weighto[i]*out[i];
00571
00572
00573 if(anninput.aTotlen>50) prob=0.;
00574
00575 return prob;
00576 }
00577
00578
00579 Float_t MadDpAnalysis::PID(Int_t event, Int_t method)
00580 {
00581
00582
00583 if(!LoadEvent(event)) return -10;
00584 Int_t track = -1;
00585 if(!LoadLargestTrackFromEvent(event,track)) return -10;
00586 if(method!=0) return -10;
00587
00588 Float_t ProbNC = 1.;
00589 Float_t ProbCC = 1.;
00590 Float_t PidCC = 0.;
00591
00592
00593 Float_t var1=eventSummary->plane.n;
00594 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00595 var1=ntpEvent->plane.end-ntpEvent->plane.beg+1;
00596 }
00597 Float_t var2=0;
00598 Float_t var3=0;
00599
00600
00601 Float_t phtrack=ntpTrack->ph.sigcor;
00602 Float_t phevent=ntpEvent->ph.sigcor;
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613 if (phevent>0) {var2=phtrack/phevent;}
00614
00615 if (ntpTrack->plane.n!=0) var3 = float(phtrack)
00616 /float(ntpTrack->plane.n);
00617
00618
00619 Float_t prob1=fLikeliHist[0]->GetBinContent(fLikeliHist[0]->FindBin(var1));
00620 Float_t prob2=fLikeliHist[2]->GetBinContent(fLikeliHist[2]->FindBin(var2));
00621 Float_t prob3=fLikeliHist[4]->GetBinContent(fLikeliHist[4]->FindBin(var3));
00622
00623 if (prob1==0) {prob1=0.0001;}
00624 if (prob2==0) {prob2=0.0001;}
00625 if (prob3==0) {prob3=0.0001;}
00626
00627 ProbCC=prob1*prob2*prob3;
00628
00629 prob1=fLikeliHist[1]->GetBinContent(fLikeliHist[1]->FindBin(var1));
00630 prob2=fLikeliHist[3]->GetBinContent(fLikeliHist[3]->FindBin(var2));
00631 prob3=fLikeliHist[5]->GetBinContent(fLikeliHist[5]->FindBin(var3));
00632
00633 if (prob1==0) {prob1=0.0001;}
00634 if (prob2==0) {prob2=0.0001;}
00635 if (prob3==0) {prob3=0.0001;}
00636
00637 ProbNC=prob1*prob2*prob3;
00638
00639 PidCC=-(sqrt(-TMath::Log(ProbCC))-sqrt(-TMath::Log(ProbNC)));
00640
00641 return PidCC;
00642 }
00643
00644
00645 Float_t MadDpAnalysis::RecoAnalysisEnergy(Int_t event){
00646
00647 if(!LoadEvent(event)) return 0;
00648 int largestTrack = -1;
00649 LoadLargestTrackFromEvent(event,largestTrack);
00650 float emu = RecoMuEnergy(0,largestTrack);
00651 float ehad = RecoShwEnergy(event);
00652 float ereco = emu + ehad;
00653 return ereco;
00654
00655 }
00656
00657 void MadDpAnalysis::MakeMyFile(std::string tag,int tarpos, int pdfbinning){
00658
00659
00660
00661
00662 std::string savename = "DPHistos_" + tag + ".root";
00663 TFile save(savename.c_str(),"RECREATE");
00664 save.cd();
00665
00666
00667
00668
00669
00670
00671 Int_t evlbins=0;
00672
00673 if (pdfbinning==0) evlbins=60;
00674 else evlbins=300;
00675
00676
00677 TH1F *evlength_nc = new TH1F("Event length - nc",
00678 "Event length - nc",
00679 evlbins,0,600);
00680 evlength_nc->SetXTitle("Event length (planes)");
00681 TH1F *evlength_cc = new TH1F("Event length - cc",
00682 "Event length - cc",
00683 evlbins,0,600);
00684 evlength_cc->SetXTitle("Event length (planes)");
00685
00686
00687
00688 TH1F *phfrac_nc = new TH1F("Track ph frac - nc",
00689 "Track ph frac - nc",
00690 50,0,1);
00691 phfrac_nc->SetXTitle("pulse height fraction");
00692
00693
00694 TH1F *phfrac_cc = new TH1F("Track ph frac - cc",
00695 "Track ph frac - cc",
00696 50,0,1);
00697 phfrac_cc->SetXTitle("pulse height fraction");
00698
00699
00700
00701 Int_t phplanebins=0;
00702 Float_t phplanemax=0;
00703
00704 if (pdfbinning==0) {
00705 phplanebins=50;
00706 phplanemax=5000;
00707 }
00708 else {
00709 phplanebins=100;
00710 phplanemax=3000;
00711 }
00712
00713 TH1F *phplane_nc = new TH1F("ph per plane (nc)",
00714 "ph per plane (nc)",
00715 phplanebins,0,phplanemax);
00716 phplane_nc->SetXTitle("pulse height per plane");
00717 TH1F *phplane_cc = new TH1F("ph per plane (cc)",
00718 "ph per plane (cc)",
00719 phplanebins,0,phplanemax);
00720 phplane_cc->SetXTitle("pulse height per plane");
00721
00722
00723 evlength_nc->SetDirectory(&save);
00724 evlength_cc->SetDirectory(&save);
00725
00726 phfrac_nc->SetDirectory(&save);
00727 phfrac_cc->SetDirectory(&save);
00728
00729 phplane_nc->SetDirectory(&save);
00730 phplane_cc->SetDirectory(&save);
00731
00732
00733
00734
00735
00736
00737 for(int i=0;i<Nentries;i++){
00738
00739 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
00740 << "% done" << std::endl;
00741 if(!this->GetEntry(i)) continue;
00742
00743 Int_t nevent = eventSummary->nevent;
00744 if(nevent==0) continue;
00745
00746 Int_t trkcount=0;
00747 Int_t shwcount=0;
00748
00749 Int_t evtmin=0;
00750 Int_t evtmax=nevent;
00751
00752 if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) {
00753
00754 Int_t maxevent=0;
00755 Float_t maxph=0;
00756 for(int j=0;j<nevent;j++){
00757 if(!LoadEvent(j)) continue;
00758 Float_t evtpheight=ntpEvent->ph.sigcor;
00759 if (evtpheight>maxph) {
00760 maxph=evtpheight;
00761 maxevent=j;
00762 }
00763 }
00764 evtmin=maxevent;
00765 evtmax=maxevent+1;
00766 }
00767
00768
00769
00770 for(int j=evtmin;j<evtmax;j++){
00771
00772
00773 if(!LoadEvent(j)) continue;
00774
00775 Int_t ntrack = 0;
00776 ntrack=ntpEvent->ntrack;
00777 Int_t nshower = 0;
00778 nshower=ntpEvent->nshower;
00779
00780 Float_t totph=0;
00781
00782
00783 for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
00784 LoadTrack(ii);
00785
00786 totph+=ntpTrack->ph.sigcor;
00787 LoadTHTrack(ii);
00788
00789 }
00790
00791 trkcount+=ntrack;
00792
00793
00794 for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
00795 LoadShower(ii);
00796
00797 totph+=ntpShower->ph.sigcor;
00798 }
00799
00800 shwcount+=nshower;
00801
00802
00803
00804 Int_t cc_nc = 0;
00805 Int_t flavour = 0;
00806 Int_t detector = -1;
00807 Int_t is_fid = 0;
00808 Double_t neu_e = -1;
00809 Double_t weight = -1;
00810
00811
00812 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00813 detector = 2;
00814
00815 if (evlength_cc->GetNbinsX()==evlbins) {
00816 if (pdfbinning==0) {
00817 evlength_cc->SetBins(50,0,500);
00818 evlength_nc->SetBins(50,0,500);
00819 }
00820 else {
00821 evlength_cc->SetBins(250,0,500);
00822 evlength_nc->SetBins(250,0,500);
00823 }
00824 }
00825 is_fid = 1;
00826 if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||
00827 (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||
00828 sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)
00829 +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5) is_fid = 0;
00830 }
00831 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00832 detector = 1;
00833
00834 if (pdfbinning==1 && evlength_cc->GetBinLowEdge(250)>200) {
00835 evlength_cc->SetBins(150,0,300);
00836 evlength_nc->SetBins(150,0,300);
00837 }
00838 if (pdfbinning==0 && evlength_cc->GetBinLowEdge(50)>200) {
00839 evlength_cc->SetBins(60,0,300);
00840 evlength_nc->SetBins(60,0,300);
00841 }
00842
00843
00844
00845 is_fid = 1;
00846 if(ntpEvent->vtx.z<1 || ntpEvent->vtx.z>5 ||
00847 sqrt(((ntpEvent->vtx.x-1.4885)*(ntpEvent->vtx.x-1.4885)) +
00848 ((ntpEvent->vtx.y-0.1397)*(ntpEvent->vtx.y-0.1397)))>1) is_fid = 0;
00849 }
00850
00851
00852
00853 Int_t mcevent=-1;
00854
00855 if(LoadTruthForRecoTH(j,mcevent)) {
00856
00857
00858 cc_nc = IAction(mcevent);
00859 flavour = Flavour(mcevent);
00860 neu_e = TrueNuEnergy(mcevent);
00861
00862
00863 if(tarpos==2) weight=1;
00864 if(tarpos==1) weight=1;
00865 if(tarpos==3) weight=1;
00866 }
00867
00868 int track_index = -1;
00869 LoadLargestTrackFromEvent(j,track_index);
00870
00871 if (track_index!=-1) {
00872
00873
00874 if (is_fid) {
00875
00876 Float_t evlength=0;
00877 Float_t phfrac=0;
00878 Float_t phplane=0;
00879
00880 evlength=ntpEvent->plane.n;
00881 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00882 evlength=ntpEvent->plane.end-ntpEvent->plane.beg+1;
00883 }
00884
00885 Float_t phtrack=ntpTrack->ph.sigcor;
00886 Float_t phevent=ntpEvent->ph.sigcor;
00887 if (phevent>0) {phfrac=phtrack/phevent;}
00888
00889 Float_t trkplane=ntpTrack->plane.n;
00890 if (trkplane>0) {phplane=phtrack/trkplane;}
00891
00892 if(flavour==2 && cc_nc==1) {
00893 evlength_cc->Fill(evlength,weight);
00894 phfrac_cc->Fill(phfrac,weight);
00895 phplane_cc->Fill(phplane,weight);
00896
00897
00898
00899 }
00900
00901 else if(cc_nc==0) {
00902 evlength_nc->Fill(evlength,weight);
00903 phfrac_nc->Fill(phfrac,weight);
00904 phplane_nc->Fill(phplane,weight);
00905 }
00906
00907
00908 }
00909 }
00910 }
00911 }
00912
00913
00914
00915 save.Write();
00916 save.Close();
00917 }
00918
00919 Float_t MadDpAnalysis::SingleTimeFrame(Int_t snarlentry,Int_t nentries){
00920
00921 Int_t tf=0,tfbef=0,tfaft=0;
00922
00923
00924 if(this->GetEntry(snarlentry)) tf = ntpHeader->GetTimeFrame();
00925 if(snarlentry<(nentries-1) && this->GetEntry(snarlentry+1)) tfaft = ntpHeader->GetTimeFrame();
00926 if(snarlentry>0 && this->GetEntry(snarlentry-1)) tfbef = ntpHeader->GetTimeFrame();
00927
00928 Float_t singletimeframe=0;
00929 if(snarlentry<(nentries-1) && snarlentry>0 ){
00930 if(tf!=tfaft && tf!=tfbef) singletimeframe=1;
00931 }
00932 if(snarlentry==(nentries-1)) {
00933 if(tf!=tfbef) singletimeframe=1;
00934 }
00935 if(snarlentry==0) {
00936 if(tf!=tfaft) singletimeframe=1;
00937 }
00938
00939 this->GetEntry(snarlentry);
00940 return singletimeframe;
00941
00942 }
00943
00944 void MadDpAnalysis::CreatePAN(std::string tag, Int_t scanfilter, Int_t mcneartype){
00945
00946
00947
00948
00949 std::string savename = "PAN_" + tag + ".root";
00950 TFile save(savename.c_str(),"RECREATE");
00951 save.cd();
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964 SpillInfo pppinfo;
00965
00966
00967
00968
00969 Float_t ann_aTotrk =0.;
00970 Float_t ann_aTotstp =0.;
00971 Float_t ann_aTotph =0.;
00972 Float_t ann_aTotlen =0.;
00973 Float_t ann_aPhperpl =0.;
00974 Float_t ann_aPhperstp =0.;
00975
00976
00977 Float_t ann_aTrkpass =0.;
00978 Float_t ann_aTrkph =0.;
00979 Float_t ann_aTrklen =0.;
00980 Float_t ann_aTrkphperpl =0.;
00981 Float_t ann_aTrkphper =0.;
00982 Float_t ann_aTrkplu =0.;
00983 Float_t ann_aTrkplv =0.;
00984 Float_t ann_aTrkstp =0.;
00985
00986 Float_t ann_aShwph =0.;
00987 Float_t ann_aShwstp =0.;
00988 Float_t ann_aShwdig =0.;
00989 Float_t ann_aShwpl =0.;
00990 Float_t ann_aShwphper =0.;
00991 Float_t ann_aShwphperpl =0.;
00992 Float_t ann_aShwphperdig =0.;
00993 Float_t ann_aShwphperstp =0.;
00994 Float_t ann_aShwplu =0.;
00995 Float_t ann_aShwplv =0.;
00996 Float_t ann_aPh3 =0.;
00997 Float_t ann_aPh6 =0.;
00998 Float_t ann_aPhlast =0.;
00999 Float_t ann_aPhcommon =0.;
01000
01001
01002
01003 Float_t true_enu;
01004 Float_t true_pxnu;
01005 Float_t true_pynu;
01006 Float_t true_pznu;
01007 Float_t true_etgt;
01008 Float_t true_pxtgt;
01009 Float_t true_pytgt;
01010 Float_t true_pztgt;
01011 Float_t true_emu;
01012 Float_t true_eshw;
01013 Float_t true_x;
01014 Float_t true_y;
01015 Float_t true_q2;
01016 Float_t true_w2;
01017 Float_t true_dircosneu;
01018 Float_t true_dircosz;
01019 Float_t true_vtxx;
01020 Float_t true_vtxy;
01021 Float_t true_vtxz;
01022
01023
01024 Int_t flavour;
01025 Int_t process;
01026 Int_t nucleus;
01027 Int_t cc_nc;
01028 Int_t initial_state;
01029 Int_t had_fs;
01030
01031
01032
01033
01034
01035 Int_t detector;
01036 Int_t run;
01037 Int_t snarl;
01038 Int_t nevent;
01039 Int_t event;
01040 Int_t subrun;
01041 Int_t trigsrc;
01042 Int_t mcevent;
01043 Int_t ntrack;
01044 Int_t nshower;
01045
01046 Int_t is_fid;
01047
01048 Int_t is_fid_dp;
01049 Int_t is_fid_ns;
01050
01051 Int_t is_cev;
01052 Int_t is_cont_trk;
01053 Int_t is_mc;
01054 Int_t pass_basic;
01055 Float_t pid0;
01056 Float_t pid1;
01057 Float_t pid2;
01058 Float_t annpid;
01059 Float_t annpid_f1p1;
01060 Float_t annpid_f1p2;
01061 Float_t annpid_f1p3;
01062 Float_t annpid_f2p1;
01063 Float_t annpid_f2p2;
01064 Float_t annpid_f2p3;
01065
01066 Int_t pass;
01067
01068 Float_t reco_enu;
01069 Float_t reco_emu;
01070 Float_t reco_eshw;
01071
01072 Int_t pass_fd_qualcuts;
01073 Int_t litag;
01074
01075
01076 Float_t reco_eshwcc;
01077 Float_t reco_eshwnc;
01078 Float_t reco_eshwccw;
01079 Float_t reco_eshwncw;
01080
01081 Float_t reco_eshw_sqrt;
01082 Float_t reco_qe_enu;
01083 Float_t reco_qe_q2;
01084 Float_t reco_y;
01085 Float_t reco_dircosneu;
01086 Float_t reco_dircosz;
01087 Float_t reco_vtxx;
01088 Float_t reco_vtxy;
01089 Float_t reco_vtxz;
01090
01091 Float_t evtlength;
01092 Float_t evtph;
01093 Float_t trklength;
01094 Float_t trkmomrange;
01095 Float_t trkqp;
01096 Float_t trkeqp;
01097 Float_t trkmombest;
01098 Float_t trkphfrac;
01099 Float_t trkphplane;
01100 Float_t trkds;
01101 Float_t rawph;
01102
01103 Float_t trkendz;
01104 Float_t trkendx;
01105 Float_t trkendy;
01106 Float_t trkendu;
01107 Float_t trkendv;
01108 Float_t trkvtxz;
01109 Float_t trkvtxx;
01110 Float_t trkvtxy;
01111 Float_t trkvtxu;
01112 Float_t trkvtxv;
01113 Float_t trkplbu;
01114 Float_t trkplbv;
01115 Float_t trkpleu;
01116 Float_t trkplev;
01117 Float_t trkple;
01118 Float_t trkplb;
01119 Int_t trkfitpass;
01120 Int_t totdigits,totstrips;
01121 Int_t trkdigits,trkstrips;
01122 Float_t trkph,trkchi2,trkndof;
01123 Int_t trkdiffuv;
01124
01125 Float_t shwph;
01126 Int_t shwdigits, shwstrips;
01127 Float_t shwvtxz;
01128 Float_t shwvtxx;
01129 Float_t shwvtxy;
01130 Float_t shwvtxu;
01131 Float_t shwvtxv;
01132
01133
01134
01135 Int_t shwncluster;
01136
01137
01138 Int_t scandecision;
01139 ScanList scan;
01140
01141
01142 Double_t evttimemin;
01143 Double_t evttimemax;
01144 Double_t snarlt;
01145 Double_t litime;
01146
01147
01148 Float_t pot,horn,tar,vvpos,hhpos,magn;
01149
01150 Float_t potdb,horndb,tardb,vvposdb,hhposdb,vvwidthdb,hhwidthdb;
01151 Double_t timedb;
01152
01153 Double_t neartdb,fartdb,neardifftdb;
01154
01155
01156 Float_t mcflux_tpx;
01157 Float_t mcflux_tpy;
01158 Float_t mcflux_tpz;
01159 Float_t mcflux_tptype;
01160
01161
01162 Float_t singletimeframe;
01163
01164 SpillInfoBlock *spinfo = new SpillInfoBlock();
01165 spinfo->Zero();
01166 singletimeframe=-999;
01167
01168 Int_t day,month,year;
01169 Int_t pass_beamcuts;
01170
01171 TTree *tree = new TTree("pan","pan");
01172
01173 tree->Branch("true_enu",&true_enu,"true_enu/F",512000);
01174 tree->Branch("true_pxnu",&true_pxnu,"true_pxnu/F",512000);
01175 tree->Branch("true_pynu",&true_pynu,"true_pynu/F",512000);
01176 tree->Branch("true_pznu",&true_pznu,"true_pznu/F",512000);
01177 tree->Branch("true_etgt",&true_etgt,"true_etgt/F",512000);
01178 tree->Branch("true_pxtgt",&true_pxtgt,"true_pxtgt/F",512000);
01179 tree->Branch("true_pytgt",&true_pytgt,"true_pytgt/F",512000);
01180 tree->Branch("true_pztgt",&true_pztgt,"true_pztgt/F",512000);
01181 tree->Branch("true_emu",&true_emu,"true_emu/F",512000);
01182 tree->Branch("true_eshw",&true_eshw,"true_eshw/F",512000);
01183 tree->Branch("true_x",&true_x,"true_x/F",512000);
01184 tree->Branch("true_y",&true_y,"true_y/F",512000);
01185 tree->Branch("true_q2",&true_q2,"true_q2/F",512000);
01186 tree->Branch("true_w2",&true_w2,"true_w2/F",512000);
01187 tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",512000);
01188 tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",512000);
01189 tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",512000);
01190 tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",512000);
01191 tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",512000);
01192
01193 tree->Branch("flavour",&flavour,"flavour/I",512000);
01194 tree->Branch("process",&process,"process/I",512000);
01195 tree->Branch("nucleus",&nucleus,"nucleus/I",512000);
01196 tree->Branch("cc_nc",&cc_nc,"cc_nc/I",512000);
01197 tree->Branch("initial_state",&initial_state,"initial_state/I",512000);
01198 tree->Branch("had_fs",&had_fs,"had_fs/I",512000);
01199
01200
01201
01202 tree->Branch("detector",&detector,"detector/I",512000);
01203 tree->Branch("run",&run,"run/I",512000);
01204 tree->Branch("snarl",&snarl,"snarl/I",512000);
01205 tree->Branch("event",&event,"event/I",512000);
01206 tree->Branch("mcevent",&mcevent,"mcevent/I",512000);
01207 tree->Branch("ntrack",&ntrack,"ntrack/I",512000);
01208 tree->Branch("nshower",&nshower,"nshower/I",512000);
01209 tree->Branch("subrun",&subrun,"subrun/I",512000);
01210 tree->Branch("trigsrc",&trigsrc,"trigsrc/I",512000);
01211 tree->Branch("litime",&litime,"litime/D",512000);
01212 tree->Branch("is_fid",&is_fid,"is_fid/I",512000);
01213 tree->Branch("is_fid_dp",&is_fid_dp,"is_fid_dp/I",512000);
01214 tree->Branch("is_fid_ns",&is_fid_ns,"is_fid_ns/I",512000);
01215
01216
01217 tree->Branch("pass_fd_qualcuts", &pass_fd_qualcuts, "pass_fd_qualcuts/I");
01218 tree->Branch("litag", &litag, "litag/I");
01219
01220 tree->Branch("is_cev",&is_cev,"is_cev/I",512000);
01221 tree->Branch("is_cont_trk",&is_cont_trk,"is_cont_trk/I",512000);
01222 tree->Branch("is_mc",&is_mc,"is_mc/I",512000);
01223 tree->Branch("pass_basic",&pass_basic,"pass_basic/I",512000);
01224 tree->Branch("pid0",&pid0,"pid0/F",512000);
01225 tree->Branch("pid1",&pid1,"pid1/F",512000);
01226 tree->Branch("pid2",&pid2,"pid2/F",512000);
01227 tree->Branch("annpid",&annpid,"annpid/F",512000);
01228 tree->Branch("annpid_f1p1",&annpid_f1p1,"annpid_f1p1/F",512000);
01229 tree->Branch("annpid_f1p2",&annpid_f1p2,"annpid_f1p2/F",512000);
01230 tree->Branch("annpid_f1p3",&annpid_f1p3,"annpid_f1p3/F",512000);
01231 tree->Branch("annpid_f2p1",&annpid_f2p1,"annpid_f2p1/F",512000);
01232 tree->Branch("annpid_f2p2",&annpid_f2p2,"annpid_f2p2/F",512000);
01233 tree->Branch("annpid_f2p3",&annpid_f2p3,"annpid_f2p3/F",512000);
01234 tree->Branch("pass",&pass,"pass/I",512000);
01235
01236 tree->Branch("reco_enu",&reco_enu,"reco_enu/F",512000);
01237 tree->Branch("reco_emu",&reco_emu,"reco_emu/F",512000);
01238 tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",512000);
01239 tree->Branch("reco_eshw_sqrt",&reco_eshw_sqrt,"reco_eshw_sqrt/F",512000);
01240 tree->Branch("reco_eshwcc", &reco_eshwcc, "reco_eshwcc/F",512000);
01241 tree->Branch("reco_eshwnc", &reco_eshwnc, "reco_eshwnc/F",512000);
01242 tree->Branch("reco_eshwccw", &reco_eshwccw, "reco_eshwccw/F",512000);
01243 tree->Branch("reco_eshwncw", &reco_eshwncw, "reco_eshwncw/F",512000);
01244 tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",512000);
01245 tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",512000);
01246 tree->Branch("reco_y",&reco_y,"reco_y/F",512000);
01247 tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",512000);
01248 tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",512000);
01249 tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",512000);
01250 tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",512000);
01251 tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",512000);
01252
01253 tree->Branch("evtlength",&evtlength,"evtlength/F",512000);
01254 tree->Branch("evtph",&evtph,"evtph/F",512000);
01255 tree->Branch("trklength",&trklength,"trklength/F",512000);
01256 tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",512000);
01257 tree->Branch("trkqp",&trkqp,"trkqp/F",512000);
01258 tree->Branch("trkeqp",&trkeqp,"trkeqp/F",512000);
01259 tree->Branch("trkmombest",&trkmombest,"trkmombest/F",512000);
01260 tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",512000);
01261 tree->Branch("trkphplane",&trkphplane,"trkphplane/F",512000);
01262 tree->Branch("rawph",&rawph,"rawph/F",512000);
01263
01264 tree->Branch("trkendz",&trkendz,"trkendz/F",512000);
01265 tree->Branch("trkendu",&trkendu,"trkendu/F",512000);
01266 tree->Branch("trkendv",&trkendv,"trkendv/F",512000);
01267 tree->Branch("trkendx",&trkendx,"trkendx/F",512000);
01268 tree->Branch("trkendy",&trkendy,"trkendy/F",512000);
01269 tree->Branch("trkplbu",&trkplbu,"trkplbu/F",512000);
01270 tree->Branch("trkplbv",&trkplbv,"trkplbv/F",512000);
01271 tree->Branch("trkpleu",&trkpleu,"trkpleu/F",512000);
01272 tree->Branch("trkplev",&trkplev,"trkplev/F",512000);
01273 tree->Branch("trkple",&trkple,"trkple/F",512000);
01274 tree->Branch("trkplb",&trkplb,"trkplb/F",512000);
01275
01276 tree->Branch("trkvtxz",&trkvtxz,"trkvtxz/F",512000);
01277 tree->Branch("trkvtxu",&trkvtxu,"trkvtxu/F",512000);
01278 tree->Branch("trkvtxv",&trkvtxv,"trkvtxv/F",512000);
01279 tree->Branch("trkvtxx",&trkvtxx,"trkvtxx/F",512000);
01280 tree->Branch("trkvtxy",&trkvtxy,"trkvtxy/F",512000);
01281 tree->Branch("trkds", &trkds,"trkds/F",512000);
01282
01283 tree->Branch("evttimemin",&evttimemin,"evttimemin/D");
01284 tree->Branch("evttimemax",&evttimemax,"evttimemax/D");
01285
01286 tree->Branch("totdigits",&totdigits,"totdigits/I");
01287 tree->Branch("totstrips",&totstrips,"totstrips/I");
01288 tree->Branch("trkdigits",&trkdigits,"trkdigits/I");
01289 tree->Branch("trkstrips",&trkstrips,"trkstrips/I");
01290 tree->Branch("trkph",&trkph,"trkph/F");
01291 tree->Branch("trkdiffuv",&trkdiffuv,"trkdiffuv/I",512000);
01292
01293 tree->Branch("trkchi2",&trkchi2,"trkchi2/F");
01294 tree->Branch("trkndof",&trkndof,"trkndof/F");
01295 tree->Branch("trkfitpass",&trkfitpass,"trkfitpass/I",512000);
01296
01297 tree->Branch("shwdigits",&shwdigits,"shwdigits/I");
01298 tree->Branch("shwstrips",&shwstrips,"shwstrips/I");
01299 tree->Branch("shwph",&shwph,"shwph/F");
01300
01301 tree->Branch("shwvtxz",&shwvtxz,"shwvtxz/F",512000);
01302 tree->Branch("shwvtxu",&shwvtxu,"shwvtxu/F",512000);
01303 tree->Branch("shwvtxv",&shwvtxv,"shwvtxv/F",512000);
01304 tree->Branch("shwvtxx",&shwvtxx,"shwvtxx/F",512000);
01305 tree->Branch("shwvtxy",&shwvtxy,"shwvtxy/F",512000);
01306 tree->Branch("shwncluster",&shwncluster,"shwncluster/I");
01307
01308
01309 tree->Branch("scandecision",&scandecision,"scandecision/I");
01310
01311
01312 tree->Branch("pot", &pot, "pot/F",512000);
01313 tree->Branch("horn", &horn, "horn/F",512000);
01314 tree->Branch("tar", &tar, "tar/F",512000);
01315 tree->Branch("vvpos", &vvpos, "vvpos/F",512000);
01316 tree->Branch("hhpos", &hhpos, "hhpos/F",512000);
01317 tree->Branch("magn", &magn, "magn/F",512000);
01318 tree->Branch("singletimeframe",&singletimeframe,"singletimeframe/F",512000);
01319
01320
01321 tree->Branch("potdb", &potdb, "potdb/F",512000);
01322 tree->Branch("horndb", &horndb, "horndb/F",512000);
01323 tree->Branch("tardb", &tardb, "tardb/F",512000);
01324 tree->Branch("vvposdb", &vvposdb, "vvposdb/F",512000);
01325 tree->Branch("hhposdb", &hhposdb, "hhposdb/F",512000);
01326 tree->Branch("vvwidthdb", &vvwidthdb, "vvwidthdb/F",512000);
01327 tree->Branch("hhwidthdb", &hhwidthdb, "hhwidthdb/F",512000);
01328 tree->Branch("timedb", &timedb, "timedb/D");
01329
01330
01331 tree->Branch("neartdb", &neartdb, "neartdb/D");
01332 tree->Branch("fartdb", &fartdb, "fartdb/D");
01333 tree->Branch("neardifftdb", &neardifftdb,"neardifftdb/D");
01334 tree->Branch("snarlt", &snarlt, "snarlt/D");
01335
01336
01337 tree->Branch("mcflux_tpx", &mcflux_tpx, "mcflux_tpx/F");
01338 tree->Branch("mcflux_tpy", &mcflux_tpy, "mcflux_tpy/F");
01339 tree->Branch("mcflux_tpz", &mcflux_tpz, "mcflux_tpz/F");
01340 tree->Branch("mcflux_tptype", &mcflux_tptype, "mcflux_tptype/F");
01341
01342
01343
01344 tree->Branch("ann_aTotrk", &ann_aTotrk, "ann_aTotrk/F");
01345 tree->Branch("ann_aTotstp", &ann_aTotstp, "ann_aTotstp/F");
01346 tree->Branch("ann_aTotph", &ann_aTotph, "ann_aTotph/F");
01347 tree->Branch("ann_aTotlen", &ann_aTotlen, "ann_aTotlen/F");
01348 tree->Branch("ann_aPhperpl", &ann_aPhperpl, "ann_aPhperpl/F");
01349 tree->Branch("ann_aPhperstp", &ann_aPhperstp, "ann_aPhperstp/F");
01350 tree->Branch("ann_aTrkpass", &ann_aTrkpass, "ann_aTrkpass/F");
01351 tree->Branch("ann_aTrkph", &ann_aTrkph, "ann_aTrkph/F");
01352 tree->Branch("ann_aTrklen", &ann_aTrklen, "ann_aTrklen/F");
01353 tree->Branch("ann_aTrkphperpl", &ann_aTrkphperpl, "ann_aTrkphperpl/F");
01354 tree->Branch("ann_aTrkphper", &ann_aTrkphper, "ann_aTrkphper/F");
01355 tree->Branch("ann_aTrkplu", &ann_aTrkplu, "ann_aTrkplu/F");
01356 tree->Branch("ann_aTrkplv", &ann_aTrkplv, "ann_aTrkplv/F");
01357 tree->Branch("ann_aTrkstp", &ann_aTrkstp, "ann_aTrkstp/F");
01358 tree->Branch("ann_aShwph", &ann_aShwph, "ann_aShwph/F");
01359 tree->Branch("ann_aShwstp", &ann_aShwstp, "ann_aShwstp/F");
01360 tree->Branch("ann_aShwdig", &ann_aShwdig, "ann_aShwdig/F");
01361 tree->Branch("ann_aShwpl", &ann_aShwpl, "ann_aShwpl/F");
01362 tree->Branch("ann_aShwphper", &ann_aShwphper, "ann_aShwphper/F");
01363 tree->Branch("ann_aShwphperpl", &ann_aShwphperpl, "ann_aShwphperpl/F");
01364 tree->Branch("ann_aShwphperdig", &ann_aShwphperdig,"ann_aShwphperdig/F");
01365 tree->Branch("ann_aShwphperstp", &ann_aShwphperstp,"ann_aShwphperstp/F");
01366 tree->Branch("ann_aShwplu", &ann_aShwplu, "ann_aShwplu/F");
01367 tree->Branch("ann_aShwplv", &ann_aShwplv, "ann_aShwplv/F");
01368 tree->Branch("ann_aPh3", &ann_aPh3, "ann_aPh3/F");
01369 tree->Branch("ann_aPh6", &ann_aPh6, "ann_aPh6/F");
01370 tree->Branch("ann_aPhlast", &ann_aPhlast, "ann_aPhlast/F");
01371 tree->Branch("ann_aPhcommon", &ann_aPhcommon, "ann_aPhcommon/F");
01372
01373 tree->Branch("day", &day, "day/I");
01374 tree->Branch("month", &month, "month/I");
01375 tree->Branch("year", &year, "year/I");
01376 tree->Branch("pass_beamcuts", &pass_beamcuts, "pass_beamcuts/I");
01377 Bool_t first_ann = true;
01378 for(int i=0;i<Nentries;i++){
01379
01380 if(i%400==0) std::cout << float(i)*100./float(Nentries)
01381 << "% done" << std::endl;
01382
01383 if(!this->GetEntry(i)) continue;
01384
01385 nevent = eventSummary->nevent;
01386 if(nevent==0) continue;
01387
01388 run = ntpHeader->GetRun();
01389 snarl = ntpHeader->GetSnarl();
01390 subrun = ntpHeader->GetSubRun();
01391 trigsrc = ntpHeader->GetTrigSrc();
01392 day = eventSummary->date.day;
01393 month = eventSummary->date.month;
01394 year = eventSummary->date.year;
01395
01396 Double_t snarltime =ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9);
01397 Int_t month = eventSummary->date.month;
01398 litime = eventSummary->litime;
01399
01400 pass_fd_qualcuts =passfail(snarltime);
01401
01402 litag=0;
01403 Int_t numshower=eventSummary->nshower;
01404 Float_t allph=eventSummary->ph.raw;
01405 Int_t plbeg=eventSummary->plane.beg;
01406 Int_t plend=eventSummary->plane.end;
01407
01408 if (numshower) {
01409 if (allph>1e6) litag+=10;
01410
01411 if (((plbeg==1 || plbeg==2) && (plend==63 || plend==64)) ||
01412 ((plbeg==65 || plbeg==66) && (plend==127 || plend==128)) ||
01413 ((plbeg==129 || plbeg==130) && (plend==191 || plend==192)) ||
01414 ((plbeg==193 || plbeg==194) && (plend==247 || plend==248))) litag++;
01415 if (((plbeg==250 || plbeg==251) && (plend==312 || plend==313)) ||
01416 ((plbeg==314 || plbeg==315) && (plend==376 || plend==377)) ||
01417 ((plbeg==378 || plbeg==379) && (plend==440 || plend==441)) ||
01418 ((plbeg==442 || plbeg==443) && (plend==484 || plend==485))) litag++;
01419 }
01420
01421
01422
01423
01424
01425 scandecision=scan.GetDecision(run,snarl);
01426
01427
01428 if (scanfilter==0 || scandecision>0) {
01429
01430
01431
01432 if(ntpHeader->GetVldContext().GetSimFlag()!=4){
01433 pppinfo.GetSpillInfo(month,year,snarltime,*spinfo);
01434
01435 pot =spinfo->GetPot();
01436 horn =spinfo->GetHorn();
01437 hhpos =spinfo->GetHpos();
01438 vvpos =spinfo->GetVpos();
01439 tar =spinfo->GetTgtpos();
01440 magn =spinfo->GetMagnet();
01441
01442 potdb = -999;
01443 horndb = -999;
01444 tardb = -999;
01445 vvposdb = -999;
01446 hhposdb = -999;
01447 vvwidthdb = -999;
01448 hhwidthdb = -999;
01449 timedb = -999;
01450 neartdb = -999;
01451 fartdb = -999;
01452 neardifftdb = -999;
01453
01454
01455 Detector::Detector_t detectort;
01456 SimFlag::SimFlag_t simflag;
01457 VldTimeStamp timestamp;
01458 VldContext cx;
01459 cx = ntpHeader->GetVldContext();
01460 detectort = ntpHeader->GetVldContext().GetDetector();
01461 simflag = ntpHeader->GetVldContext().GetSimFlag();
01462 timestamp = ntpHeader->GetVldContext().GetTimeStamp();
01463
01464
01465 BDSpillAccessor &spillAccess=BDSpillAccessor::Get();
01466 const BeamMonSpill *spillInfoDB = spillAccess.LoadSpill(timestamp);
01467 if(spillInfoDB) {
01468 if(spillInfoDB->fTortgt !=0.) potdb = spillInfoDB->fTortgt;
01469 else if(spillInfoDB->fTrtgtd !=0.) potdb = spillInfoDB->fTrtgtd;
01470 else if(spillInfoDB->fTor101 !=0.) potdb = spillInfoDB->fTor101;
01471 else potdb = -999;
01472
01473 horndb = spillInfoDB->fHornCur;
01474 tardb = (int)spillInfoDB->GetStatusBits().target_in;
01475 vvposdb = spillInfoDB->fTargBpmY[0];
01476 hhposdb = spillInfoDB->fTargBpmX[0];
01477 vvwidthdb = spillInfoDB->fProfWidY;
01478 hhwidthdb = spillInfoDB->fProfWidX;
01479 timedb = spillInfoDB->SpillTime();
01480
01481
01482 if (potdb>0.5 && fabs(snarltime-timedb)<1 &&
01483 (hhposdb<0 && hhposdb>-0.002) &&
01484 (vvposdb>0 && vvposdb<0.002) &&
01485 (hhwidthdb>0.0001 && hhwidthdb<0.0022) &&
01486 (vvwidthdb>0.0001 && vvwidthdb<0.005) &&
01487 (horndb<-155 && horndb>-200)) pass_beamcuts=1;
01488
01489 }
01490
01491
01492
01493 SpillTimeFinder& spillfinder = SpillTimeFinder::Instance();
01494 fartdb = spillfinder.GetTimeOfNearestSpill(cx);
01495 const SpillTimeND& spnd = spillfinder.GetNearestSpill(cx);
01496 neartdb = spnd.GetTimeStamp().GetSec()+(spnd.GetTimeStamp().GetNanoSec()/1e9);
01497 neardifftdb = spillfinder.GetTimeToNearestSpill(cx);
01498 snarlt = ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9);
01499
01500 }
01501
01502
01503
01504 Int_t trkcount=0;
01505 Int_t shwcount=0;
01506 Float_t totph=0;
01507
01508 Int_t evtmin=0;
01509 Int_t evtmax=nevent;
01510
01511 if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) {
01512
01513
01514 Int_t maxevent=0;
01515 Float_t maxph=0;
01516 for(int j=0;j<nevent;j++){
01517 if(!LoadEvent(j)) continue;
01518
01519
01520
01521 Float_t evtpheight=ntpEvent->ph.sigcor;
01522 if (evtpheight>maxph) {
01523 maxph=evtpheight;
01524 maxevent=j;
01525 }
01526 }
01527
01528 evtmin=maxevent;
01529 evtmax=maxevent+1;
01530 }
01531
01532
01533
01534
01535 for(int j=evtmin;j<evtmax;j++){
01536
01537 if(!LoadEvent(j)) continue;
01538
01539 totph = 0;
01540
01541
01542 event = j;
01543 ntrack = ntpEvent->ntrack;
01544 nshower = ntpEvent->nshower;
01545 totdigits=ntpEvent->ndigit;
01546 totstrips=ntpEvent->nstrip;
01547
01548
01549
01550 mcflux_tpx=0; mcflux_tpy=0; mcflux_tpz=0; mcflux_tptype=0;
01551 true_enu = 0; true_emu = 0; true_eshw = 0;
01552 true_pxnu = 0; true_pynu = 0; true_pznu = 0;
01553 true_x = 0; true_y = 0; true_q2 = 0; true_w2 = 0;
01554 true_dircosneu = 0; true_dircosz = 0;
01555 true_vtxx = 0; true_vtxy = 0; true_vtxz = 0;
01556 flavour = 0; process = 0; nucleus = 0; cc_nc = 0;
01557 initial_state = 0; had_fs = 0;
01558 true_etgt = 0; true_pxtgt = 0; true_pytgt = 0; true_pztgt = 0;
01559 detector = -1; mcevent = -1;
01560 is_fid = 0; is_cev = 0; is_mc = 0; pass_basic = 0;
01561 pid0 = -999; pid1 = -999; pid2 = -999; pass = 0;
01562 reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_eshw_sqrt = 0;
01563 reco_qe_enu = 0; reco_qe_q2 = 0; reco_y = 0; reco_dircosneu = 0;
01564 reco_dircosz = 0;
01565 reco_vtxx = -999; reco_vtxy = -999; reco_vtxz = -999;
01566 evtlength = -1; trklength = -1; trkphfrac = -1; trkphplane = -1;
01567 trkmomrange = -999; trkqp = -999; trkeqp = -999;
01568 trkendz=100.; trkendu=100.; trkendv=100.; trkendx=100.; trkendy=100.;
01569 trkendz=999; trkendu=999; trkendv=999; trkendx=999;
01570 trkendy=-999;trkvtxz=-999; trkvtxu=-999; trkvtxv=-999; trkvtxx=-999; trkvtxy=-999;
01571 trkds=-1; evttimemin=-1; evttimemax=-1;
01572 evtph=0.;
01573 annpid=-999.;
01574 annpid_f1p1=-999.;annpid_f1p2=-999.;annpid_f1p3=-999.;
01575 annpid_f2p1=-999.;annpid_f2p2=-999.;annpid_f2p3=-999.;
01576 reco_eshwcc=0.;reco_eshwnc=0.;reco_eshwccw=0.;reco_eshwncw=0.;
01577 is_fid_dp=0;is_fid_ns=0;
01578 is_cont_trk=0; trkmombest=-999.;
01579
01580 trkfitpass=0;
01581 trkdiffuv=0;
01582 trkdigits=0; trkstrips=0; trkph=0; trkndof=0; trkchi2=0;
01583 shwdigits=0; shwstrips=0; shwph=0;
01584 shwvtxu=0;
01585 shwvtxv=0;
01586 shwvtxx=0;
01587 shwvtxy=0;
01588 shwvtxz=0;
01589 shwncluster=0;
01590
01591
01592 ann_aTotrk =0.;
01593 ann_aTotstp =0.;
01594 ann_aTotph =0.;
01595 ann_aTotlen =0.;
01596 ann_aPhperpl =0.;
01597 ann_aPhperstp =0.;
01598
01599
01600 ann_aTrkpass =0.;
01601 ann_aTrkph =0.;
01602 ann_aTrklen =0.;
01603 ann_aTrkphperpl =0.;
01604 ann_aTrkphper =0.;
01605 ann_aTrkplu =0.;
01606 ann_aTrkplv =0.;
01607 ann_aTrkstp =0.;
01608
01609 ann_aShwph =0.;
01610 ann_aShwstp =0.;
01611 ann_aShwdig =0.;
01612 ann_aShwpl =0.;
01613 ann_aShwphper =0.;
01614 ann_aShwphperpl =0.;
01615 ann_aShwphperdig =0.;
01616 ann_aShwphperstp =0.;
01617 ann_aShwplu =0.;
01618 ann_aShwplv =0.;
01619 ann_aPh3 =0.;
01620 ann_aPh6 =0.;
01621 ann_aPhlast =0.;
01622 ann_aPhcommon =0.;
01623
01624 rawph=ntpEvent->ph.raw;
01625
01626 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01627 detector = 2;
01628 totph=eventSummary->ph.sigcor;
01629
01630 is_fid = 1;
01631 if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
01632 }
01633 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01634 detector = 1;
01635
01636
01637 for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
01638 LoadTrack(ii);
01639 totph+=ntpTrack->ph.sigcor;
01640 }
01641 trkcount+=ntrack;
01642
01643 for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
01644 LoadShower(ii);
01645 totph+=ntpShower->ph.sigcor;
01646 }
01647 shwcount+=nshower;
01648
01649
01650 is_fid = 1;
01651 if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
01652 }
01653
01654 if (detector==2) is_fid_dp=PassFidNew(ntpEvent->index);
01655 else is_fid_dp=PassFid(ntpEvent->index);
01656
01657 if (detector==2) is_fid_ns=PassFidNewN(ntpEvent->index);
01658 else is_fid_ns=PassFid(ntpEvent->index);
01659
01660
01661 if(ntpHeader->GetVldContext().GetSimFlag()==4){
01662 if(LoadTruthForRecoTH(j,mcevent)) {
01663 Float_t *vtx = TrueVtx(mcevent);
01664 Float_t *nu3mom = TrueNuMom(mcevent);
01665 Float_t *targ4vec = Target4Vector(mcevent);
01666
01667 is_mc = 1;
01668 nucleus = Nucleus(mcevent);
01669 flavour = Flavour(mcevent);
01670 initial_state = Initial_state(mcevent);
01671 had_fs = HadronicFinalState(mcevent);
01672 process = IResonance(mcevent);
01673 cc_nc = IAction(mcevent);
01674 true_enu = TrueNuEnergy(mcevent);
01675 true_pxnu = nu3mom[0];
01676 true_pynu = nu3mom[1];
01677 true_pznu = nu3mom[2];
01678 true_emu = TrueMuEnergy(mcevent);
01679 true_eshw = TrueShwEnergy(mcevent);
01680 true_x = X(mcevent);
01681 true_y = Y(mcevent);
01682 true_q2 = Q2(mcevent);
01683 true_w2 = W2(mcevent);
01684 true_dircosz = TrueMuDCosZ(mcevent);
01685 true_dircosneu = TrueMuDCosNeu(mcevent);
01686 true_vtxx = vtx[0];
01687 true_vtxy = vtx[1];
01688 true_vtxz = vtx[2];
01689 true_etgt = targ4vec[3];
01690 true_pxtgt = targ4vec[0];
01691 true_pytgt = targ4vec[1];
01692 true_pztgt = targ4vec[2];
01693
01694
01695
01696
01697
01698
01699
01700 if (LoadTruth(mcevent)) {
01701 mcflux_tpx=ntpTruth->flux.tpx;
01702 mcflux_tpy=ntpTruth->flux.tpy;
01703 mcflux_tpz=ntpTruth->flux.tpz;
01704 mcflux_tptype=ntpTruth->flux.tptype;
01705 }
01706
01707 delete [] vtx;
01708 delete [] nu3mom;
01709 delete [] targ4vec;
01710 }
01711 }
01712
01713 if(PassBasicCuts(event)) {pass_basic=1;}
01714
01715
01716 reco_vtxx = ntpEvent->vtx.x;
01717 reco_vtxy = ntpEvent->vtx.y;
01718 reco_vtxz = ntpEvent->vtx.z;
01719 evtlength = ntpEvent->plane.end-ntpEvent->plane.beg+1;
01720 evtph = ntpEvent->ph.sigcor;
01721
01722
01723
01724
01725 AnnInputBlock anninput=AnnVar(event);
01726
01727 ann_aTotrk = anninput.aTotrk;
01728 ann_aTotstp = anninput.aTotstp;
01729 ann_aTotph = anninput.aTotph;
01730 ann_aTotlen = anninput.aTotlen;
01731 ann_aPhperpl = anninput.aPhperpl;
01732 ann_aPhperstp = anninput.aPhperstp;
01733
01734
01735 ann_aTrkpass = anninput.aTrkpass;
01736 ann_aTrkph = anninput.aTrkph;
01737 ann_aTrklen = anninput.aTrklen;
01738 ann_aTrkphperpl = anninput.aTrkphperpl;
01739 ann_aTrkphper = anninput.aTrkphper;
01740 ann_aTrkplu = anninput.aTrkplu;
01741 ann_aTrkplv = anninput.aTrkplv;
01742 ann_aTrkstp = anninput.aTrkstp;
01743
01744 ann_aShwph = anninput.aShwph;
01745 ann_aShwstp = anninput.aShwstp;
01746 ann_aShwdig = anninput.aShwdig;
01747 ann_aShwpl = anninput.aShwpl;
01748 ann_aShwphper = anninput.aShwphper;
01749 ann_aShwphperpl = anninput.aShwphperpl;
01750 ann_aShwphperdig = anninput.aShwphperdig ;
01751 ann_aShwphperstp = anninput.aShwphperstp;
01752 ann_aShwplu = anninput.aShwplu;
01753 ann_aShwplv = anninput.aShwplv;
01754 ann_aPh3 = anninput.aPh3;
01755 ann_aPh6 = anninput.aPh6;
01756 ann_aPhlast = anninput.aPhlast;
01757 ann_aPhcommon = anninput.aPhcommon;
01758
01759
01760
01761
01762
01763
01764 evttimemin = anninput.aTimemin-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9;
01765 evttimemax = anninput.aTimemax-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9;
01766
01767
01768
01769
01770
01771 Int_t target=1;
01772 if(ntpHeader->GetVldContext().GetSimFlag()!=4){
01773
01774
01775
01776
01777
01778 if(mcneartype==1) target=1;
01779 if(mcneartype==2) target=2;
01780 if(mcneartype==3) target=3;
01781 }
01782 if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==1){
01783 if(mcneartype==1) target=1;
01784 if(mcneartype==2) target=2;
01785 if(mcneartype==3) target=3;
01786 }
01787 if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==2){
01788 target=0;
01789 }
01790
01791
01792 Int_t fid =1;
01793 Int_t prior = 1;
01794
01795 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01796 if(!PassFid(event)) annpid=-999;
01797 else {
01798 annpid=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
01799 first_ann=false;
01800 }
01801 }
01802 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01803 if(!is_fid_ns){
01804 annpid_f2p1=-999;
01805 annpid_f2p2=-999;
01806 annpid_f2p3=-999;
01807 }
01808 if(is_fid_ns){
01809 fid=2;
01810 prior=1;
01811 annpid_f2p1=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
01812 prior=2;
01813 annpid_f2p2=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
01814 prior=3;
01815 annpid_f2p3=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
01816 first_ann=false;
01817 }
01818 if(!is_fid_dp){
01819 annpid_f1p1=-999;
01820 annpid_f1p2=-999;
01821 annpid_f1p3=-999;
01822 }
01823 if(is_fid_dp){
01824 fid=1;
01825 prior=1;
01826 annpid_f1p1=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
01827 prior=2;
01828 annpid_f1p2=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
01829 prior=3;
01830 annpid_f1p3=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
01831 first_ann=false;
01832 }
01833 }
01834
01835
01836
01837 bool isdata=true;
01838 if(ntpHeader->GetVldContext().GetSimFlag()==4) isdata=false;
01839
01840
01841 int track_index = -1;
01842 LoadLargestTrackFromEvent(j,track_index);
01843 int shower_index = -1;
01844
01845 LoadShower_Jim(j,shower_index,detector);
01846 if(shower_index==-1) nshower=0;
01847 if(shower_index!=-1){
01848
01849 shwdigits = ntpShower->ndigit;
01850 shwstrips = ntpShower->nstrip;
01851 shwph = ntpShower->ph.sigcor;
01852 shwvtxu = ntpShower->vtx.u;
01853 shwvtxv = ntpShower->vtx.v;
01854 shwvtxx = ntpShower->vtx.x;
01855 shwvtxy = ntpShower->vtx.y;
01856 shwvtxz = ntpShower->vtx.z;
01857 shwncluster = ntpShower->ncluster;
01858
01859
01860 Int_t shwtype=0;
01861 const Int_t shw_correction_mode=1;
01862 shwtype=-1;
01863 reco_eshw = RecoShwEnergy(shower_index,shwtype,
01864 detector,isdata,
01865 shw_correction_mode);
01866 shwtype=0;
01867 reco_eshwcc = RecoShwEnergy(shower_index,shwtype,
01868 detector,isdata,
01869 shw_correction_mode);
01870 shwtype=2;
01871 reco_eshwnc = RecoShwEnergy(shower_index,shwtype,
01872 detector,isdata,
01873 shw_correction_mode);
01874 shwtype=1;
01875 reco_eshwccw = RecoShwEnergy(shower_index,shwtype,
01876 detector,isdata,
01877 shw_correction_mode);
01878 shwtype=3;
01879 reco_eshwncw = RecoShwEnergy(shower_index,shwtype,
01880 detector,isdata,
01881 shw_correction_mode);
01882
01883 reco_eshw_sqrt = RecoShwEnergySqrt(event);
01884
01885 }
01886
01887
01888 if(track_index!=-1){
01889
01890 reco_emu = RecoMuEnergy(0,track_index,isdata,
01891 ntpHeader->GetVldContext().GetDetector());
01892 trklength = ntpTrack->plane.end-ntpTrack->plane.beg+1;
01893 trkmomrange = ntpTrack->momentum.range;
01894 trkqp = ntpTrack->momentum.qp;
01895 trkeqp = ntpTrack->momentum.eqp;
01896 trkendz = ntpTrack->end.z;
01897 trkendu = ntpTrack->end.u;
01898 trkendv = ntpTrack->end.v;
01899 trkendx = ntpTrack->end.x;
01900 trkendy = ntpTrack->end.y;
01901
01902 trkvtxz = ntpTrack->vtx.z;
01903 trkvtxu = ntpTrack->vtx.u;
01904 trkvtxv = ntpTrack->vtx.v;
01905 trkvtxx = ntpTrack->vtx.x;
01906 trkvtxy = ntpTrack->vtx.y;
01907 trkds = ntpTrack->ds;
01908 trkplbu = ntpTrack->plane.begu;
01909 trkplbv = ntpTrack->plane.begv;
01910 trkpleu = ntpTrack->plane.endu;
01911 trkplev = ntpTrack->plane.endv;
01912 trkple = ntpTrack->plane.end;
01913 trkplb = ntpTrack->plane.beg;
01914 trkfitpass = ntpTrack->fit.pass;
01915 trkdiffuv = abs(ntpTrack->plane.nu-ntpTrack->plane.nv);
01916 trkdigits = ntpTrack->ndigit;
01917 trkstrips = ntpTrack->nstrip;
01918 trkph = ntpTrack->ph.sigcor;
01919 trkndof = ntpTrack->fit.ndof;
01920 trkchi2 = ntpTrack->fit.chi2;
01921
01922 Float_t phtrack =ntpTrack->ph.sigcor;
01923 Float_t phevent =ntpEvent->ph.sigcor;
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934 if (phevent>0) {trkphfrac=phtrack/phevent;}
01935 if(ntpTrack->plane.n>0) {
01936 trkphplane = ntpTrack->ph.sigcor/ntpTrack->plane.n;
01937 }
01938
01939 is_cont_trk = ntpTrack->contained;
01940 trkmombest =ntpTrack->momentum.best;
01941 }
01942 else {
01943 trklength = 0;
01944 trkmomrange = 0;
01945 trkqp = 0;
01946 trkeqp = 0;
01947 trkphfrac = 0;
01948 trkphplane = 0;
01949 }
01950
01951 reco_enu = reco_emu + reco_eshwcc;
01952 reco_qe_enu = RecoQENuEnergy(track_index);
01953 if (reco_enu>0) {
01954 reco_y = reco_eshwcc/reco_enu;
01955 }
01956 reco_qe_q2 = RecoQEQ2(track_index);
01957 reco_dircosz = RecoMuDCosZ(track_index);
01958 reco_dircosneu = RecoMuDCosNeu(track_index);
01959 is_cev = IsFidAll(track_index);
01960
01961
01962
01963
01964 if(is_fid_dp){
01965 pid0 = PID(event,0);
01966 pid1 = PID(event,1);
01967 pid2 = PID(event,2);
01968 if(PassAnalysisCuts(event)) pass = 1;
01969 }
01970 else{
01971 pid0 = -999.;
01972 pid1 = -999.;
01973 pid2 = -999.;
01974 pass = 0;
01975 }
01976 tree->Fill();
01977
01978 }
01979 }
01980 }
01981
01982
01983 delete spinfo;
01984
01985 save.cd();
01986 tree->Write();
01987 save.Write();
01988 save.Close();
01989 }
01990
01991
01992
01993 void MadDpAnalysis::ReadPIDFile(std::string tag){
01994
01995 fLikeliFile = new TFile(tag.c_str(),"READ");
01996
01997 if(fLikeliFile->IsOpen() && !fLikeliFile->IsZombie()) {
01998
01999 fLikeliFile->cd();
02000
02001 fLikeliHist[0] = (TH1F*) fLikeliFile->Get("Event length - cc");
02002 fLikeliHist[1] = (TH1F*) fLikeliFile->Get("Event length - nc");
02003 fLikeliHist[2] = (TH1F*) fLikeliFile->Get("Track ph frac - cc");
02004 fLikeliHist[3] = (TH1F*) fLikeliFile->Get("Track ph frac - nc");
02005 fLikeliHist[4] = (TH1F*) fLikeliFile->Get("ph per plane (cc)");
02006 fLikeliHist[5] = (TH1F*) fLikeliFile->Get("ph per plane (nc)");
02007
02008 for(int i=0;i<6;i++){
02009 float integ = fLikeliHist[i]->Integral();
02010 fLikeliHist[i]->Scale(1./integ);
02011 }
02012 }
02013 else fLikeliFile = NULL;
02014 }
02015
02016
02017 bool MadDpAnalysis::InFidVol(const Detector::Detector_t& det,
02018 const Float_t& x, const Float_t& y,
02019 const Float_t& z){
02020
02021
02022
02023
02024 bool is_fid=false;
02025
02026 if(det==Detector::kFar){
02027
02028
02029 is_fid = true;
02030 if(z<0.5 || z>29.4 ||
02031 (z<16.5&&z>14.5) ||
02032 sqrt((x*x)
02033 +(y*y))>3.5) is_fid = false;
02034 }
02035 else if(det==Detector::kNear){
02036
02037
02038 is_fid = true;
02039 if(z<1 || z>5 ||
02040 sqrt(((x-1.4885)*(x-1.4885)) +
02041 ((y-0.1397)*(y-0.1397)))>1) is_fid = false;
02042 }
02043
02044 return is_fid;
02045
02046 }
02047
02048 #endif // #ifdef maddpanalysis_cxx