loon: MNtpModule.cxx Source File

00001 #include "TGraph.h"
00002 #include "TMath.h"
00003 //#include "TClonesArray.h"
00004 
00005 #include "MNtpModule.h"
00006 #include "MNtp.h"
00007 #include "MBSpillAccessor.h"
00008 #include "MBSpill.h"
00009 
00010 #include "RecoBase/PropagationVelocity.h"
00011 #include "Record/RecRecordImp.h"
00012 #include "Record/RecCandHeader.h"
00013 
00014 #include "DataUtil/EnergyCorrections.h"
00015 using namespace EnergyCorrections;
00016 
00017 #include "Validity/VldContext.h"
00018 #include "Validity/VldTimeStamp.h"
00019 
00020 #include "UgliGeometry/UgliGeomHandle.h"
00021 #include "UgliGeometry/UgliStripHandle.h"
00022 
00023 #include "Conventions/PlaneView.h"
00024 #include "Conventions/ReleaseType.h"
00025 
00026 #include "StandardNtuple/NtpStRecord.h"
00027 #include "CandNtupleSR/NtpSRStrip.h"
00028 #include "CandNtupleSR/NtpSREvent.h"
00029 #include "CandNtupleSR/NtpSRTrack.h"
00030 #include "CandNtupleSR/NtpSRShower.h"
00031 
00032 #include "MCNtuple/NtpMCTruth.h"
00033 #include "MCNtuple/NtpMCStdHep.h"
00034 #include "TruthHelperNtuple/NtpTHTrack.h"
00035 
00036 #include "AstroUtil/AstCoordinate.h"
00037 
00038 #include "NueAna/NueAnaTools/SntpHelpers.h"
00039 
00040 #include "MessageService/MsgService.h"
00041 #include "MinosObjectMap/MomNavigator.h"
00042 #include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
00043 
00044 #include "VertexFinder/NtpVtxFinder/NtpVtxFinder.h"
00045 
00046 JOBMODULE(MNtpModule, "MNtpModule",
00047           "a job module to fill the information for MiniBooNE neutrino study");
00048 CVSID("$Id: MNtpModule.cxx,v 1.21 2008/10/04 23:59:24 tjyang Exp $");
00049 //......................................................................
00050 
00051 //zenith and azimuth angles of the expected MB neutrino direction
00052 const float theta_mb = 83.5;
00053 const float phi_mb = 172.8;
00054 
00055 MNtpModule::MNtpModule()
00056   : fFile("mntp.root")
00057   , fN(0)
00058   , fM(0)
00059   , IsTriggerOn(false)
00060   , numfiles(0)
00061 {
00062   fMNtp = new MNtp();
00063   tfit_dt_ds_pos = new TF1("tfit_dt_ds_pos","(1./299792458.0)*x+[0]",0,30) ;
00064   tfit_dt_ds_neg = new TF1("tfit_dt_ds_neg","(-1./299792458.0)*x+[0]",0,30) ;
00065 }
00066 //......................................................................
00067 
00068 MNtpModule::~MNtpModule()
00069 {
00070 }
00071 
00072 //......................................................................
00073 
00074 void MNtpModule::BeginJob()
00075 {
00076   fNtpFile = new TFile(fFile.c_str(),"RECREATE");  
00077   fNtuple = new TTree("mini","Analysis Tree");
00078   fNtuple->Branch("mntp","MNtp",&fMNtp,64000,2);
00079 
00080 
00081 }
00082 
00083 //......................................................................
00084 
00085 void MNtpModule::EndJob()
00086 {
00087   fNtuple->Write();
00088   fNtpFile->Close();
00089 }
00090 
00091 //......................................................................
00092 
00093 void MNtpModule::BeginFile()
00094 {
00095   numfiles++;
00096   cout <<"=====================================================" <<endl;
00097   cout << "This is the NO. " << numfiles << " ntuples." << endl;  
00098 }
00099 
00100 //......................................................................
00101 
00102 JobCResult MNtpModule::Ana(const MomNavigator* mom)
00103 {
00104   //reset everything at the beginning of each snarl
00105   fMNtp->ResetAll();
00106 
00107   //get run no,subrun no and snarl no
00108   RecRecordImp<RecCandHeader> *rr = 
00109     dynamic_cast<RecRecordImp<RecCandHeader>*>
00110     (mom->GetFragment("RecRecordImp<RecCandHeader>"));
00111   if (rr) {
00112     fMNtp->run    = rr->GetHeader().GetRun();
00113     fMNtp->subrun = rr->GetHeader().GetSubRun();
00114     fMNtp->snarl  = rr->GetHeader().GetSnarl();
00115     fDetectorType = rr->GetHeader().GetVldContext().GetDetector();
00116     fSimFlag = rr->GetHeader().GetVldContext().GetSimFlag();
00117   }
00118   else {
00119     MSG("MNtpModule", Msg::kWarning) << "Could not get RecCandHeader from MOM" << endl;
00120   }
00121 
00122 
00123   //get NtpStRecord
00124   NtpStRecord *record = dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord"));
00125   if (!record) {
00126     MSG("MNtpModule", Msg::kError) << "Could not get NtpStRecord from MOM" << endl;
00127     return JobCResult::kFailed;
00128   }
00129 
00130   //get no. of events
00131   int nevt = record->evthdr.nevent;
00132   if (!nevt) {
00133     MSG("MNtpModule", Msg::kDebug) << "There is no event in this snarl"<<endl;
00134     return JobCResult::kPassed;
00135   }
00136 
00137   fMNtp->nevt = nevt;
00138 
00139   //get VldContext, UgliGeomHandle and VldTimeStamp
00140   VldContext vldc = *(record->GetVldContext());
00141   VldTimeStamp vldts = vldc.GetTimeStamp();
00142   //UgliGeomHandle ugh(vldc);
00143 
00144   fDetector = record->GetHeader().GetVldContext().GetDetector();
00145   fBeam = BeamType::kLE ;
00146 
00147   fMNtp->sec = (int)vldts.GetSec();
00148   fMNtp->nsec = (int)vldts.GetNanoSec();
00149 
00150   fMNtp->day = (int)record->evthdr.date.day;
00151   fMNtp->month = (int)record->evthdr.date.month;
00152 
00153   //get trigger source:
00154   fMNtp->trgsrc = record->GetHeader().GetTrigSrc();
00155 
00156   //get run type:
00157   fMNtp->runtype = record->GetHeader().GetRunType();
00158 
00159   string relName = record->GetTitle();
00160   string mcinfo = "";
00161   if (fSimFlag==SimFlag::kMC){
00162     mcinfo = "Carrot";
00163     string temp = record->mchdr.geninfo.codename;
00164     if (temp.size() != 0) mcinfo = temp;
00165     if (mcinfo == "development") mcinfo = "daikon_04";
00166   }
00167   int release = ReleaseType::MakeReleaseType(relName, mcinfo);
00168   string title = ReleaseType::AsString(release);
00169   static bool checkrel = true;
00170   if (checkrel){
00171     cout<<"Release "<<release<<": "<<title<<endl;
00172     checkrel = false;
00173   }
00174   // ENERGY CORRECTIONS - 1) set reco version
00175   if(ReleaseType::IsCedar(release))
00176     EnergyCorrections::SetCorrectionVersion(EnergyCorrections::kCedar);
00177   if(ReleaseType::IsBirch(release))
00178     EnergyCorrections::SetCorrectionVersion(EnergyCorrections::kBirch);
00179 
00180   // ENERGY CORRECTIONS - 2) set release type for MC or data
00181 
00182   // MC or data
00183   ReleaseType::Release_t reltype=release;
00184 
00185 
00186   // ENERGY CORRECTIONS - 3) set correction version - use default
00187 
00188   EnergyCorrections::WhichCorrection_t corrver=EnergyCorrections::kDefault;
00189 
00190   
00191   
00192   if (fSimFlag==SimFlag::kData && !IsTriggerOn){
00193     //get MinoBooNE spill
00194     MBSpillAccessor& spillAccess = MBSpillAccessor::Get();
00195     
00196     const MBSpill *spillInfo = spillAccess.LoadSpill(vldc.GetTimeStamp());
00197     if (spillInfo) {
00198       fMNtp->dt = double(spillInfo->GetTimeStamp()-vldc.GetTimeStamp());
00199       fMNtp->mbsec = spillInfo->GetSec();
00200       fMNtp->mbnsec = spillInfo->GetNanoSec();
00201       fMNtp->mbpot = spillInfo->GetPOT();
00202       fMNtp->mbhorncurr = spillInfo->GetHorn_current();
00203     }
00204   }
00205 
00206   for(int ievt = 0; ievt<nevt; ievt++){//loop over all events
00207     fMNtp->Reset();  //reset all the variables excepts the header information
00208     fMNtp->evtno = ievt;
00209     NtpSREvent *event = SntpHelpers::GetEvent(ievt,record);
00210     if (!event) continue;
00211 
00212     if (IsTriggerOn){
00213       fN = 10 ;
00214       fM = 12 ;
00215       fMNtp->MNTrig = PlanTriggerSim(event,record) ;
00216 
00217       fN = 10 ;
00218       fM = 600 ;
00219       fMNtp->PlaneTrig = PlanTriggerSim(event,record) ;
00220 
00221       fN = 4 ;
00222       fM = 5 ;
00223       fMNtp->MN45Trig = PlanTriggerSim(event,record) ; 
00224     }
00225 
00226     //index will be -1 if no track/shower present in event
00227     //if there is a track, find the primary shower close to the track vertex
00228     //if there is no track, find the biggest shower
00229     //pull out from Mad
00230     int track_index   = -1;
00231     int shower_index  = -1;
00232     if(ReleaseType::IsBirch(release)){
00233       if (LoadLargestTrackFromEvent(event,record,track_index)){
00234         if(!LoadShowerAtTrackVertex(event,record,track_index,shower_index)){
00235           LoadLargestShowerFromEvent(event,record,shower_index);
00236         }
00237       }
00238       else LoadLargestShowerFromEvent(event,record,shower_index);
00239     }
00240     else if(ReleaseType::IsCedar(release)){
00241       if (event->ntrack){
00242         track_index = event->trk[0];
00243       }
00244       if (event->nshower){
00245         shower_index = event->shw[0];
00246       }
00247     }
00248 
00249     if (track_index!=-1){//found track
00250       NtpSRTrack *track = SntpHelpers::GetTrack(track_index,record);
00251       if (track) {
00252         this->FillTrkInfo(fMNtp,track,record);
00253         fMNtp->trkp_rangecor = FullyCorrectMomentumFromRange(fMNtp->trkp_range,vldc,reltype,corrver);
00254         fMNtp->trkp_fitcor = FullyCorrectSignedMomentumFromCurvature(fMNtp->trkp_fit,vldc,reltype,corrver);      
00255         fMNtp->recoPmu = 0.;
00256         if (fMNtp->trkp_rangecor>0){
00257           fMNtp->recoPmu = fMNtp->trkp_rangecor;
00258         }
00259         if (fMNtp->trkpassfit&&fMNtp->trkexit&&TMath::Abs(fMNtp->trkeqp/fMNtp->trkqp)<0.3){
00260           fMNtp->recoPmu = TMath::Abs(fMNtp->trkp_fitcor);
00261         }
00262         
00263         fMNtp->recoEmu = sqrt(pow(fMNtp->recoPmu,2)+0.105658*0.105658);
00264       }
00265     }
00266 
00267     if (shower_index!=-1){//found shower
00268       NtpSRShower *shower = SntpHelpers::GetShower(shower_index,record);
00269       if (shower) {
00270         this->FillShwInfo(fMNtp,shower);
00271         fMNtp->shwlinCCcor = FullyCorrectShowerEnergy(fMNtp->shwlinCCgev,CandShowerHandle::kCC,vldc,reltype,corrver);
00272         fMNtp->shwwtCCcor = FullyCorrectShowerEnergy(fMNtp->shwwtCCgev,CandShowerHandle::kWtCC,vldc,reltype,corrver);
00273         fMNtp->recoEshw = 0;
00274         if (fMNtp->shwlinCCcor>0) fMNtp->recoEshw = fMNtp->shwlinCCcor;
00275       }
00276     }
00277     else {
00278       fMNtp->shwgev = 0;
00279       fMNtp->shwlinCCgev = 0;
00280       fMNtp->shwwtCCgev = 0;
00281       fMNtp->shwlinCCcor = 0;
00282       fMNtp->shwwtCCcor = 0;
00283       fMNtp->shwlinNCgev = 0;
00284       fMNtp->shwwtNCgev = 0;
00285       fMNtp->shwEMgev = 0;
00286       fMNtp->shwtotgev = 0;
00287       fMNtp->recoEshw = 0;
00288     }
00289 
00290     //calculate david's cc pid
00291     if (track_index!=-1){
00292       NtpSRTrack *track = SntpHelpers::GetTrack(track_index,record);
00293       if(dpid.ChoosePDFs(fDetector,fBeam))  // david pid 
00294          fMNtp->dav_cc_pid = dpid.CalcPID(track,event,
00295                                   &(record->evthdr),fDetector,0);
00296     }  
00297 
00298     //fill event information
00299     this->FillEvtInfo(fMNtp,event);
00300     //correcting vtx
00301     if(ReleaseType::IsCedar(release)){
00302       NtpVtxFinder vtxf;
00303       vtxf.SetTargetEvent(ievt,record);
00304       if (vtxf.FindVertex()>0){
00305         fMNtp->evtvtxx = vtxf.VtxX();
00306         fMNtp->evtvtxy = vtxf.VtxY();
00307         fMNtp->evtvtxz = vtxf.VtxZ();
00308         fMNtp->evtvtxu = vtxf.VtxU();
00309         fMNtp->evtvtxv = vtxf.VtxV();
00310         if (fMNtp->evtvtxz>0.5&&fMNtp->evtvtxz<6&&
00311             pow(fMNtp->evtvtxx-1.4885,2)+pow(fMNtp->evtvtxy-0.1397,2)<1){
00312           fMNtp->fid = 1;
00313         }//fiducial volume cut
00314         else fMNtp->fid = 0;
00315       }
00316     }
00317 
00318 
00319    //calculate event direction: track+shower
00320     if (track_index!=-1&&shower_index!=-1){
00321       NtpSRTrack *track = SntpHelpers::GetTrack(track_index,record);
00322       NtpSRShower *shower = SntpHelpers::GetShower(shower_index,record);
00323       float trkdx = track->vtx.dcosx;
00324       float trkdy = track->vtx.dcosy;
00325       float trkdz = track->vtx.dcosz;
00326       float shwdx = shower->vtx.dcosx;
00327       float shwdy = shower->vtx.dcosy;
00328       float shwdz = shower->vtx.dcosz;
00329       float trkmom = track->momentum.range;
00330       float shwmom = shower->ph.gev;
00331       if (shower->shwph.linCCgev>0) shwmom = shower->shwph.linCCgev;
00332       float px = trkmom*trkdx + shwmom*shwdx;
00333       float py = trkmom*trkdy + shwmom*shwdy;
00334       float pz = trkmom*trkdz + shwmom*shwdz;      
00335 
00336       if (px||py||pz){
00337         float norm = sqrt(px*px+py*py+pz*pz);
00338         float dcosx = -px/norm;
00339         float dcosy = -py/norm;
00340         float dcosz = -pz/norm;
00341         double altitude, evtazimuth_tmp;
00342         AstUtil::LocalToHorizon(dcosx,dcosy,dcosz,Detector::kNear,altitude,evtazimuth_tmp);
00343         fMNtp->evtzenith = 90. - altitude;
00344         fMNtp->evtazimuth = evtazimuth_tmp;
00345         fMNtp->evtcosa = CalCosine(fMNtp->evtzenith,fMNtp->evtazimuth,theta_mb,phi_mb);
00346       }
00347     }
00348     else if (track_index!=-1){
00349       fMNtp->evtzenith = fMNtp->trkzenith;
00350       fMNtp->evtazimuth = fMNtp->trkazimuth;
00351       fMNtp->evtcosa = CalCosine(fMNtp->evtzenith,fMNtp->evtazimuth,theta_mb,phi_mb);
00352 
00353     }
00354     
00355     int nstp = event->nstrip; //no. of strips
00356       
00357     // the following is to check the timing of the earliest hit of an event is consistent with
00358     // the snarl time (just for MB events)
00359     if (fSimFlag==SimFlag::kData && !IsTriggerOn){
00360       //find the earliest time in the event
00361       //because Peter worried there might be a difference between 
00362       //the snarl timestamp and the earliest hit time.
00363       VldTimeStamp t0(2020,1,1,0,0,0);
00364       for (int istp = 0; istp<nstp; istp++){
00365         NtpSRStrip *strip = SntpHelpers::GetStrip(event->stp[istp],record);
00366         //cout<<strip->time1<<" "<<VldTimeStamp(strip->time1)<<" "<<t0<<endl;
00367         if (fMNtp->sec+strip->time1<double(t0)){
00368           t0 = VldTimeStamp(fMNtp->sec+strip->time1);
00369         }
00370       }
00371       
00372       fMNtp->evtsec = (int)t0.GetSec();
00373       fMNtp->evtnsec = (int)t0.GetNanoSec();
00374 
00375       MBSpillAccessor& spillAccess = MBSpillAccessor::Get();
00376 
00377       const MBSpill *spillInfo0 = spillAccess.LoadSpill(t0);
00378       if (spillInfo0) {
00379         fMNtp->dt_evt = double(spillInfo0->GetTimeStamp()-t0);
00380       }
00381 //      if (spillInfo&&spillInfo0&&spillInfo->GetTimeStamp()!=spillInfo0->GetTimeStamp()){
00382 //      MSG("MNtpModule", Msg::kWarning) <<"Different time stamps between snarl and earliest hit "<<vldts<<" "<<t0<<endl;
00383 //      }
00384     }
00385   
00386     //find the highest strip number in fu,pu and v views.
00387     //this is not as useful as trkendu,trkendv and trkendy
00388     for (int istp = 0; istp < nstp; istp++){
00389       NtpSRStrip *strip = SntpHelpers::GetStrip(event->stp[istp],record);
00390       if(strip->planeview==PlaneView::kU){
00391         if((strip->plane-1)%5==0){//fully instrumented
00392           if (strip->strip>fMNtp->maxuf) fMNtp->maxuf = strip->strip;
00393         }
00394         else {//partially instrumented
00395           if (strip->strip>fMNtp->maxup) fMNtp->maxup = strip->strip;
00396         }
00397       }
00398       else if (strip->planeview==PlaneView::kV){
00399         if (strip->strip>fMNtp->maxv) fMNtp->maxv = strip->strip;
00400       }
00401     }
00402     
00403     //fill mc information
00404     if (fSimFlag==SimFlag::kMC){//mc
00405       int index = SntpHelpers::GetEvent2MCIndex(fMNtp->evtno,record);
00406       NtpMCTruth *mctruth = SntpHelpers::GetMCTruth(index,record);
00407       TClonesArray& stdhepArray = *(record->stdhep);
00408       if (mctruth) {
00409         fMNtp->nuenergy = mctruth->p4neu[3];
00410         fMNtp->nuPx = mctruth->p4neu[0];
00411         fMNtp->nuPy = mctruth->p4neu[1];
00412         fMNtp->nuPz = mctruth->p4neu[2];
00413         fMNtp->tarE = mctruth->p4tgt[3];
00414         fMNtp->tarPx = mctruth->p4tgt[0];
00415         fMNtp->tarPy = mctruth->p4tgt[1];
00416         fMNtp->tarPz = mctruth->p4tgt[2];
00417         fMNtp->ires = mctruth->iresonance;
00418         fMNtp->inu = mctruth->inu;
00419         fMNtp->iact = mctruth->iaction;
00420         fMNtp->initial_state = this->Initial_state(mctruth,stdhepArray);
00421         fMNtp->nucleus = this->Nucleus(mctruth);
00422         fMNtp->had_fs = this->HadronicFinalState(mctruth,stdhepArray);
00423         fMNtp->ptype = mctruth->flux.ptype;
00424         fMNtp->tptype = mctruth->flux.tptype;
00425         fMNtp->pdpx = mctruth->flux.pdpx;
00426         fMNtp->pdpy = mctruth->flux.pdpy;
00427         fMNtp->pdpz = mctruth->flux.pdpz;
00428         fMNtp->vx = mctruth->flux.vx;
00429         fMNtp->vy = mctruth->flux.vy;
00430         fMNtp->vz = mctruth->flux.vz;
00431         fMNtp->tpx = mctruth->flux.tpx;
00432         fMNtp->tpy = mctruth->flux.tpy;
00433         fMNtp->tpz = mctruth->flux.tpz;
00434         fMNtp->tvx = mctruth->flux.tvx;
00435         fMNtp->tvy = mctruth->flux.tvy;
00436         fMNtp->tvz = mctruth->flux.tvz;
00437         fMNtp->x = mctruth->x;
00438         fMNtp->y = mctruth->y;
00439         fMNtp->Q2 = mctruth->q2;
00440         fMNtp->W2 = mctruth->w2;
00441         fMNtp->showerE = mctruth->p4shw[3];
00442         if (pow(mctruth->p4mu1[0],2)+pow(mctruth->p4mu1[1],2)+pow(mctruth->p4mu1[2],2)>0) fMNtp->muonP = sqrt(pow(mctruth->p4mu1[0],2)+pow(mctruth->p4mu1[1],2)+pow(mctruth->p4mu1[2],2));
00443       }
00444       if (track_index!=-1){
00445         NtpTHTrack *thtrack = dynamic_cast<NtpTHTrack *>((*record->thtrk)[track_index]);
00446         if (thtrack){
00447           NtpMCStdHep *stdhep = dynamic_cast<NtpMCStdHep *>((*record->stdhep)[thtrack->trkstdhep]);
00448           if (stdhep){
00449             fMNtp->trkpid = stdhep->IdHEP;
00450           }
00451         }
00452       }
00453     }    
00454 
00455     fNtuple->Fill();
00456 
00457   }
00458   return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00459 }
00460 
00461 //......................................................................
00462 
00463 int MNtpModule::Initial_state(NtpMCTruth *mctruth,
00464                               TClonesArray &stdhepArray){
00465   //copied from MadQuantities::Initial_state
00466   Int_t initial_state=0;
00467   Int_t nStdHep = stdhepArray.GetEntries();
00468 
00469   int protneut = -1;  // 0 = neutron, 1 = proton, 2 = N, 3 = A
00470   int nubarnu = 0;    // +1 = neutrino, -1 = antineutrino
00471 
00472   for(int i=0;i<nStdHep;i++){
00473  
00474     NtpMCStdHep *ntpStdHep = dynamic_cast<NtpMCStdHep *>(stdhepArray[i]);
00475 
00476     if(ntpStdHep->mc==mctruth->index){
00477 
00478       if(ntpStdHep->IstHEP==0){  //initial state particle
00479         if(abs(ntpStdHep->IdHEP)==12 ||
00480            abs(ntpStdHep->IdHEP)==14 ||
00481            abs(ntpStdHep->IdHEP)==16){   //(anti)neutrino
00482           nubarnu = ntpStdHep->IdHEP/abs(ntpStdHep->IdHEP);  //get sign
00483         }
00484       }
00485       if(ntpStdHep->IstHEP==11){    //target nucleon
00486         if(ntpStdHep->IdHEP==2212) protneut = 1;  //proton
00487         else if(ntpStdHep->IdHEP==2112) protneut = 0;  //neutron
00488         else if(abs(ntpStdHep->IdHEP)>1000000000) protneut = 2;  //nucleus
00489         else protneut = 3; //atom - probably never get here since IdHEP A>N?
00490       }
00491     }
00492   }
00493 
00494   if(protneut==1 && nubarnu==1)  initial_state=1;  //p + v
00495   if(protneut==0 && nubarnu==1)  initial_state=2;  //n + v
00496   if(protneut==1 && nubarnu==-1) initial_state=3;  //p + vbar
00497   if(protneut==0 && nubarnu==-1) initial_state=4;  //n + vbar
00498   if(protneut==2 && nubarnu==1)  initial_state=5;  //N + v
00499   if(protneut==3 && nubarnu==1)  initial_state=6;  //A + v
00500   if(protneut==2 && nubarnu==-1) initial_state=7;  //N + vbar
00501   if(protneut==3 && nubarnu==-1) initial_state=8;  //A + vbar
00502 
00503   return initial_state;
00504 }
00505 
00506 //......................................................................
00507 
00508 int MNtpModule::Nucleus(NtpMCTruth *mctruth){
00509 
00510   Int_t z=int(mctruth->z);
00511   Int_t a=int(mctruth->a);
00512   Int_t nucleus=1;
00513 
00514   switch (z) {
00515     //  case 1:
00516     //nucleus=0;   // free nucleon
00517     //break;
00518   case 1:
00519     switch (a) {
00520     case 1:
00521       nucleus=256;   // hydrogen1
00522       break;
00523     case 2:
00524       nucleus=257;   // hydrogen2
00525       break;
00526     case 3:
00527       nucleus=258;   // hydrogen2
00528       break;
00529     default:
00530       nucleus=256;   // hydrogen1
00531       break;
00532     }
00533     break;
00534   case 6:
00535     switch (a) {
00536     case 11:
00537       nucleus=274; // carbon11   
00538       break;
00539     case 12:
00540       nucleus=275; // carbon12
00541       break;
00542     case 13:
00543       nucleus=276; // carbon13
00544       break;
00545     case 14:
00546       nucleus=277; // carbon14
00547       break;
00548     default:
00549       nucleus=275; // carbon12
00550       break;
00551     }
00552     break;
00553   case 7:
00554     switch (a) {
00555     case 13:
00556       nucleus=278; // nitrogen13   
00557       break;
00558     case 14:
00559       nucleus=279; // nitrogen14
00560       break;
00561     case 15:
00562       nucleus=280; // nitrogen15
00563       break;
00564     case 16:
00565       nucleus=281; // nitrogen16
00566       break;
00567     case 17:
00568       nucleus=282; // nitrogen17
00569       break;
00570     default:
00571       nucleus=279; // nitrogen14
00572       break;
00573     }
00574     break;
00575   case 8:
00576     switch (a) {
00577     case 15:
00578       nucleus=283;   // oxygen15
00579       break;
00580     case 16:
00581       nucleus=284;   // oxygen16
00582       break;
00583     case 17:
00584       nucleus=285;   // oxygen17
00585       break;
00586     case 18:
00587       nucleus=286;   // oxygen18
00588       break;
00589     default:
00590       nucleus=284;   // oxygen16
00591       break;
00592     }
00593     break;
00594   case 13:
00595     switch (a) {
00596     case 26:
00597       nucleus=303;   // aluminium26
00598       break;
00599     case 27:
00600       nucleus=304;   // aluminium27
00601       break;
00602     case 28:
00603       nucleus=305;   // aluminium28
00604       break;
00605     case 29:
00606       nucleus=306;   // aluminium29
00607       break;
00608     default:
00609       nucleus=304;   // aluminium27
00610       break;
00611     }
00612     break;
00613   case 14:
00614     switch (a) {
00615     case 27:
00616       nucleus=307;   // silicon27
00617       break;
00618     case 28:
00619       nucleus=308;   // silicon28
00620       break;
00621     case 29:
00622       nucleus=309;   // silicon29
00623       break;
00624     case 30:
00625       nucleus=310;   // silicon30
00626       break;
00627     default:         
00628       nucleus=308;   // silicon28
00629       break;
00630     }
00631     break;
00632   case 15:
00633     switch (a) {
00634     case 30:
00635       nucleus=311;   //phosphorus30
00636       break;
00637     case 31:
00638       nucleus=312;   //phosphorus31
00639       break;
00640     case 32:
00641       nucleus=313;   //phosphorus32
00642       break;
00643     case 33:
00644       nucleus=314;   //phosphorus33
00645       break;
00646     default:
00647       nucleus=312;   //phosphorus31
00648       break;
00649     }
00650     break;
00651   case 16:
00652     switch (a) {
00653     case 31:
00654       nucleus=315;   //sulphur31
00655       break;
00656     case 32:
00657       nucleus=316;   //sulphur32
00658       break;
00659     case 33:
00660       nucleus=317;   //sulphur33
00661       break;
00662     case 34:
00663       nucleus=318;   //sulphur34
00664       break;
00665     case 35:
00666       nucleus=319;   //sulphur35
00667       break;
00668     case 36:
00669       nucleus=320;   //sulphur36
00670       break;
00671     default:
00672       nucleus=316;   //sulphur32
00673       break;
00674     }
00675     break;
00676   case 22:
00677     switch (a) {
00678     case 45:
00679       nucleus=347;   //titanium45
00680       break;
00681     case 46:
00682       nucleus=348;   //titanium46
00683       break;
00684     case 47:
00685       nucleus=349;   //titanium47
00686       break;
00687     case 48:
00688       nucleus=350;   //titanium48
00689       break;
00690     case 49:
00691       nucleus=351;   //titanium49
00692       break;
00693     case 50:
00694       nucleus=352;   //titanium50
00695       break;
00696     default:
00697       nucleus=350;   //titanium48
00698       break;
00699     }
00700     break;
00701   case 23:
00702     switch (a) {
00703     case 49:
00704       nucleus=353;   //vanadium49
00705       break;
00706     case 50:
00707       nucleus=354;   //vanadium50
00708       break;
00709     case 51:
00710       nucleus=355;   //vanadium51
00711       break;
00712     case 52:
00713       nucleus=356;   //vanadium52
00714       break;
00715     case 53:
00716       nucleus=357;   //vanadium53
00717       break;
00718     default:
00719       nucleus=355;   //vanadium51
00720       break;
00721     }
00722     break;
00723   case 24:
00724     switch (a) {
00725     case 49:
00726       nucleus=358;   //chromium49
00727       break;
00728     case 50:
00729       nucleus=359;   //chromium50
00730       break;
00731     case 51:
00732       nucleus=360;   //chromium51
00733       break;
00734     case 52:
00735       nucleus=361;   //chromium52
00736       break;
00737     case 53:
00738       nucleus=362;   //chromium53
00739       break;
00740     case 54:
00741       nucleus=363;   //chromium54
00742       break;
00743     default:
00744       nucleus=361;   //chromium52
00745       break;
00746     }
00747     break;
00748   case 25:
00749     switch (a) {
00750     case 53:
00751       nucleus=364;   //manganese53
00752       break;
00753     case 54:
00754       nucleus=365;   //manganese54
00755       break;    
00756     case 55:
00757       nucleus=366;   //manganese55
00758       break;    
00759     case 56:
00760       nucleus=367;   //manganese56
00761       break;    
00762     case 57:
00763       nucleus=368;   //manganese57
00764       break;    
00765     default:
00766       nucleus=366;   //manganese55
00767       break;
00768     }
00769     break;
00770   case 26:
00771     switch (a) {
00772     case 53:
00773       nucleus=369;   //iron53
00774       break;
00775     case 54:
00776       nucleus=370;   //iron54
00777       break;
00778     case 55:
00779       nucleus=371;   //iron55
00780       break;
00781     case 56:
00782       nucleus=372;   //iron56
00783       break;
00784     case 57:
00785       nucleus=373;   //iron57
00786       break;
00787     case 58:
00788       nucleus=374;   //iron58
00789       break;
00790     default:
00791       nucleus=372;   //iron56
00792       break;
00793     }
00794     break;
00795   case 28:
00796     switch (a) {
00797     case 57:
00798       nucleus=382;   //nickel57
00799       break;
00800     case 58:
00801       nucleus=383;   //nickel58
00802       break;
00803     case 59:
00804       nucleus=384;   //nickel59
00805       break;
00806     case 60:
00807       nucleus=385;   //nickel60
00808       break;
00809     case 61:
00810       nucleus=386;   //nickel61
00811       break;
00812     case 62:
00813       nucleus=387;   //nickel62
00814       break;
00815     case 63:
00816       nucleus=388;   //nickel63
00817       break;
00818     case 64:
00819       nucleus=389;   //nickel64
00820       break;
00821     default:
00822       nucleus=383;   //nickel58
00823       break;
00824     }
00825     break;
00826   case 29:
00827     switch (a) {
00828     case 62:
00829       nucleus=390;   //copper62
00830       break;
00831     case 63:
00832       nucleus=391;   //copper63
00833       break;
00834     case 64:
00835       nucleus=392;   //copper64
00836       break;
00837     case 65:
00838       nucleus=393;   //copper65
00839       break;
00840     case 66:
00841       nucleus=394;   //copper66
00842       break;
00843     case 67:
00844       nucleus=395;   //copper67
00845       break;
00846     default:
00847       nucleus=392;   //copper64
00848       break;
00849     }
00850     break;
00851   default:
00852     nucleus=1;   // unknown
00853     break;
00854   }
00855 
00856   return nucleus;
00857 }
00858 
00859 //......................................................................
00860 
00861 int MNtpModule::HadronicFinalState(NtpMCTruth *mctruth,
00862                                    TClonesArray &stdhepArray){
00863   Int_t hfs=0;
00864   Int_t proc = mctruth->iresonance;
00865   Int_t nStdHep = stdhepArray.GetEntries();
00866   if(proc==1002){
00867     for(int i=0;i<nStdHep;i++){
00868 
00869       NtpMCStdHep *ntpStdHep = dynamic_cast<NtpMCStdHep *>(stdhepArray[i]);
00870 
00871       if(ntpStdHep->mc==mctruth->index && ntpStdHep->IstHEP==3 &&
00872          !(abs(ntpStdHep->IdHEP)==15)){  //not a tau lepton
00873         hfs = ntpStdHep->IdHEP;
00874         break;
00875       }
00876     }
00877     hfs = hfs%1000;
00878   }
00879   else {
00880     for(int i=0;i<nStdHep;i++){
00881       NtpMCStdHep *ntpStdHep = dynamic_cast<NtpMCStdHep *>(stdhepArray[i]);
00882       if(ntpStdHep->mc==mctruth->index && ntpStdHep->IstHEP==3){
00883         hfs = ntpStdHep->IdHEP;
00884         break;
00885       }
00886     }
00887     hfs = hfs%1000;
00888   }
00889   return hfs;
00890 }
00891 
00892 //......................................................................
00893 
00894 bool MNtpModule::PlanTriggerSim(NtpSREvent* event,
00895                                 NtpStRecord* record)
00896 {
00897   int nstp = event->nstrip; //no. of strips
00898 
00899   std::vector<Int_t> planeMap(600);
00900   for(UInt_t i = 0;i<planeMap.size();i++) planeMap[i]=0;
00901 
00902   Int_t totplanes = 0;
00903   Int_t firstPlane = 9999;
00904   Int_t lastPlane = -9999;
00905 
00906   for (int istp = 0; istp<nstp; istp++){
00907     NtpSRStrip *strip = SntpHelpers::GetStrip(event->stp[istp],record);
00908     Int_t plane = strip->plane ;
00909 
00910     if(plane > (Int_t)planeMap.size()) continue;
00911     if(strip->ph1.pe == 0 || strip->ndigit==0) continue ;
00912     
00913     // Add plane to map.
00914     if(planeMap[plane] ==0){
00915       totplanes++;
00916       if(plane<firstPlane) firstPlane = plane;
00917       if(plane>lastPlane) lastPlane = plane;
00918     }
00919     
00920     planeMap[plane]++;
00921   }
00922 
00923   if(totplanes < fN) return false;
00924   
00925   for(int i = firstPlane; i <= lastPlane; i++) {
00926     Int_t nPlane = 0;
00927     
00928     for (Int_t j=0; (j<fM) && (i+j<=lastPlane); j++) {
00929       if ( planeMap[i+j]>0 ) nPlane++;
00930     }
00931     if (nPlane>=fN) {
00932       return true;
00933     }
00934     if (i+fM>lastPlane) break; // We're done... stop looking.
00935   }
00936   
00937   return false;
00938 }
00939 
00940 //......................................................................
00941 
00942 bool MNtpModule::LoadLargestTrackFromEvent(NtpSREvent* event,
00943                                       NtpStRecord* record,
00944                                       int& trkidx)
00945 {
00946   bool status = false;
00947   float longest_trk = 0;
00948   for (int i=0; i<event->ntrack; i++){
00949     NtpSRTrack *track = SntpHelpers::GetTrack(event->trk[i],record);
00950     if (!track) continue;
00951     int trklen = abs(track->plane.end-track->plane.beg)+1;
00952     if (trklen>longest_trk){
00953       trkidx = event->trk[i];
00954       longest_trk = trklen;
00955     }
00956   }
00957   if (longest_trk>0){
00958     status = true;
00959   }
00960   return status;
00961 }
00962 
00963 //.....................................................................
00964 
00965 bool MNtpModule::LoadShowerAtTrackVertex(NtpSREvent* event,
00966                                     NtpStRecord* record,
00967                                     int  trkidx,
00968                                     int& shwidx)
00969 {
00970   bool status = false;
00971   NtpSRTrack *track = SntpHelpers::GetTrack(trkidx,record);
00972   if (!track) return false;
00973   float closest_ph = 0.;
00974   float closest_vtx = 1000.;
00975   const float close_enough = 7*0.06; //about 1 interaction length
00976   for (int i=0; i<event->nshower; i++){
00977     NtpSRShower *shower = SntpHelpers::GetShower(event->shw[i],record);
00978     if (!shower) continue;
00979     if (fMNtp->shwtotgev<=0) fMNtp->shwtotgev = shower->ph.gev;
00980     else fMNtp->shwtotgev += shower->ph.gev;
00981     float vtx_sep = TMath::Abs(shower->vtx.z-track->vtx.z);
00982 
00983     //if this is the closest shower
00984     if (vtx_sep<closest_vtx){
00985       shwidx = event->shw[i];
00986       closest_vtx = vtx_sep;
00987       if (shower->shwph.linCCgev>0) closest_ph = shower->shwph.linCCgev;
00988       else closest_ph = shower->ph.gev;
00989     }
00990     // if this isn't the closest shower but it's within close_enough m of the
00991     // track vertex and has a larger pulseheight
00992     else if(vtx_sep<close_enough && shower->ph.gev>closest_ph){
00993       shwidx = event->shw[i];
00994       closest_vtx = vtx_sep;
00995 
00996       if (shower->shwph.linCCgev>0) closest_ph = shower->shwph.linCCgev;
00997       else closest_ph = shower->ph.gev;
00998     }
00999   }
01000   if (closest_vtx<close_enough){
01001     status = true;
01002   }
01003   else shwidx = -1;
01004   return status;
01005 }
01006 
01007 //.........................................................................
01008 
01009 bool MNtpModule::LoadLargestShowerFromEvent(NtpSREvent* event,
01010                                        NtpStRecord* record,
01011                                        int& shwidx)
01012 {
01013   bool status = false;
01014   float largest_ph = 0;
01015   for (int i=0; i<event->nshower; i++){
01016     NtpSRShower *shower = SntpHelpers::GetShower(event->shw[i],record);
01017     if (!shower) continue;
01018     if (shower->ph.gev>largest_ph){
01019       shwidx = event->shw[i];
01020       largest_ph = shower->ph.gev;
01021     }
01022     if (fMNtp->shwtotgev<=0) fMNtp->shwtotgev = shower->ph.gev;
01023     else fMNtp->shwtotgev += shower->ph.gev;
01024   }
01025   if (largest_ph>0){
01026     status = true;
01027   }
01028   return status;
01029 }
01030 
01031 //........................................................................
01032 
01033 void MNtpModule::FillEvtInfo(MNtp* fMNtp, NtpSREvent* event)
01034 {
01035   fMNtp->ntrk = event->ntrack;
01036   fMNtp->nshw = event->nshower;
01037   fMNtp->evtvtxx = event->vtx.x;
01038   fMNtp->evtvtxy = event->vtx.y;
01039   fMNtp->evtvtxz = event->vtx.z;
01040   fMNtp->evtvtxu = event->vtx.u;
01041   fMNtp->evtvtxv = event->vtx.v;
01042 
01043   if (event->vtx.z>0.5&&event->vtx.z<6&&
01044       pow(event->vtx.x-1.4885,2)+pow(event->vtx.y-0.1397,2)<1){
01045     fMNtp->fid = 1;
01046   }//fiducial volume cut
01047   else fMNtp->fid = 0;
01048 
01049   fMNtp->recoEnu = fMNtp->recoEmu;
01050   if (fMNtp->recoEshw>0) fMNtp->recoEnu += fMNtp->recoEshw;
01051 
01052   if (fMNtp->recoPmu>0){
01053     fMNtp->recoQ2 = 2*fMNtp->recoEnu*fMNtp->recoEmu-2*fMNtp->recoPmu*fMNtp->recoEnu*fMNtp->trkcosa - 0.105658*0.105658;
01054     float M = (0.9383+0.9396)/2;
01055     fMNtp->recoW2 = M*M+2*M*fMNtp->recoEshw-fMNtp->recoQ2;
01056   }
01057 }
01058 
01059 //.......................................................................
01060 
01061 void MNtpModule::FillTrkInfo(MNtp* fMNtp, NtpSRTrack* track, NtpStRecord* record)
01062 {
01063   fMNtp->trkplanes = track->plane.n;
01064   fMNtp->trklength = abs(track->plane.end-track->plane.beg)+1;
01065   fMNtp->trkp_range = track->momentum.range;
01066 //  bool isdata=false;
01067 //  if (fSimFlag == SimFlag::kData) isdata=true;
01068 //  fMNtp->trkp_rangecor = CorrectMomentumFromRange(fMNtp->trkp_range, isdata, fDetectorType);
01069   if (track->fit.pass){
01070     fMNtp->trkpassfit = 1;
01071   }
01072   else fMNtp->trkpassfit = 0;
01073   fMNtp->trkqp = track->momentum.qp;
01074   fMNtp->trkeqp = track->momentum.eqp;
01075   if (fMNtp->trkqp){
01076     fMNtp->trkp_fit = 1/fMNtp->trkqp;
01077     //fMNtp->trkp_fitcor = CorrectSignedMomentumFromCurvature(fMNtp->trkp_fit,isdata, fDetectorType);
01078   }
01079   fMNtp->trkvtxx = track->vtx.x;
01080   fMNtp->trkvtxy = track->vtx.y;
01081   fMNtp->trkvtxz = track->vtx.z;
01082   fMNtp->trkvtxu = track->vtx.u;
01083   fMNtp->trkvtxv = track->vtx.v;
01084   fMNtp->trkendx = track->end.x;
01085   fMNtp->trkendy = track->end.y;
01086   fMNtp->trkendz = track->end.z;
01087   fMNtp->trkendu = track->end.u;
01088   fMNtp->trkendv = track->end.v;
01089   fMNtp->trkzenith = track->cr.zenith;
01090   fMNtp->trkazimuth = track->cr.azimuth;
01091   fMNtp->trkcosa = CalCosine(fMNtp->trkzenith,fMNtp->trkazimuth,theta_mb,phi_mb);
01092   fMNtp->trkfiddr = track->fidall.dr;
01093   //determine wheather the track exits the detector
01094   //doc-1547 p6
01095   bool out = true;
01096   //"pitt" fiducial volume
01097   if (track->end.u>0.3&&track->end.u<1.8
01098       &&track->end.v>-1.8&&track->end.v<-0.3&&track->end.x<2.4
01099       &&pow(track->end.x,2)+pow(track->end.y,2)>0.64) out = false;
01100   if ((track->end.z < 7&&out)
01101       ||(track->end.z > 7&&pow(track->end.x-0.8,2)/(1.7*1.7)+pow(track->end.y,2)/(1.4*1.4)>1&&pow(track->end.x,2)+pow(track->end.y,2)>1)
01102       ||(track->end.z>15.6)){
01103     fMNtp->trkexit = 1;
01104   }
01105   else fMNtp->trkexit = 0;
01106   //cout<<out<<" "<<track->end.z<<" "<<pow(track->end.x-0.8,2)/(1.7*1.7)+pow(track->end.y,2)/(1.4*1.4)<<" "<<pow(track->end.x,2)+pow(track->end.y,2)<<" "<<fMNtp->trkexit<<endl;
01107 
01108   //calculate rms from dt_ds fitting
01109   //get VldContext, UgliGeomHandle and VldTimeStamp
01110   VldContext vldc = *(record->GetVldContext());
01111   UgliGeomHandle ugh(vldc);  vector<double> spathLength;
01112   vector<double> st0;
01113   int nstp = track->nstrip;
01114   for (int i = 0; i<nstp; i++){
01115     spathLength.push_back(track->ds-track->stpds[i]);
01116     NtpSRStrip *strip = SntpHelpers::GetStrip(track->stp[i],record);
01117     PlexStripEndId seid(Detector::kNear,strip->plane,strip->strip,StripEnd::kWest);
01118     UgliStripHandle stripHandle = ugh.GetStripHandle(seid);
01119     float halfLength = stripHandle.GetHalfLength();
01120     const TVector3 ghitxyz(track->stpx[i],track->stpy[i],track->stpz[i]);
01121     TVector3 lhitxyz = stripHandle.GlobalToLocal(ghitxyz);
01122     float fiberDist = (halfLength - lhitxyz.x() + stripHandle.ClearFiber(StripEnd::kWest) + stripHandle.WlsPigtail(StripEnd::kWest));
01123     //st0.push_back(track->stpt1[i]-fiberDist/PropagationVelocity::Velocity());
01124     st0.push_back(strip->time1-fiberDist/PropagationVelocity::Velocity());
01125   }
01126   TGraph *gr = new TGraph(st0.size(),&spathLength[0],&st0[0]);
01127   tfit_dt_ds_neg -> SetParameter(0,0.0) ;
01128   gr->Fit("tfit_dt_ds_neg","QR") ;
01129   
01130   tfit_dt_ds_pos -> SetParameter(0,0.0) ;
01131   gr->Fit("tfit_dt_ds_pos","QR") ;
01132   
01133   double trms1 = 0;
01134   double trms2 = 0;
01135   
01136   for (unsigned i = 0; i<st0.size(); i++){
01137     double x,y;
01138     gr->GetPoint(i,x,y);
01139     trms1 += pow(y-tfit_dt_ds_pos->Eval(x),2);
01140     trms2 += pow(y-tfit_dt_ds_neg->Eval(x),2);
01141   }
01142   trms1/=st0.size();
01143   trms2/=st0.size();
01144   trms1=sqrt(trms1)*1e9;
01145   trms2=sqrt(trms2)*1e9;
01146   fMNtp->trktrms_pos = trms1;
01147   fMNtp->trktrms_neg = trms2;
01148 }
01149 
01150 //.......................................................................
01151 
01152 void MNtpModule::FillShwInfo(MNtp* fMNtp, NtpSRShower* shower)
01153 {
01154   fMNtp->shwvtxx = shower->vtx.x;
01155   fMNtp->shwvtxy = shower->vtx.y;
01156   fMNtp->shwvtxz = shower->vtx.z;
01157   fMNtp->shwvtxu = shower->vtx.u;
01158   fMNtp->shwvtxv = shower->vtx.v;
01159   fMNtp->shwgev = shower->ph.gev;
01160   fMNtp->shwlinCCgev = shower->shwph.linCCgev;
01161   fMNtp->shwwtCCgev = shower->shwph.wtCCgev;
01162   fMNtp->shwlinNCgev = shower->shwph.linNCgev;
01163   fMNtp->shwwtNCgev = shower->shwph.wtNCgev;
01164   fMNtp->shwEMgev = shower->shwph.EMgev;
01165   //make corrections for shower energy
01166 //  fMNtp->shwlinCCcor = CorrectShowerEnergy(fMNtp->shwlinCCgev,fDetectorType,CandShowerHandle::kCC); 
01167 //  fMNtp->shwwtCCcor = CorrectShowerEnergy(fMNtp->shwwtCCgev,fDetectorType,CandShowerHandle::kWtCC); 
01168 
01169 }
01170 
01171 //.......................................................................
01172 
01173 float MNtpModule::CalCosine(float theta1,float phi1, float theta2, float phi2)
01174 {
01175   float t1 = theta1*3.1416/180;
01176   float t2 = theta2*3.1416/180;
01177   float p1 = phi1*3.1416/180;
01178   float p2 = phi2*3.1416/180;
01179   return TMath::Sin(t1)*TMath::Sin(t2)*TMath::Cos(p1-p2)+TMath::Cos(t1)*TMath::Cos(t2);
01180 }
01181 
01182 const Registry& MNtpModule::DefaultConfig() const
01183 {
01187   static Registry r; // Default configuration for module
01188 
01189   // Set name of config
01190   std::string name = this->GetName();
01191   name += ".config.default";
01192   r.SetName(name.c_str());
01193 
01194   // Set values in configuration
01195   r.UnLockValues();
01196   r.Set("File","mntp.root");
01197   r.Set("TriggerOn",0) ;
01198   r.LockValues();
01199 
01200   return r;
01201 }
01202 
01203 //......................................................................
01204 
01205 void MNtpModule::Config(const Registry& r)
01206 {
01210   fFile = r.GetCharString("File");
01211   IsTriggerOn = (bool)r.GetInt("TriggerOn") ;
01212 }
01213