loon: MadDpAnalysis.cxx Source File

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   //  if(PassBasicCuts(event)) return true;
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   // for events with no track, use event vertex
00107 
00108   if (track==-1 && (ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>28.0 ||   //ends
00109            (ntpEvent->vtx.z<16.2 && ntpEvent->vtx.z>14.3) ||  //between SMs
00110            (pow(ntpEvent->vtx.x,2)+           //radial cut
00111             pow(ntpEvent->vtx.y,2))>14)) return false;
00112 
00113   // otherwise, use track vertex
00114 
00115   if (track!=-1 && (ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>28.0 ||   //ends
00116            (ntpTrack->vtx.z<16.2 && ntpTrack->vtx.z>14.3) ||  //between SMs
00117            (pow(ntpTrack->vtx.x,2)+           //radial cut
00118             pow(ntpTrack->vtx.y,2))>14)) return false;  
00119 
00120   return true; 
00121 
00122 }
00123 
00124  // new fiducial region - used by DP FD event analysis
00125 
00126 Bool_t MadDpAnalysis::PassFidNewN(Int_t event)  // NEW !!
00127 { 
00128 
00129   if(!LoadEvent(event)) return false;
00130 
00131   Int_t track = -1;
00132   LoadLargestTrackFromEvent(event,track);
00133 
00134   // for events with no track, use event vertex
00135 
00136   if (track==-1 && (ntpEvent->vtx.z<1.0 || ntpEvent->vtx.z>27.0 ||   //ends
00137            (ntpEvent->vtx.z<17.0 && ntpEvent->vtx.z>14.0) ||  //between SMs
00138            (pow(ntpEvent->vtx.x,2)+           //radial cut
00139             pow(ntpEvent->vtx.y,2))>(3.5*3.5))) return false;
00140 
00141   // otherwise, use track vertex
00142 
00143   if (track!=-1 && (ntpTrack->vtx.z<1.0 || ntpTrack->vtx.z>27.0 ||   //ends
00144            (ntpTrack->vtx.z<17. && ntpTrack->vtx.z>14.0) ||  //between SMs
00145            (pow(ntpTrack->vtx.x,2)+           //radial cut
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    //   Int_t track = -1;
00157    // LoadLargestTrackFromEvent(event,track);  
00158    // if(!PassTrackCuts(track)) return false;  
00159    //   if(!IsFidVtxEvt(ntpEvent->index)) return false;  
00160 //  if(!IsFidVtx(track)) return false;
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 // ANN VARIABLES
00171 //
00172 AnnInputBlock MadDpAnalysis::AnnVar(Int_t event) 
00173 {
00174 AnnInputBlock anninput;
00175   // Initialize
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 // Calculate variables needed for ANN (and a few more)
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 //  LoadLargestShowerFromEvent(event,shower_index);
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; // prepei na to alla3w
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 // loop over strips to compute ph3 ph6 phlast phcommon
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  // TIMING INFO // 
00271 //  Find max min time of events by loading strips for ND 
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        // shower strips
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        // tracks strips
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       } // end if strip != -1 (!!!)
00363      } // end of looping over event strips....    
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    // DET IS THE DETECTOR 1 = NEAR 2 = FAR
00397    // TAR IS THE ENERGY   1 = LE 2 = PME 3 = PHE
00398    // FID IS THE FIDUCIAL REGION FOR FAR , 2 = DAVIDS  1 = MINE (MORE STRICT)
00399    // PRIOR IS THE FAR TRAINED METHOD    , 1 = NO OSCILLATIONS , 2 = DM2 0.0025 SINTHETA2 = 0.95 
00400    //                                      3 = DM2 0.002 SINTHETA2 = 0.95
00401       
00402       
00403 //      int    inneuron  = 13;
00404       int    inneuron  = 7;
00405       int    hidneuron = 15;
00406       double var;
00407      
00408       double out[15]; // input neurons
00409       double rin[15]; // first hidden layer
00410      
00411 //      double weight1[13][15];// weights first layer
00412       double weight1[7][15];
00413       double constant1[15];   // constant first layer
00414       
00415       double weighto[15];// weights second  (output) layer
00416       double constanto[1];   // constant second (output) layer  
00417       double prob;
00418       
00419 // Initialize    
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 //   INPUT VARIABLES
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    // out[1]  =  (anninput.aTotph/50000.); 
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       out[0]  = anninput.aTrkphperpl/2000.;
00486       out[1]  = anninput.aTotrk*anninput.aTrkpass;
00487       out[2]  = (anninput.aTotstp/100.);
00488       out[3]  = (anninput.aTotph/50000.);   
00489       out[4]  = (anninput.aTotlen/40);
00490       out[5]  = (anninput.aPhperpl/5000.);      
00491       out[6]  = (anninput.aPhperstp/1000.);
00492       out[10] = ((anninput.aTrkplv-anninput.aShwplv)+20.)/40.;
00493       out[11] = ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.;
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      //out[1]  =  (anninput.aTotph/100000.); 
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        out[0]  = anninput.aTrkphperpl/5000.;
00510        out[1]  = anninput.aTotrk*anninput.aTrkpass;
00511        out[2]  = (anninput.aTotstp/100.);
00512        out[3]  = (anninput.aTotph/100000.);   
00513        out[4]  = (anninput.aTotlen/40);
00514        out[5]  = (anninput.aPhperpl/10000.);    
00515        out[6]  = (anninput.aPhperstp/2000.);  
00516        out[10] = ((anninput.aTrkplv-anninput.aShwplv)+10.)/20.;
00517        out[11] = ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.;
00518       */
00519       }
00520       /*
00521       out[7]  = (anninput.aPh3/(anninput.aTotph+1.));
00522       out[8]  = (anninput.aPh6/(anninput.aTotph+1.));
00523       out[9]  = (anninput.aPhlast/(anninput.aTotph+1.));
00524       out[12] = (anninput.aShwphperdig/800.);
00525       */
00526 // INPUT FILE
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 // end of reading the weights
00556    
00557        // first layer         
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        // output                        
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 //SP likelihood analysis
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   /* In case of Birch
00605   // apply 1.8% correction to pulse height variables
00606   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar && 
00607      ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kData  ){
00608     phtrack=phtrack*1.018;
00609     phevent=phevent*1.018;
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   // pdf binning: 0 - old binning, max 50 bins, 1 - new finer binning
00661 
00662   std::string savename = "DPHistos_" + tag + ".root";
00663   TFile save(savename.c_str(),"RECREATE"); 
00664   save.cd();
00665   
00666   
00667 
00668 
00669   //#declare lots of histos:
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   //  gDirectory->ls();
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       // EVENT LOOP - fill tree once for each reconstructed event
00769 
00770       for(int j=evtmin;j<evtmax;j++){       
00771 
00772 
00773       if(!LoadEvent(j)) continue; //no event found
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       //get detector type:
00812       if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00813         detector = 2;
00814         //fiducial criteria on vtx for far detector
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 ||   //ends
00827            (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||  //between SMs
00828            sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
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         //fiducial criteria on vtx for near detector
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       //check for a corresponding mc event
00852  
00853       Int_t mcevent=-1;
00854 
00855       if(LoadTruthForRecoTH(j,mcevent)) {
00856         //true info for tree:
00857 
00858         cc_nc             = IAction(mcevent);
00859         flavour           = Flavour(mcevent);
00860         neu_e             = TrueNuEnergy(mcevent);  
00861                           
00862    // tarpos   1 = LE 2 = ME 3 = HE     
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         //      if (is_fid && trkpass) {
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           //      cout<<"I've filled a histogram, where it is, I don't know"<<endl;
00897           //      cout<<phplane_cc->GetEntries()<<endl;
00898           //      gDirectory->ls();
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         //      std::cout<<"evlenght "<<evlength<<" weight "<<weight<<std::endl;
00907 
00908         }
00909       }
00910     }
00911   }
00912   //#close file  
00913   //  std::cout<<"writing hist "<<endl;
00914   //  save.ls();
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 // 1 = le-10 2 = pme 3 = phe there is no other way in an MC file to now what target positions it
00947 // was...!
00948   
00949   std::string savename = "PAN_" + tag + ".root";
00950   TFile save(savename.c_str(),"RECREATE"); 
00951   save.cd();
00952   /*
00953   GnumiInterface gnumi;
00954   if(!gnumi.Status()) {
00955     cout << "Error MadAnalysis::CreatePAN() - Flux files not open."
00956          << " Will not be filling NuParent info"
00957          << endl;
00958     cout << "Set environmental variable $GNUMIAUX to point to the "
00959          << "directory containing the gnumi files to fix this!"
00960          << endl;
00961 
00962   }
00963 */
00964  SpillInfo pppinfo;
00965  
00966   
00967   //ann quantities
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   //PAN Quantities 
01002   //Truth:
01003   Float_t true_enu;       // true neutrino energy (GeV)
01004   Float_t true_pxnu;      // true neutrino momentum-x (GeV)
01005   Float_t true_pynu;      // true neutrino momentum-y (GeV)
01006   Float_t true_pznu;      // true neutrino momentum-z (GeV)
01007   Float_t true_etgt;       // true target energy (GeV)
01008   Float_t true_pxtgt;      // true target momentum-x (GeV)
01009   Float_t true_pytgt;      // true target momentum-y (GeV)
01010   Float_t true_pztgt;      // true target momentum-z (GeV)
01011   Float_t true_emu;       // true muon energy
01012   Float_t true_eshw;      // true shower energy
01013   Float_t true_x;         // true x
01014   Float_t true_y;         // true y
01015   Float_t true_q2;        // true q2
01016   Float_t true_w2;        // true w2
01017   Float_t true_dircosneu; // true muon dircos w.r.t neutrino
01018   Float_t true_dircosz;   // track z-direction cosine
01019   Float_t true_vtxx;      // true x vtx
01020   Float_t true_vtxy;      // true y vtx
01021   Float_t true_vtxz;      // true z vtx
01022 
01023   //For Neugen:
01024   Int_t flavour;          // true flavour: 1-e 2-mu 3-tau
01025   Int_t process;          // process: 1001-QEL 1002-RES 1003-DIS 1004-COH
01026   Int_t nucleus;          // target nucleus: 274-C 284-O 372-Fe
01027   Int_t cc_nc;            // cc/nc flag: 1-cc 2-nc
01028   Int_t initial_state;    // initial state: 1-vp 2-vn 3-vbarp 4-vbarn ...
01029   Int_t had_fs;           // hadronic final state, number between 200-300
01030   
01031   //For Beam Reweighting:
01032  // NuParent *nuparent = new NuParent();
01033 
01034   //Reco Quantities
01035   Int_t detector;         // Near = 1, Far = 2;
01036   Int_t run;              // run number
01037   Int_t snarl;            // snarl number
01038   Int_t nevent;           // number of events in snarl
01039   Int_t event;            // event index
01040   Int_t subrun;           // subrun number     
01041   Int_t trigsrc;          // trigger source    
01042   Int_t mcevent;          // mc event index
01043   Int_t ntrack;           // number of tracks in event
01044   Int_t nshower;          // number of showers in event
01045 
01046   Int_t is_fid;           // pass fid cut. true = 1, false = 0
01047   
01048   Int_t is_fid_dp;        // pass dp analysis fid cut. true = 1, false = 0
01049   Int_t is_fid_ns;        // pass  nd analysis fid cut true=1, false=0 // !! NEW
01050   
01051   Int_t is_cev;           // fully contained. true = 1, false = 0 
01052   Int_t is_cont_trk;      // new containemnt support by IsContained() for the largest track
01053   Int_t is_mc;            // is a corresponding mc event found
01054   Int_t pass_basic;       // Pass basic cuts. true = 1, false = 0
01055   Float_t pid0;           // pid parameter (usual method)
01056   Float_t pid1;           // pid parameter (alternative method 1)
01057   Float_t pid2;           // pid parameter (alternative method 2)
01058   Float_t annpid;         // ANN PID
01059   Float_t annpid_f1p1;    // ANN PID FAR FIDUCIAL DP PRIOR 1 // UN-OSICLLATED 
01060   Float_t annpid_f1p2;    // ANN PID FAR FIDUCIAL DP PRIOR 2 // OSCILLATED WITH DM2 0.0025 SIN2THETA 0.95
01061   Float_t annpid_f1p3;    // ANN PID FAR FIDUCIAL DP PRIOR 3 // OSCILLATED WITH DM2 0.0020 SIN2THETA 0.95
01062   Float_t annpid_f2p1;    // ANN PID FAR FIDUCIAL NS PRIOR 1  UN-OSICLLATED 
01063   Float_t annpid_f2p2;    // ANN PID FAR FIDUCIAL NS PRIOR 2 // OSCILLATED WITH DM2 0.0025 SIN2THETA 0.95
01064   Float_t annpid_f2p3;    // ANN PID FAR FIDUCIAL NS PRIOR 3 // OSCILLATED WITH DM2 0.0020 SIN2THETA 0.95
01065   
01066   Int_t pass;             // pass analysis cuts. true = 1, false = 0
01067 
01068   Float_t reco_enu;       // reco neutrino energy
01069   Float_t reco_emu;       // reco muon energy
01070   Float_t reco_eshw;      // reco shower energy  (shw.ph.gev/1.23)
01071   
01072   Int_t pass_fd_qualcuts;  // pass FD quality cuts (HV trips etc)
01073   Int_t litag;  
01074 
01075 
01076   Float_t  reco_eshwcc;    // linear cc version
01077   Float_t  reco_eshwnc;    // linear nc version 
01078   Float_t  reco_eshwccw;   // weighted cc version 
01079   Float_t  reco_eshwncw;   // weighted nc version 
01080   
01081   Float_t reco_eshw_sqrt; // reco shower energy using sqrt method
01082   Float_t reco_qe_enu;    // reco qe neutrino energy
01083   Float_t reco_qe_q2;     // reco qe q2
01084   Float_t reco_y;         // reco y
01085   Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino
01086   Float_t reco_dircosz;   // reco track vtx z-dircos
01087   Float_t reco_vtxx;      // reco x vtx
01088   Float_t reco_vtxy;      // reco y vtx
01089   Float_t reco_vtxz;      // reco z vtx
01090 
01091   Float_t evtlength;      // reco event length
01092   Float_t evtph;          // reco event ph
01093   Float_t trklength;      // reco track length
01094   Float_t trkmomrange;    // reco track momentum from range
01095   Float_t trkqp;          // reco track q/p
01096   Float_t trkeqp;         // reco track q/p error
01097   Float_t trkmombest;     // reco best track momentum based on new containment support
01098   Float_t trkphfrac;      // reco pulse height fraction in track
01099   Float_t trkphplane;     // reco average track pulse height/plane
01100   Float_t trkds;          // track DS 
01101   Float_t rawph;          // Raw ph
01102 // Additional track info
01103   Float_t trkendz;  // track end Z position
01104   Float_t trkendx;  // track end X position
01105   Float_t trkendy;  // track end Y position
01106   Float_t trkendu;  // track end U position
01107   Float_t trkendv;  // track end V position
01108   Float_t trkvtxz;  // track begin Z position
01109   Float_t trkvtxx;  // track begin X position
01110   Float_t trkvtxy;  // track begin Y position
01111   Float_t trkvtxu;  // track begin U position
01112   Float_t trkvtxv;  // track begin V position
01113   Float_t trkplbu;  // track begin plane u 
01114   Float_t trkplbv;  // track begin plane v 
01115   Float_t trkpleu;  // track end   plane u 
01116   Float_t trkplev;  // track end   plane v 
01117   Float_t trkple;   // track end   plane  
01118   Float_t trkplb;   // track begin   plane 
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    // shower variables
01125   Float_t shwph;
01126   Int_t   shwdigits, shwstrips;
01127   Float_t shwvtxz;  // shower begin Z position
01128   Float_t shwvtxx;  // shower begin X position
01129   Float_t shwvtxy;  // shower begin Y position
01130   Float_t shwvtxu;  // shower begin U position
01131   Float_t shwvtxv;  // shower begin V position
01132  
01133   // cluster variables
01134 
01135    Int_t shwncluster;
01136   
01137   // SCAN DECISION   // <- define an instance of ScanList class here
01138   Int_t scandecision;
01139   ScanList scan;
01140   
01141 // EVENT TIMING
01142    Double_t evttimemin;  //MIN STRIP TIME OF EVENT
01143    Double_t evttimemax;  //MAX STRIP TIME OF EVENT
01144    Double_t snarlt;
01145    Double_t litime; 
01146  
01147 // BEAM INFO
01148    Float_t  pot,horn,tar,vvpos,hhpos,magn;    
01149 // Beam INFO DATABASE
01150    Float_t potdb,horndb,tardb,vvposdb,hhposdb,vvwidthdb,hhwidthdb;
01151    Double_t timedb;
01152 // TIME INFO DATABASE
01153    Double_t neartdb,fartdb,neardifftdb;
01154 
01155 // mc flux info for xf,pt weighting
01156    Float_t mcflux_tpx;
01157    Float_t mcflux_tpy;
01158    Float_t mcflux_tpz;
01159    Float_t mcflux_tptype;
01160 
01161 // IS SNARL CUT IN PIECES 
01162    Float_t singletimeframe;
01163 // Beam Info Block    
01164    SpillInfoBlock *spinfo = new SpillInfoBlock();  
01165    spinfo->Zero(); 
01166    singletimeframe=-999;
01167  //
01168   Int_t day,month,year;
01169   Int_t pass_beamcuts;     // pass default beam quality cuts   
01170   //Trees
01171   TTree *tree = new TTree("pan","pan");
01172   //Truth
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   //Neugen
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   //NuParent
01200   //tree->Branch("nuparent","NuParent",&nuparent,8000,1);
01201   //Reco
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   // SCAN DECISIONS 
01309   tree->Branch("scandecision",&scandecision,"scandecision/I");
01310 
01311   //BEAM INFO 
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   // BEAM INFO DB
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   // TIME DB  
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  // mc flux info
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 // ann vars
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     // SCAN INFO  
01424     
01425     scandecision=scan.GetDecision(run,snarl);
01426 
01427     // only pass through scanned snarls if scanfilter=1
01428     if (scanfilter==0 || scandecision>0) {  
01429 
01430       // Get Beam Monitorint Info if NOT MC: 
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         // BEAM INFO DB
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           // apply beam quality cuts
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         // TIME DB
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       //singletimeframe=SingleTimeFrame(i,Nentries);
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         //      cout << "MULTIPLE events!!\n";
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            //  cout << "evtlength=" << ntpEvent->plane.n << " totph=" 
01519            //   << ntpEvent->ph.sigcor <<  "ntrack=" << ntpEvent->ntrack
01520            //   << " nshower=" << ntpEvent->nshower << endl;
01521            Float_t evtpheight=ntpEvent->ph.sigcor;
01522            if (evtpheight>maxph) {
01523              maxph=evtpheight;
01524              maxevent=j;
01525            }
01526          }
01527          //      cout << "BIGGEST event=" << maxevent << endl;
01528          evtmin=maxevent;
01529          evtmax=maxevent+1;
01530       }
01531 
01532 
01533       // EVENT LOOP - fill tree once for each reconstructed event
01534 
01535       for(int j=evtmin;j<evtmax;j++){ 
01536 
01537         if(!LoadEvent(j)) continue; //no event found
01538 
01539         totph = 0;
01540 
01541         //set event index:
01542         event = j;
01543         ntrack = ntpEvent->ntrack;
01544         nshower = ntpEvent->nshower;
01545         totdigits=ntpEvent->ndigit;
01546         totstrips=ntpEvent->nstrip;
01547                 
01548         // Initialize everything
01549         //nuparent->Zero();  
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         //get detector type:
01626         if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
01627           detector = 2;
01628           totph=eventSummary->ph.sigcor;
01629           //fiducial criteria on vtx for far detector
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           //get total charge in trk/shw
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           //fiducial criteria on vtx for near detector
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         //check for a corresponding mc event     
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           //true info for tree:
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           if(gnumi.Status()) {
01696             if(isST) gnumi.GetParent(strecord,mcevent,*nuparent);
01697             else     gnumi.GetParent(mcrecord,mcevent,*nuparent);
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;} // make this reflect pass/fail flag
01714         
01715         //reco info for tree:
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         // different definition for neardet - accounts for 
01723         // uninstrumented planes
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         // unused:      Double_t trgtime=eventSummary->trigtime;
01760         //      evttimemin = anninput.aTimemin-trgtime;
01761         //      evttimemax = anninput.aTimemax-trgtime;
01762  
01763         // cedar - use snarl time rather than trigger time
01764         evttimemin = anninput.aTimemin-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9;
01765         evttimemax = anninput.aTimemax-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9;
01766         
01767 
01768         // ANN PID // 
01769         
01770      
01771         Int_t target=1;
01772         if(ntpHeader->GetVldContext().GetSimFlag()!=4){
01773           //if(tar<5000.  && tar>3500.)    target=1;
01774           //if(tar>5000.  && tar<40000.)   target=2;
01775           //if(tar>40000. && tar<100000.)  target=3; // Used to work just fine for the pre shutdown data
01776           // after the shutdown this variable was not properly filled from ACNET data. in fact it is not
01777           // recoverable since it is corrupted in ACNET data.
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; // FOR THE MOMENT IN THE ABSENCE OF ANY FAR PME PHE MC 
01789         }                       
01790         // SET FIDUCIAL AND TARGET DEFAULT IS DAVIDS FIDUCIAL AND NO OSCILLATION TRAINING
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          // MAK March 8, 2005
01835          // data or MC?
01836          // needed by energy corrections below
01837          bool isdata=true;
01838          if(ntpHeader->GetVldContext().GetSimFlag()==4) isdata=false;
01839          
01840         //index will -1 if no track/shower present in event
01841         int track_index   = -1;
01842         LoadLargestTrackFromEvent(j,track_index);
01843         int shower_index  = -1;
01844         //LoadLargestShowerFromEvent(j,shower_index); // Change
01845         LoadShower_Jim(j,shower_index,detector);
01846         if(shower_index==-1) nshower=0;
01847         if(shower_index!=-1){ //check shower is present before using ntpShower
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            // MAK March 8, 2005
01860            Int_t shwtype=0; 
01861            const Int_t shw_correction_mode=1; // Niki's re-tuning
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){ //check track is present before using ntpTrack
01889            // MAK March 8, 2005    
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           /* In case of Birch
01926           // apply 1.8% correction to pulse height variables
01927           if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar && 
01928             ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kData  ){
01929             phtrack=phtrack*1.018;
01930             phevent=phevent*1.018;
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         // only calculate likelihood PID for events that pass event 
01962         // fiducial cuts and track quality cuts 
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       } // end of EVENT LOOP
01979     } //  end of SCANINFO loop 
01980   } // end of NTUPLE ENTRIES loop
01981 
01982  // delete nuparent;
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()) { //if file is found
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   // Mike Kordosky: 11 August 2005
02021   // Defines a fiducial volume used to pass events which fill the PDFs
02022   // Transcribed from MakeMyFile() above
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 ||   //ends
02031        (z<16.5&&z>14.5) ||  //between SMs
02032        sqrt((x*x)           //radial cut
02033             +(y*y))>3.5) is_fid = false;
02034   }
02035   else if(det==Detector::kNear){
02036     
02037     //fiducial criteria on vtx for near detector
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