00001 #ifndef madpidanalysis_cxx
00002 #define madpidanalysis_cxx
00003
00004 #include <iostream>
00005 #include <cmath>
00006 #include <map>
00007 #include "Mad/MadPIDAnalysis.h"
00008 #include "TClonesArray.h"
00009 #include "TMath.h"
00010 #include "Validity/VldContext.h"
00011 #include "Conventions/Detector.h"
00012
00013 #include "Mad/MadAnalysis.h"
00014 #include "MCReweight/MCReweight.h"
00015 #include "MCReweight/GnumiInterface.h"
00016 #include "MCReweight/NuParent.h"
00017 #include "Mad/SpillInfo.h"
00018 #include "Mad/SpillInfoBlock.h"
00019 #include "BeamDataUtil/BDSpillAccessor.h"
00020 #include "BeamDataUtil/BeamMonSpill.h"
00021 #include "SpillTiming/SpillTimeFinder.h"
00022 #include "Riostream.h"
00023 #include "Mad/ScanList.h"
00024 #include "Mad/fddataquality.h"
00025 #include "DcsUser/Dcs_Mag_Near.h"
00026
00027 #include "BeamDataUtil/BMSpillAna.h"
00028
00029 #include "CandNtupleSR/NtpSRDataQuality.h"
00030 #include "CandNtupleSR/NtpSRDetStatus.h"
00031
00032 #include "DcsUser/HvStatusFinder.h"
00033 #include "DcsUser/CoilTools.h"
00034 #include "SpillTiming/SpillServerMonFinder.h"
00035
00036 #include "DataUtil/DataQualDB.h"
00037
00038 #include "NtupleUtils/LISieve.h"
00039 #include "Conventions/ReleaseType.h"
00040
00041
00042 using namespace std;
00043
00044
00045
00046
00047
00048
00049
00050 MadPIDAnalysis::MadPIDAnalysis(TChain *chainSR,TChain *chainMC,
00051 TChain *chainTH,TChain *chainEM)
00052 {
00053
00054 if(!chainSR) {
00055 record = 0;
00056 strecord = 0;
00057 emrecord = 0;
00058 mcrecord = 0;
00059 threcord = 0;
00060 Clear();
00061 whichSource = -1;
00062 isMC = true;
00063 isTH = true;
00064 isEM = true;
00065 Nentries = -1;
00066 return;
00067 }
00068
00069 InitChain(chainSR,chainMC,chainTH,chainEM);
00070 whichSource = 0;
00071 fLikeliFile = NULL;
00072 for(int i=0;i<6;i++) fLikeliHist[i] = NULL;
00073
00074 potall=0.;
00075 potcut=0.;
00076
00077 }
00078
00079
00080 MadPIDAnalysis::MadPIDAnalysis(JobC *j,string path,int entries)
00081 {
00082 record = 0;
00083 strecord = 0;
00084 emrecord = 0;
00085 mcrecord = 0;
00086 threcord = 0;
00087 Clear();
00088 isMC = true;
00089 isTH = true;
00090 isEM = true;
00091 Nentries = entries;
00092 whichSource = 1;
00093 jcPath = path;
00094 fJC = j;
00095 fLikeliFile = NULL;
00096 for(int i=0;i<6;i++) fLikeliHist[i] = NULL;
00097 }
00098
00099
00100 MadPIDAnalysis::~MadPIDAnalysis()
00101 {
00102
00103
00104 if(fLikeliFile) fLikeliFile->Close();
00105
00106 cout << "pot (all) =" << potall << endl;
00107 cout << "pot (beamcuts)=" << potcut << endl;
00108
00109 }
00110
00111
00112 Bool_t MadPIDAnalysis::PassTruthSignal(Int_t mcevent)
00113 {
00114 if(!LoadTruth(mcevent)) return false;
00115 if(abs(ntpTruth->inu)==14&&ntpTruth->iaction==1) return true;
00116 return false;
00117 }
00118
00119
00120 Bool_t MadPIDAnalysis::PassAnalysisCuts(Int_t event)
00121 {
00122 if(!LoadEvent(event)) return false;
00123
00124 Float_t pidcut=-0.4;
00125 if(ntpHeader->GetVldContext().GetDetector()==1) {pidcut=-0.1;}
00126
00127 if(!(PID(event,0)>pidcut)) return false;
00128 return true;
00129 }
00130
00131
00132 Bool_t MadPIDAnalysis::PassFidNew(Int_t event)
00133 {
00134
00135 if(!LoadEvent(event)) return false;
00136
00137 Int_t track = -1;
00138 LoadLargestTrackFromEvent(event,track);
00139
00140
00141
00142 if (track==-1 && (ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>28.0 ||
00143 (ntpEvent->vtx.z<16.2 && ntpEvent->vtx.z>14.3) ||
00144 (pow(ntpEvent->vtx.x,2)+
00145 pow(ntpEvent->vtx.y,2))>14)) return false;
00146
00147
00148
00149 if (track!=-1 && (ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>28.0 ||
00150 (ntpTrack->vtx.z<16.2 && ntpTrack->vtx.z>14.3) ||
00151 (pow(ntpTrack->vtx.x,2)+
00152 pow(ntpTrack->vtx.y,2))>14)) return false;
00153
00154 return true;
00155
00156 }
00157
00158
00159
00160 Bool_t MadPIDAnalysis::PassFidNewN(Int_t event)
00161 {
00162
00163
00164 if(!LoadEvent(event)) return false;
00165
00166 Int_t track = -1;
00167 LoadLargestTrackFromEvent(event,track);
00168
00169
00170
00171 if (track==-1 && (ntpEvent->vtx.z<1.0 || ntpEvent->vtx.z>27.0 ||
00172 (ntpEvent->vtx.z<17.0 && ntpEvent->vtx.z>14.0) ||
00173 (pow(ntpEvent->vtx.x,2)+
00174 pow(ntpEvent->vtx.y,2))>(3.5*3.5))) return false;
00175
00176
00177
00178 if (track!=-1 && (ntpTrack->vtx.z<1.0 || ntpTrack->vtx.z>27.0 ||
00179 (ntpTrack->vtx.z<17. && ntpTrack->vtx.z>14.0) ||
00180 (pow(ntpTrack->vtx.x,2)+
00181 pow(ntpTrack->vtx.y,2))>(3.5*3.5))) return false;
00182
00183 return true;
00184
00185 }
00186
00187
00188 Bool_t MadPIDAnalysis::PassBasicCuts(Int_t event)
00189 {
00190 if(!LoadEvent(event)) return false;
00191 if(ntpEvent->ntrack==0) return false;
00192
00193
00194
00195
00196
00197 return true;
00198 }
00199 Bool_t MadPIDAnalysis::PassFid(Int_t event)
00200 {
00201 if(!LoadEvent(event)) return false;
00202 if(!IsFidVtxEvt(ntpEvent->index)) return false;
00203 return true;
00204 }
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214 Bool_t MadPIDAnalysis::PassTimingCuts(Double_t timediff)
00215 {
00216
00217
00218 if (timediff<-20e-6) return false;
00219 if (timediff>30e-6) return false;
00220
00221
00222 return true;
00223
00224 }
00225
00226
00227
00228
00229 Bool_t MadPIDAnalysis::PassFidND(Int_t event)
00230 {
00231 if(!LoadEvent(event)) return false;
00232
00233
00234
00235 if(ntpEvent->ntrack==0) return false;
00236
00237 Int_t track = -1;
00238 LoadLargestTrackFromEvent(event,track);
00239
00240
00241
00242 if (track!=-1 && (ntpTrack->vtx.z<1.0 || ntpTrack->vtx.z>5.0 ||
00243
00244 (pow(ntpTrack->vtx.x-1.4828,2)+
00245 pow(ntpTrack->vtx.y-0.2384,2))>(1.0))) return false;
00246
00247
00248
00249
00250 return true;
00251
00252 }
00253
00254
00255
00256
00257 AnnInputBlock MadPIDAnalysis::AnnVar(Int_t event)
00258 {
00259
00260 AnnInputBlock anninput;
00261
00262
00263
00264 anninput.aTotrk =0.;
00265 anninput.aTotstp =0.;
00266 anninput.aTotph =0.;
00267 anninput.aTotlen =0.;
00268 anninput.aPhperpl =0.;
00269 anninput.aPhperstp =0.;
00270
00271
00272 anninput.aTrkpass =0.;
00273 anninput.aTrkph =0.;
00274 anninput.aTrklen =0.;
00275 anninput.aTrkphperpl =0.;
00276 anninput.aTrkphper =0.;
00277 anninput.aTrkplu =0.;
00278 anninput.aTrkplv =0.;
00279 anninput.aTrkstp =0.;
00280
00281 anninput.aShwph =0.;
00282 anninput.aShwstp =0.;
00283 anninput.aShwdig =0.;
00284 anninput.aShwpl =0.;
00285 anninput.aShwphper =0.;
00286 anninput.aShwphperpl =0.;
00287 anninput.aShwphperdig =0.;
00288 anninput.aShwphperstp =0.;
00289 anninput.aShwplu =0.;
00290 anninput.aShwplv =0.;
00291 anninput.aShwstpu =0.;
00292 anninput.aShwstpv =0.;
00293 anninput.aPh3 =0.;
00294 anninput.aPh6 =0.;
00295 anninput.aPhlast =0.;
00296 anninput.aPhcommon =0.;
00297 anninput.aTimemax =0.;
00298 anninput.aTimemin =0.;
00299
00300
00301
00302
00303 Int_t detector = 0;
00304
00305 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) detector=1;
00306 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar) detector=2;
00307
00308
00309 LoadEvent(event) ;
00310 int track_index = -1;
00311 LoadLargestTrackFromEvent(event,track_index);
00312 int shower_index = -1;
00313
00314 LoadShower_Jim(event,shower_index,detector);
00315
00316 anninput.aTotrk =ntpEvent->ntrack;
00317 anninput.aTotstp =ntpEvent->nstrip;
00318 anninput.aTotph =ntpEvent->ph.sigcor;
00319 anninput.aTotlen =ntpEvent->plane.end-ntpEvent->plane.beg+1;
00320
00321 if (ntpEvent->plane.n>0) anninput.aPhperpl =ntpEvent->ph.sigcor/ntpEvent->plane.n;
00322 if (ntpEvent->nstrip>0) anninput.aPhperstp =ntpEvent->ph.sigcor/ntpEvent->nstrip;
00323
00324
00325 if(track_index!=-1) {
00326 anninput.aTrkpass =ntpTrack->fit.pass;
00327 anninput.aTrkph =ntpTrack->ph.sigcor;
00328 anninput.aTrklen =ntpTrack->plane.ntrklike;
00329
00330 if(ntpTrack->plane.ntrklike>0) anninput.aTrkphperpl =ntpTrack->ph.sigcor/ntpTrack->plane.ntrklike;
00331 if(ntpEvent->ph.sigcor>0) anninput.aTrkphper =ntpTrack->ph.sigcor/ntpEvent->ph.sigcor;
00332
00333 anninput.aTrkplu =ntpTrack->plane.nu;
00334 anninput.aTrkplv =ntpTrack->plane.nv;
00335 anninput.aTrkstp =ntpTrack->nstrip;
00336 }
00337
00338 if(shower_index!=-1) {
00339 anninput.aShwph =ntpShower->ph.sigcor;
00340 anninput.aShwstp =ntpShower->nstrip;
00341 anninput.aShwdig =ntpShower->ndigit;
00342 anninput.aShwpl =ntpShower->plane.end-ntpShower->plane.beg+1;
00343
00344 if (ntpEvent->ph.sigcor>0) anninput.aShwphper =ntpShower->ph.sigcor/ntpEvent->ph.sigcor;
00345 anninput.aShwphperpl =ntpShower->ph.sigcor/(ntpShower->plane.end-ntpShower->plane.beg+1);
00346
00347 if (ntpShower->ndigit>0) anninput.aShwphperdig =ntpShower->ph.sigcor/ntpShower->ndigit;
00348 if (ntpShower->nstrip>0) anninput.aShwphperstp =ntpShower->ph.sigcor/ntpShower->nstrip;
00349
00350 anninput.aShwplu =ntpShower->plane.nu;
00351 anninput.aShwplv =ntpShower->plane.nv;
00352 }
00353
00354
00355 Int_t ssind,ssplind;
00356 Double_t ssphtot;
00357 Bool_t foundshw,foundtrk;
00358 Int_t planes=ntpEvent->plane.beg;
00359
00360 Double_t ph3,ph6,phlast,phcommon,phnowere,hitcommon,hitnowere;
00361 ph3=0.;ph6=0.;phlast=0.;phcommon=0.;phnowere=0.;hitcommon=0.;hitnowere=0.;
00362
00363
00364
00365 Double_t timemin=9.*10e30;
00366 Double_t timemax=-9999999;
00367
00368 Float_t shwstpu=0;
00369 Float_t shwstpv=0;
00370
00371 if(shower_index!=-1) {
00372 for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
00373 Int_t stp_indexs=((ntpShower->stp)[jj]);
00374 if(stp_indexs!=-1){
00375 LoadStrip(stp_indexs);
00376 if((ntpStrip->plane)%2==1) shwstpu=shwstpu+1;
00377 if((ntpStrip->plane)%2==0) shwstpv=shwstpv+1;
00378 }
00379 }
00380 }
00381
00382
00383
00384 for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
00385 Int_t stp_index=((ntpEvent->stp)[evss]);
00386 if(stp_index!=-1){
00387 LoadStrip(stp_index);
00388
00389
00390 if(detector==1){
00391 Double_t striptime=ntpStrip->time1;
00392 if(striptime<=timemin) timemin=striptime;
00393 if(striptime>=timemax) timemax=striptime;
00394 }
00395
00396 if(detector==2){
00397 Double_t striptime1=ntpStrip->time1;
00398 Double_t striptime0=ntpStrip->time0;
00399 Double_t striptime=0;
00400 if(striptime1>0 && striptime0<0) striptime=striptime1;
00401 if(striptime0>0 && striptime1<0) striptime=striptime0;
00402 if(striptime0>0 && striptime1>0) striptime=(striptime0+striptime1)/2.;
00403 if(striptime<=timemin) timemin=striptime;
00404 if(striptime>=timemax) timemax=striptime;
00405
00406 }
00407
00408 ssind=ntpStrip->strip;
00409 ssplind=ntpStrip->plane;
00410 ssphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
00411
00412 foundshw=false;
00413 foundtrk=false;
00414
00415 Double_t phstrips=0;
00416 Double_t phstript=0;
00417
00418
00419 Int_t sshwind,sshwplind;
00420 Double_t sshwphtot;
00421
00422 if(shower_index!=-1) {
00423 for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
00424 Int_t stp_indexs=((ntpShower->stp)[jj]);
00425 if(stp_indexs!=-1){
00426 LoadStrip(stp_indexs);
00427
00428 if(!foundshw){
00429 sshwind=ntpStrip->strip;
00430 sshwplind=ntpStrip->plane;
00431 sshwphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
00432 if(sshwind==ssind && sshwplind==ssplind) {
00433 foundshw=true;
00434 phstrips=sshwphtot;
00435 }
00436 }
00437 }
00438 }
00439 }
00440
00441
00442 Int_t strkind,strkplind;
00443 Double_t strkphtot;
00444
00445 if(track_index!=-1) {
00446 for(Int_t jj=0;jj<ntpTrack->nstrip;jj++){
00447 Int_t stp_indext=((ntpTrack->stp)[jj]);
00448
00449 if(stp_indext!=-1){
00450 LoadStrip(stp_indext);
00451 if(!foundtrk){
00452 strkind=ntpStrip->strip;
00453 strkplind=ntpStrip->plane;
00454 strkphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
00455 if(strkind==ssind && strkplind==ssplind) {
00456 foundtrk=true;
00457 phstript=strkphtot;
00458 }
00459 }
00460 }
00461 }
00462 }
00463
00464 if(foundshw && foundtrk) {
00465 hitcommon=hitcommon+1;
00466 phcommon=phcommon+phstrips+phstript;
00467 }
00468
00469 if(!foundshw && ! foundtrk) {
00470 hitnowere=hitnowere+1;
00471 phnowere=phnowere+ssphtot;
00472 }
00473
00474 if(ssplind>=planes && ssplind<=(planes+3)){
00475 ph3=ph3+ssphtot;
00476 }
00477 else if(ssplind>(planes+3) && ssplind<=(planes+6)){
00478 ph6=ph6+ssphtot;
00479 }
00480 else {
00481 phlast=phlast+ssphtot;
00482 }
00483
00484 }
00485 }
00486
00487 anninput.aTimemin = timemin;
00488 anninput.aTimemax = timemax;
00489
00490
00491
00492 anninput.aPh3 =ph3;
00493 anninput.aPh6 =ph6;
00494 anninput.aPhlast =phlast;
00495 anninput.aPhcommon =phcommon;
00496 anninput.aShwstpu =shwstpu;
00497 anninput.aShwstpv =shwstpv;
00498
00499
00500 return anninput;
00501
00502 }
00503
00504
00505
00506 Double_t MadPIDAnalysis::Sigmoid(Double_t x)
00507 {
00508 double sig;
00509 if(x>37.){
00510 sig = 1.;
00511 }
00512 else if(x<-37.){
00513 sig = 0.;
00514 }
00515 else {
00516 sig = 1./(1.+exp(-x));
00517 }
00518 return sig;
00519 }
00520
00521
00522
00523 Double_t MadPIDAnalysis::GetAnnPid(AnnInputBlock anninput,Int_t det, Int_t tar,Int_t fid, Int_t prior, Bool_t first_ann)
00524 {
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535 int inneuron = 7;
00536 int hidneuron = 15;
00537 double var;
00538
00539 double out[15];
00540 double rin[15];
00541
00542
00543 double weight1[7][15];
00544 double constant1[15];
00545
00546 double weighto[15];
00547 double constanto[1];
00548 double prob;
00549
00550
00551
00552 for(Int_t i=0;i<hidneuron;i++) {
00553 out[i]=0.;
00554 rin[i]=0.;
00555 constant1[i]=0.;
00556 }
00557
00558 for(Int_t i=0;i<inneuron;i++) {
00559 for(Int_t j=0;j<hidneuron;j++){
00560 weight1[i][j]=0.;
00561 }
00562 }
00563
00564 constanto[0]=0.;
00565 for(Int_t i=0;i<hidneuron;i++){
00566 weighto[i]=0.;
00567 }
00568 Int_t all = inneuron*hidneuron+hidneuron+hidneuron+1;
00569 Int_t all1 = inneuron*hidneuron+hidneuron;
00570 Int_t n1 =0;
00571 Int_t nw1 =0;
00572 Int_t no =0;
00573 Int_t nwo =0;
00574 ifstream weightfile;
00575
00576
00577 if(first_ann){
00578
00579 if(det==2){
00580 if(prior==1){
00581 if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00582 if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00583 }
00584 if(prior==2){
00585 if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00586 if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00587 }
00588 if(prior==3){
00589 if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00590 if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat");
00591 }
00592 }
00593
00594 if(det==1){
00595 if(tar==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearle_cedar_daikon7v2.dat");
00596 if(tar==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearpme_cedar_daikon7v2.dat");
00597 if(tar==3) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearphe_cedar_daikon7v2.dat");
00598 }
00599 for(Int_t k=0;k<all;k++){
00600 weightfile >> var ;
00601 Ann_weight_vector[k]=var;
00602 }
00603 weightfile.close();
00604 }
00605 if(det==2){
00606
00607 out[0] = (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/2000.;
00608 out[1] = (anninput.aPhperpl/10000.);
00609
00610 out[2] = (anninput.aTotlen/40);
00611 out[3] = (anninput.aPh3/(anninput.aTotph+1.));
00612 out[4] = (anninput.aPh6/(anninput.aTotph+1.));
00613 out[5] = (anninput.aPhlast/(anninput.aTotph+1.));
00614 out[6] = ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.;
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626 }
00627
00628
00629 if(det==1){
00630
00631 out[0] = (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/5000.;
00632 out[1] = (anninput.aPhperpl/10000.);
00633
00634 out[2] = (anninput.aTotlen/40);
00635 out[3] = (anninput.aPh3/(anninput.aTotph+1.));
00636 out[4] = (anninput.aPh6/(anninput.aTotph+1.));
00637 out[5] = (anninput.aPhlast/(anninput.aTotph+1.));
00638 out[6] = ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.;
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650 }
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660 for(Int_t k=0;k<all;k++){
00661 if(k<all1){
00662 if(k%(inneuron+1)==0){
00663 nw1=0;
00664 constant1[n1]=Ann_weight_vector[k];
00665 n1=n1+1;
00666 }
00667 else {
00668 weight1[nw1][n1-1]=Ann_weight_vector[k];
00669 nw1=nw1+1;
00670 }
00671 }
00672 else if(k>=all1){
00673 if((k-all1)%(hidneuron+1)==0){
00674 nwo=0;
00675 constanto[no]=Ann_weight_vector[k];
00676 no=no+1;
00677 }
00678 else {
00679 weighto[nwo]=Ann_weight_vector[k];
00680 nwo=nwo+1;
00681 }
00682 }
00683 }
00684
00685
00686
00687
00688
00689 for(Int_t i=0;i<hidneuron;i++){
00690 rin[i]=constant1[i];
00691 for(Int_t j=0;j<inneuron;j++){
00692 rin[i]=rin[i]+weight1[j][i]*out[j];
00693 }
00694 }
00695
00696 for(Int_t i=0;i<hidneuron;i++){
00697 out[i]=Sigmoid(rin[i]);
00698 }
00699
00700 prob=constanto[0];
00701 for(Int_t i=0;i<hidneuron;i++) prob=prob+weighto[i]*out[i];
00702
00703
00704 if(anninput.aTotlen>50) prob=0.;
00705
00706 return prob;
00707 }
00708
00709
00710
00711 Float_t MadPIDAnalysis::PID(Int_t event, Int_t method)
00712 {
00713
00714
00715 if(!LoadEvent(event)) return -10;
00716 Int_t track = -1;
00717 if(!LoadLargestTrackFromEvent(event,track)) return -10;
00718 if(method!=0) return -10;
00719
00720 Float_t ProbNC = 1.;
00721 Float_t ProbCC = 1.;
00722 Float_t PidCC = 0.;
00723
00724
00725 Float_t var1=eventSummary->plane.n;
00726 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00727 var1=ntpEvent->plane.end-ntpEvent->plane.beg+1;
00728 }
00729 Float_t var2=0;
00730 Float_t var3=0;
00731
00732
00733 Float_t phtrack=ntpTrack->ph.sigcor;
00734 Float_t phevent=ntpEvent->ph.sigcor;
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745 if (phevent>0) {var2=phtrack/phevent;}
00746
00747 if (ntpTrack->plane.n!=0) var3 = float(phtrack)
00748 /float(ntpTrack->plane.n);
00749
00750
00751 Float_t prob1=fLikeliHist[0]->GetBinContent(fLikeliHist[0]->FindBin(var1));
00752 Float_t prob2=fLikeliHist[2]->GetBinContent(fLikeliHist[2]->FindBin(var2));
00753 Float_t prob3=fLikeliHist[4]->GetBinContent(fLikeliHist[4]->FindBin(var3));
00754
00755 if (prob1==0) {prob1=0.0001;}
00756 if (prob2==0) {prob2=0.0001;}
00757 if (prob3==0) {prob3=0.0001;}
00758
00759 ProbCC=prob1*prob2*prob3;
00760
00761 prob1=fLikeliHist[1]->GetBinContent(fLikeliHist[1]->FindBin(var1));
00762 prob2=fLikeliHist[3]->GetBinContent(fLikeliHist[3]->FindBin(var2));
00763 prob3=fLikeliHist[5]->GetBinContent(fLikeliHist[5]->FindBin(var3));
00764
00765 if (prob1==0) {prob1=0.0001;}
00766 if (prob2==0) {prob2=0.0001;}
00767 if (prob3==0) {prob3=0.0001;}
00768
00769 ProbNC=prob1*prob2*prob3;
00770
00771 PidCC=-(sqrt(-TMath::Log(ProbCC))-sqrt(-TMath::Log(ProbNC)));
00772
00773 return PidCC;
00774 }
00775
00776
00777 Float_t MadPIDAnalysis::RecoAnalysisEnergy(Int_t event){
00778
00779 if(!LoadEvent(event)) return 0;
00780 int largestTrack = -1;
00781 LoadLargestTrackFromEvent(event,largestTrack);
00782 float emu = RecoMuEnergy(0,largestTrack);
00783 float ehad = RecoShwEnergy(event);
00784 float ereco = emu + ehad;
00785 return ereco;
00786
00787 }
00788
00789 void MadPIDAnalysis::MakeMyFile(std::string tag,int tarpos, int pdfbinning){
00790
00791
00792
00793
00794
00795 std::string savename = "DPHistos_" + tag + ".root";
00796 TFile save(savename.c_str(),"RECREATE");
00797 save.cd();
00798
00799
00800
00801
00802
00803 Int_t evlbins=0;
00804
00805 if (pdfbinning==0) evlbins=60;
00806 else evlbins=300;
00807
00808
00809 TH1F *evlength_nc = new TH1F("Event length - nc",
00810 "Event length - nc",
00811 evlbins,0,600);
00812 evlength_nc->SetXTitle("Event length (planes)");
00813 TH1F *evlength_cc = new TH1F("Event length - cc",
00814 "Event length - cc",
00815 evlbins,0,600);
00816 evlength_cc->SetXTitle("Event length (planes)");
00817
00818
00819
00820 TH1F *phfrac_nc = new TH1F("Track ph frac - nc",
00821 "Track ph frac - nc",
00822 50,0,1);
00823 phfrac_nc->SetXTitle("pulse height fraction");
00824
00825
00826 TH1F *phfrac_cc = new TH1F("Track ph frac - cc",
00827 "Track ph frac - cc",
00828 50,0,1);
00829 phfrac_cc->SetXTitle("pulse height fraction");
00830
00831
00832
00833 Int_t phplanebins=0;
00834 Float_t phplanemax=0;
00835
00836 if (pdfbinning==0) {
00837 phplanebins=50;
00838 phplanemax=5000;
00839 }
00840 else {
00841 phplanebins=100;
00842 phplanemax=3000;
00843 }
00844
00845 TH1F *phplane_nc = new TH1F("ph per plane (nc)",
00846 "ph per plane (nc)",
00847 phplanebins,0,phplanemax);
00848 phplane_nc->SetXTitle("pulse height per plane");
00849 TH1F *phplane_cc = new TH1F("ph per plane (cc)",
00850 "ph per plane (cc)",
00851 phplanebins,0,phplanemax);
00852 phplane_cc->SetXTitle("pulse height per plane");
00853
00854
00855 evlength_nc->SetDirectory(&save);
00856 evlength_cc->SetDirectory(&save);
00857
00858 phfrac_nc->SetDirectory(&save);
00859 phfrac_cc->SetDirectory(&save);
00860
00861 phplane_nc->SetDirectory(&save);
00862 phplane_cc->SetDirectory(&save);
00863
00864
00865
00866
00867
00868
00869
00870 for(int i=0;i<Nentries;i++){
00871
00872 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
00873 << "% done" << std::endl;
00874 if(!this->GetEntry(i)) continue;
00875
00876 Int_t nevent = eventSummary->nevent;
00877 if(nevent==0) continue;
00878
00879 Int_t trkcount=0;
00880 Int_t shwcount=0;
00881
00882 Int_t evtmin=0;
00883 Int_t evtmax=nevent;
00884
00885 if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) {
00886
00887 Int_t maxevent=0;
00888 Float_t maxph=0;
00889 for(int j=0;j<nevent;j++){
00890 if(!LoadEvent(j)) continue;
00891 Float_t evtpheight=ntpEvent->ph.sigcor;
00892 if (evtpheight>maxph) {
00893 maxph=evtpheight;
00894 maxevent=j;
00895 }
00896 }
00897 evtmin=maxevent;
00898 evtmax=maxevent+1;
00899 }
00900
00901
00902
00903 for(int j=evtmin;j<evtmax;j++){
00904
00905
00906 if(!LoadEvent(j)) continue;
00907
00908 Int_t ntrack = 0;
00909 ntrack=ntpEvent->ntrack;
00910 Int_t nshower = 0;
00911 nshower=ntpEvent->nshower;
00912
00913 Float_t totph=0;
00914
00915
00916 for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
00917 LoadTrack(ii);
00918
00919 totph+=ntpTrack->ph.sigcor;
00920 LoadTHTrack(ii);
00921
00922 }
00923
00924 trkcount+=ntrack;
00925
00926
00927 for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
00928 LoadShower(ii);
00929
00930 totph+=ntpShower->ph.sigcor;
00931 }
00932
00933 shwcount+=nshower;
00934
00935
00936
00937 Int_t cc_nc = 0;
00938 Int_t flavour = 0;
00939 Int_t detector = -1;
00940 Int_t is_fid = 0;
00941 Double_t neu_e = -1;
00942 Double_t weight = -1;
00943
00944
00945
00946 is_fid = IsFid_2008(ntpEvent->index);
00947
00948
00949
00950
00951 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00952 detector = 2;
00953
00954 if (evlength_cc->GetNbinsX()==evlbins) {
00955 if (pdfbinning==0) {
00956 evlength_cc->SetBins(50,0,500);
00957 evlength_nc->SetBins(50,0,500);
00958 }
00959 else {
00960 evlength_cc->SetBins(250,0,500);
00961 evlength_nc->SetBins(250,0,500);
00962 }
00963 }
00964
00965 }
00966 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00967 detector = 1;
00968
00969 if (pdfbinning==1 && evlength_cc->GetBinLowEdge(250)>200) {
00970 evlength_cc->SetBins(150,0,300);
00971 evlength_nc->SetBins(150,0,300);
00972 }
00973 if (pdfbinning==0 && evlength_cc->GetBinLowEdge(50)>200) {
00974 evlength_cc->SetBins(60,0,300);
00975 evlength_nc->SetBins(60,0,300);
00976 }
00977
00978
00979 }
00980
00981
00982
00983 Int_t mcevent=-1;
00984
00985 if(LoadTruthForRecoTH(j,mcevent)) {
00986
00987
00988 cc_nc = IAction(mcevent);
00989 flavour = Flavour(mcevent);
00990 neu_e = TrueNuEnergy(mcevent);
00991
00992
00993 if(tarpos==2) weight=1;
00994 if(tarpos==1) weight=1;
00995 if(tarpos==3) weight=1;
00996 }
00997
00998 int track_index = -1;
00999 LoadLargestTrackFromEvent(j,track_index);
01000
01001 if (track_index!=-1) {
01002
01003
01004 Int_t trkpass=ntpTrack->fit.pass;
01005 Int_t nd_reclaim=0;
01006
01007 Int_t trkplbu = ntpTrack->plane.begu;
01008 Int_t trkplbv = ntpTrack->plane.begv;
01009 Int_t trkpleu = ntpTrack->plane.endu;
01010 Int_t trkplev = ntpTrack->plane.endv;
01011
01012
01013 if (abs(trkplbu-trkplbv)<=5 && abs(trkpleu-trkplev)<=40 && (trkpleu<270&&trkplev<270)) nd_reclaim=1;
01014
01015
01016
01017 if (detector==1 && trkpass==0 && nd_reclaim==1) trkpass=1;
01018
01019
01020 if (is_fid && trkpass) {
01021
01022
01023 Float_t evlength=0;
01024 Float_t phfrac=0;
01025 Float_t phplane=0;
01026
01027 evlength=ntpEvent->plane.n;
01028 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
01029 evlength=ntpEvent->plane.end-ntpEvent->plane.beg+1;
01030 }
01031
01032 Float_t phtrack=ntpTrack->ph.sigcor;
01033 Float_t phevent=ntpEvent->ph.sigcor;
01034 if (phevent>0) {phfrac=phtrack/phevent;}
01035
01036 Float_t trkplane=ntpTrack->plane.n;
01037 if (trkplane>0) {phplane=phtrack/trkplane;}
01038
01039 if(flavour==2 && cc_nc==1) {
01040 evlength_cc->Fill(evlength,weight);
01041 phfrac_cc->Fill(phfrac,weight);
01042 phplane_cc->Fill(phplane,weight);
01043 }
01044
01045 else if(cc_nc==0) {
01046 evlength_nc->Fill(evlength,weight);
01047 phfrac_nc->Fill(phfrac,weight);
01048 phplane_nc->Fill(phplane,weight);
01049 }
01050
01051
01052 }
01053 }
01054 }
01055 }
01056
01057
01058
01059 save.Write();
01060 save.Close();
01061 }
01062
01063
01064 Float_t MadPIDAnalysis::SingleTimeFrame(Int_t snarlentry,Int_t nentries){
01065
01066 Int_t tf=0,tfbef=0,tfaft=0;
01067
01068
01069 if(this->GetEntry(snarlentry)) tf = ntpHeader->GetTimeFrame();
01070 if(snarlentry<(nentries-1) && this->GetEntry(snarlentry+1)) tfaft = ntpHeader->GetTimeFrame();
01071 if(snarlentry>0 && this->GetEntry(snarlentry-1)) tfbef = ntpHeader->GetTimeFrame();
01072
01073 Float_t singletimeframe=0;
01074 if(snarlentry<(nentries-1) && snarlentry>0 ){
01075 if(tf!=tfaft && tf!=tfbef) singletimeframe=1;
01076 }
01077 if(snarlentry==(nentries-1)) {
01078 if(tf!=tfbef) singletimeframe=1;
01079 }
01080 if(snarlentry==0) {
01081 if(tf!=tfaft) singletimeframe=1;
01082 }
01083
01084 this->GetEntry(snarlentry);
01085 return singletimeframe;
01086
01087 }
01088
01089
01090 void MadPIDAnalysis::CreatePAN(std::string tag, Int_t scanfilter, Int_t mcneartype, EnergyCorrections::WhichCorrection_t corrver){
01091
01092
01093
01094
01095 std::string savename = "PAN_" + tag + ".root";
01096 TFile save(savename.c_str(),"RECREATE");
01097 save.cd();
01098
01099
01100 SpillInfo pppinfo;
01101
01102
01103
01104
01105 NtpSRDataQuality *ntpDataQual;
01106 NtpSRDetStatus *ntpDetStatus;
01107
01108
01109
01110
01111 BMSpillAna bm;
01112 bm.UseDatabaseCuts();
01113
01114
01115
01116
01117 Float_t ann_aTotrk =0.;
01118 Float_t ann_aTotstp =0.;
01119 Float_t ann_aTotph =0.;
01120 Float_t ann_aTotlen =0.;
01121 Float_t ann_aPhperpl =0.;
01122 Float_t ann_aPhperstp =0.;
01123
01124
01125 Float_t ann_aTrkpass =0.;
01126 Float_t ann_aTrkph =0.;
01127 Float_t ann_aTrklen =0.;
01128 Float_t ann_aTrkphperpl =0.;
01129 Float_t ann_aTrkphper =0.;
01130 Float_t ann_aTrkplu =0.;
01131 Float_t ann_aTrkplv =0.;
01132 Float_t ann_aTrkstp =0.;
01133
01134 Float_t ann_aShwph =0.;
01135 Float_t ann_aShwstp =0.;
01136 Float_t ann_aShwdig =0.;
01137 Float_t ann_aShwpl =0.;
01138 Float_t ann_aShwphper =0.;
01139 Float_t ann_aShwphperpl =0.;
01140 Float_t ann_aShwphperdig =0.;
01141 Float_t ann_aShwphperstp =0.;
01142 Float_t ann_aShwplu =0.;
01143 Float_t ann_aShwplv =0.;
01144 Float_t ann_aShwstpu =0.;
01145 Float_t ann_aShwstpv =0.;
01146 Float_t ann_aPh3 =0.;
01147 Float_t ann_aPh6 =0.;
01148 Float_t ann_aPhlast =0.;
01149 Float_t ann_aPhcommon =0.;
01150
01151
01152
01153 Float_t true_enu;
01154 Float_t true_pxnu;
01155 Float_t true_pynu;
01156 Float_t true_pznu;
01157 Float_t true_etgt;
01158 Float_t true_pxtgt;
01159 Float_t true_pytgt;
01160 Float_t true_pztgt;
01161 Float_t true_emu;
01162 Float_t true_eshw;
01163 Float_t true_x;
01164 Float_t true_y;
01165 Float_t true_q2;
01166 Float_t true_w2;
01167 Float_t true_dircosneu;
01168 Float_t true_dircosz;
01169 Float_t true_vtxx;
01170 Float_t true_vtxy;
01171 Float_t true_vtxz;
01172
01173
01174 Int_t flavour;
01175 Int_t process;
01176 Int_t nucleus;
01177 Int_t cc_nc;
01178 Int_t initial_state;
01179 Int_t had_fs;
01180
01181
01182
01183
01184
01185 Int_t detector;
01186 Int_t run;
01187 Int_t snarl;
01188 Int_t nevent;
01189 Int_t event;
01190 Int_t subrun;
01191 Int_t trigsrc;
01192 Int_t mcevent;
01193 Int_t ntrack;
01194 Int_t nshower;
01195
01196 Int_t is_fid;
01197 Int_t is_fid_dp;
01198 Int_t is_fid_ns;
01199 Int_t is_fid_ab;
01200 Int_t is_fid_2008;
01201
01202
01203 Int_t is_cev;
01204 Int_t is_cont_trk;
01205 Int_t is_mc;
01206 Int_t pass_basic;
01207 Int_t pass_cccuts;
01208 Int_t pass_cccuts_2007;
01209 Int_t pass_cccuts_2008;
01210
01211 Float_t pid0;
01212 Float_t pid1;
01213 Float_t pid2;
01214 Float_t annpid;
01215 Float_t annpid_f1p1;
01216 Float_t annpid_f1p2;
01217 Float_t annpid_f1p3;
01218 Float_t annpid_f2p1;
01219 Float_t annpid_f2p2;
01220 Float_t annpid_f2p3;
01221
01222 Int_t pass;
01223 Int_t passcut;
01224
01225 Float_t reco_enu;
01226 Float_t reco_emu;
01227 Float_t reco_eshw;
01228
01229 Int_t pass_fd_qualcuts;
01230 Int_t litag;
01231 Int_t lisieve;
01232
01233
01234
01235 Float_t reco_eshwcc;
01236 Float_t reco_eshwnc;
01237 Float_t reco_eshwccw;
01238 Float_t reco_eshwncw;
01239
01240 Float_t reco_eshwcc_uncorr;
01241 Float_t reco_emu_uncorr;
01242
01243
01244 Float_t reco_eshw_sqrt;
01245 Float_t reco_qe_enu;
01246 Float_t reco_qe_q2;
01247 Float_t reco_y;
01248 Float_t reco_dircosneu;
01249 Float_t reco_dircosz;
01250 Float_t reco_vtxx;
01251 Float_t reco_vtxy;
01252 Float_t reco_vtxz;
01253
01254 Float_t evtlength;
01255 Float_t evtph;
01256 Float_t trklength;
01257 Float_t trkmomrange;
01258 Float_t trkqp;
01259 Float_t trkeqp;
01260 Float_t trkmombest;
01261 Float_t trkphfrac;
01262 Float_t trkphplane;
01263 Float_t trkds;
01264 Float_t rawph;
01265
01266 Float_t trkendz;
01267 Float_t trkendx;
01268 Float_t trkendy;
01269 Float_t trkendu;
01270 Float_t trkendv;
01271 Float_t trkvtxz;
01272 Float_t trkvtxx;
01273 Float_t trkvtxy;
01274 Float_t trkvtxu;
01275 Float_t trkvtxv;
01276 Float_t trkplbu;
01277 Float_t trkplbv;
01278 Float_t trkpleu;
01279 Float_t trkplev;
01280 Float_t trkple;
01281 Float_t trkplb;
01282 Float_t trkcosx;
01283 Float_t trkcosy;
01284 Float_t trkcosz;
01285 Float_t evanglet;
01286 Int_t trkfitpass;
01287 Int_t totdigits,totstrips;
01288 Int_t trkdigits,trkstrips;
01289 Float_t trkph,trkchi2,trkndof;
01290 Int_t trkdiffuv;
01291
01292
01293 Float_t shwph;
01294 Int_t shwdigits, shwstrips;
01295 Float_t shwvtxz;
01296 Float_t shwvtxx;
01297 Float_t shwvtxy;
01298 Float_t shwvtxu;
01299 Float_t shwvtxv;
01300
01301
01302
01303 Int_t shwncluster;
01304
01305
01306 Int_t scandecision;
01307 ScanList scan;
01308
01309
01310 Double_t evttimemin;
01311 Double_t evttimemax;
01312 Double_t snarlt;
01313 Double_t litime;
01314
01315
01316 Float_t pot,horn,tar,vvpos,hhpos,magn;
01317
01318 Float_t potdb,horndb,tardb,vvposdb,hhposdb,vvwidthdb,hhwidthdb;
01319 Double_t timedb;
01320
01321 Double_t neartdb,fartdb,neardifftdb;
01322
01323
01324 Float_t mcflux_tpx;
01325 Float_t mcflux_tpy;
01326 Float_t mcflux_tpz;
01327 Float_t mcflux_tptype;
01328 Float_t mcflux_nenergyN;
01329 Float_t mcflux_nenergyF;
01330 Float_t mcflux_nwtN;
01331 Float_t mcflux_nwtF;
01332 Float_t mcflux_nimpwt;
01333
01334 Int_t mcflux_fluxrun;
01335 Int_t mcflux_fluxevtno;
01336
01337
01338
01339 Int_t inu;
01340
01341
01342
01343
01344 Float_t abid_costheta;
01345 Float_t abid_eventenergy;
01346 Int_t abid_trackcharge;
01347 Int_t abid_trackenergy;
01348 Int_t abid_tracklikeplanes;
01349 Float_t abid_trackphfrac;
01350 Float_t abid_trackphmean;
01351 Float_t abid_trackqpsigmaqp;
01352 Float_t abid_y;
01353
01354
01355
01356
01357
01358 Float_t roid_numplanes;
01359 Float_t roid_meansigcor;
01360 Float_t roid_lowphhiph;
01361 Float_t roid_trkphfracclose;
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373 Int_t pass_fd_qualcuts_db;
01374
01375 Int_t dataqual_cratemask;
01376 Int_t dataqual_spilltype;
01377 Int_t dataqual_spillerror;
01378 Int_t dataqual_spillstatus;
01379 Int_t dataqual_coldchips;
01380
01381 Float_t dataqual_coilcurr_sm1;
01382 Float_t dataqual_coilcurr_sm2;
01383
01384
01385 Int_t alec_hvstatus;
01386 Int_t alec_coilstatus;
01387 Int_t alec_spillstatus;
01388
01389
01390
01391
01392
01393 Int_t pass_fd_timing;
01394
01395 Int_t nd_reclaim;
01396
01397
01398 Int_t coilok;
01399 Int_t coilpol;
01400
01401 Int_t numplanes;
01402
01403
01404 Float_t evtcom;
01405 Float_t evtslc;
01406 Float_t evtpur;
01407
01408
01409
01410
01411
01412 Float_t singletimeframe;
01413
01414 SpillInfoBlock *spinfo = new SpillInfoBlock();
01415 spinfo->Zero();
01416 singletimeframe=-999;
01417
01418 Int_t day,month,year;
01419 Int_t pass_beamcuts;
01420
01421
01422
01423
01424
01425 TTree *tree = new TTree("pan","pan");
01426
01427
01428 tree->Branch("true_enu",&true_enu,"true_enu/F",512000);
01429 tree->Branch("true_pxnu",&true_pxnu,"true_pxnu/F",512000);
01430 tree->Branch("true_pynu",&true_pynu,"true_pynu/F",512000);
01431 tree->Branch("true_pznu",&true_pznu,"true_pznu/F",512000);
01432 tree->Branch("true_etgt",&true_etgt,"true_etgt/F",512000);
01433 tree->Branch("true_pxtgt",&true_pxtgt,"true_pxtgt/F",512000);
01434 tree->Branch("true_pytgt",&true_pytgt,"true_pytgt/F",512000);
01435 tree->Branch("true_pztgt",&true_pztgt,"true_pztgt/F",512000);
01436 tree->Branch("true_emu",&true_emu,"true_emu/F",512000);
01437 tree->Branch("true_eshw",&true_eshw,"true_eshw/F",512000);
01438 tree->Branch("true_x",&true_x,"true_x/F",512000);
01439 tree->Branch("true_y",&true_y,"true_y/F",512000);
01440 tree->Branch("true_q2",&true_q2,"true_q2/F",512000);
01441 tree->Branch("true_w2",&true_w2,"true_w2/F",512000);
01442 tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",512000);
01443 tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",512000);
01444 tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",512000);
01445 tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",512000);
01446 tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",512000);
01447
01448 tree->Branch("flavour",&flavour,"flavour/I",512000);
01449 tree->Branch("process",&process,"process/I",512000);
01450 tree->Branch("nucleus",&nucleus,"nucleus/I",512000);
01451 tree->Branch("cc_nc",&cc_nc,"cc_nc/I",512000);
01452 tree->Branch("initial_state",&initial_state,"initial_state/I",512000);
01453 tree->Branch("had_fs",&had_fs,"had_fs/I",512000);
01454
01455
01456
01457 tree->Branch("detector",&detector,"detector/I",512000);
01458 tree->Branch("run",&run,"run/I",512000);
01459 tree->Branch("snarl",&snarl,"snarl/I",512000);
01460 tree->Branch("event",&event,"event/I",512000);
01461 tree->Branch("mcevent",&mcevent,"mcevent/I",512000);
01462 tree->Branch("ntrack",&ntrack,"ntrack/I",512000);
01463 tree->Branch("nshower",&nshower,"nshower/I",512000);
01464 tree->Branch("subrun",&subrun,"subrun/I",512000);
01465 tree->Branch("trigsrc",&trigsrc,"trigsrc/I",512000);
01466 tree->Branch("litime",&litime,"litime/D",512000);
01467 tree->Branch("is_fid",&is_fid,"is_fid/I",512000);
01468 tree->Branch("is_fid_dp",&is_fid_dp,"is_fid_dp/I",512000);
01469 tree->Branch("is_fid_ns",&is_fid_ns,"is_fid_ns/I",512000);
01470
01471
01472 tree->Branch("pass_fd_qualcuts", &pass_fd_qualcuts, "pass_fd_qualcuts/I");
01473 tree->Branch("litag", &litag, "litag/I");
01474 tree->Branch("lisieve", &lisieve, "lisieve/I");
01475
01476
01477 tree->Branch("is_cev",&is_cev,"is_cev/I",512000);
01478 tree->Branch("is_cont_trk",&is_cont_trk,"is_cont_trk/I",512000);
01479 tree->Branch("is_mc",&is_mc,"is_mc/I",512000);
01480 tree->Branch("pass_basic",&pass_basic,"pass_basic/I",512000);
01481 tree->Branch("pass_cccuts",&pass_cccuts,"pass_cccuts/I",512000);
01482 tree->Branch("pass_cccuts_2007",&pass_cccuts_2007,"pass_cccuts_2007/I",512000);
01483 tree->Branch("pass_cccuts_2008",&pass_cccuts_2008,"pass_cccuts_2008/I",512000);
01484 tree->Branch("pid0",&pid0,"pid0/F",512000);
01485 tree->Branch("pid1",&pid1,"pid1/F",512000);
01486 tree->Branch("pid2",&pid2,"pid2/F",512000);
01487 tree->Branch("annpid",&annpid,"annpid/F",512000);
01488 tree->Branch("annpid_f1p1",&annpid_f1p1,"annpid_f1p1/F",512000);
01489 tree->Branch("annpid_f1p2",&annpid_f1p2,"annpid_f1p2/F",512000);
01490 tree->Branch("annpid_f1p3",&annpid_f1p3,"annpid_f1p3/F",512000);
01491 tree->Branch("annpid_f2p1",&annpid_f2p1,"annpid_f2p1/F",512000);
01492 tree->Branch("annpid_f2p2",&annpid_f2p2,"annpid_f2p2/F",512000);
01493 tree->Branch("annpid_f2p3",&annpid_f2p3,"annpid_f2p3/F",512000);
01494 tree->Branch("pass",&pass,"pass/I",512000);
01495 tree->Branch("passcut", &passcut, "passcut/I");
01496
01497 tree->Branch("reco_enu",&reco_enu,"reco_enu/F",512000);
01498 tree->Branch("reco_emu",&reco_emu,"reco_emu/F",512000);
01499 tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",512000);
01500 tree->Branch("reco_eshw_sqrt",&reco_eshw_sqrt,"reco_eshw_sqrt/F",512000);
01501 tree->Branch("reco_eshwcc", &reco_eshwcc, "reco_eshwcc/F",512000);
01502 tree->Branch("reco_eshwcc_uncorr", &reco_eshwcc_uncorr, "reco_eshwcc_uncorr/F",512000);
01503 tree->Branch("reco_emu_uncorr", &reco_emu_uncorr, "reco_emu_uncorr/F",512000);
01504 tree->Branch("reco_eshwnc", &reco_eshwnc, "reco_eshwnc/F",512000);
01505 tree->Branch("reco_eshwccw", &reco_eshwccw, "reco_eshwccw/F",512000);
01506 tree->Branch("reco_eshwncw", &reco_eshwncw, "reco_eshwncw/F",512000);
01507 tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",512000);
01508 tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",512000);
01509 tree->Branch("reco_y",&reco_y,"reco_y/F",512000);
01510 tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",512000);
01511 tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",512000);
01512 tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",512000);
01513 tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",512000);
01514 tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",512000);
01515
01516 tree->Branch("evtlength",&evtlength,"evtlength/F",512000);
01517 tree->Branch("evtph",&evtph,"evtph/F",512000);
01518 tree->Branch("trklength",&trklength,"trklength/F",512000);
01519 tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",512000);
01520 tree->Branch("trkqp",&trkqp,"trkqp/F",512000);
01521 tree->Branch("trkeqp",&trkeqp,"trkeqp/F",512000);
01522 tree->Branch("trkmombest",&trkmombest,"trkmombest/F",512000);
01523 tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",512000);
01524 tree->Branch("trkphplane",&trkphplane,"trkphplane/F",512000);
01525 tree->Branch("rawph",&rawph,"rawph/F",512000);
01526 tree->Branch("numplanes",&numplanes,"numplanes/I",512000);
01527
01528 tree->Branch("trkendz",&trkendz,"trkendz/F",512000);
01529 tree->Branch("trkendu",&trkendu,"trkendu/F",512000);
01530 tree->Branch("trkendv",&trkendv,"trkendv/F",512000);
01531 tree->Branch("trkendx",&trkendx,"trkendx/F",512000);
01532 tree->Branch("trkendy",&trkendy,"trkendy/F",512000);
01533 tree->Branch("trkplbu",&trkplbu,"trkplbu/F",512000);
01534 tree->Branch("trkplbv",&trkplbv,"trkplbv/F",512000);
01535 tree->Branch("trkpleu",&trkpleu,"trkpleu/F",512000);
01536 tree->Branch("trkplev",&trkplev,"trkplev/F",512000);
01537 tree->Branch("trkple",&trkple,"trkple/F",512000);
01538 tree->Branch("trkplb",&trkplb,"trkplb/F",512000);
01539 tree->Branch("trkcosx",&trkcosx,"trkcosx/F",512000);
01540 tree->Branch("trkcosy",&trkcosy,"trkcosy/F",512000);
01541 tree->Branch("trkcosz",&trkcosz,"trkcosz/F",512000);
01542
01543 tree->Branch("trkvtxz",&trkvtxz,"trkvtxz/F",512000);
01544 tree->Branch("trkvtxu",&trkvtxu,"trkvtxu/F",512000);
01545 tree->Branch("trkvtxv",&trkvtxv,"trkvtxv/F",512000);
01546 tree->Branch("trkvtxx",&trkvtxx,"trkvtxx/F",512000);
01547 tree->Branch("trkvtxy",&trkvtxy,"trkvtxy/F",512000);
01548 tree->Branch("trkds", &trkds,"trkds/F",512000);
01549 tree->Branch("evanglet",&evanglet,"evanglet/F",512000);
01550
01551 tree->Branch("evttimemin",&evttimemin,"evttimemin/D");
01552 tree->Branch("evttimemax",&evttimemax,"evttimemax/D");
01553
01554 tree->Branch("totdigits",&totdigits,"totdigits/I");
01555 tree->Branch("totstrips",&totstrips,"totstrips/I");
01556 tree->Branch("trkdigits",&trkdigits,"trkdigits/I");
01557 tree->Branch("trkstrips",&trkstrips,"trkstrips/I");
01558 tree->Branch("trkph",&trkph,"trkph/F");
01559 tree->Branch("trkdiffuv",&trkdiffuv,"trkdiffuv/I",512000);
01560
01561 tree->Branch("trkchi2",&trkchi2,"trkchi2/F");
01562 tree->Branch("trkndof",&trkndof,"trkndof/F");
01563 tree->Branch("trkfitpass",&trkfitpass,"trkfitpass/I",512000);
01564
01565 tree->Branch("shwdigits",&shwdigits,"shwdigits/I");
01566 tree->Branch("shwstrips",&shwstrips,"shwstrips/I");
01567 tree->Branch("shwph",&shwph,"shwph/F");
01568
01569 tree->Branch("shwvtxz",&shwvtxz,"shwvtxz/F",512000);
01570 tree->Branch("shwvtxu",&shwvtxu,"shwvtxu/F",512000);
01571 tree->Branch("shwvtxv",&shwvtxv,"shwvtxv/F",512000);
01572 tree->Branch("shwvtxx",&shwvtxx,"shwvtxx/F",512000);
01573 tree->Branch("shwvtxy",&shwvtxy,"shwvtxy/F",512000);
01574 tree->Branch("shwncluster",&shwncluster,"shwncluster/I");
01575
01576
01577 tree->Branch("scandecision",&scandecision,"scandecision/I");
01578
01579
01580 tree->Branch("pot", &pot, "pot/F",512000);
01581 tree->Branch("horn", &horn, "horn/F",512000);
01582 tree->Branch("tar", &tar, "tar/F",512000);
01583 tree->Branch("vvpos", &vvpos, "vvpos/F",512000);
01584 tree->Branch("hhpos", &hhpos, "hhpos/F",512000);
01585 tree->Branch("magn", &magn, "magn/F",512000);
01586 tree->Branch("singletimeframe",&singletimeframe,"singletimeframe/F",512000);
01587
01588
01589
01590
01591 tree->Branch("potdb", &potdb, "potdb/F",512000);
01592 tree->Branch("horndb", &horndb, "horndb/F",512000);
01593 tree->Branch("tardb", &tardb, "tardb/F",512000);
01594 tree->Branch("vvposdb", &vvposdb, "vvposdb/F",512000);
01595 tree->Branch("hhposdb", &hhposdb, "hhposdb/F",512000);
01596 tree->Branch("vvwidthdb", &vvwidthdb, "vvwidthdb/F",512000);
01597 tree->Branch("hhwidthdb", &hhwidthdb, "hhwidthdb/F",512000);
01598 tree->Branch("timedb", &timedb, "timedb/D");
01599
01600
01601 tree->Branch("neartdb", &neartdb, "neartdb/D");
01602 tree->Branch("fartdb", &fartdb, "fartdb/D");
01603 tree->Branch("neardifftdb", &neardifftdb,"neardifftdb/D");
01604 tree->Branch("snarlt", &snarlt, "snarlt/D");
01605
01606
01607 tree->Branch("mcflux_tpx", &mcflux_tpx, "mcflux_tpx/F");
01608 tree->Branch("mcflux_tpy", &mcflux_tpy, "mcflux_tpy/F");
01609 tree->Branch("mcflux_tpz", &mcflux_tpz, "mcflux_tpz/F");
01610 tree->Branch("mcflux_tptype", &mcflux_tptype, "mcflux_tptype/F");
01611 tree->Branch("mcflux_nenergyN", &mcflux_nenergyN, "mcflux_nenergyN/F");
01612 tree->Branch("mcflux_nenergyF", &mcflux_nenergyF, "mcflux_nenergyF/F");
01613 tree->Branch("mcflux_nwtN", &mcflux_nwtN, "mcflux_nwtN/F");
01614 tree->Branch("mcflux_nwtF", &mcflux_nwtF, "mcflux_nwtF/F");
01615 tree->Branch("mcflux_nimpwt", &mcflux_nimpwt, "mcflux_nimpwt/F");
01616 tree->Branch("mcflux_fluxrun", &mcflux_fluxrun, "mcflux_fluxrun/I");
01617 tree->Branch("mcflux_fluxevtno", &mcflux_fluxevtno, "mcflux_fluxevtno/I");
01618
01619
01620
01621 tree->Branch("ann_aTotrk", &ann_aTotrk, "ann_aTotrk/F");
01622 tree->Branch("ann_aTotstp", &ann_aTotstp, "ann_aTotstp/F");
01623 tree->Branch("ann_aTotph", &ann_aTotph, "ann_aTotph/F");
01624 tree->Branch("ann_aTotlen", &ann_aTotlen, "ann_aTotlen/F");
01625 tree->Branch("ann_aPhperpl", &ann_aPhperpl, "ann_aPhperpl/F");
01626 tree->Branch("ann_aPhperstp", &ann_aPhperstp, "ann_aPhperstp/F");
01627 tree->Branch("ann_aTrkpass", &ann_aTrkpass, "ann_aTrkpass/F");
01628 tree->Branch("ann_aTrkph", &ann_aTrkph, "ann_aTrkph/F");
01629 tree->Branch("ann_aTrklen", &ann_aTrklen, "ann_aTrklen/F");
01630 tree->Branch("ann_aTrkphperpl", &ann_aTrkphperpl, "ann_aTrkphperpl/F");
01631 tree->Branch("ann_aTrkphper", &ann_aTrkphper, "ann_aTrkphper/F");
01632 tree->Branch("ann_aTrkplu", &ann_aTrkplu, "ann_aTrkplu/F");
01633 tree->Branch("ann_aTrkplv", &ann_aTrkplv, "ann_aTrkplv/F");
01634 tree->Branch("ann_aTrkstp", &ann_aTrkstp, "ann_aTrkstp/F");
01635 tree->Branch("ann_aShwph", &ann_aShwph, "ann_aShwph/F");
01636 tree->Branch("ann_aShwstp", &ann_aShwstp, "ann_aShwstp/F");
01637 tree->Branch("ann_aShwdig", &ann_aShwdig, "ann_aShwdig/F");
01638 tree->Branch("ann_aShwpl", &ann_aShwpl, "ann_aShwpl/F");
01639 tree->Branch("ann_aShwphper", &ann_aShwphper, "ann_aShwphper/F");
01640 tree->Branch("ann_aShwphperpl", &ann_aShwphperpl, "ann_aShwphperpl/F");
01641 tree->Branch("ann_aShwphperdig", &ann_aShwphperdig,"ann_aShwphperdig/F");
01642 tree->Branch("ann_aShwphperstp", &ann_aShwphperstp,"ann_aShwphperstp/F");
01643 tree->Branch("ann_aShwplu", &ann_aShwplu, "ann_aShwplu/F");
01644 tree->Branch("ann_aShwplv", &ann_aShwplv, "ann_aShwplv/F");
01645 tree->Branch("ann_aShwstpu", &ann_aShwstpu, "ann_aShwstpu/F");
01646 tree->Branch("ann_aShwstpv", &ann_aShwstpv, "ann_aShwstpv/F");
01647 tree->Branch("ann_aPh3", &ann_aPh3, "ann_aPh3/F");
01648 tree->Branch("ann_aPh6", &ann_aPh6, "ann_aPh6/F");
01649 tree->Branch("ann_aPhlast", &ann_aPhlast, "ann_aPhlast/F");
01650 tree->Branch("ann_aPhcommon", &ann_aPhcommon, "ann_aPhcommon/F");
01651
01652 tree->Branch("day", &day, "day/I");
01653 tree->Branch("month", &month, "month/I");
01654 tree->Branch("year", &year, "year/I");
01655 tree->Branch("pass_beamcuts", &pass_beamcuts, "pass_beamcuts/I");
01656
01657
01658
01659
01660 tree->Branch("abid_costheta", &abid_costheta, "abid_costheta/F");
01661 tree->Branch("abid_eventenergy", &abid_eventenergy, "abid_eventenergy/F");
01662 tree->Branch("abid_trackcharge", &abid_trackcharge, "abid_trackcharge/I");
01663 tree->Branch("abid_trackenergy", &abid_trackenergy, "abid_trackenergy/i");
01664 tree->Branch("abid_tracklikeplanes", &abid_tracklikeplanes, "abid_tracklikeplanes/I");
01665 tree->Branch("abid_trackphfrac", &abid_trackphfrac, "abid_trackphfrac/F");
01666 tree->Branch("abid_trackphmean", &abid_trackphmean, "abid_trackphmean/F");
01667 tree->Branch("abid_trackqpsigmaqp", &abid_trackqpsigmaqp, "abid_trackqpsigmaqp/F");
01668 tree->Branch("abid_y", &abid_y, "abid_y/F");
01669
01670
01671
01672
01673 tree->Branch("roid_numplanes", &roid_numplanes, "roid_numplanes/F");
01674 tree->Branch("roid_meansigcor", &roid_meansigcor, "roid_meansigcor/F");
01675 tree->Branch("roid_lowphhiph", &roid_lowphhiph, "roid_lowphhiph/F");
01676 tree->Branch("roid_trkphfracclose", &roid_trkphfracclose, "roid_trkphfracclose/F");
01677
01678
01679
01680
01681
01682
01683
01684
01685 tree->Branch("inu", &inu, "inu/I");
01686
01687
01688
01689
01690
01691 tree->Branch("pass_fd_qualcuts_db", &pass_fd_qualcuts_db, "pass_fd_qualcuts_db/I");
01692
01693 tree->Branch("dataqual_cratemask", &dataqual_cratemask, "dataqual_cratemask/I");
01694 tree->Branch("dataqual_spilltype", &dataqual_spilltype, "dataqual_spilltype/I");
01695 tree->Branch("dataqual_spillerror", &dataqual_spillerror, "dataqual_spillerror/I");
01696 tree->Branch("dataqual_spillstatus", &dataqual_spillstatus, "dataqual_spillstatus/I");
01697 tree->Branch("dataqual_coldchips", &dataqual_coldchips, "dataqual_coldchips/I");
01698 tree->Branch("dataqual_coilcurr_sm1", &dataqual_coilcurr_sm1, "dataqual_coilcurr_sm1/F");
01699 tree->Branch("dataqual_coilcurr_sm2", &dataqual_coilcurr_sm2, "dataqual_coilcurr_sm2/F");
01700
01701
01702 tree->Branch("alec_hvstatus", &alec_hvstatus, "alec_hvstatus/I");
01703 tree->Branch("alec_coilstatus", &alec_coilstatus, "alec_coilstatus/I");
01704 tree->Branch("alec_spillstatus", &alec_spillstatus, "alec_spillstatus/I");
01705
01706
01707
01708
01709 tree->Branch("is_fid_ab",&is_fid_ab,"is_fid_ab/I",512000);
01710 tree->Branch("is_fid_2008",&is_fid_2008,"is_fid_2008/I",512000);
01711
01712
01713
01714 tree->Branch("pass_fd_timing",&pass_fd_timing,"pass_fd_timing/I",512000);
01715
01716
01717 tree->Branch("nd_reclaim",&nd_reclaim,"nd_reclaim/I",512000);
01718
01719 tree->Branch("coilok",&coilok,"coilok/I",512000);
01720 tree->Branch("coilpol",&coilpol,"coilpol/I",512000);
01721
01722 tree->Branch("evtcom", &evtcom,"evtcom/F");
01723 tree->Branch("evtslc", &evtslc,"evtslc/F");
01724 tree->Branch("evtpur", &evtpur,"evtpur/F");
01725
01726
01727
01728
01729
01730 TTree *potcount = new TTree("potcount","potcount");
01731 potcount->Branch("potall",&potall,"potall/F",512000);
01732 potcount->Branch("potcut",&potcut,"potcut/F",512000);
01733
01734
01735
01736
01737
01738 Bool_t first_ann = true;
01739
01740
01741
01742 for(int i=0;i<Nentries;i++){
01743
01744 if(i%400==0) std::cout << float(i)*100./float(Nentries)
01745 << "% done" << std::endl;
01746
01747 if(!this->GetEntry(i)) continue;
01748
01749
01750 nevent = eventSummary->nevent;
01751
01752
01753
01754
01755
01756 run = ntpHeader->GetRun();
01757 snarl = ntpHeader->GetSnarl();
01758 subrun = ntpHeader->GetSubRun();
01759 trigsrc = ntpHeader->GetTrigSrc();
01760 day = eventSummary->date.day;
01761 month = eventSummary->date.month;
01762 year = eventSummary->date.year;
01763
01764
01765 Double_t snarltime =ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9);
01766
01767 litime = eventSummary->litime;
01768
01769
01770
01771
01772
01773 pass_fd_qualcuts =passfail(snarltime);
01774
01775
01776
01777
01778
01779
01780 ReleaseType::Release_t reltype=strecord->GetRelease();
01781
01782
01783
01784
01785
01786
01787 ntpDataQual = &(strecord->dataquality);
01788 ntpDetStatus = &(strecord->detstatus);
01789
01790
01791 pass_fd_qualcuts_db=1;
01792
01793 dataqual_cratemask=0;
01794 dataqual_spilltype=0;
01795 dataqual_spillerror=0;
01796 dataqual_spillstatus=0;
01797 dataqual_coldchips=0;
01798 dataqual_coilcurr_sm1=0;
01799 dataqual_coilcurr_sm2=0;
01800
01801
01802
01803 alec_hvstatus=1;
01804 alec_coilstatus=1;
01805 alec_spillstatus=1;
01806
01807 coilok=0;
01808 coilpol=0;
01809
01810
01811
01812
01813
01814 if(ntpHeader->GetVldContext().GetSimFlag()==1 && ntpHeader->GetVldContext().GetDetector()==2) {
01815
01816
01817
01818 HvStatusFinder& hv = HvStatusFinder::Instance();
01819 HvStatus::HvStatus_t hv_ok=hv.GetHvStatus(ntpHeader->GetVldContext(),60,1);
01820
01821 SpillServerMonFinder& smon=SpillServerMonFinder::Instance();
01822 const SpillServerMon& spill_near=smon.GetNearestSpill(ntpHeader->GetVldContext());
01823 VldTimeStamp dt=spill_near.GetSpillTime()-ntpHeader->GetVldContext().GetTimeStamp();
01824
01825 Int_t dt_sec=abs(dt.GetSec());
01826 Int_t gps_error=spill_near.GetSpillTimeError();
01827
01828 Bool_t coil_ok=CoilTools::IsOK(ntpHeader->GetVldContext());
01829
01830
01831 if(!HvStatus::Good(hv_ok)) alec_hvstatus=0;
01832 if(!coil_ok) alec_coilstatus=0;
01833 if (dt_sec<360 && gps_error>1000) alec_spillstatus=0;
01834
01835
01836 dataqual_cratemask=ntpDataQual->cratemask;
01837 dataqual_spilltype=ntpDataQual->spilltype;
01838 dataqual_spillerror=ntpDataQual->spilltimeerror;
01839 dataqual_spillstatus=ntpDataQual->spillstatus;
01840 dataqual_coldchips=ntpDataQual->coldchips;
01841 dataqual_coilcurr_sm1=ntpDetStatus->coilcurrent1;
01842 dataqual_coilcurr_sm2=ntpDetStatus->coilcurrent2;
01843
01844
01845
01846 if (!DataUtil::IsGoodFDData(strecord)) pass_fd_qualcuts_db=0;
01847
01848
01849 }
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859 Bool_t ro_pass=false;
01860 ro_pass=ro_interface.FillSnarl(strecord);
01861
01862
01863
01864
01865
01866 litag=0;
01867
01868 Int_t numshower=eventSummary->nshower;
01869 Float_t allph=eventSummary->ph.raw;
01870 Int_t plbeg=eventSummary->plane.beg;
01871 Int_t plend=eventSummary->plane.end;
01872
01873 if (numshower) {
01874 if (allph>1e6) litag+=10;
01875
01876 if (((plbeg==1 || plbeg==2) && (plend==63 || plend==64)) ||
01877 ((plbeg==65 || plbeg==66) && (plend==127 || plend==128)) ||
01878 ((plbeg==129 || plbeg==130) && (plend==191 || plend==192)) ||
01879 ((plbeg==193 || plbeg==194) && (plend==247 || plend==248))) litag++;
01880 if (((plbeg==250 || plbeg==251) && (plend==312 || plend==313)) ||
01881 ((plbeg==314 || plbeg==315) && (plend==376 || plend==377)) ||
01882 ((plbeg==378 || plbeg==379) && (plend==440 || plend==441)) ||
01883 ((plbeg==442 || plbeg==443) && (plend==484 || plend==485))) litag++;
01884 }
01885
01886
01887
01888
01889
01890 lisieve=0;
01891
01892 if (ntpHeader->GetVldContext().GetDetector()==2) {
01893
01894 if (LISieve::IsLI(*strecord)) lisieve=1;
01895
01896 }
01897
01898
01899
01900
01901 magn=0;
01902
01903
01904
01905
01906
01907
01908 scandecision=scan.GetDecision(run,snarl);
01909
01910
01911 if (scanfilter==0 || scandecision>0) {
01912
01913
01914
01915 if(ntpHeader->GetVldContext().GetSimFlag()!=4){
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933 potdb = -999;
01934 horndb = -999;
01935 tardb = -999;
01936 vvposdb = -999;
01937 hhposdb = -999;
01938 vvwidthdb = -999;
01939 hhwidthdb = -999;
01940 timedb = -999;
01941 neartdb = -999;
01942 fartdb = -999;
01943 neardifftdb = -999;
01944
01945
01946
01947 Detector::Detector_t detectort;
01948 SimFlag::SimFlag_t simflag;
01949 VldTimeStamp timestamp;
01950 VldContext cx;
01951 cx = ntpHeader->GetVldContext();
01952 detectort = ntpHeader->GetVldContext().GetDetector();
01953 simflag = ntpHeader->GetVldContext().GetSimFlag();
01954 timestamp = ntpHeader->GetVldContext().GetTimeStamp();
01955
01956 pass_beamcuts=0;
01957
01958
01959
01960
01961 Bool_t coil_ok=CoilTools::IsOK(cx);
01962 Bool_t coil_pol=CoilTools::IsReverse(cx);
01963
01964 coilpol=1;
01965
01966 if (coil_pol) coilpol=-1;
01967 if (coil_ok) coilok=1;
01968
01969
01970
01971 const Dcs_Mag_Near* magnear =
01972 CoilTools::Instance().GetMagNear(cx);
01973
01974 if (magnear) {
01975 magn=magnear->GetCurrent();
01976 }
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991 BDSpillAccessor &spillAccess=BDSpillAccessor::Get();
01992 const BeamMonSpill *spillInfoDB = spillAccess.LoadSpill(timestamp);
01993 if(spillInfoDB) {
01994 if(spillInfoDB->fTortgt !=0.) potdb = spillInfoDB->fTortgt;
01995 else if(spillInfoDB->fTrtgtd !=0.) potdb = spillInfoDB->fTrtgtd;
01996 else if(spillInfoDB->fTor101 !=0.) potdb = spillInfoDB->fTor101;
01997 else potdb = -999;
01998
01999 horndb = spillInfoDB->fHornCur;
02000 tardb = (int)spillInfoDB->GetStatusBits().target_in;
02001 vvposdb = spillInfoDB->fTargBpmY[0];
02002 hhposdb = spillInfoDB->fTargBpmX[0];
02003 vvwidthdb = spillInfoDB->fProfWidY;
02004 hhwidthdb = spillInfoDB->fProfWidX;
02005 timedb = spillInfoDB->SpillTime();
02006
02007
02008
02009
02010 bm.SetSpill(*spillInfoDB);
02011 bm.SetSnarlTime(timestamp);
02012
02013 if (bm.SelectSpill()==1) pass_beamcuts=1;
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026 if (potdb>0.2 && fabs(snarltime-timedb)<1) potall=potall+potdb;
02027
02028
02029
02030
02031
02032
02033 if (pass_beamcuts && horndb<-155 && coilok==1 && coilpol==1) potcut=potcut+potdb;
02034
02035 }
02036
02037
02038
02039
02040 SpillTimeFinder& spillfinder = SpillTimeFinder::Instance();
02041 fartdb = spillfinder.GetTimeOfNearestSpill(cx);
02042 const SpillTimeND& spnd = spillfinder.GetNearestSpill(cx);
02043 neartdb = spnd.GetTimeStamp().GetSec()+(spnd.GetTimeStamp().GetNanoSec()/1e9);
02044 neardifftdb = spillfinder.GetTimeToNearestSpill(cx);
02045 snarlt = ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9);
02046
02047 }
02048
02049
02050
02051 Int_t trkcount=0;
02052 Int_t shwcount=0;
02053 Float_t totph=0;
02054
02055 Int_t evtmin=0;
02056 Int_t evtmax=nevent;
02057
02058 if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) {
02059
02060
02061 Int_t maxevent=0;
02062 Float_t maxph=0;
02063 for(int j=0;j<nevent;j++){
02064 if(!LoadEvent(j)) continue;
02065
02066
02067
02068 Float_t evtpheight=ntpEvent->ph.sigcor;
02069 if (evtpheight>maxph) {
02070 maxph=evtpheight;
02071 maxevent=j;
02072 }
02073 }
02074
02075 evtmin=maxevent;
02076 evtmax=maxevent+1;
02077 }
02078
02079
02080
02081
02082
02083 for (int j=evtmin;j<evtmax;j++){
02084
02085 if(!LoadEvent(j)) continue;
02086
02087 totph = 0;
02088
02089
02090 event = j;
02091 ntrack = ntpEvent->ntrack;
02092 nshower = ntpEvent->nshower;
02093 totdigits=ntpEvent->ndigit;
02094 totstrips=ntpEvent->nstrip;
02095
02096
02097
02098
02099 mcflux_tpx=0; mcflux_tpy=0; mcflux_tpz=0; mcflux_tptype=0;
02100 mcflux_nenergyN=0; mcflux_nenergyF=0; mcflux_nwtN=0; mcflux_nwtF=0;
02101 mcflux_nimpwt=0; mcflux_fluxrun=0; mcflux_fluxevtno=0;
02102 true_enu = 0; true_emu = 0; true_eshw = 0;
02103 true_pxnu = 0; true_pynu = 0; true_pznu = 0;
02104 true_x = 0; true_y = 0; true_q2 = 0; true_w2 = 0;
02105 true_dircosneu = 0; true_dircosz = 0;
02106 true_vtxx = 0; true_vtxy = 0; true_vtxz = 0;
02107 flavour = 0; process = 0; nucleus = 0; cc_nc = 0;
02108 initial_state = 0; had_fs = 0;
02109 true_etgt = 0; true_pxtgt = 0; true_pytgt = 0; true_pztgt = 0;
02110 detector = -1; mcevent = -1;
02111 is_fid = 0; is_cev = 0; is_mc = 0; pass_basic = 0;
02112 pid0 = -999; pid1 = -999; pid2 = -999; pass = 0; passcut=0;
02113 reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_eshw_sqrt = 0;
02114
02115 reco_eshwcc_uncorr=0;
02116 reco_emu_uncorr=0;
02117
02118 reco_qe_enu = 0; reco_qe_q2 = 0; reco_y = 0; reco_dircosneu = 0;
02119 reco_dircosz = 0;
02120 reco_vtxx = -999; reco_vtxy = -999; reco_vtxz = -999;
02121 evtlength = -1; trklength = -1; trkphfrac = -1; trkphplane = -1;
02122 trkmomrange = -999; trkqp = -999; trkeqp = -999;
02123 trkendz=100.; trkendu=100.; trkendv=100.; trkendx=100.; trkendy=100.;
02124 trkcosx=0.; trkcosy=0.; trkcosz=0.; evanglet=0.;
02125 trkendz=999; trkendu=999; trkendv=999; trkendx=999;
02126 trkendy=-999;trkvtxz=-999; trkvtxu=-999; trkvtxv=-999; trkvtxx=-999; trkvtxy=-999;
02127 trkds=-1; evttimemin=-1; evttimemax=-1;
02128 evtph=0.;
02129 annpid=-999.;
02130 annpid_f1p1=-999.;annpid_f1p2=-999.;annpid_f1p3=-999.;
02131 annpid_f2p1=-999.;annpid_f2p2=-999.;annpid_f2p3=-999.;
02132 reco_eshwcc=0.;reco_eshwnc=0.;reco_eshwccw=0.;reco_eshwncw=0.;
02133 is_fid_dp=0;is_fid_ns=0;
02134 is_cont_trk=0; trkmombest=-999.;
02135
02136 trkfitpass=0;
02137 trkdiffuv=0;
02138 trkdigits=0; trkstrips=0; trkph=0; trkndof=0; trkchi2=0;
02139 shwdigits=0; shwstrips=0; shwph=0;
02140 shwvtxu=0;
02141 shwvtxv=0;
02142 shwvtxx=0;
02143 shwvtxy=0;
02144 shwvtxz=0;
02145 shwncluster=0;
02146 numplanes=0;
02147
02148 ann_aTotrk =0.;
02149 ann_aTotstp =0.;
02150 ann_aTotph =0.;
02151 ann_aTotlen =0.;
02152 ann_aPhperpl =0.;
02153 ann_aPhperstp =0.;
02154
02155
02156 ann_aTrkpass =0.;
02157 ann_aTrkph =0.;
02158 ann_aTrklen =0.;
02159 ann_aTrkphperpl =0.;
02160 ann_aTrkphper =0.;
02161 ann_aTrkplu =0.;
02162 ann_aTrkplv =0.;
02163 ann_aTrkstp =0.;
02164
02165 ann_aShwph =0.;
02166 ann_aShwstp =0.;
02167 ann_aShwdig =0.;
02168 ann_aShwpl =0.;
02169 ann_aShwphper =0.;
02170 ann_aShwphperpl =0.;
02171 ann_aShwphperdig =0.;
02172 ann_aShwphperstp =0.;
02173 ann_aShwplu =0.;
02174 ann_aShwplv =0.;
02175 ann_aShwstpu =0.;
02176 ann_aShwstpv =0.;
02177 ann_aPh3 =0.;
02178 ann_aPh6 =0.;
02179 ann_aPhlast =0.;
02180 ann_aPhcommon =0.;
02181
02182
02183 abid_costheta =0.;
02184 abid_eventenergy =0.;
02185 abid_trackcharge =0;
02186 abid_trackenergy =0;
02187 abid_tracklikeplanes =0;
02188 abid_trackphfrac =0.;
02189 abid_trackphmean =0.;
02190 abid_trackqpsigmaqp =0.;
02191 abid_y =0.;
02192
02193
02194
02195 roid_numplanes =0;
02196 roid_meansigcor =0.;
02197 roid_lowphhiph =0.;
02198 roid_trkphfracclose =0.;
02199
02200
02201
02202
02203
02204
02205 inu=0;
02206 pass_cccuts=0;
02207 pass_cccuts_2007=0;
02208 pass_cccuts_2008=0;
02209
02210 is_fid_ab=0;
02211 is_fid_2008=0;
02212
02213 nd_reclaim=0;
02214
02215
02216 evtcom=0.; evtpur=0.; evtslc=0.;
02217
02218
02219
02220
02221 rawph=ntpEvent->ph.raw;
02222
02223
02224 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
02225 detector = 2;
02226 totph=eventSummary->ph.sigcor;
02227
02228 is_fid = 1;
02229 if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
02230 }
02231 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
02232 detector = 1;
02233
02234
02235 for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) {
02236 LoadTrack(ii);
02237 totph+=ntpTrack->ph.sigcor;
02238 }
02239 trkcount+=ntrack;
02240
02241 for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) {
02242 LoadShower(ii);
02243 totph+=ntpShower->ph.sigcor;
02244 }
02245 shwcount+=nshower;
02246
02247
02248 is_fid = 1;
02249
02250
02251 if(!PassFidND(ntpEvent->index)) is_fid = 0;
02252 }
02253
02254
02255
02256
02257 if (detector==2) is_fid_dp=PassFidNew(ntpEvent->index);
02258 else is_fid_dp=PassFidND(ntpEvent->index);
02259
02260 if (detector==2) is_fid_ns=PassFidNewN(ntpEvent->index);
02261 else is_fid_ns=PassFidND(ntpEvent->index);
02262
02263
02264 if (detector==2) is_fid_ab=IsFidFD_AB(ntpEvent->index);
02265 else is_fid_ab=PassFidND(ntpEvent->index);
02266
02267
02268
02269
02270
02271 is_fid_2008=IsFid_2008(ntpEvent->index);
02272
02273
02274
02275
02276
02277 if(ntpHeader->GetVldContext().GetSimFlag()==4){
02278 if(LoadTruthForRecoTH(j,mcevent)) {
02279
02280 if(LoadTHEvent(j)) {
02281 evtcom = (ntpTHEvent->completeall);
02282 evtslc = (ntpTHEvent->completeslc);
02283 evtpur = (ntpTHEvent->purity);
02284 }
02285
02286 Float_t *vtx = TrueVtx(mcevent);
02287 Float_t *nu3mom = TrueNuMom(mcevent);
02288 Float_t *targ4vec = Target4Vector(mcevent);
02289
02290 is_mc = 1;
02291 nucleus = Nucleus(mcevent);
02292 flavour = Flavour(mcevent);
02293 initial_state = Initial_state(mcevent);
02294 had_fs = HadronicFinalState(mcevent);
02295 process = IResonance(mcevent);
02296 cc_nc = IAction(mcevent);
02297 true_enu = TrueNuEnergy(mcevent);
02298 true_pxnu = nu3mom[0];
02299 true_pynu = nu3mom[1];
02300 true_pznu = nu3mom[2];
02301 true_emu = TrueMuEnergy(mcevent);
02302 true_eshw = TrueShwEnergy(mcevent);
02303 true_x = X(mcevent);
02304 true_y = Y(mcevent);
02305 true_q2 = Q2(mcevent);
02306 true_w2 = W2(mcevent);
02307 true_dircosz = TrueMuDCosZ(mcevent);
02308 true_dircosneu = TrueMuDCosNeu(mcevent);
02309 true_vtxx = vtx[0];
02310 true_vtxy = vtx[1];
02311 true_vtxz = vtx[2];
02312 true_etgt = targ4vec[3];
02313 true_pxtgt = targ4vec[0];
02314 true_pytgt = targ4vec[1];
02315 true_pztgt = targ4vec[2];
02316 inu = INu(mcevent);
02317
02318 if (LoadTruth(mcevent)) {
02319 mcflux_tpx=ntpTruth->flux.tpx;
02320 mcflux_tpy=ntpTruth->flux.tpy;
02321 mcflux_tpz=ntpTruth->flux.tpz;
02322 mcflux_tptype=ntpTruth->flux.tptype;
02323 mcflux_nenergyN=ntpTruth->flux.nenergynear;
02324 mcflux_nenergyF=ntpTruth->flux.nenergyfar;
02325 mcflux_nwtN=ntpTruth->flux.nwtnear;
02326 mcflux_nwtF=ntpTruth->flux.nwtfar;
02327 mcflux_nimpwt=ntpTruth->flux.nimpwt;
02328 mcflux_fluxrun=ntpTruth->flux.fluxrun;
02329 mcflux_fluxevtno=ntpTruth->flux.fluxevtno;
02330 }
02331
02332 delete [] vtx;
02333 delete [] nu3mom;
02334 delete [] targ4vec;
02335 }
02336 }
02337
02338 if(PassBasicCuts(event)) {pass_basic=1;}
02339
02340
02341
02342
02343 reco_vtxx = ntpEvent->vtx.x;
02344 reco_vtxy = ntpEvent->vtx.y;
02345 reco_vtxz = ntpEvent->vtx.z;
02346 evtlength = ntpEvent->plane.end-ntpEvent->plane.beg+1;
02347 evtph = ntpEvent->ph.sigcor;
02348
02349 numplanes = ntpEvent->plane.n;
02350
02351
02352
02353
02354 AnnInputBlock anninput=AnnVar(event);
02355
02356 ann_aTotrk = anninput.aTotrk;
02357 ann_aTotstp = anninput.aTotstp;
02358 ann_aTotph = anninput.aTotph;
02359 ann_aTotlen = anninput.aTotlen;
02360 ann_aPhperpl = anninput.aPhperpl;
02361 ann_aPhperstp = anninput.aPhperstp;
02362
02363
02364 ann_aTrkpass = anninput.aTrkpass;
02365 ann_aTrkph = anninput.aTrkph;
02366 ann_aTrklen = anninput.aTrklen;
02367 ann_aTrkphperpl = anninput.aTrkphperpl;
02368 ann_aTrkphper = anninput.aTrkphper;
02369 ann_aTrkplu = anninput.aTrkplu;
02370 ann_aTrkplv = anninput.aTrkplv;
02371 ann_aTrkstp = anninput.aTrkstp;
02372
02373 ann_aShwph = anninput.aShwph;
02374 ann_aShwstp = anninput.aShwstp;
02375 ann_aShwdig = anninput.aShwdig;
02376 ann_aShwpl = anninput.aShwpl;
02377 ann_aShwphper = anninput.aShwphper;
02378 ann_aShwphperpl = anninput.aShwphperpl;
02379 ann_aShwphperdig = anninput.aShwphperdig ;
02380 ann_aShwphperstp = anninput.aShwphperstp;
02381 ann_aShwplu = anninput.aShwplu;
02382 ann_aShwplv = anninput.aShwplv;
02383 ann_aShwstpu = anninput.aShwstpu;
02384 ann_aShwstpv = anninput.aShwstpv;
02385 ann_aPh3 = anninput.aPh3;
02386 ann_aPh6 = anninput.aPh6;
02387 ann_aPhlast = anninput.aPhlast;
02388 ann_aPhcommon = anninput.aPhcommon;
02389
02390
02391
02392
02393
02394
02395
02396 evttimemin = anninput.aTimemin-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9;
02397 evttimemax = anninput.aTimemax-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9;
02398
02399
02400
02401 pass_fd_timing=1;
02402
02403 if (detector==2 && ntpHeader->GetVldContext().GetSimFlag()==1) {
02404 pass_fd_timing=PassTimingCuts(snarlt-fartdb+evttimemin);
02405 }
02406
02407
02408
02409
02410
02411
02412 Int_t target=1;
02413 if(ntpHeader->GetVldContext().GetSimFlag()!=4){
02414
02415
02416
02417
02418
02419 if(mcneartype==1) target=1;
02420 if(mcneartype==2) target=2;
02421 if(mcneartype==3) target=3;
02422 }
02423 if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==1){
02424 if(mcneartype==1) target=1;
02425 if(mcneartype==2) target=2;
02426 if(mcneartype==3) target=3;
02427 }
02428 if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==2){
02429 target=0;
02430 }
02431
02432
02433 Int_t fid =1;
02434 Int_t prior = 1;
02435
02436 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
02437 if(!PassFidND(event) && !PassFid(event) && is_fid_2008==0) annpid=-999;
02438 else {
02439 annpid=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
02440 first_ann=false;
02441 }
02442 }
02443 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
02444 if(!is_fid_ns){
02445 annpid_f2p1=-999;
02446 annpid_f2p2=-999;
02447 annpid_f2p3=-999;
02448 }
02449 if(is_fid_ns){
02450 fid=2;
02451 prior=1;
02452 annpid_f2p1=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
02453 prior=2;
02454 annpid_f2p2=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
02455 prior=3;
02456 annpid_f2p3=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
02457 first_ann=false;
02458 }
02459
02460
02461 if(!is_fid_2008){
02462
02463 annpid_f1p1=-999;
02464 annpid_f1p2=-999;
02465 annpid_f1p3=-999;
02466 }
02467
02468
02469
02470 if(is_fid_2008){
02471
02472 fid=1;
02473 prior=1;
02474 annpid_f1p1=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
02475 prior=2;
02476 annpid_f1p2=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
02477 prior=3;
02478 annpid_f1p3=GetAnnPid(anninput,detector,target,fid,prior,first_ann);
02479 first_ann=false;
02480 }
02481 }
02482
02483
02484
02485 bool isdata=true;
02486 if(ntpHeader->GetVldContext().GetSimFlag()==4) isdata=false;
02487
02488
02489 int track_index = -1;
02490 LoadLargestTrackFromEvent(j,track_index);
02491 int shower_index = -1;
02492
02493 LoadShower_Jim(j,shower_index,detector);
02494 if(shower_index==-1) nshower=0;
02495 if(shower_index!=-1){
02496
02497 shwdigits = ntpShower->ndigit;
02498 shwstrips = ntpShower->nstrip;
02499 shwph = ntpShower->ph.sigcor;
02500 shwvtxu = ntpShower->vtx.u;
02501 shwvtxv = ntpShower->vtx.v;
02502 shwvtxx = ntpShower->vtx.x;
02503 shwvtxy = ntpShower->vtx.y;
02504 shwvtxz = ntpShower->vtx.z;
02505 shwncluster = ntpShower->ncluster;
02506
02507 Int_t shwtype=0;
02508
02509
02510
02511
02512
02513
02514
02515
02516
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526
02527
02528
02529
02530
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543 shwtype=-1;
02544 reco_eshw = RecoShwEnergyNew(shower_index,shwtype,
02545 ntpHeader->GetVldContext(),
02546 reltype,corrver);
02547 shwtype=0;
02548 reco_eshwcc = RecoShwEnergyNew(shower_index,shwtype,
02549 ntpHeader->GetVldContext(),
02550 reltype,corrver);
02551
02552
02553
02554 shwtype=2;
02555 reco_eshwnc = ntpShower->shwph.linNCgev;
02556
02557 shwtype=1;
02558 reco_eshwccw = ntpShower->shwph.wtCCgev;
02559
02560 shwtype=3;
02561 reco_eshwncw = ntpShower->shwph.wtNCgev;
02562
02563
02564
02565
02566 reco_eshwcc_uncorr = ntpShower->shwph.linCCgev;
02567
02568
02569
02570
02571 reco_eshw_sqrt = RecoShwEnergySqrt(event);
02572
02573 }
02574
02575
02576 if(track_index!=-1){
02577
02578
02579
02580
02581
02582
02583 reco_emu = RecoMuEnergyNew(0,track_index,ntpHeader->GetVldContext(),
02584 reltype,corrver);
02585
02586
02587
02588
02589
02590 const float mumass=0.10566;
02591
02592 Float_t mr=ntpTrack->momentum.range;
02593
02594 float mc=0.0;
02595 if(ntpTrack->momentum.qp!=0) mc = 1.0/(ntpTrack->momentum.qp);
02596
02597 if(IsFidAll(track_index) || ntpTrack->momentum.qp==0) {
02598 reco_emu_uncorr=sqrt(mr*mr+ mumass*mumass);
02599 }
02600 else {
02601 if(ntpTrack->momentum.qp == 0.0) reco_emu_uncorr=10000.0;
02602 else reco_emu_uncorr=sqrt(mc*mc+ mumass*mumass);
02603 }
02604
02605
02606
02607 trklength = ntpTrack->plane.end-ntpTrack->plane.beg+1;
02608 trkmomrange = ntpTrack->momentum.range;
02609 trkqp = ntpTrack->momentum.qp;
02610 trkeqp = ntpTrack->momentum.eqp;
02611 trkendz = ntpTrack->end.z;
02612 trkendu = ntpTrack->end.u;
02613 trkendv = ntpTrack->end.v;
02614 trkendx = ntpTrack->end.x;
02615 trkendy = ntpTrack->end.y;
02616
02617 trkvtxz = ntpTrack->vtx.z;
02618 trkvtxu = ntpTrack->vtx.u;
02619 trkvtxv = ntpTrack->vtx.v;
02620 trkvtxx = ntpTrack->vtx.x;
02621 trkvtxy = ntpTrack->vtx.y;
02622 trkds = ntpTrack->ds;
02623 trkplbu = ntpTrack->plane.begu;
02624 trkplbv = ntpTrack->plane.begv;
02625 trkpleu = ntpTrack->plane.endu;
02626 trkplev = ntpTrack->plane.endv;
02627 trkple = ntpTrack->plane.end;
02628 trkplb = ntpTrack->plane.beg;
02629 trkfitpass = ntpTrack->fit.pass;
02630 trkdiffuv = abs(ntpTrack->plane.nu-ntpTrack->plane.nv);
02631 trkdigits = ntpTrack->ndigit;
02632 trkstrips = ntpTrack->nstrip;
02633 trkph = ntpTrack->ph.sigcor;
02634 trkndof = ntpTrack->fit.ndof;
02635 trkchi2 = ntpTrack->fit.chi2;
02636 trkcosx = ntpTrack->vtx.dcosx;
02637 trkcosy = ntpTrack->vtx.dcosy;
02638 trkcosz = ntpTrack->vtx.dcosz;
02639
02640
02641
02642 Float_t pbeamy,pbeamz,pbeamx;
02643
02644 pbeamy=0; pbeamz=0; pbeamx=0;
02645
02646 if(detector==2) {
02647 pbeamx=0.;
02648 pbeamy=sin(-3.33*3.14/180.);
02649 pbeamz=cos(-3.33*3.14/180.);
02650 }
02651
02652 if(detector==1) {
02653 pbeamx=0.;
02654 pbeamy=sin(3.33*3.14/180.);
02655 pbeamz=cos(3.33*3.14/180.);
02656 }
02657
02658 evanglet = trkcosx*pbeamx+trkcosy*pbeamy+trkcosz*pbeamz;
02659
02660
02661 Float_t phtrack =ntpTrack->ph.sigcor;
02662 Float_t phevent =ntpEvent->ph.sigcor;
02663
02664
02665
02666
02667 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear && trkfitpass==0) {
02668
02669 if (abs(trkplbu-trkplbv)<=5 && abs(trkpleu-trkplev)<=40 && (trkpleu<270&&trkplev<270)) nd_reclaim=1;
02670
02671 }
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683 if (phevent>0) {trkphfrac=phtrack/phevent;}
02684 if(ntpTrack->plane.n>0) {
02685 trkphplane = ntpTrack->ph.sigcor/ntpTrack->plane.n;
02686 }
02687
02688 is_cont_trk = ntpTrack->contained;
02689 trkmombest =ntpTrack->momentum.best;
02690 }
02691 else {
02692 trklength = 0;
02693 trkmomrange = 0;
02694 trkqp = 0;
02695 trkeqp = 0;
02696 trkphfrac = 0;
02697 trkphplane = 0;
02698 }
02699
02700 reco_enu = reco_emu + reco_eshwcc;
02701 reco_qe_enu = RecoQENuEnergy(track_index);
02702 if (reco_enu>0) {
02703 reco_y = reco_eshwcc/reco_enu;
02704 }
02705 reco_qe_q2 = RecoQEQ2(track_index);
02706 reco_dircosz = RecoMuDCosZ(track_index);
02707 reco_dircosneu = RecoMuDCosNeu(track_index);
02708 is_cev = IsFidAll(track_index);
02709
02710
02711
02712
02713
02714 if (detector==2 && ntrack>0 && reco_dircosneu>0.6 && is_fid_dp==1) pass_cccuts=1;
02715 if (detector==2 && ntrack>0 && reco_dircosneu>0.6 && is_fid_ab==1) pass_cccuts_2007=1;
02716 if (detector==2 && ntrack>0 && reco_dircosneu>0.6 && is_fid_2008==1) pass_cccuts_2008=1;
02717
02718
02719
02720
02721
02722 if(detector==2){
02723
02724 Bool_t general=false;
02725 Bool_t tr1=false;
02726 Bool_t tr2=false;
02727 Bool_t tr3=false;
02728 Bool_t tr4=false;
02729 Bool_t tr5=false;
02730 Bool_t shw1=false;
02731 Bool_t shw2=false;
02732 Bool_t shw3=false;
02733
02734 Float_t phtot = ntpEvent->ph.sigcor;
02735
02736 passcut=0;
02737
02738
02739
02740 if((ann_aPhperstp<2000.&& phtot>2000 && evtlength>2 && (evttimemax-evttimemin)*1e6>0.01 && fabs(litime-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9)>0.5 && nevent<=2)) general = true;
02741
02742 if(general){
02743 if(track_index!=-1){
02744
02745
02746
02747
02748 if(!((trkcosx*trkcosx+trkcosy*trkcosy)>0.6 && (trkvtxx*trkvtxx+trkvtxy*trkvtxy)>13)) tr1=true;
02749 if(!((trkcosx*trkcosx+trkcosy*trkcosy)>0.3 && (trkvtxx*trkvtxx+trkvtxy*trkvtxy)>15)) tr2=true;
02750 if(!((trkcosx*trkcosx+trkcosy*trkcosy)>0.9 && (trkvtxx*trkvtxx+trkvtxy*trkvtxy)>10)) tr3=true;
02751
02752
02753 if(trkplb>2 && trkplb!=250 && trkplb!=251 && fabs(trkcosy)<0.7 && fabs(trkvtxy-trkendy)<5.5&&fabs(trkcosz)>0.5) tr4=true;
02754
02755 if (trkcosx==0) {tr5=true;}
02756 else {
02757 if (fabs(trkcosx)*(trkcosy/trkcosx)<0.8) {tr5=true;}
02758 }
02759 if(tr1&&tr2&&tr3&&tr4&&tr5) passcut=1;
02760
02761 }
02762
02763 if(track_index==-1){
02764
02765 Double_t shwvx=ntpShower->vtx.x;
02766 Double_t shwvy=ntpShower->vtx.y;
02767 Double_t shwvz=ntpShower->vtx.z;
02768
02769 Double_t shwcosx=ntpShower->vtx.dcosx;
02770 Double_t shwcosy=ntpShower->vtx.dcosy;
02771 Double_t shwcosz=ntpShower->vtx.dcosz;
02772 Double_t shwplb =ntpShower->plane.beg;
02773 Double_t shwph =ntpShower->ph.sigcor;
02774
02775 Double_t eex=ntpEvent->end.x;
02776 Double_t eey=ntpEvent->end.y;
02777 Double_t evx=ntpEvent->vtx.x;
02778 Double_t evy=ntpEvent->vtx.y;
02779
02780 Double_t shwstpphtemp;
02781 Double_t shwphmax=-1;
02782 Double_t shwphmin=10e9;
02783 Double_t shphpl[500];
02784
02785 for(Int_t spl=0;spl<500;spl++) shphpl[spl]=0;
02786
02787 Double_t shwstpu=0;
02788 Double_t shwstpv=0;
02789 Double_t shwphmaxpl;
02790 Double_t shphplmax;
02791 Double_t shphplmaxv;
02792 Double_t shwphminpl;
02793
02794 for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
02795
02796 Int_t stp_indexs=((ntpShower->stp)[jj]);
02797
02798 if (stp_indexs==-1) continue;
02799
02800 LoadStrip(stp_indexs);
02801
02802 if((ntpStrip->plane)%2==1) shwstpu=shwstpu+1;
02803 if((ntpStrip->plane)%2==0) shwstpv=shwstpv+1;
02804
02805 shwstpphtemp=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
02806
02807 if(shwstpphtemp>shwphmax) {
02808 shwphmax =shwstpphtemp;
02809 shwphmaxpl =ntpStrip->plane-ntpShower->plane.beg;
02810 }
02811
02812 if(shwstpphtemp<shwphmin) {
02813 shwphmin =shwstpphtemp;
02814 shwphminpl =ntpStrip->plane-ntpShower->plane.beg;
02815 }
02816
02817 for(Int_t spl=ntpShower->plane.beg;spl<=ntpShower->plane.end;spl++){
02818 if(spl==ntpStrip->plane){
02819 shphpl[spl] =shphpl[spl]+shwstpphtemp;
02820 }
02821 }
02822 }
02823
02824 shphplmax=-1;
02825 shphplmaxv=-1;
02826
02827 for(Int_t spl=ntpShower->plane.beg;spl<=ntpShower->plane.end;spl++){
02828
02829 if(shphpl[spl]>shphplmaxv){
02830 shphplmax =spl;
02831 shphplmaxv =shphpl[spl];
02832 }
02833
02834 }
02835
02836 shphplmax=shphplmax-ntpShower->plane.beg;
02837
02838 Double_t diff1=evx-eex;
02839 Double_t diff2=evy-eey;
02840 Double_t diff3=shwstpu-shwstpv;
02841
02842 if((shwvx*shwvx+shwvy*shwvy)<19. && shwvz>0.2 && shwvz<29.5 && fabs(diff1)<2. && fabs(diff2)<2. && fabs(shwvx)<3.5 &&fabs(shwvy)<3.5) shw1= true;
02843
02844 if((shwcosx*shwcosx+shwcosy*shwcosy)<0.8 && shwcosz>0.6 && shwplb<480) shw2=true;
02845
02846 if(fabs(diff3)<30. && (shphplmaxv/shwph)<0.9 && (shwphmax/shwph)<0.8) shw3=true;
02847
02848 if(shw1&&shw2&&shw3) passcut=1;
02849
02850 }
02851 }
02852 }
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863 if(ntrack){
02864
02865
02866
02867
02868 pid0 = PID(event,0);
02869
02870
02871
02872
02873
02874 pid1=-1.;
02875
02876
02877
02878
02879 pid1 = andy_id.CalcPID(ntpEvent,strecord);
02880
02881 abid_costheta =andy_id.GetRecoCosTheta(ntpEvent,strecord);
02882 abid_eventenergy =andy_id.GetRecoE(ntpEvent,strecord);
02883 abid_trackcharge =andy_id.GetTrackEMCharge(ntpEvent,strecord);
02884 abid_trackenergy =andy_id.GetTrackPlanes(ntpEvent,strecord);
02885 abid_tracklikeplanes =andy_id.GetTrackLikePlanes(ntpEvent,strecord);
02886 abid_trackphfrac =andy_id.GetTrackPHfrac(ntpEvent,strecord);
02887 abid_trackphmean =andy_id.GetTrackPHmean(ntpEvent,strecord);
02888 abid_trackqpsigmaqp =andy_id.GetTrackQPsigmaQP(ntpEvent,strecord);
02889 abid_y =andy_id.GetRecoY(ntpEvent,strecord);
02890
02891
02892
02893
02894
02895
02896 if (ro_pass) {
02897
02898
02899 pid2=ro_interface.GetVar("knn_pid",ntpEvent);
02900
02901 roid_numplanes = ro_interface.GetVar("knn_01",ntpEvent);
02902 roid_meansigcor = ro_interface.GetVar("knn_10",ntpEvent);
02903 roid_lowphhiph = ro_interface.GetVar("knn_20",ntpEvent);
02904 roid_trkphfracclose = ro_interface.GetVar("knn_40",ntpEvent);
02905
02906
02907 }
02908
02909
02910
02911 if(PassAnalysisCuts(event)) pass = 1;
02912 }
02913 else{
02914 pid0 = -999.;
02915 pid1 = -999.;
02916 pid2 = -999.;
02917 pass = 0;
02918 }
02919 cout << "\t\tFilling tree\n";
02920 tree->Fill();
02921
02922 }
02923 }
02924 }
02925
02926
02927 delete spinfo;
02928
02929 potcount->Fill();
02930
02931 save.cd();
02932 tree->Write();
02933 potcount->Write();
02934 save.Write();
02935 save.Close();
02936
02937 }
02938
02939
02940
02941 void MadPIDAnalysis::ReadPIDFile(std::string tag){
02942
02943 fLikeliFile = new TFile(tag.c_str(),"READ");
02944
02945 if(fLikeliFile->IsOpen() && !fLikeliFile->IsZombie()) {
02946
02947 fLikeliFile->cd();
02948
02949 fLikeliHist[0] = (TH1F*) fLikeliFile->Get("Event length - cc");
02950 fLikeliHist[1] = (TH1F*) fLikeliFile->Get("Event length - nc");
02951 fLikeliHist[2] = (TH1F*) fLikeliFile->Get("Track ph frac - cc");
02952 fLikeliHist[3] = (TH1F*) fLikeliFile->Get("Track ph frac - nc");
02953 fLikeliHist[4] = (TH1F*) fLikeliFile->Get("ph per plane (cc)");
02954 fLikeliHist[5] = (TH1F*) fLikeliFile->Get("ph per plane (nc)");
02955
02956 for(int i=0;i<6;i++){
02957 float integ = fLikeliHist[i]->Integral();
02958 fLikeliHist[i]->Scale(1./integ);
02959 }
02960 }
02961 else fLikeliFile = NULL;
02962 }
02963
02964
02965
02966 void MadPIDAnalysis::ReadAbPIDFile(std::string tag){
02967
02968
02969 cout << "Reading AbID pdfs from " << tag.c_str() << endl;
02970 andy_id.ReadPDFs(tag.c_str());
02971
02972 }
02973
02974
02975
02976 void MadPIDAnalysis::ConfigureRoID(std::string keyfile, std::string tag){
02977
02978
02979
02980 ro_registry.UnLockValues();
02981
02982
02983 ro_registry.Set("InterfaceConfigPath", keyfile.c_str());
02984 ro_registry.Set("FillkNNFilePath", tag.c_str());
02985
02986
02987
02988
02989 ro_interface.Config(ro_registry);
02990
02991 }
02992
02993
02994
02995
02996
02997 bool MadPIDAnalysis::InFidVol(const Detector::Detector_t& det,
02998 const Float_t& x, const Float_t& y,
02999 const Float_t& z){
03000
03001
03002
03003
03004 bool is_fid=false;
03005
03006 if(det==Detector::kFar){
03007
03008
03009 is_fid = true;
03010 if(z<0.5 || z>29.4 ||
03011 (z<16.5&&z>14.5) ||
03012 sqrt((x*x)
03013 +(y*y))>3.5) is_fid = false;
03014 }
03015 else if(det==Detector::kNear){
03016
03017
03018 is_fid = true;
03019 if(z<1 || z>5 ||
03020 sqrt(((x-1.4828)*(x-1.4828)) +
03021 ((y-0.2384)*(y-0.2384)))>1) is_fid = false;
03022 }
03023
03024 return is_fid;
03025
03026 }
03027
03028
03029 void MadPIDAnalysis::MakeAbIDFile(std::string tag){
03030
03031
03032
03033
03034 std::string savename = "ABHistos_" + tag + ".root";
03035
03036
03037 for(int i=0;i<Nentries;i++){
03038
03039
03040 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
03041 << "% done" << std::endl;
03042 if(!this->GetEntry(i)) continue;
03043
03044 Int_t nevent = eventSummary->nevent;
03045 if(nevent==0) continue;
03046
03047
03048 Int_t evtmin=0;
03049 Int_t evtmax=nevent;
03050
03051
03052
03053
03054 if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) {
03055
03056 Int_t maxevent=0;
03057 Float_t maxph=0;
03058 for(int j=0;j<nevent;j++){
03059 if(!LoadEvent(j)) continue;
03060 Float_t evtpheight=ntpEvent->ph.sigcor;
03061 if (evtpheight>maxph) {
03062 maxph=evtpheight;
03063 maxevent=j;
03064 }
03065 }
03066 evtmin=maxevent;
03067 evtmax=maxevent+1;
03068 }
03069
03070
03071
03072 for(int j=evtmin;j<evtmax;j++){
03073
03074
03075 if(!LoadEvent(j)) continue;
03076
03077
03078 Int_t ccnc = -1;
03079
03080
03081
03082 Int_t mcevent=-1;
03083
03084 if(LoadTruthForRecoTH(j,mcevent)) {
03085
03086
03087 Int_t iact = IAction(mcevent);
03088 Int_t inu = INu(mcevent);
03089
03090
03091 if(abs(inu)==14 && iact==1) ccnc=1;
03092 else ccnc=0;
03093
03094 }
03095
03096
03097 if (andy_id.PassCuts(ntpEvent,strecord)) {
03098 andy_id.FillPDFs(ntpEvent,strecord,ccnc);
03099 }
03100
03101 }
03102
03103 }
03104
03105
03106
03107
03108 andy_id.NormalizePDFs();
03109 andy_id.WritePDFs(savename.c_str());
03110
03111
03112 }
03113
03114
03115
03116
03117 #endif // #ifdef madpidanalysis_cxx
