00001
00002
00003
00004
00005
00007
00008 #include <set>
00009 #include <cmath>
00010
00011 #include "TClonesArray.h"
00012 #include "TH2F.h"
00013
00014 #include "BeamDataUtil/BDSpillAccessor.h"
00015 #include "Conventions/BeamType.h"
00016 #include "DataUtil/EnergyCorrections.h"
00017 #include "MessageService/MsgService.h"
00018 #include "MessageService/MsgFormat.h"
00019 #include "MCNtuple/NtpMCStdHep.h"
00020 #include "MCNtuple/NtpMCTruth.h"
00021 #include "CandNtupleSR/NtpSREvent.h"
00022 #include "CandNtupleSR/NtpSRShower.h"
00023 #include "CandNtupleSR/NtpSRStrip.h"
00024 #include "CandNtupleSR/NtpSRTrack.h"
00025 #include "StandardNtuple/NtpStRecord.h"
00026 #include "TruthHelperNtuple/NtpTHEvent.h"
00027 #include "TruthHelperNtuple/NtpTHShower.h"
00028 #include "TruthHelperNtuple/NtpTHTrack.h"
00029 #include "DataUtil/PlaneOutline.h"
00030 #include "Conventions/ReleaseType.h"
00031 #include "UgliGeometry/UgliGeomHandle.h"
00032
00033 #include "NtupleUtils/LISieve.h"
00034 #include "NtupleUtils/NuEvent.h"
00035 #include "NtupleUtils/NuLibrary.h"
00036 #include "NtupleUtils/NuReco.h"
00037
00038 using std::endl;
00039 using std::cout;
00040 using std::map;
00041 using std::vector;
00042
00043 namespace EnergyCorrections {}
00044 using namespace EnergyCorrections;
00045
00046 CVSID("$Id: NuReco.cxx,v 1.75 2008/05/16 17:31:15 nickd Exp $");
00047
00048
00049
00050 NuReco::NuReco()
00051 {
00052 MSG("NuReco",Msg::kDebug)
00053 <<"Running NuReco Constructor..."<<endl;
00054
00055
00056 MSG("NuReco",Msg::kDebug)
00057 <<"Finished NuReco Constructor"<<endl;
00058 }
00059
00060
00061 NuReco::~NuReco()
00062 {
00063 MSG("NuReco",Msg::kDebug)
00064 <<"Running NuReco Destructor..."<<endl;
00065
00066
00067 MSG("NuReco",Msg::kDebug)
00068 <<"Finished NuReco Destructor"<<endl;
00069 }
00070
00071
00072
00073 Double_t NuReco::GetEvtEnergy(NuEvent& nu) const
00074 {
00075 Msg::LogLevel_t logLevel=Msg::kDebug;
00076
00077 MAXMSG("NuReco",logLevel,200)
00078 <<"GetEvtEnergy: evt="<<nu.evt<<", trks="<<nu.ntrk
00079 <<", shws="<<nu.nshw<<endl;
00080
00081
00082 this->SetBestTrk(nu);
00083
00084 this->SetBestShw(nu);
00085
00086 this->CalcEvtVariables(nu);
00087 this->CalcTrkVariables(nu);
00088 this->CalcTrkReclamation(nu);
00089
00090 this->CalcShwVariables(nu);
00091
00092
00093 nu.energy=nu.trkEn+nu.shwEn;
00094
00095 nu.energyNoRw=nu.energy;
00096 this->CalcKinematicVariables(nu);
00097
00098 MAXMSG("NuReco",Msg::kInfo,10)
00099 <<endl
00100 <<"GetEvtEnergy:"
00101 <<" nshw="<<nu.nshw
00102 <<", ntrk="<<nu.ntrk
00103 <<", primshw="<<nu.primshw
00104 <<", primtrk="<<nu.primtrk
00105 <<endl
00106 <<"shwExists1,2,3="<<nu.shwExists1<<","
00107 <<nu.shwExists2<<","<<nu.shwExists3
00108 <<", trkExists1,2,3="<<nu.trkExists1<<","
00109 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
00110
00111 MAXMSG("NuReco",Msg::kInfo,10)
00112 <<"E="<<nu.energy<<", trkEn="<<nu.trkEn<<", shwEn="<<nu.shwEn
00113 <<endl
00114 <<"dircosnu="<<nu.dirCosNu
00115 <<", y="<<nu.y<<", q2="<<nu.q2
00116 <<", w2="<<nu.w2<<", x="<<nu.x<<endl;
00117
00118
00119 if (nu.ntrk>3 || nu.nshw>7) {
00120 MAXMSG("NuReco",Msg::kWarning,3)
00121 <<endl<<"Lots of trks/shws"
00122 <<", ntrk="<<nu.ntrk<<", nshw="<<nu.nshw
00123 <<", primtrk="<<nu.primtrk<<", primshw="<<nu.primshw
00124 <<", trkEn="<<nu.trkEn<<", shwEn="<<nu.shwEn
00125 <<endl
00126 <<"shwEn1="<<nu.shwEnCor1<<", shw2="<<nu.shwEnCor2
00127 <<", shw3="<<nu.shwEnCor3<<", shw4="<<nu.shwEnCor4
00128 <<", shw5="<<nu.shwEnCor5
00129
00130 <<endl
00131 <<"trkCv1="<<nu.trkEnCorCurv1<<", trkCv2="<<nu.trkEnCorCurv2
00132 <<", trkCv3="<<nu.trkEnCorCurv3<<endl
00133 <<"trkRg1="<<nu.trkEnCorRange1<<", trkRg2="<<nu.trkEnCorRange2
00134 <<", trkRg3="<<nu.trkEnCorRange3<<endl;
00135 if (nu.ntrk>3) {
00136 MAXMSG("NuReco",Msg::kError,300)
00137 <<endl<<"Lots of trks/shws"
00138 <<", ntrk="<<nu.ntrk<<", nshw="<<nu.nshw
00139 <<", primtrk="<<nu.primtrk<<", primshw="<<nu.primshw
00140 <<", trkEn="<<nu.trkEn<<", shwEn="<<nu.shwEn
00141 <<endl
00142 <<"shwEn1="<<nu.shwEnCor1<<", shw2="<<nu.shwEnCor2
00143 <<", shw3="<<nu.shwEnCor3<<", shw4="<<nu.shwEnCor4
00144 <<", shw5="<<nu.shwEnCor5
00145
00146 <<endl
00147 <<"trkCv1="<<nu.trkEnCorCurv1<<", trkCv2="<<nu.trkEnCorCurv2
00148 <<", trkCv3="<<nu.trkEnCorCurv3<<endl
00149 <<"trkRg1="<<nu.trkEnCorRange1<<", trkRg2="<<nu.trkEnCorRange2
00150 <<", trkRg3="<<nu.trkEnCorRange3<<endl;
00151 }
00152 }
00153
00154 return nu.energy;
00155 }
00156
00157
00158
00159 void NuReco::CalcKinematicVariables(NuEvent& nu) const
00160 {
00161 nu.dirCosNu=this->GetDirCosNu(nu);
00162
00164
00165
00166
00167
00168
00169
00170
00171 const Double_t M=(0.93827 + 0.93957)/2.0;
00172
00173
00174 if (nu.trkEn+nu.shwEn) nu.y=nu.shwEn/(nu.trkEn+nu.shwEn);
00175 nu.q2=2*nu.energy*nu.trkEn*(1.0-nu.dirCosNu);
00176 nu.w2=M*M-nu.q2+2*M*nu.shwEn;
00177 if (nu.shwEn) nu.x=nu.q2/(2*M*nu.shwEn);
00178
00179 }
00180
00181
00182
00183 void NuReco::CalcEvtVariables(NuEvent& nu) const
00184 {
00185
00186 nu.rEvtVtx=this->GetRadius(nu.xEvtVtx,nu.yEvtVtx,nu);
00187 nu.rEvtEnd=this->GetRadius(nu.xEvtEnd,nu.yEvtEnd,nu);
00188 nu.evtVtxUVDiffPl=nu.planeEvtBegu-nu.planeEvtBegv;
00189 nu.distToEdgeEvtVtx=this->GetSmallestDeepDistToEdge(nu);
00190 }
00191
00192
00193
00194 void NuReco::CalcTrkVariables(NuEvent& nu) const
00195 {
00196 if (nu.trkExists) {
00197 nu.rTrkVtx=this->GetRadius(nu.xTrkVtx,nu.yTrkVtx,nu);
00198 nu.rTrkEnd=this->GetRadius(nu.xTrkEnd,nu.yTrkEnd,nu);
00199
00200
00201 if (nu.qp) nu.sigqp_qp=nu.sigqp/nu.qp;
00202 if (nu.ndof) {
00203 nu.chi2PerNdof=nu.chi2/nu.ndof;
00204 nu.prob=TMath::Prob(nu.chi2,static_cast<Int_t>(nu.ndof));
00205 }
00206
00207
00208 nu.containmentFlagCC0093Std=this->GetContainmentFlagCC0093Std(nu);
00209 nu.containmentFlagCC0250Std=this->GetContainmentFlagCC0250Std(nu);
00210 nu.containmentFlagPitt=this->GetContainmentFlagPitt(nu);
00211 nu.containmentFlag=this->GetContainmentFlag(nu);
00212
00213
00214 nu.trkEnRange=this->GetTrackEnergyFromRange(nu);
00215 nu.trkEnCurv=this->GetTrackEnergyFromCurv(nu);
00216 nu.usedRange=this->UseRange(nu);
00217 nu.usedCurv=this->UseCurvature(nu);
00218 nu.trkEn=this->GetTrackEnergy(nu);
00219 nu.trkEnNoRw=nu.trkEn;
00220 }
00221 else {
00222 MAXMSG("NuReco",Msg::kWarning,10)
00223 <<"No track, nu.primtrk="<<nu.primtrk<<endl;
00224 nu.rTrkVtx=-999;
00225 nu.rTrkEnd=-999;
00226
00227 nu.sigqp_qp=-999;
00228 nu.chi2PerNdof=-1;
00229 nu.prob=-1;
00230
00231 nu.containmentFlagCC0093Std=-1;
00232 nu.containmentFlagPitt=-1;
00233 nu.containmentFlag=-1;
00234
00235 nu.trkEnRange=0;
00236 nu.trkEnCurv=0;
00237 nu.usedRange=0;
00238 nu.usedCurv=0;
00239 nu.trkEn=0;
00240 nu.trkEnNoRw=nu.trkEn;
00241 }
00242 }
00243
00244
00245
00246 void NuReco::CalcTrkReclamation(NuEvent& nu) const
00247 {
00248 if (nu.trkExists && nu.trkEnCurv==0) {
00249
00250
00251 if (nu.detector==Detector::kFar) {
00252 MAXMSG("NuReco",Msg::kDebug,100)
00253 <<"Looks like track reclamation for the FD..."
00254 <<", usedCurv="<<nu.usedCurv<<", trkEnCurv="<<nu.trkEnCurv
00255 <<", trkEnRange="<<nu.trkEnRange
00256 <<", trkfitpass="<<nu.trkfitpass
00257 <<endl;
00258 }
00259
00260 MAXMSG("NuReco",Msg::kDebug,10)
00261 <<endl<<"Setting curv=range"
00262 <<", usedCurv="<<nu.usedCurv<<", trkEnCurv="<<nu.trkEnCurv
00263 <<", trkEnRange="<<nu.trkEnRange<<", trkfitpass="<<nu.trkfitpass
00264 <<endl;
00265 nu.trkEnCurv=nu.trkEnRange;
00266 nu.trkEn=nu.trkEnRange;
00267 nu.trkEnNoRw=nu.trkEn;
00268 }
00269 else if (!nu.trkExists && nu.trkEnCurv==0) {
00270 MAXMSG("NuReco",Msg::kWarning,10)
00271 <<"nu.trkExists="<<nu.trkExists<<", nu.trkEnCurv="<<nu.trkEnCurv
00272 <<endl;
00273 }
00274 }
00275
00276
00277
00278 void NuReco::CalcShwVariables(NuEvent& nu) const
00279 {
00280 if (nu.shwExists) {
00281 nu.shwEn=this->GetShowerEnergy(nu);
00282 nu.shwEnNoRw=nu.shwEn;
00283 }
00284 else {
00285 nu.shwEn=0;
00286 nu.shwEnNoRw=nu.shwEn;
00287 }
00288 }
00289
00290
00291
00292 void NuReco::CalcExtraTruthVariables(NuEvent& nu) const
00293 {
00294
00295 const VldContext& vc=this->GetVldContext(nu);
00296 UgliGeomHandle ugh(vc);
00297
00298
00299 const NuLibrary& lib=NuLibrary::Instance();
00300
00301
00302 TVector3 uvz=lib.reco.xyz2uvz(nu.vtxxMC,nu.vtxyMC,nu.vtxzMC,vc);
00303
00304
00305 nu.vtxuMC=uvz.X();
00306 nu.vtxvMC=uvz.Y();
00307
00308
00309 nu.planeTrkVtxMC=ugh.GetPlaneIdFromZ(nu.vtxzMC).GetPlane();
00310 nu.rTrkVtxMC=lib.reco.GetRadius(nu.vtxxMC,nu.vtxyMC,nu);
00311 }
00312
00313
00314
00315 void NuReco::ApplyReweights(NuEvent& nu) const
00316 {
00317 if (nu.simFlag==SimFlag::kData) {
00318 MAXMSG("NuReco",Msg::kInfo,1)
00319 <<"Not doing ApplyReweights for simFlag="<<nu.simFlag<<endl;
00320
00321 nu.trkEnRw=nu.trkEnNoRw;
00322 nu.shwEnRw=nu.shwEnNoRw;
00323 nu.energyRw=nu.trkEnRw+nu.shwEnRw;
00324 return;
00325 }
00326
00327
00328 nu.trkEnRw=nu.trkEnNoRw*nu.trkEnWeight;
00329 nu.shwEnRw=nu.shwEnNoRw*nu.shwEnWeight;
00330 nu.energyRw=nu.trkEnRw+nu.shwEnRw;
00331
00332
00333 if (nu.applyEnergyShifts) {
00334 MAXMSG("NuReco",Msg::kInfo,1)
00335 <<"Applying SKZP energy shifts"<<endl;
00336 nu.trkEn=nu.trkEnRw;
00337 nu.shwEn=nu.shwEnRw;
00338 nu.energy=nu.energyRw;
00339 }
00340 else {
00341 MAXMSG("NuReco",Msg::kInfo,1)
00342 <<"NOT applying SKZP energy shifts"<<endl;
00343 }
00344
00345
00346
00347 if (nu.anaVersion==NuCuts::kJJH1 ||
00348 nu.anaVersion==NuCuts::kJJE1 ||
00349 nu.anaVersion==NuCuts::kFullDST || NuCuts::IsNMBPQ(nu) ) {
00350 MAXMSG("NuReco",Msg::kInfo,1)
00351 <<"Setting the detector weight for NuMuBar selection"<<endl;
00352 nu.detectorWeight=nu.detectorWeightNMB;
00353 }
00354 else {
00355 MAXMSG("NuReco",Msg::kInfo,1)
00356 <<"Setting the detector weight for NuMu selection"<<endl;
00357 nu.detectorWeight=nu.detectorWeightNM;
00358 }
00359
00360
00361
00362 nu.rw=1;
00363
00364
00365 if (nu.applyBeamWeight) {
00366 MAXMSG("NuReco",Msg::kInfo,1)
00367 <<"Applying beam weight"<<endl;
00368 nu.rw*=nu.beamWeight;
00369 }
00370 else {
00371 MAXMSG("NuReco",Msg::kInfo,1)
00372 <<"NOT applying beam weight"<<endl;
00373 }
00374
00375
00376 if (nu.applyDetectorWeight) {
00377 MAXMSG("NuReco",Msg::kInfo,1)
00378 <<"Applying detector weight"<<endl;
00379 nu.rw*=nu.detectorWeight;
00380 }
00381 else {
00382 MAXMSG("NuReco",Msg::kInfo,1)
00383 <<"NOT applying detector weight"<<endl;
00384 }
00385
00386
00387 if (nu.applyGeneratorWeight) {
00388 MAXMSG("NuReco",Msg::kInfo,1)
00389 <<"Applying generator weight"<<endl;
00390 nu.rw*=nu.generatorWeight;
00391 }
00392 else {
00393 MAXMSG("NuReco",Msg::kInfo,1)
00394 <<"NOT applying generator weight"<<endl;
00395 }
00396
00397
00398 nu.rwActual=nu.rw;
00399
00401
00402 nu.fluxErr=nu.fluxErrTotalErrorAfterTune;
00403 }
00404
00405
00406
00407 Int_t NuReco::FDRCBoundary(const NuEvent& nu) const
00408 {
00410 Int_t litag=0;
00411
00412 if(nu.detector!=Detector::kFar) return litag;
00413
00414
00415
00416
00417
00418 Int_t numshower=nu.nshw;
00419 Float_t allph=nu.rawPhEvt;
00420
00421
00422 Int_t plbeg=nu.planeEvtHdrBeg;
00423 Int_t plend=nu.planeEvtHdrEnd;
00424
00425 if (numshower) {
00426 if (allph>1e6) litag+=10;
00427
00428 if (((plbeg==1 || plbeg==2) && (plend==63 || plend==64)) ||
00429 ((plbeg==65 || plbeg==66) && (plend==127 || plend==128)) ||
00430 ((plbeg==129 || plbeg==130) && (plend==191 || plend==192)) ||
00431 ((plbeg==193 || plbeg==194) && (plend==247 || plend==248))) litag++;
00432 if (((plbeg==250 || plbeg==251) && (plend==312 || plend==313)) ||
00433 ((plbeg==314 || plbeg==315) && (plend==376 || plend==377)) ||
00434 ((plbeg==378 || plbeg==379) && (plend==440 || plend==441)) ||
00435 ((plbeg==442 || plbeg==443) && (plend==484 || plend==485))) litag++;
00436 }
00437 return litag;
00438 }
00439
00440
00441
00442 Double_t NuReco::GetTrackEnergy(const NuEvent& nu) const
00443 {
00444
00445
00446 if (this->UseRange(nu)) {
00447 return this->GetTrackEnergyFromRange(nu);
00448 }
00449
00450 else if (this->UseCurvature(nu)){
00451 return this->GetTrackEnergyFromCurv(nu);
00452 }
00453 else {
00454 cout<<"Ahhhh, can't use range or curvature"<<endl;
00455 return -1;
00456 }
00457 }
00458
00459
00460
00461 Double_t NuReco::GetTrackEnergyFromRange(const NuEvent& nu) const
00462 {
00463
00464 if (true) {
00465 MAXMSG("NuReco",Msg::kInfo,1)
00466 <<"GetTrackEnergyFromRange: (re)doing EnergyCorrection"<<endl;
00467 Double_t trkEn=this->GetTrackEnergyFromRangeCor
00468 (nu.trkMomentumRange,nu);
00469 MAXMSG("NuReco",Msg::kDebug,3000)
00470 <<"Range: trkEn="<<trkEn<<", trkEnCorRange="<<nu.trkEnCorRange
00471 <<", trkMomentumRange="<<nu.trkMomentumRange<<endl;
00472 return trkEn;
00473 }
00474 else {
00475 MAXMSG("NuReco",Msg::kInfo,1)
00476 <<"GetTrackEnergyFromRange: NOT (re)doing EnergyCorrection"<<endl;
00477 return nu.trkEnRange;
00478 }
00479 }
00480
00481
00482
00483 void NuReco::CapTrackMomentum(Double_t& trkMom) const
00484 {
00485
00486 if (trkMom>(1000*Munits::GeV)) trkMom=1000*Munits::GeV;
00487 else if (trkMom<(-1000*Munits::GeV)) trkMom=-1000*Munits::GeV;
00488 }
00489
00490
00491
00492 Double_t NuReco::GetTrackEnergyFromRangeCor(Double_t trkMomentumRange,
00493 const NuEvent& nu) const
00494 {
00495
00496 if (trkMomentumRange<0) return -1;
00497
00498 this->CapTrackMomentum(trkMomentumRange);
00499
00500 Double_t mr=0;
00501 VldContext vc=this->GetVldContext(nu);
00502 mr=EnergyCorrections::FullyCorrectMomentumFromRange
00503 (trkMomentumRange,vc,nu.releaseType);
00504
00505 return sqrt(pow(mr,2)+pow(0.105658357*(Munits::GeV),2));
00506 }
00507
00508
00509
00510 Double_t NuReco::GetTrackEnergyFromCurv(const NuEvent& nu) const
00511 {
00512
00513 if (true) {
00514 MAXMSG("NuReco",Msg::kInfo,1)
00515 <<"GetTrackEnergyFromCurv: (re)doing EnergyCorrection"<<endl;
00516 Double_t trkEn=this->GetTrackEnergyFromCurvCor(nu.qp,nu);
00517 Double_t p=-1;
00518 if (nu.qp) p=1./nu.qp;
00519 MAXMSG("NuReco",Msg::kDebug,3000)
00520 <<"Curv: trkEn="<<trkEn<<", trkEnCorCurv="<<nu.trkEnCorCurv
00521 <<", 1/qp="<<p<<endl;
00522 return trkEn;
00523 }
00524 else {
00525 MAXMSG("NuReco",Msg::kInfo,1)
00526 <<"GetTrackEnergyFromCurv: NOT (re)doing EnergyCorrection"<<endl;
00527 return nu.trkEnCurv;
00528 }
00529 }
00530
00531
00532
00533 Double_t NuReco::GetTrackEnergyFromCurvCor(Double_t qp,
00534 const NuEvent& nu) const
00535 {
00539 if (qp!=0){
00540 Double_t p=1./qp;
00541 this->CapTrackMomentum(p);
00542 VldContext vc=this->GetVldContext(nu);
00543 p=EnergyCorrections::FullyCorrectSignedMomentumFromCurvature
00544 (p,vc,nu.releaseType);
00545 return sqrt(pow(p,2)+pow(0.105658357*(Munits::GeV),2));
00546 }
00547 else return 0.;
00548 }
00549
00550
00551
00552 Double_t NuReco::GetShowerEnergy(const NuEvent& nu) const
00553 {
00554
00555
00556
00557
00558
00559
00560 if (true) {
00561 MAXMSG("NuReco",Msg::kInfo,1)
00562 <<"GetShowerEnergy: (re)doing EnergyCorrection"<<endl;
00563 Double_t shwEnCor=this->GetShowerEnergyCor(nu.shwEnNoCor,nu);
00564 MAXMSG("NuReco",Msg::kDebug,300)
00565 <<"Redo energy cor: shwEnNoCor="<<nu.shwEnNoCor
00566 <<", shwEnCor="<<shwEnCor<<endl;
00567 return shwEnCor;
00568 }
00569 else {
00570 MAXMSG("NuReco",Msg::kInfo,1)
00571 <<"GetShowerEnergy: NOT (re)doing EnergyCorrection"<<endl;
00572 return nu.shwEnCor;
00573 }
00574 }
00575
00576
00577
00578 Double_t NuReco::GetShowerEnergyCor(Double_t shwEn,
00579 const NuEvent& nu) const
00580 {
00581
00582 if (shwEn<0) return -1;
00583
00584
00585 VldContext vc=this->GetVldContext(nu);
00586 return EnergyCorrections::FullyCorrectShowerEnergy
00587 (shwEn,CandShowerHandle::kCC,vc,nu.releaseType);
00588 }
00589
00590
00591
00592 const NtpSRTrack* NuReco::GetLongestTrackTLikePlanes
00593 (const NtpStRecord& ntp,const NtpSREvent& evt) const
00594 {
00595 TClonesArray& trkTca=(*ntp.trk);
00596
00597
00598 MAXMSG("NuReco",Msg::kInfo,1)
00599 <<"Setting best track as longest track (track-like planes)"<<endl;
00600 Int_t trkToUse=0;
00601 Double_t trkLength=0;
00602 if (evt.ntrack>1){
00603 for(Int_t itrk=0;itrk<evt.ntrack;itrk++){
00604 const NtpSRTrack& trk=
00605 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
00606
00607 MAXMSG("NuReco",Msg::kDebug,200)
00608 <<"trk="<<itrk<<"/"<<evt.ntrack
00609 <<", pass="<<trk.fit.pass
00610 <<", ntrklike="<<trk.plane.ntrklike
00611 <<endl;
00612
00613
00614 if (trk.plane.ntrklike>trkLength) {
00615 trkLength=trk.plane.ntrklike;
00616 trkToUse=itrk;
00617 }
00618 }
00619 MAXMSG("NuReco",Msg::kDebug,200)
00620 <<"selected trkToUse="<<trkToUse<<endl;
00621 const NtpSRTrack& trk=
00622 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[trkToUse]]);
00623 return &trk;
00624 }
00625 else if (evt.ntrack==1){
00626 const NtpSRTrack& trk=
00627 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
00628 return &trk;
00629 }
00630 else return 0;
00631 }
00632
00633
00634
00635 Int_t NuReco::GetLongestTrack(const NuEvent& nu) const
00636 {
00637 if (nu.ntrk<=0) return -1;
00638 else if (nu.ntrk==1) return 1;
00639 else if (nu.ntrk==2) {
00640 Int_t trkIndex=1;
00641 Int_t longestTrack=nu.trkLength1;
00642 if (nu.trkLength2>longestTrack) {
00643 trkIndex=2;
00644 longestTrack=nu.trkLength2;
00645 }
00646 return trkIndex;
00647 }
00648 else {
00649 Int_t trkIndex=1;
00650 Int_t longestTrack=nu.trkLength1;
00651 if (nu.trkLength2>longestTrack) {
00652 trkIndex=2;
00653 longestTrack=nu.trkLength2;
00654 }
00655 if (nu.trkLength3>longestTrack) {
00656 trkIndex=3;
00657 longestTrack=nu.trkLength3;
00658 }
00659 return trkIndex;
00660 }
00661 }
00662
00663
00664
00665 const NtpSRTrack* NuReco::GetLongestTrack(const NtpStRecord& ntp,
00666 const NtpSREvent& evt) const
00667 {
00668 TClonesArray& trkTca=(*ntp.trk);
00669
00670
00671 MAXMSG("NuReco",Msg::kInfo,1)
00672 <<"Setting best track as longest track (end-beg+1)"<<endl;
00673 Int_t trkToUse=0;
00674 Double_t trkLength=0;
00675 if (evt.ntrack>1){
00676 for(Int_t itrk=0;itrk<evt.ntrack;itrk++){
00677 const NtpSRTrack& trk=
00678 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
00679
00680 Int_t absLength=abs(trk.plane.end-trk.plane.beg+1);
00681
00682
00683 if (absLength>trkLength) {
00684 trkLength=absLength;
00685 trkToUse=itrk;
00686 }
00687
00688 MAXMSG("NuReco",Msg::kDebug,200)
00689 <<"itrk="<<itrk<<" of "<<evt.ntrack
00690 <<", trkToUse="<<trkToUse
00691 <<", pass="<<trk.fit.pass
00692 <<", ntrklike="<<trk.plane.ntrklike
00693 <<", absLength="<<absLength
00694 <<", longest="<<trkLength
00695 <<endl;
00696 }
00697 MAXMSG("NuReco",Msg::kDebug,200)
00698 <<"Finished loop, selected trkToUse="<<trkToUse<<endl;
00699 const NtpSRTrack& trk=
00700 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[trkToUse]]);
00701 return &trk;
00702 }
00703 else if (evt.ntrack==1){
00704 const NtpSRTrack& trk=
00705 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
00706 return &trk;
00707 }
00708 else return 0;
00709 }
00710
00711
00712
00713 void NuReco::SetBestShw(NuEvent& nu) const
00714 {
00715
00716 if (nu.nshw<=0) {
00717 nu.shwExists=false;
00718 return;
00719 }
00720
00721
00722 Int_t bestShower=this->GetBestShower(nu);
00723
00724
00725 if (bestShower<1) {
00726 MAXMSG("NuReco",Msg::kDebug,3)
00727 <<"No best shower, run="<<nu.run<<", snarl="<<nu.snarl
00728 <<", primshw="<<nu.primshw
00729 <<", nshw="<<nu.nshw
00730 <<", shwExists1,2,3="<<nu.shwExists1<<","
00731 <<nu.shwExists2<<","<<nu.shwExists3<<endl;
00732
00733
00734 nu.shwExists=false;
00735 nu.ndigitShw=-1;
00736 nu.nstripShw=-1;
00737 nu.shwEnCor=-1;
00738 nu.shwEnNoCor=-1;
00739 nu.shwEnMip=-1;
00740 nu.planeShwBeg=-1;
00741 nu.planeShwEnd=-1;
00742 nu.xShwVtx=-999;
00743 nu.yShwVtx=-999;
00744 nu.zShwVtx=-999;
00745 return;
00746 }
00747 else if (bestShower==1) {
00748 MAXMSG("NuReco",Msg::kDebug,100)
00749 <<"Best shower=first, primshw="<<nu.primshw
00750 <<", nshw="<<nu.nshw
00751 <<", shwExists1,2,3="<<nu.shwExists1<<","
00752 <<nu.shwExists2<<","<<nu.shwExists3<<endl;
00753 nu.shwExists=nu.shwExists1;
00754 nu.ndigitShw=nu.ndigitShw1;
00755 nu.nstripShw=nu.nstripShw1;
00756 nu.shwEnCor=nu.shwEnCor1;
00757 nu.shwEnMip=nu.shwEnMip1;
00758 nu.shwEnNoCor=nu.shwEnNoCor1;
00759 nu.planeShwBeg=nu.planeShwBeg1;
00760 nu.planeShwEnd=nu.planeShwEnd1;
00761 nu.xShwVtx=nu.xShwVtx1;
00762 nu.yShwVtx=nu.yShwVtx1;
00763 nu.zShwVtx=nu.zShwVtx1;
00764 }
00765 else if (bestShower==2) {
00766 MAXMSG("NuReco",Msg::kWarning,3)
00767 <<"Best shower=second, primshw="<<nu.primshw
00768 <<", nshw="<<nu.nshw
00769 <<", shwExists1,2,3="<<nu.shwExists1<<","
00770 <<nu.shwExists2<<","<<nu.shwExists3<<endl;
00771 nu.shwExists=nu.shwExists2;
00772 nu.ndigitShw=nu.ndigitShw2;
00773 nu.nstripShw=nu.nstripShw2;
00774 nu.shwEnCor=nu.shwEnCor2;
00775 nu.shwEnMip=nu.shwEnMip2;
00776 nu.shwEnNoCor=nu.shwEnNoCor2;
00777 nu.planeShwBeg=nu.planeShwBeg2;
00778 nu.planeShwEnd=nu.planeShwEnd2;
00779 nu.xShwVtx=nu.xShwVtx2;
00780 nu.yShwVtx=nu.yShwVtx2;
00781 nu.zShwVtx=nu.zShwVtx2;
00782 }
00783 else if (bestShower==3) {
00784 MAXMSG("NuReco",Msg::kWarning,3)
00785 <<"Best shower=third, primshw="<<nu.primshw
00786 <<", nshw="<<nu.nshw
00787 <<", shwExists1,2,3="<<nu.shwExists1<<","
00788 <<nu.shwExists2<<","<<nu.shwExists3<<endl;
00789 nu.shwExists=nu.shwExists3;
00790 nu.ndigitShw=nu.ndigitShw3;
00791 nu.nstripShw=nu.nstripShw3;
00792 nu.shwEnCor=nu.shwEnCor3;
00793 nu.shwEnMip=nu.shwEnMip3;
00794 nu.shwEnNoCor=nu.shwEnNoCor3;
00795 nu.planeShwBeg=nu.planeShwBeg3;
00796 nu.planeShwEnd=nu.planeShwEnd3;
00797 nu.xShwVtx=nu.xShwVtx3;
00798 nu.yShwVtx=nu.yShwVtx3;
00799 nu.zShwVtx=nu.zShwVtx3;
00800 }
00801 else cout<<"Ahhhhhhhhhhhh"<<endl;
00802 }
00803
00804
00805
00806 void NuReco::SetBestTrk(NuEvent& nu) const
00807 {
00808
00809 if (nu.ntrk<=0) {
00810 nu.trkExists=false;
00811 return;
00812 }
00813
00814
00815 Int_t bestTrack=this->GetBestTrack(nu);
00816
00817
00818 if (bestTrack<1) {
00819 MAXMSG("NuReco",Msg::kInfo,3)
00820 <<"No best trk, primtrk="<<nu.primtrk
00821 <<", ntrk="<<nu.ntrk
00822 <<", trkExists1,2,3="<<nu.trkExists1<<","
00823 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
00824
00825
00826 nu.trkExists=false;
00827 nu.trkIndex=-1;
00828 nu.ndigitTrk=-1;
00829 nu.nstripTrk=-1;
00830 nu.trkEnCorRange=-1;
00831 nu.trkEnCorCurv=-1;
00832 nu.trkMomentumRange=-1;
00833 nu.containedTrk=0;
00834 nu.trkfitpass=-1;
00835 nu.trkvtxdcosz=-999;
00836 nu.trkvtxdcosy=-999;
00837 nu.trknplane=-999;
00838 nu.charge=0;
00839 nu.qp=-999;
00840 nu.sigqp=-1;
00841 nu.chi2=-1;
00842 nu.ndof=0;
00843 nu.qpFraction=-1;
00844 nu.trkVtxUVDiffPl=-999;
00845 nu.trkLength=-1;
00846 nu.planeTrkNu=-1;
00847 nu.planeTrkNv=-1;
00848 nu.ntrklike=-1;
00849 nu.trkphsigcor=-1;
00850 nu.trkphsigmap=-1;
00851
00852 nu.jitter=-1;
00853 nu.jPID=-999;
00854 nu.majC=0;
00855
00856
00857
00858 nu.smoothMajC=-999;
00859
00860
00861
00862 nu.xTrkVtx=-999;
00863 nu.yTrkVtx=-999;
00864 nu.zTrkVtx=-999;
00865 nu.uTrkVtx=-999;
00866 nu.vTrkVtx=-999;
00867 nu.planeTrkVtx=-1;
00868 nu.planeTrkBeg=-1;
00869 nu.planeTrkBegu=-1;
00870 nu.planeTrkBegv=-1;
00871
00872 nu.xTrkEnd=-999;
00873 nu.yTrkEnd=-999;
00874 nu.zTrkEnd=-999;
00875 nu.uTrkEnd=-999;
00876 nu.vTrkEnd=-999;
00877 nu.planeTrkEnd=-1;
00878 nu.planeTrkEndu=-1;
00879 nu.planeTrkEndv=-1;
00880
00881 nu.drTrkFidall=-999;
00882 nu.dzTrkFidall=-999;
00883 nu.drTrkFidvtx=-999;
00884 nu.dzTrkFidvtx=-999;
00885 nu.drTrkFidend=-999;
00886 nu.dzTrkFidend=-999;
00887 return;
00888 }
00889 else if (bestTrack==1) {
00890 MAXMSG("NuReco",Msg::kDebug,3)
00891 <<"Best trk=first, primtrk="<<nu.primtrk
00892 <<", ntrk="<<nu.ntrk
00893 <<", trkExists1,2,3="<<nu.trkExists1<<","
00894 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
00895 nu.trkExists=nu.trkExists1;
00896 nu.trkIndex=nu.trkIndex1;
00897 nu.ndigitTrk=nu.ndigitTrk1;
00898 nu.nstripTrk=nu.nstripTrk1;
00899 nu.trkEnCorRange=nu.trkEnCorRange1;
00900 nu.trkEnCorCurv=nu.trkEnCorCurv1;
00901 nu.trkMomentumRange=nu.trkMomentumRange1;
00902 nu.containedTrk=nu.containedTrk1;
00903 nu.trkfitpass=nu.trkfitpass1;
00904 nu.trkvtxdcosz=nu.trkvtxdcosz1;
00905 nu.trkvtxdcosy=nu.trkvtxdcosy1;
00906 nu.trknplane=nu.trknplane1;
00907 nu.charge=nu.charge1;
00908 nu.qp=nu.qp1;
00909 nu.sigqp=nu.sigqp1;
00910 nu.chi2=nu.chi21;
00911 nu.ndof=nu.ndof1;
00912 nu.qpFraction=nu.qpFraction1;
00913 nu.trkVtxUVDiffPl=nu.trkVtxUVDiffPl1;
00914 nu.trkLength=nu.trkLength1;
00915 nu.planeTrkNu=nu.planeTrkNu1;
00916 nu.planeTrkNv=nu.planeTrkNv1;
00917 nu.ntrklike=nu.ntrklike1;
00918 nu.trkphsigcor=nu.trkphsigcor1;
00919 nu.trkphsigmap=nu.trkphsigmap1;
00920
00921
00922
00923 nu.xTrkVtx=nu.xTrkVtx1;
00924 nu.yTrkVtx=nu.yTrkVtx1;
00925 nu.zTrkVtx=nu.zTrkVtx1;
00926 nu.uTrkVtx=nu.uTrkVtx1;
00927 nu.vTrkVtx=nu.vTrkVtx1;
00928 nu.planeTrkVtx=nu.planeTrkVtx1;
00929 nu.planeTrkBeg=nu.planeTrkBeg1;
00930 nu.planeTrkBegu=nu.planeTrkBegu1;
00931 nu.planeTrkBegv=nu.planeTrkBegv1;
00932
00933 nu.xTrkEnd=nu.xTrkEnd1;
00934 nu.yTrkEnd=nu.yTrkEnd1;
00935 nu.zTrkEnd=nu.zTrkEnd1;
00936 nu.uTrkEnd=nu.uTrkEnd1;
00937 nu.vTrkEnd=nu.vTrkEnd1;
00938 nu.planeTrkEnd=nu.planeTrkEnd1;
00939 nu.planeTrkEndu=nu.planeTrkEndu1;
00940 nu.planeTrkEndv=nu.planeTrkEndv1;
00941
00942 nu.drTrkFidall=nu.drTrkFidall1;
00943 nu.dzTrkFidall=nu.dzTrkFidall1;
00944 nu.drTrkFidvtx=nu.drTrkFidvtx1;
00945 nu.dzTrkFidvtx=nu.dzTrkFidvtx1;
00946 nu.drTrkFidend=nu.drTrkFidend1;
00947 nu.dzTrkFidend=nu.dzTrkFidend1;
00948 }
00949 else if (bestTrack==2) {
00950 MAXMSG("NuReco",Msg::kInfo,3)
00951 <<"Best trk=second, primtrk="<<nu.primtrk
00952 <<", ntrk="<<nu.ntrk
00953 <<", trkExists1,2,3="<<nu.trkExists1<<","
00954 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
00955
00956 nu.trkExists=nu.trkExists2;
00957 nu.trkIndex=nu.trkIndex2;
00958 nu.ndigitTrk=nu.ndigitTrk2;
00959 nu.nstripTrk=nu.nstripTrk2;
00960 nu.trkEnCorRange=nu.trkEnCorRange2;
00961 nu.trkEnCorCurv=nu.trkEnCorCurv2;
00962 nu.trkMomentumRange=nu.trkMomentumRange2;
00963 nu.containedTrk=nu.containedTrk2;
00964 nu.trkfitpass=nu.trkfitpass2;
00965 nu.trkvtxdcosz=nu.trkvtxdcosz2;
00966 nu.trkvtxdcosy=nu.trkvtxdcosy2;
00967 nu.trknplane=nu.trknplane2;
00968 nu.charge=nu.charge2;
00969 nu.qp=nu.qp2;
00970 nu.sigqp=nu.sigqp2;
00971 nu.chi2=nu.chi22;
00972 nu.ndof=nu.ndof2;
00973 nu.qpFraction=nu.qpFraction2;
00974 nu.trkVtxUVDiffPl=nu.trkVtxUVDiffPl2;
00975 nu.trkLength=nu.trkLength2;
00976 nu.planeTrkNu=nu.planeTrkNu2;
00977 nu.planeTrkNv=nu.planeTrkNv2;
00978 nu.ntrklike=nu.ntrklike2;
00979 nu.trkphsigcor=nu.trkphsigcor2;
00980 nu.trkphsigmap=nu.trkphsigmap2;
00981
00982
00983
00984 nu.xTrkVtx=nu.xTrkVtx2;
00985 nu.yTrkVtx=nu.yTrkVtx2;
00986 nu.zTrkVtx=nu.zTrkVtx2;
00987 nu.uTrkVtx=nu.uTrkVtx2;
00988 nu.vTrkVtx=nu.vTrkVtx2;
00989 nu.planeTrkVtx=nu.planeTrkVtx2;
00990 nu.planeTrkBeg=nu.planeTrkBeg2;
00991 nu.planeTrkBegu=nu.planeTrkBegu2;
00992 nu.planeTrkBegv=nu.planeTrkBegv2;
00993
00994 nu.xTrkEnd=nu.xTrkEnd2;
00995 nu.yTrkEnd=nu.yTrkEnd2;
00996 nu.zTrkEnd=nu.zTrkEnd2;
00997 nu.uTrkEnd=nu.uTrkEnd2;
00998 nu.vTrkEnd=nu.vTrkEnd2;
00999 nu.planeTrkEnd=nu.planeTrkEnd2;
01000 nu.planeTrkEndu=nu.planeTrkEndu2;
01001 nu.planeTrkEndv=nu.planeTrkEndv2;
01002
01003 nu.drTrkFidall=nu.drTrkFidall2;
01004 nu.dzTrkFidall=nu.dzTrkFidall2;
01005 nu.drTrkFidvtx=nu.drTrkFidvtx2;
01006 nu.dzTrkFidvtx=nu.dzTrkFidvtx2;
01007 nu.drTrkFidend=nu.drTrkFidend2;
01008 nu.dzTrkFidend=nu.dzTrkFidend2;
01009 }
01010 else if (bestTrack==3) {
01011 MAXMSG("NuReco",Msg::kInfo,3)
01012 <<"Best trk=third, primtrk="<<nu.primtrk
01013 <<", ntrk="<<nu.ntrk
01014 <<", trkExists1,2,3="<<nu.trkExists1<<","
01015 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01016
01017 nu.trkExists=nu.trkExists3;
01018 nu.trkIndex=nu.trkIndex3;
01019 nu.ndigitTrk=nu.ndigitTrk3;
01020 nu.nstripTrk=nu.nstripTrk3;
01021 nu.trkEnCorRange=nu.trkEnCorRange3;
01022 nu.trkEnCorCurv=nu.trkEnCorCurv3;
01023 nu.trkMomentumRange=nu.trkMomentumRange3;
01024 nu.containedTrk=nu.containedTrk3;
01025 nu.trkfitpass=nu.trkfitpass3;
01026 nu.trkvtxdcosz=nu.trkvtxdcosz3;
01027 nu.trkvtxdcosy=nu.trkvtxdcosy3;
01028 nu.trknplane=nu.trknplane3;
01029 nu.charge=nu.charge3;
01030 nu.qp=nu.qp3;
01031 nu.sigqp=nu.sigqp3;
01032 nu.chi2=nu.chi23;
01033 nu.ndof=nu.ndof3;
01034 nu.qpFraction=nu.qpFraction3;
01035 nu.trkVtxUVDiffPl=nu.trkVtxUVDiffPl3;
01036 nu.trkLength=nu.trkLength3;
01037 nu.planeTrkNu=nu.planeTrkNu3;
01038 nu.planeTrkNv=nu.planeTrkNv3;
01039 nu.ntrklike=nu.ntrklike3;
01040 nu.trkphsigcor=nu.trkphsigcor3;
01041 nu.trkphsigmap=nu.trkphsigmap3;
01042
01043
01044
01045 nu.xTrkVtx=nu.xTrkVtx3;
01046 nu.yTrkVtx=nu.yTrkVtx3;
01047 nu.zTrkVtx=nu.zTrkVtx3;
01048 nu.uTrkVtx=nu.uTrkVtx3;
01049 nu.vTrkVtx=nu.vTrkVtx3;
01050 nu.planeTrkVtx=nu.planeTrkVtx3;
01051 nu.planeTrkBeg=nu.planeTrkBeg3;
01052 nu.planeTrkBegu=nu.planeTrkBegu3;
01053 nu.planeTrkBegv=nu.planeTrkBegv3;
01054
01055 nu.xTrkEnd=nu.xTrkEnd3;
01056 nu.yTrkEnd=nu.yTrkEnd3;
01057 nu.zTrkEnd=nu.zTrkEnd3;
01058 nu.uTrkEnd=nu.uTrkEnd3;
01059 nu.vTrkEnd=nu.vTrkEnd3;
01060 nu.planeTrkEnd=nu.planeTrkEnd3;
01061 nu.planeTrkEndu=nu.planeTrkEndu3;
01062 nu.planeTrkEndv=nu.planeTrkEndv3;
01063
01064 nu.drTrkFidall=nu.drTrkFidall3;
01065 nu.dzTrkFidall=nu.dzTrkFidall3;
01066 nu.drTrkFidvtx=nu.drTrkFidvtx3;
01067 nu.dzTrkFidvtx=nu.dzTrkFidvtx3;
01068 nu.drTrkFidend=nu.drTrkFidend3;
01069 nu.dzTrkFidend=nu.dzTrkFidend3;
01070 }
01071 else cout<<"Ahhhhhhhhhhhh"<<endl;
01072
01073
01074
01075
01076
01077
01078
01079 this->SetBestTrkMajorityCurvature(nu);
01080
01081 this->SetBestTrkSAFit(nu);
01082 }
01083
01084
01085
01086 void NuReco::SetBestTrkMajorityCurvature(NuEvent& nu) const
01087 {
01088
01089 const NuLibrary& lib=NuLibrary::Instance();
01090
01091
01092 Int_t bestTrack=lib.reco.GetBestTrack(nu);
01093
01094
01095 if (bestTrack<1) {
01096 MAXMSG("NuReco",Msg::kInfo,3)
01097 <<"SetBestTrkMajorityCurvature: No best trk, primtrk="
01098 <<nu.primtrk
01099 <<", ntrk="<<nu.ntrk
01100 <<", trkExists1,2,3="<<nu.trkExists1<<","
01101 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01102 return;
01103 }
01104 else if (bestTrack==1) {
01105 MAXMSG("NuReco",Msg::kDebug,3)
01106 <<"SetBestTrkMajorityCurvature: Best trk=first, primtrk="
01107 <<nu.primtrk
01108 <<", ntrk="<<nu.ntrk
01109 <<", trkExists1,2,3="<<nu.trkExists1<<","
01110 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01111
01112 nu.jitter=nu.jitter1;
01113 nu.jPID=nu.jPID1;
01114 nu.majC=nu.majC1;
01115
01116
01117
01118 nu.smoothMajC=nu.smoothMajC1;
01119
01120
01121 }
01122 else if (bestTrack==2) {
01123 MAXMSG("NuReco",Msg::kDebug,3)
01124 <<"SetBestTrkMajorityCurvature: Best trk=second, primtrk="
01125 <<nu.primtrk
01126 <<", ntrk="<<nu.ntrk
01127 <<", trkExists1,2,3="<<nu.trkExists1<<","
01128 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01129
01130 nu.jitter=nu.jitter2;
01131 nu.jPID=nu.jPID2;
01132 nu.majC=nu.majC2;
01133
01134
01135
01136 nu.smoothMajC=nu.smoothMajC2;
01137
01138
01139 }
01140 else if (bestTrack==3) {
01141 MAXMSG("NuReco",Msg::kDebug,3)
01142 <<"SetBestTrkMajorityCurvature: Best trk=third, primtrk="
01143 <<nu.primtrk
01144 <<", ntrk="<<nu.ntrk
01145 <<", trkExists1,2,3="<<nu.trkExists1<<","
01146 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01147
01148 nu.jitter=nu.jitter3;
01149 nu.jPID=nu.jPID3;
01150 nu.majC=nu.majC3;
01151
01152
01153
01154 nu.smoothMajC=nu.smoothMajC3;
01155
01156
01157 }
01158 else cout<<"Ahhhhhhhhhhhh"<<endl;
01159 }
01160
01161
01162
01163 void NuReco::SetBestTrkSAFit(NuEvent& nu) const
01164 {
01165
01166 const NuLibrary& lib=NuLibrary::Instance();
01167
01168
01169 Int_t bestTrack=lib.reco.GetBestTrack(nu);
01170
01171
01172 if (bestTrack<1) {
01173 MAXMSG("NuReco",Msg::kInfo,3)
01174 <<"SetBestTrkSAFit: No best trk, primtrk="
01175 <<nu.primtrk
01176 <<", ntrk="<<nu.ntrk
01177 <<", trkExists1,2,3="<<nu.trkExists1<<","
01178 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01179 return;
01180 }
01181 else if (bestTrack==1) {
01182 MAXMSG("NuReco",Msg::kDebug,3)
01183 <<"SetBestTrkSAFit: Best trk=first, primtrk="
01184 <<nu.primtrk
01185 <<", ntrk="<<nu.ntrk
01186 <<", trkExists1,2,3="<<nu.trkExists1<<","
01187 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01188
01189 nu.trkfitpassSA=nu.trkfitpassSA1;
01190 nu.trkvtxdcoszSA=nu.trkvtxdcoszSA1;
01191 nu.chargeSA=nu.chargeSA1;
01192 nu.qpSA=nu.qpSA1;
01193 nu.sigqpSA=nu.sigqpSA1;
01194 nu.chi2SA=nu.chi2SA1;
01195 nu.ndofSA=nu.ndofSA1;
01196 nu.probSA=nu.probSA1;
01197 nu.xTrkVtxSA=nu.xTrkVtxSA1;
01198 nu.yTrkVtxSA=nu.yTrkVtxSA1;
01199 nu.zTrkVtxSA=nu.zTrkVtxSA1;
01200 nu.uTrkVtxSA=nu.uTrkVtxSA1;
01201 nu.vTrkVtxSA=nu.vTrkVtxSA1;
01202 }
01203 else if (bestTrack==2) {
01204 MAXMSG("NuReco",Msg::kDebug,3)
01205 <<"SetBestTrkSAFit: Best trk=second, primtrk="
01206 <<nu.primtrk
01207 <<", ntrk="<<nu.ntrk
01208 <<", trkExists1,2,3="<<nu.trkExists1<<","
01209 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01210
01211 nu.trkfitpassSA=nu.trkfitpassSA2;
01212 nu.trkvtxdcoszSA=nu.trkvtxdcoszSA2;
01213 nu.chargeSA=nu.chargeSA2;
01214 nu.qpSA=nu.qpSA2;
01215 nu.sigqpSA=nu.sigqpSA2;
01216 nu.chi2SA=nu.chi2SA2;
01217 nu.ndofSA=nu.ndofSA2;
01218 nu.probSA=nu.probSA2;
01219 nu.xTrkVtxSA=nu.xTrkVtxSA2;
01220 nu.yTrkVtxSA=nu.yTrkVtxSA2;
01221 nu.zTrkVtxSA=nu.zTrkVtxSA2;
01222 nu.uTrkVtxSA=nu.uTrkVtxSA2;
01223 nu.vTrkVtxSA=nu.vTrkVtxSA2;
01224 }
01225 else if (bestTrack==3) {
01226 MAXMSG("NuReco",Msg::kDebug,3)
01227 <<"SetBestTrkSAFit: Best trk=third, primtrk="
01228 <<nu.primtrk
01229 <<", ntrk="<<nu.ntrk
01230 <<", trkExists1,2,3="<<nu.trkExists1<<","
01231 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01232
01233 nu.trkfitpassSA=nu.trkfitpassSA3;
01234 nu.trkvtxdcoszSA=nu.trkvtxdcoszSA3;
01235 nu.chargeSA=nu.chargeSA3;
01236 nu.qpSA=nu.qpSA3;
01237 nu.sigqpSA=nu.sigqpSA3;
01238 nu.chi2SA=nu.chi2SA3;
01239 nu.ndofSA=nu.ndofSA3;
01240 nu.probSA=nu.probSA3;
01241 nu.xTrkVtxSA=nu.xTrkVtxSA3;
01242 nu.yTrkVtxSA=nu.yTrkVtxSA3;
01243 nu.zTrkVtxSA=nu.zTrkVtxSA3;
01244 nu.uTrkVtxSA=nu.uTrkVtxSA3;
01245 nu.vTrkVtxSA=nu.vTrkVtxSA3;
01246 }
01247 else cout<<"Ahhhhhhhhhhhh"<<endl;
01248 }
01249
01250
01251
01252 Int_t NuReco::GetTrackCharge(const NtpSRTrack& trk) const
01253 {
01254 Int_t charge=this->GetTrackCharge(trk.momentum.qp);
01255 if (trk.momentum.qp==0) {
01256 MAXMSG("NuReco",Msg::kInfo,3)
01257 <<"Note: assigning charge to be -1 since trk.momentum.qp="
01258 <<trk.momentum.qp<<", trkfitpass="<<trk.fit.pass<<endl;
01259 }
01260 return charge;
01261 }
01262
01263
01264
01265 Int_t NuReco::GetTrackCharge(Double_t qp) const
01266 {
01267 Int_t charge=qp<0 ? -1 : 1;
01268 if (qp==0) {
01269 charge=-1;
01270 }
01271
01272 return charge;
01273 }
01274
01275
01276
01277 Int_t NuReco::GetPrimaryTrack(const NuEvent& nu) const
01278 {
01279 if (nu.trkExists1) return 1;
01280 else return 0;
01281 }
01282
01283
01284
01285 Int_t NuReco::GetPrimaryShowerCCStd(const NuEvent& nu) const
01286 {
01289 if (nu.nshw<=0) {
01290 MAXMSG("NuReco",Msg::kDebug,3)
01291 <<"GetPrimaryShower: No showers so return 0"<<endl;
01292 return 0;
01293 }
01294
01295 if (nu.primshw<0) {
01296
01297 return 0;
01298 }
01299
01300
01301 if (!nu.trkExists) return 1;
01302
01303 MAXMSG("NuReco",Msg::kInfo,1)
01304 <<"GetPrimaryShw: require trk vtx within 0.5 m || >2 GeV"<<endl;
01305
01306
01307
01308
01309 if (fabs(nu.zTrkVtx-nu.zShwVtx1)<(0.5*Munits::m) ||
01310 (fabs(nu.zTrkVtx-nu.zShwVtx1)>=(0.5*Munits::m) &&
01311 nu.shwEnCor1>(2*Munits::GeV))) {
01312 return 1;
01313 }
01314 else {
01315 return 0;
01316 }
01317 }
01318
01319
01320
01321 Int_t NuReco::GetPrimaryShowerStdReco(const NuEvent& nu) const
01322 {
01323 if (nu.nshw<=0) {
01324 MAXMSG("NuReco",Msg::kInfo,3)
01325 <<"GetPrimaryShower: No showers so return 0"<<endl;
01326 return 0;
01327 }
01328
01329 MAXMSG("NuReco",Msg::kInfo,1)
01330 <<"GetPrimaryShw: using evt.primshw"<<endl;
01331 if (nu.primshw<0) return 0;
01332 else {
01333 if (nu.primshw==0 && nu.shwExists1) return 1;
01334 else if (nu.primshw==1 && nu.shwExists2) {
01335 MAXMSG("NuReco",Msg::kWarning,300)
01336 <<"nu.primshw="<<nu.primshw<<endl;
01337 return 2;
01338 }
01339 else if (nu.primshw==2 && nu.shwExists3) {
01340 MAXMSG("NuReco",Msg::kWarning,300)
01341 <<"nu.primshw="<<nu.primshw<<endl;
01342 return 3;
01343 }
01344 else {
01345 cout<<"ahhhhhhhh nu.primshw="<<nu.primshw<<endl;
01346 return 0;
01347 }
01348 }
01349 }
01350
01351
01352
01353 Int_t NuReco::GetBestTrack(const NuEvent& nu) const
01354 {
01355
01356
01359
01360 if (nu.anaVersion==NuCuts::kJJH1 ||
01361 nu.anaVersion==NuCuts::kJJE1 || NuCuts::IsNMB(nu) ) {
01362 MAXMSG("NuReco",Msg::kInfo,1)
01363 <<"Setting best track as longest track"<<endl;
01364 return this->GetLongestTrack(nu);
01365 }
01366 else if (nu.anaVersion==NuCuts::kCC0093Std) {
01367 MAXMSG("NuReco",Msg::kInfo,1)
01368 <<"Setting best track as longest track"<<endl;
01369 return this->GetLongestTrack(nu);
01370 }
01371 else if (nu.anaVersion==NuCuts::kCC0250Std ||
01372 nu.anaVersion==NuCuts::kCC0325Std) {
01373 MAXMSG("NuReco",Msg::kInfo,1)
01374 <<"Setting best track as longest track"<<endl;
01375 return this->GetLongestTrack(nu);
01376 }
01377 else {
01378 MAXMSG("NuReco",Msg::kInfo,1)
01379 <<"Setting best track as primary track (from reconstruction)"
01380 <<endl;
01381 return this->GetPrimaryTrack(nu);
01382 }
01383 }
01384
01385
01386
01387 const NtpSRTrack* NuReco::GetTrackWithIndexX(const NtpStRecord& ntp,
01388 const NtpSREvent& evt,
01389 Int_t index) const
01390 {
01391 TClonesArray& trkTca=(*ntp.trk);
01392
01393 if (evt.ntrack>=index+1 && index>=0){
01394 const NtpSRTrack* trk=
01395 dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[index]]);
01396 return trk;
01397 }
01398 else {
01399 return 0;
01400 }
01401 }
01402
01403
01404
01405 const NtpSRShower* NuReco::GetShowerWithIndexX(const NtpStRecord& ntp,
01406 const NtpSREvent& evt,
01407 Int_t index) const
01408 {
01409 TClonesArray& shwTca=(*ntp.shw);
01410
01411 if (evt.nshower>=index+1 && index>=0){
01412 const NtpSRShower* shw=
01413 dynamic_cast<NtpSRShower*>(shwTca[evt.shw[index]]);
01414 return shw;
01415 }
01416 else {
01417 return 0;
01418 }
01419 }
01420
01421
01422
01423 Int_t NuReco::GetBestShower(const NuEvent& nu) const
01424 {
01427
01428 if (nu.anaVersion==NuCuts::kJJH1 ||
01429 nu.anaVersion==NuCuts::kJJE1 || NuCuts::IsNMBPQ(nu) ) {
01430 MAXMSG("NuReco",Msg::kInfo,1)
01431 <<"Setting best shower as primary shower (from reconstruction)"
01432 <<endl;
01433 return this->GetPrimaryShowerStdReco(nu);
01434 }
01435 else if (nu.anaVersion==NuCuts::kCC0093Std) {
01436 MAXMSG("NuReco",Msg::kInfo,1)
01437 <<"Setting best shower according to CC std"<<endl;
01438 return this->GetPrimaryShowerCCStd(nu);
01439 }
01440 else if (nu.anaVersion==NuCuts::kCC0250Std ||
01441 nu.anaVersion==NuCuts::kCC0325Std || NuCuts::IsNMBNQ(nu) ) {
01442 MAXMSG("NuReco",Msg::kInfo,1)
01443 <<"Setting best shower according to CC std"<<endl;
01444 return this->GetPrimaryShowerCCStd(nu);
01445 }
01446 else {
01447 MAXMSG("NuReco",Msg::kInfo,1)
01448 <<"Setting best shower as primary shower (from reconstruction)"
01449 <<endl;
01450 return this->GetPrimaryShowerStdReco(nu);
01451 }
01452 }
01453
01454
01455
01456 Bool_t NuReco::UseRange(const NuEvent& nu) const
01457 {
01458 if (nu.containmentFlag==1 ||
01459 nu.containmentFlag==3) return true;
01460 else return false;
01461 }
01462
01463
01464
01465 Bool_t NuReco::UseCurvature(const NuEvent& nu) const
01466 {
01467 if (nu.containmentFlag==2 ||
01468 nu.containmentFlag==4) return true;
01469 else return false;
01470 }
01471
01472
01473
01474 VldContext NuReco::GetVldContext(const NuEvent& nu) const
01475 {
01476
01477 VldTimeStamp time(nu.timeSec,nu.timeNanoSec);
01478
01479
01480 VldContext vc(static_cast<Detector::Detector_t>(nu.detector),
01481 static_cast<SimFlag::SimFlag_t>(nu.simFlag),time);
01482 return vc;
01483 }
01484
01485
01486
01487 void NuReco::GetEvtsPerSlc(const NtpStRecord& ntp,
01488 std::map<Int_t,Int_t>& evtsPerSlc) const
01489 {
01490 TClonesArray& evtTca=(*ntp.evt);
01491 const Int_t numEvts=evtTca.GetEntriesFast();
01492
01493
01494 for (Int_t ntpevt=0;ntpevt<numEvts;++ntpevt) {
01495 const NtpSREvent* pevt=
01496 dynamic_cast<NtpSREvent*>(evtTca[ntpevt]);
01497 const NtpSREvent& evt=(*pevt);
01498
01499
01500 evtsPerSlc[evt.slc]++;
01501 }
01502 }
01503
01504
01505
01506 void NuReco::GetEvtDeltaTs(const NtpStRecord& ntp,
01507 std::map<Int_t,Double_t>& deltaTs) const
01508 {
01509 Msg::LogLevel_t logLevel=Msg::kDebug;
01510
01511 const TClonesArray& evtTca=(*ntp.evt);
01512 const Int_t numEvts=evtTca.GetEntriesFast();
01513
01514 map<Int_t,Double_t> evtMedianTimes;
01515 multiset<Double_t> allMedianTimes;
01516 MsgFormat ffmt("%9.11f");
01517
01518
01519 for (Int_t i=0;i<numEvts;i++) {
01520 const NtpSREvent& evt=
01521 *dynamic_cast<NtpSREvent*>(evtTca[i]);
01522 Double_t medTime=this->GetEvtMedianTime(ntp,evt);
01523 evtMedianTimes[i]=medTime;
01524 allMedianTimes.insert(medTime);
01525 }
01526
01527
01528 typedef map<Int_t,Double_t>::iterator evtTimesIt;
01529 for (evtTimesIt it=evtMedianTimes.begin();
01530 it!=evtMedianTimes.end();++it){
01531
01532
01533 multiset<Double_t>::iterator tIt=allMedianTimes.find(it->second);
01534 if (tIt!=allMedianTimes.end()){
01535
01536 multiset<Double_t>::iterator tBeforeIt=tIt;
01537 --tBeforeIt;
01538 multiset<Double_t>::iterator tAfterIt=tIt;
01539 ++tAfterIt;
01540
01541 Double_t deltaAfter=-1;
01542 if (tAfterIt!=allMedianTimes.end()){
01543 MAXMSG("NuReco",logLevel,500)
01544 <<"t="<<ffmt(*tIt)
01545 <<", tAfterIt="<<ffmt(*tAfterIt)<<endl;
01546 deltaAfter=fabs((*tAfterIt)-(*tIt));
01547 }
01548
01549 Double_t deltaBefore=-1;
01550 if (tBeforeIt!=allMedianTimes.end()){
01551 MAXMSG("NuReco",logLevel,500)
01552 <<"t="<<ffmt(*tIt)
01553 <<", tBeforeIt="<<ffmt(*tBeforeIt)<<endl;
01554 deltaBefore=fabs((*tBeforeIt)-(*tIt));
01555 }
01556
01557 Double_t smallestDelta=deltaAfter;
01558 if (deltaAfter!=-1 && deltaBefore!=-1){
01559 if (deltaBefore<deltaAfter) smallestDelta=deltaBefore;
01560 }
01561 if (deltaBefore==-1){
01562 smallestDelta=deltaAfter;
01563 }
01564 if (deltaAfter==-1){
01565 smallestDelta=deltaBefore;
01566 }
01567
01568 MAXMSG("NuReco",logLevel,500)
01569 <<"smallestDelta="<<ffmt(smallestDelta)
01570 <<", deltaBefore="<<ffmt(deltaBefore)
01571 <<", deltaAfter="<<ffmt(deltaAfter)
01572 <<endl;
01573
01574
01575
01576 if (smallestDelta>=0) deltaTs[it->first]=smallestDelta;
01577 else cout<<"Ahhh bad delta T, dT="<<smallestDelta<<endl;
01578 }
01579 else cout<<"Ahh, not found time"<<endl;
01580 }
01581 }
01582
01583
01584
01585 Double_t NuReco::GetTrkQPFraction(const NtpSRTrack& trk) const
01586 {
01587 Int_t nPositive=0;
01588 Int_t nNegative=0;
01589 Int_t nZero=0;
01590
01591 for (Int_t i=0;i<trk.nstrip;i++){
01592
01593 Double_t stpfitqp=trk.stpfitqp[i];
01594
01595
01596 if (stpfitqp>0) nPositive++;
01597 else if (stpfitqp<0) nNegative++;
01598 else nZero++;
01599
01600 MAXMSG("NuReco",Msg::kDebug,1000)
01601 <<"stpfitqp="<<stpfitqp<<endl;
01602 }
01603
01604 Double_t qpFraction=-1;
01605 if (trk.nstrip>0) qpFraction=1.*nPositive/trk.nstrip;
01606
01607 if (trk.nstrip!=nPositive+nNegative) {
01608 MAXMSG("NuReco",Msg::kDebug,100)
01609 <<"??? stpfitqp: nPos="<<nPositive<<", nNeg="<<nNegative
01610 <<", nZero="<<nZero<<", sumN+P="<<nPositive+nNegative
01611 <<", trkn="<<trk.nstrip
01612 <<", qpFraction="<<qpFraction
01613 <<endl;
01614 }
01615 else {
01616 MAXMSG("NuReco",Msg::kDebug,100)
01617 <<"stpfitqp: nPos="<<nPositive<<", nNeg="<<nNegative
01618 <<", nZero="<<nZero<<", sumN+P="<<nPositive+nNegative
01619 <<", trkn="<<trk.nstrip
01620 <<", qpFraction="<<qpFraction
01621 <<endl;
01622 }
01623
01624 return qpFraction;
01625 }
01626
01627
01628
01629 Bool_t NuReco::IsTrkWithMuonFromNotNuNotPiNotKaonNotCharm
01630 (const NtpStRecord& ntp,const NtpSRTrack& trk) const
01631 {
01633 return false;
01634
01635 TClonesArray& thtrkTca=(*ntp.thtrk);
01636 const Int_t numthtrks=thtrkTca.GetEntriesFast();
01637 if (numthtrks<=0) {
01638 MAXMSG("NuReco",Msg::kDebug,1)
01639 <<"No THTracks, so can't do IsMuonFromNotNuNotPiNotCharm..."
01640 <<endl;
01641 return false;
01642 }
01643 MAXMSG("NuReco",Msg::kDebug,1)
01644 <<"Found THTrack, IsMuonFromNotNuNotPiNotCharm..."<<endl;
01645 const NtpTHTrack& thtrk=
01646 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
01647
01648
01649
01650 TClonesArray& mcTca=(*ntp.mc);
01651 const NtpMCTruth& mc=
01652 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
01653
01654 TClonesArray& hepTca=(*ntp.stdhep);
01655
01656
01657 Int_t muonEvent=-1;
01658 Bool_t isFromOther=false;
01659
01660 Int_t parent0=-999;
01661 Int_t parent1=-999;
01662
01663 Int_t parent0b=-999;
01664 Int_t parent1b=-999;
01665
01666 Int_t parent0c=-999;
01667 Int_t parent1c=-999;
01668
01669
01670 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
01671 const NtpMCStdHep& stdhep=
01672 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
01673
01674 if (abs(stdhep.IdHEP)==13){
01675 muonEvent=stdhep.mc;
01676 MAXMSG("NuReco",Msg::kDebug,3000)
01677 <<"Found muon, mc="<<muonEvent
01678 <<", ihep="<<ihep
01679 <<", id="<<stdhep.IdHEP
01680 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
01681 <<", parent=["<<stdhep.parent[0]
01682 <<","<<stdhep.parent[1]<<"]"
01683 <<", E="<<stdhep.p4[3]
01684 <<endl;
01685
01686 if (stdhep.parent[0]-stdhep.parent[1]!=0) {
01687 MSG("NuReco",Msg::kWarning)
01688 <<"Ahhhhhhhh, multiple parents of a muon"
01689 <<", parent=["<<parent0<<","<<parent1<<"]"<<endl;
01690 this->PrintStdhep(ntp,mc);
01691 }
01692
01693 if (parent0==-999) {
01694 parent0=stdhep.parent[0];
01695 parent1=stdhep.parent[1];
01696 }
01697 else if (parent0!=-999 && parent0b==-999) {
01698 MAXMSG("NuReco",Msg::kDebug,3000)
01699 <<"Found second muon, 1st parent=["
01700 <<parent0<<","<<parent1<<"]"
01701 <<", 2nd parent=["<<stdhep.parent[0]
01702 <<","<<stdhep.parent[1]<<"]"
01703 <<endl;
01704 parent0b=stdhep.parent[0];
01705 parent1b=stdhep.parent[1];
01706 }
01707 else if (parent0b!=-999 && parent0c==-999) {
01708 MAXMSG("NuReco",Msg::kDebug,3000)
01709 <<endl
01710 <<"Found third muon, parentb=["
01711 <<parent0b<<","<<parent1b<<"]"
01712 <<", new parent=["<<stdhep.parent[0]
01713 <<","<<stdhep.parent[1]<<"]"
01714 <<endl;
01715 parent0c=stdhep.parent[0];
01716 parent1c=stdhep.parent[1];
01717
01718 }
01719 else if (parent0c!=-999) {
01720 MAXMSG("NuReco",Msg::kDebug,3000)
01721 <<endl
01722 <<"Ahhh, already found >=3 muons, parentb=["
01723 <<parent0b<<","<<parent1b<<"]"
01724 <<", new parent=["<<stdhep.parent[0]
01725 <<","<<stdhep.parent[1]<<"]"
01726 <<endl;
01727
01728 }
01729 }
01730 }
01731
01732
01733
01734
01735 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
01736 const NtpMCStdHep& stdhep=
01737 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
01738
01739 if (parent0 ==(Int_t)stdhep.index || parent1 ==(Int_t)stdhep.index ||
01740 parent0b==(Int_t)stdhep.index || parent1b==(Int_t)stdhep.index) {
01741
01742 if (abs(stdhep.IdHEP)!=14) {
01743 MAXMSG("NuReco",Msg::kDebug,3000)
01744 <<endl<<endl<<endl
01745 <<"**** Found non-nu muon parent..."
01746 <<", i="<<stdhep.index
01747 <<", id="<<stdhep.IdHEP
01748 <<", parent0="<<parent0
01749 <<", parent1="<<parent1
01750 <<", parent0b="<<parent0b
01751 <<", parent1b="<<parent1b
01752 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
01753 <<endl<<endl<<endl;
01754 if (!((abs(stdhep.IdHEP)>200 && abs(stdhep.IdHEP)<400) ||
01755 (abs(stdhep.IdHEP)==411 || abs(stdhep.IdHEP)==421 ||
01756 abs(stdhep.IdHEP)==431 || abs(stdhep.IdHEP)==4122))){
01757 isFromOther=true;
01758 MAXMSG("NuReco",Msg::kInfo,3000)
01759 <<"From a ?? ?? !!!"
01760 <<" m="<<stdhep.mass
01761 <<", E="<<stdhep.p4[3]
01762 <<", id="<<stdhep.IdHEP
01763 <<endl<<endl;
01764
01765 this->PrintStdhep(ntp,mc);
01766 break;
01767 }
01768 }
01769 }
01770 }
01771
01772 if (isFromOther) return true;
01773 else return false;
01774 }
01775
01776
01777
01778 Bool_t NuReco::IsTrkWithMuonFromPionDecay(const NtpStRecord& ntp,
01779 const NtpSRTrack& trk) const
01780 {
01783
01784 TClonesArray& thtrkTca=(*ntp.thtrk);
01785 const Int_t numthtrks=thtrkTca.GetEntriesFast();
01786 if (numthtrks<=0) {
01787 MAXMSG("NuReco",Msg::kDebug,1)
01788 <<"No THTracks, so can't do IsTrkWithMuonFromPionDecay..."<<endl;
01789 return false;
01790 }
01791 MAXMSG("NuReco",Msg::kDebug,1)
01792 <<"Found THTrack, IsTrkWithMuonFromPionDecay..."<<endl;
01793 const NtpTHTrack& thtrk=
01794 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
01795
01796
01797
01798 TClonesArray& mcTca=(*ntp.mc);
01799 const NtpMCTruth& mc=
01800 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
01801
01802 TClonesArray& hepTca=(*ntp.stdhep);
01803
01804
01805 Int_t muonEvent=-1;
01806 Bool_t isFromPion=false;
01807
01808 Int_t parent0=-999;
01809 Int_t parent1=-999;
01810
01811 Int_t parent0b=-999;
01812 Int_t parent1b=-999;
01813
01814 Int_t parent0c=-999;
01815 Int_t parent1c=-999;
01816
01817
01818 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
01819 const NtpMCStdHep& stdhep=
01820 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
01821
01822 if (abs(stdhep.IdHEP)==13){
01823 muonEvent=stdhep.mc;
01824 MAXMSG("NuReco",Msg::kDebug,3000)
01825 <<"Found muon, mc="<<muonEvent
01826 <<", ihep="<<ihep
01827 <<", id="<<stdhep.IdHEP
01828 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
01829 <<", parent=["<<stdhep.parent[0]
01830 <<","<<stdhep.parent[1]<<"]"
01831 <<", E="<<stdhep.p4[3]
01832 <<endl;
01833
01834 if (stdhep.parent[0]-stdhep.parent[1]!=0) {
01835 MSG("NuReco",Msg::kWarning)
01836 <<"Ahhhhhhhh, multiple parents of a muon"
01837 <<", parent=["<<parent0<<","<<parent1<<"]"<<endl;
01838 this->PrintStdhep(ntp,mc);
01839 }
01840
01841 if (parent0==-999) {
01842 parent0=stdhep.parent[0];
01843 parent1=stdhep.parent[1];
01844 }
01845 else if (parent0!=-999 && parent0b==-999) {
01846 MAXMSG("NuReco",Msg::kDebug,3000)
01847 <<"Found second muon, 1st parent=["
01848 <<parent0<<","<<parent1<<"]"
01849 <<", 2nd parent=["<<stdhep.parent[0]
01850 <<","<<stdhep.parent[1]<<"]"
01851 <<endl;
01852 parent0b=stdhep.parent[0];
01853 parent1b=stdhep.parent[1];
01854 }
01855 else if (parent0b!=-999 && parent0c==-999) {
01856 MAXMSG("NuReco",Msg::kDebug,3000)
01857 <<endl
01858 <<"Found third muon, parentb=["
01859 <<parent0b<<","<<parent1b<<"]"
01860 <<", new parent=["<<stdhep.parent[0]
01861 <<","<<stdhep.parent[1]<<"]"
01862 <<endl;
01863 parent0c=stdhep.parent[0];
01864 parent1c=stdhep.parent[1];
01865 }
01866 else if (parent0c!=-999) {
01867 MAXMSG("NuReco",Msg::kDebug,3000)
01868 <<endl
01869 <<"Ahhh, already found >=3 muons, parentb=["
01870 <<parent0b<<","<<parent1b<<"]"
01871 <<", new parent=["<<stdhep.parent[0]
01872 <<","<<stdhep.parent[1]<<"]"
01873 <<endl;
01874
01875 }
01876 }
01877 }
01878
01879
01880
01881
01882 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
01883 const NtpMCStdHep& stdhep=
01884 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
01885
01886 if (parent0 ==(Int_t)stdhep.index || parent1 ==(Int_t)stdhep.index ||
01887 parent0b==(Int_t)stdhep.index || parent1b==(Int_t)stdhep.index) {
01888
01889 if (abs(stdhep.IdHEP)!=14) {
01890 MAXMSG("NuReco",Msg::kDebug,3000)
01891 <<endl
01892 <<"**** Found non-nu muon parent..."
01893 <<", i="<<stdhep.index
01894 <<", id="<<stdhep.IdHEP
01895 <<", parent0="<<parent0
01896 <<", parent1="<<parent1
01897 <<", parent0b="<<parent0b
01898 <<", parent1b="<<parent1b
01899 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
01900 <<endl;
01901 if (abs(stdhep.IdHEP)>200 && abs(stdhep.IdHEP)<300){
01902 isFromPion=true;
01903 MAXMSG("NuReco",Msg::kDebug,3000)
01904 <<"From a pion!!!"
01905 <<" m="<<stdhep.mass
01906 <<", E="<<stdhep.p4[3]<<endl<<endl;
01907 break;
01908 }
01909 }
01910 }
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939 }
01940
01941 if (isFromPion) return true;
01942 else return false;
01943 }
01944
01945
01946
01947 Bool_t NuReco::IsTrkWithMuonFromKaonDecay(const NtpStRecord& ntp,
01948 const NtpSRTrack& trk) const
01949 {
01952
01953 TClonesArray& thtrkTca=(*ntp.thtrk);
01954 const Int_t numthtrks=thtrkTca.GetEntriesFast();
01955 if (numthtrks<=0) {
01956 MAXMSG("NuReco",Msg::kDebug,1)
01957 <<"No THTracks, so can't do IsTrkWithMuonFromKaonDecay..."<<endl;
01958 return false;
01959 }
01960 MAXMSG("NuReco",Msg::kDebug,1)
01961 <<"Found THTrack, IsTrkWithMuonFromKaonDecay..."<<endl;
01962 const NtpTHTrack& thtrk=
01963 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
01964
01965
01966
01967 TClonesArray& mcTca=(*ntp.mc);
01968 const NtpMCTruth& mc=
01969 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
01970
01971 TClonesArray& hepTca=(*ntp.stdhep);
01972
01973
01974 Int_t muonEvent=-1;
01975 Bool_t isFromKaon=false;
01976
01977 Int_t parent0=-999;
01978 Int_t parent1=-999;
01979
01980 Int_t parent0b=-999;
01981 Int_t parent1b=-999;
01982
01983 Int_t parent0c=-999;
01984 Int_t parent1c=-999;
01985
01986
01987 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
01988 const NtpMCStdHep& stdhep=
01989 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
01990
01991 if (abs(stdhep.IdHEP)==13){
01992 muonEvent=stdhep.mc;
01993 MAXMSG("NuReco",Msg::kDebug,3000)
01994 <<"Found muon, mc="<<muonEvent
01995 <<", ihep="<<ihep
01996 <<", id="<<stdhep.IdHEP
01997 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
01998 <<", parent=["<<stdhep.parent[0]
01999 <<","<<stdhep.parent[1]<<"]"
02000 <<", E="<<stdhep.p4[3]
02001 <<endl;
02002
02003 if (stdhep.parent[0]-stdhep.parent[1]!=0) {
02004 MSG("NuReco",Msg::kWarning)
02005 <<"Ahhhhhhhh, multiple parents of a muon"
02006 <<", parent=["<<parent0<<","<<parent1<<"]"<<endl;
02007 this->PrintStdhep(ntp,mc);
02008 }
02009
02010 if (parent0==-999) {
02011 parent0=stdhep.parent[0];
02012 parent1=stdhep.parent[1];
02013 }
02014 else if (parent0!=-999 && parent0b==-999) {
02015 MAXMSG("NuReco",Msg::kDebug,3000)
02016 <<"Found second muon, 1st parent=["
02017 <<parent0<<","<<parent1<<"]"
02018 <<", 2nd parent=["<<stdhep.parent[0]
02019 <<","<<stdhep.parent[1]<<"]"
02020 <<endl;
02021 parent0b=stdhep.parent[0];
02022 parent1b=stdhep.parent[1];
02023 }
02024 else if (parent0b!=-999 && parent0c==-999) {
02025 MAXMSG("NuReco",Msg::kDebug,3000)
02026 <<endl
02027 <<"Found third muon, parentb=["
02028 <<parent0b<<","<<parent1b<<"]"
02029 <<", new parent=["<<stdhep.parent[0]
02030 <<","<<stdhep.parent[1]<<"]"
02031 <<endl;
02032 parent0c=stdhep.parent[0];
02033 parent1c=stdhep.parent[1];
02034 }
02035 else if (parent0c!=-999) {
02036 MAXMSG("NuReco",Msg::kDebug,3000)
02037 <<endl
02038 <<"Ahhh, already found >=3 muons, parentb=["
02039 <<parent0b<<","<<parent1b<<"]"
02040 <<", new parent=["<<stdhep.parent[0]
02041 <<","<<stdhep.parent[1]<<"]"
02042 <<endl;
02043
02044 }
02045 }
02046 }
02047
02048
02049
02050
02051 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02052 const NtpMCStdHep& stdhep=
02053 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02054
02055 if (parent0 ==(Int_t)stdhep.index || parent1 ==(Int_t)stdhep.index ||
02056 parent0b==(Int_t)stdhep.index || parent1b==(Int_t)stdhep.index) {
02057
02058 if (abs(stdhep.IdHEP)!=14) {
02059 MAXMSG("NuReco",Msg::kDebug,3000)
02060 <<endl
02061 <<"**** Found non-nu muon parent..."
02062 <<", i="<<stdhep.index
02063 <<", id="<<stdhep.IdHEP
02064 <<", parent0="<<parent0
02065 <<", parent1="<<parent1
02066 <<", parent0b="<<parent0b
02067 <<", parent1b="<<parent1b
02068 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02069 <<endl;
02070 if (abs(stdhep.IdHEP)>300 && abs(stdhep.IdHEP)<400){
02071 isFromKaon=true;
02072 MAXMSG("NuReco",Msg::kDebug,3000)
02073 <<"From a kaon!!!"
02074 <<" m="<<stdhep.mass
02075 <<", E="<<stdhep.p4[3]<<endl<<endl;
02076 break;
02077 }
02078 }
02079 }
02080 }
02081
02082 if (isFromKaon) return true;
02083 else return false;
02084 }
02085
02086
02087
02088 Bool_t NuReco::IsTrkWithDimuon(const NtpStRecord& ntp,
02089 const NtpSRTrack& trk) const
02090 {
02091 TClonesArray& thtrkTca=(*ntp.thtrk);
02092 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02093 if (numthtrks<=0) {
02094 MAXMSG("NuReco",Msg::kDebug,1)
02095 <<"No THTracks, so can't do IsDimuon..."<<endl;
02096 return false;
02097 }
02098 MAXMSG("NuReco",Msg::kDebug,1)
02099 <<"Found THTrack, IsDimuon..."<<endl;
02100 const NtpTHTrack& thtrk=
02101 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02102
02103
02104
02105 TClonesArray& mcTca=(*ntp.mc);
02106 const NtpMCTruth& mc=
02107 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
02108
02109 TClonesArray& hepTca=(*ntp.stdhep);
02110
02111
02112 Int_t charmEvent=-1;
02113 Bool_t isDimuon=false;
02114
02115 Int_t child0=-999;
02116 Int_t child1=-999;
02117
02118
02119 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02120 const NtpMCStdHep& stdhep=
02121 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02122
02123 if (abs(stdhep.IdHEP)==411 || abs(stdhep.IdHEP)==421 ||
02124 abs(stdhep.IdHEP)==431 || abs(stdhep.IdHEP)==4122){
02125 charmEvent=stdhep.mc;
02126 MAXMSG("NuReco",Msg::kDebug,3000)
02127 <<"Found Charm particle, mc="<<charmEvent
02128 <<", ihep="<<ihep
02129 <<", id="<<stdhep.IdHEP
02130 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02131 <<endl;
02132
02133 child0=stdhep.child[0];
02134 child1=stdhep.child[1];
02135 }
02136 else if ((abs(stdhep.IdHEP)>4000 && abs(stdhep.IdHEP)<5000) ||
02137 (abs(stdhep.IdHEP)>400 && abs(stdhep.IdHEP)<500)) {
02138 MAXMSG("NuReco",Msg::kInfo,3000)
02139 <<endl<<endl
02140 <<"What is this particle??? Charm?, mc="<<charmEvent
02141 <<", ihep="<<ihep
02142 <<", id="<<stdhep.IdHEP
02143 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02144 <<endl<<endl<<endl;
02145 }
02146
02147 if (child0==(Int_t)stdhep.index || child1==(Int_t)stdhep.index) {
02148 MAXMSG("NuReco",Msg::kDebug,3000)
02149 <<"Found Charm child..."
02150 <<", i="<<stdhep.index
02151 <<", child0="<<child0
02152 <<", child1="<<child1
02153 <<", id="<<stdhep.IdHEP
02154 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02155 <<endl;
02156 if (abs(stdhep.IdHEP)==13){
02157 isDimuon=true;
02158 MAXMSG("NuReco",Msg::kDebug,3000)
02159 <<"It's a muon!!!"
02160 <<" m="<<stdhep.mass
02161 <<", E="<<stdhep.p4[3]<<endl<<endl;
02162 break;
02163 }
02164 }
02165 }
02166
02167 if (isDimuon) return true;
02168 else return false;
02169 }
02170
02171
02172
02173 Bool_t NuReco::PrintCharm(const NtpStRecord& ntp,
02174 const NtpSREvent& evt) const
02175 {
02176 TClonesArray& thevtTca=(*ntp.thevt);
02177 const Int_t numthevts=thevtTca.GetEntriesFast();
02178 if (numthevts<=0) {
02179 MAXMSG("NuReco",Msg::kDebug,1)
02180 <<"No THEvents, so can't GetTruthInfo..."<<endl;
02181 return false;
02182 }
02183 MAXMSG("NuReco",Msg::kDebug,1)
02184 <<"Found THEvent, GetTruthInfo..."<<endl;
02185 const NtpTHEvent& thevt=
02186 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
02187
02188
02189
02190 TClonesArray& mcTca=(*ntp.mc);
02191 const NtpMCTruth& mc=
02192 *(dynamic_cast<NtpMCTruth*>(mcTca[thevt.neumc]));
02193
02194 TClonesArray& hepTca=(*ntp.stdhep);
02195
02196
02197 Int_t charmEvent=-1;
02198
02199 Int_t child0=-999;
02200 Int_t child1=-999;
02201
02202
02203 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02204 const NtpMCStdHep& stdhep=
02205 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02206
02207 if (abs(stdhep.IdHEP)==411 || abs(stdhep.IdHEP)==421 ||
02208 abs(stdhep.IdHEP)==431 || abs(stdhep.IdHEP)==4122){
02209 charmEvent=stdhep.mc;
02210 MAXMSG("NuReco",Msg::kDebug,3000)
02211 <<"Found Charm particle, mc="<<charmEvent
02212 <<", ihep="<<ihep
02213 <<", id="<<stdhep.IdHEP
02214 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02215 <<endl;
02216
02217 child0=stdhep.child[0];
02218 child1=stdhep.child[1];
02219 }
02220
02221 if (child0==(Int_t)stdhep.index || child1==(Int_t)stdhep.index) {
02222 MAXMSG("NuReco",Msg::kDebug,3000)
02223 <<"Found Charm child..."
02224 <<", i="<<stdhep.index
02225 <<", child0="<<child0
02226 <<", child1="<<child1
02227 <<", id="<<stdhep.IdHEP
02228 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02229 <<endl;
02230 if (abs(stdhep.IdHEP)==13){
02231 MAXMSG("NuReco",Msg::kDebug,3000)
02232 <<"It's a muon!!!"
02233 <<" m="<<stdhep.mass
02234 <<", E="<<stdhep.p4[3]<<endl<<endl;
02235 }
02236 }
02237
02238 }
02239
02240
02241 if (charmEvent!=-1){
02242 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02243 const NtpMCStdHep& stdhep=
02244 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02245
02246 MAXMSG("NuReco",Msg::kInfo,3000)
02247 <<"ihep="<<ihep
02248 <<", mc="<<stdhep.mc
02249 <<", id="<<stdhep.IdHEP
02250 <<", Ist="<<stdhep.IstHEP
02251 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02252 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02253 <<", m="<<stdhep.mass
02254 <<", E="<<stdhep.p4[3]
02255 <<endl;
02256 }
02257 return true;
02258 }
02259 else return false;
02260 }
02261
02262
02263
02264 void NuReco::PrintStdhep(const NtpStRecord& ntp,
02265 const NtpMCTruth& mc) const
02266 {
02267 TClonesArray& hepTca=(*ntp.stdhep);
02268
02269
02270 MAXMSG("NuReco",Msg::kInfo,3000)
02271 <<"PrintStdhep: mc.stdhep=["<<mc.stdhep[0]
02272 <<","<<mc.stdhep[1]<<"]"<<endl;
02273
02274
02275 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02276 const NtpMCStdHep& stdhep=
02277 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02278
02279 MAXMSG("NuReco",Msg::kInfo,3000)
02280 <<"ihep="<<ihep
02281 <<", mc="<<stdhep.mc
02282 <<", id="<<stdhep.IdHEP
02283 <<", Ist="<<stdhep.IstHEP
02284 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02285 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02286 <<", m="<<stdhep.mass
02287 <<", E="<<stdhep.p4[3]
02288 <<endl;
02289 }
02290 }
02291
02292
02293
02294 Double_t NuReco::GetEvtMedianTime(const NtpStRecord& ntp,
02295 const NtpSREvent& evt) const
02296 {
02297
02298
02299 multiset<Double_t> times;
02300 const TClonesArray& stpTca=(*ntp.stp);
02301
02302
02303 MAXMSG("NuReco",Msg::kDebug,200)
02304 <<"evt.nstrip="<<evt.nstrip<<endl;
02305 for (Int_t i=0;i<evt.nstrip;i++) {
02306 const NtpSRStrip& stp=
02307 *(dynamic_cast<NtpSRStrip*>(stpTca[evt.stp[i]]));
02308
02309 Double_t time=stp.time0;
02310 if (time<0 || time>1) time=stp.time1;
02311
02312 if (time>0 && time<=1) {
02313 times.insert(time);
02314 }
02315
02316
02317 MAXMSG("NuReco",Msg::kDebug,500)
02318 <<" Strip index="<<stp.index<<", stp="<<stp.strip
02319 <<", t="<<time
02320 <<", pl="<<stp.plane
02321 <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
02322 }
02323
02324
02325 multiset<Double_t>::const_iterator it=times.begin();
02326 advance(it,times.size()/2);
02327 Double_t medianTime=*it;
02328 MAXMSG("NuReco",Msg::kDebug,100)
02329 <<"Median time="<<medianTime<<endl;
02330
02331 return medianTime;
02332 }
02333
02334
02335
02336 void NuReco::GetBestTruthIndex(const NtpStRecord& ntp,
02337 const NtpSREvent& evt,NuEvent& nu) const
02338 {
02339 TClonesArray& thevtTca=(*ntp.thevt);
02340 const Int_t numthevts=thevtTca.GetEntriesFast();
02341 if (numthevts<=0) {
02342 MAXMSG("NuReco",Msg::kDebug,1)
02343 <<"No THEvents, so can't GetBestTruthIndex..."<<endl;
02344 return;
02345 }
02346 MAXMSG("NuReco",Msg::kDebug,1)
02347 <<"Found THEvent, GetBestTruthIndex..."<<endl;
02348
02349 Int_t mcTrk=-1;
02350 Int_t mcShw=-1;
02351 Int_t mcEvt=-1;
02352
02353 const NtpTHEvent& thevt=
02354 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
02355 mcEvt=thevt.neumc;
02356
02357
02358
02359
02360
02361 TClonesArray& thshwTca=(*ntp.thshw);
02362 const Int_t numthshws=thshwTca.GetEntriesFast();
02363 if (numthshws>0) {
02364 Int_t bestShower=this->GetBestShower(nu);
02365 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,
02366 bestShower-1);
02367 if (pshw) {
02368 const NtpSRShower& shw=*pshw;
02369 const NtpTHShower& thshw=
02370 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
02371 mcShw=thshw.neumc;
02372 }
02373 }
02374
02375
02376 TClonesArray& thtrkTca=(*ntp.thtrk);
02377 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02378 if (numthtrks>0) {
02379 Int_t bestTrack=this->GetBestTrack(nu);
02380 const NtpSRTrack* ptrk=this->GetTrackWithIndexX(ntp,evt,
02381 bestTrack-1);
02382 if (ptrk) {
02383 const NtpSRTrack& trk=*ptrk;
02384 const NtpTHTrack& thtrk=
02385 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02386 mcTrk=thtrk.neumc;
02387 }
02388 }
02389
02390
02391 nu.mcTrk=mcTrk;
02392 nu.mcShw=mcShw;
02393 nu.mcEvt=mcEvt;
02394 this->ChooseTruthIndexToUse(nu);
02395 }
02396
02397
02398
02399 void NuReco::GetPrimaryTruthIndex(const NtpStRecord& ntp,
02400 const NtpSREvent& evt,NuEvent& nu) const
02401 {
02402 TClonesArray& thevtTca=(*ntp.thevt);
02403 const Int_t numthevts=thevtTca.GetEntriesFast();
02404 if (numthevts<=0) {
02405 MAXMSG("NuReco",Msg::kDebug,1)
02406 <<"No THEvents, so can't GetPrimaryTruthIndex..."<<endl;
02407 return;
02408 }
02409 MAXMSG("NuReco",Msg::kDebug,1)
02410 <<"Found THEvent, GetPrimaryTruthIndex..."<<endl;
02411
02412 Int_t mcTrk=-1;
02413 Int_t mcShw=-1;
02414 Int_t mcEvt=-1;
02415
02416 const NtpTHEvent& thevt=
02417 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
02418 mcEvt=thevt.neumc;
02419
02420
02421
02422
02423
02424 TClonesArray& thshwTca=(*ntp.thshw);
02425 const Int_t numthshws=thshwTca.GetEntriesFast();
02426 if (numthshws>0) {
02427 Int_t bestShower=1;
02428 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,
02429 bestShower-1);
02430 if (pshw) {
02431 const NtpSRShower& shw=*pshw;
02432 const NtpTHShower& thshw=
02433 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
02434 mcShw=thshw.neumc;
02435 }
02436 }
02437
02438
02439 TClonesArray& thtrkTca=(*ntp.thtrk);
02440 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02441 if (numthtrks>0) {
02442 Int_t bestTrack=1;
02443 const NtpSRTrack* ptrk=this->GetTrackWithIndexX(ntp,evt,
02444 bestTrack-1);
02445 if (ptrk) {
02446 const NtpSRTrack& trk=*ptrk;
02447 const NtpTHTrack& thtrk=
02448 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02449 mcTrk=thtrk.neumc;
02450 }
02451 }
02452
02453
02454 nu.mcTrk=mcTrk;
02455 nu.mcShw=mcShw;
02456 nu.mcEvt=mcEvt;
02457 this->ChooseTruthIndexToUse(nu);
02458 }
02459
02460
02461
02462 void NuReco::ChooseTruthIndexToUse(NuEvent& nu) const
02463 {
02464
02465
02466
02467
02468
02469
02470
02471
02472
02473
02474
02475 if (nu.mcTrk>=0) nu.mc=nu.mcTrk;
02476 else nu.mc=nu.mcShw;
02477 }
02478
02479
02480
02481 void NuReco::GetTruthInfo(const NtpStRecord& ntp,
02482 const NtpSREvent& evt,NuEvent& nu) const
02483 {
02484
02485
02486 this->GetPrimaryTruthIndex(ntp,evt,nu);
02487
02488
02489
02490
02491 if (nu.mc<0) {
02492 MAXMSG("NuReco",Msg::kInfo,1)
02493 <<"Not Getting TruthInfo for simFlag="<<nu.simFlag<<endl;
02494 return;
02495 }
02496 MAXMSG("NuReco",Msg::kInfo,1)
02497 <<"Getting TruthInfo for simFlag="<<nu.simFlag<<endl;
02498
02499
02500
02501 TClonesArray& mcTca=(*ntp.mc);
02502 const NtpMCTruth& mc=
02503 *(dynamic_cast<NtpMCTruth*>(mcTca[nu.mc]));
02504
02505
02506 this->GetTruthInfo(ntp,mc,nu);
02507
02508 if (mc.iaction==1){
02509 MAXMSG("NuReco",Msg::kDebug,300)
02510 <<"y="<<mc.y<<endl
02511 <<"p4neu[3]="<<mc.p4neu[3]<<", p4tgt[3]="<<mc.p4tgt[3]
02512 <<", sum="<<mc.p4neu[3]+mc.p4tgt[3]
02513 <<endl
02514 <<"p4mu1[3]="<<mc.p4mu1[3]<<", mc.p4shw[3]="<<mc.p4shw[3]
02515 <<", sum="<<fabs(mc.p4mu1[3])+mc.p4shw[3]
02516 <<", labY="<<mc.p4shw[3]/(fabs(mc.p4mu1[3])+mc.p4shw[3])
02517 <<endl;
02518 }
02519 else if (mc.iaction==0){
02520 MAXMSG("NuReco",Msg::kDebug,300)
02521 <<"NC: y="<<mc.y<<endl
02522 <<"p4neu[3]="<<mc.p4neu[3]
02523 <<", p4tgt[3]="<<mc.p4tgt[3]
02524 <<endl
02525 <<"p4mu1[3]="<<mc.p4mu1[3]<<", p4shw[3]="<<mc.p4shw[3]
02526 <<endl;
02527 }
02528 }
02529
02530
02531
02532 void NuReco::GetTruthInfo(const NtpStRecord& ntp,
02533 const NtpMCTruth& mc,NuEvent& nu) const
02534 {
02535
02536 Double_t trueEn=mc.y*mc.p4neu[3];
02537 if (mc.iaction==1) trueEn=mc.p4neu[3];
02538
02539 Double_t muEn=(1-mc.y)*mc.p4neu[3];
02540 if (mc.iaction==0) muEn=0;
02541 Double_t hadEn=mc.y*mc.p4neu[3];
02542
02543 nu.energyMC=trueEn;
02544 nu.neuEnMC=mc.p4neu[3];
02545 nu.neuPxMC=mc.p4neu[0];
02546 nu.neuPyMC=mc.p4neu[1];
02547 nu.neuPzMC=mc.p4neu[2];
02548 nu.tgtEnMC=mc.p4tgt[3];
02549 nu.tgtPxMC=mc.p4tgt[0];
02550 nu.tgtPyMC=mc.p4tgt[1];
02551 nu.tgtPzMC=mc.p4tgt[2];
02552 nu.zMC=static_cast<Int_t>(mc.z);
02553 nu.aMC=static_cast<Int_t>(mc.a);
02554 nu.yMC=mc.y;
02555 nu.y2MC=mc.p4shw[3]/(fabs(mc.p4mu1[3])+mc.p4shw[3]);
02556 nu.xMC=mc.x;
02557 nu.q2MC=mc.q2;
02558 nu.w2MC=mc.w2;
02559 nu.trkEnMC=mc.p4mu1[3];
02560 nu.trkEn2MC=muEn;
02561 nu.shwEnMC=mc.p4shw[3];
02562 nu.shwEn2MC=hadEn;
02563 nu.iaction=mc.iaction;
02564 nu.iresonance=mc.iresonance;
02565 nu.inu=mc.inu;
02566 nu.inunoosc=mc.inunoosc;
02567 nu.itg=mc.itg;
02568
02569 nu.vtxxMC=mc.vtxx;
02570 nu.vtxyMC=mc.vtxy;
02571 nu.vtxzMC=mc.vtxz;
02572
02573 nu.tpx=mc.flux.tpx;
02574 nu.tpy=mc.flux.tpy;
02575 nu.tpz=mc.flux.tpz;
02576 nu.tptype=mc.flux.tptype;
02577 nu.ppvx=mc.flux.ppvx;
02578 nu.ppvy=mc.flux.ppvy;
02579 nu.ppvz=mc.flux.ppvz;
02580 nu.ptype=mc.flux.ptype;
02581 nu.tgen=mc.flux.tgen;
02582
02583 this->CalcExtraTruthVariables(nu);
02584
02585
02586
02587 nu.nucleusMC=this->GetNucleus(ntp,nu);
02588 nu.initialStateMC=this->GetInitialState(ntp,nu);
02589 nu.hadronicFinalStateMC=this->GetHadronicFinalState(ntp,nu);
02590 }
02591
02592
02593
02594 void NuReco::PrintTrueEnergy(const NtpStRecord& ntp,
02595 const NtpSREvent& evt,Double_t recoEn) const
02596 {
02597 TClonesArray& thevtTca=(*ntp.thevt);
02598 const Int_t numthevts=thevtTca.GetEntriesFast();
02599 if (numthevts<=0) {
02600 MAXMSG("NuReco",Msg::kInfo,1)
02601 <<"No THEvents, so can't print truth information..."<<endl;
02602 MAXMSG("NuReco",Msg::kInfo,20)
02603 <<"recoEn="<<recoEn<<endl;
02604 return;
02605 }
02606 MAXMSG("NuReco",Msg::kInfo,1)
02607 <<"Found THEvent, printing info..."<<endl;
02608 const NtpTHEvent& thevt=
02609 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
02610
02611
02612 TClonesArray& mcTca=(*ntp.mc);
02613 const NtpMCTruth& mc=
02614 *(dynamic_cast<NtpMCTruth*>(mcTca[thevt.neumc]));
02615 MAXMSG("NuReco",Msg::kDebug,100)
02616 <<"Found mc="<<mc.index<<", thevt.index="<<thevt.index<<endl;
02617
02618 Double_t nuEn=mc.p4neu[3];
02619
02620 Double_t trueEn=mc.y*mc.p4neu[3];
02621 if (mc.iaction==1) trueEn=mc.p4neu[3];
02622
02623 Double_t muEn=(1-mc.y)*mc.p4neu[3];
02624 Double_t hadEn=mc.y*mc.p4neu[3];
02625
02626 string siaction="NC";
02627 if (mc.iaction==1) siaction="CC";
02628
02629 if (mc.iaction==1) {
02630 MAXMSG("NuReco",Msg::kInfo,100)
02631 <<"RecoEn="<<recoEn
02632 <<", "<<siaction
02633 <<" truth: En="<<trueEn<<", mu="<<muEn<<", had="<<hadEn
02634 <<", (y="<<mc.y<<")"
02635 <<endl;
02636 }
02637 else {
02638 MAXMSG("NuReco",Msg::kInfo,20)
02639 <<"RecoEn="<<recoEn
02640 <<", "<<siaction
02641 <<" truth: had="<<hadEn
02642 <<", (Nu="<<nuEn<<", y="<<mc.y<<")"
02643 <<endl;
02644 }
02645 }
02646
02647
02648
02649
02650
02651
02652
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686
02687
02688
02689
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706
02707 Int_t NuReco::GetContainmentFlag(const NuEvent& nu) const
02708 {
02709 Int_t flag=-1;
02710
02711
02712 if (nu.anaVersion==NuCuts::kCC0093Std) {
02713 MAXMSG("NuReco",Msg::kInfo,1)
02714 <<"Using CC0093Std containment criteria"<<endl;
02715 flag=nu.containmentFlagCC0093Std;
02716 }
02717 else if (nu.anaVersion==NuCuts::kCC0250Std) {
02718 MAXMSG("NuReco",Msg::kInfo,1)
02719 <<"Using CC0250Std containment criteria"<<endl;
02720 flag=nu.containmentFlagCC0250Std;
02721 }
02722 else if (nu.anaVersion==NuCuts::kCC0325Std || NuCuts::IsNMBNQ(nu)) {
02723 MAXMSG("NuReco",Msg::kInfo,1)
02724 <<"Using CC0325Std hybrid containment criteria"<<endl;
02725 flag=nu.containmentFlagCC0250Std;
02726
02727
02728 if (nu.detector==Detector::kNear) {
02729 if (nu.xTrkEnd>1.3*Munits::m) {
02730 flag=this->GetContainmentFlagStdReco(nu);
02731 }
02732 }
02733 }
02734 else {
02735 if (ReleaseType::IsBirch(nu.recoVersion)) {
02736 MAXMSG("NuReco",Msg::kInfo,1)
02737 <<"Using Pitt containment criteria"<<endl;
02738
02739 flag=nu.containmentFlagPitt;
02740 }
02741 else {
02742 MAXMSG("NuReco",Msg::kInfo,1)
02743 <<"Using std reco containment criteria"<<endl;
02744 flag=this->GetContainmentFlagStdReco(nu);
02745
02746
02747
02748 }
02749 }
02750
02751
02752 if (!(flag==1 || flag==2 || flag==3 || flag==4)) {
02753 MSG("NuReco",Msg::kWarning)<<"flag="<<flag<<endl;
02754 }
02755
02756
02757 return flag;
02758 }
02759
02760
02761
02762 Int_t NuReco::GetContainmentFlagStdReco(const NuEvent& nu) const
02763 {
02764 Int_t flag=-1;
02765
02766
02767
02768
02769
02770
02771 if (nu.detector==Detector::kNear) {
02772 if (nu.containedTrk) {
02773
02774 if (nu.zTrkEnd<7.0) flag=1;
02775 else flag=3;
02776 }
02777 else {
02778
02779 if (nu.zTrkEnd<7.0) flag=2;
02780 else flag=4;
02781 }
02782 }
02783 else if (nu.detector==Detector::kFar) {
02784
02785 if (nu.containedTrk) flag=1;
02786 else flag=2;
02787 }
02788 else cout<<"Ahh, detector="<<nu.detector<<endl;
02789
02790 return flag;
02791 }
02792
02793
02794
02795 Int_t NuReco::GetContainmentFlagCC0250Std(const NuEvent& nu) const
02796 {
02797 const Float_t x=nu.xTrkEnd;
02798 const Float_t y=nu.yTrkEnd;
02799 const Float_t z=nu.zTrkEnd;
02800 const Float_t r=sqrt(x*x+y*y);
02801
02802
02803
02804 const Float_t xMin=-1.65*(Munits::m);
02805 const Float_t xMax=+2.7*(Munits::m);
02806 const Float_t yMin=-1.65*(Munits::m);
02807 const Float_t yMax=+1.65*(Munits::m);
02808
02809 const Float_t uMax=+3.55*(Munits::m);
02810
02811
02812 const Float_t zMax=+15*(Munits::m);
02813
02814
02815 const Float_t planeMin=4;
02816 const Float_t planeMax=475;
02817 const Float_t rMax=TMath::Sqrt(14)*(Munits::m);
02818
02819 Bool_t contained=false;
02820 Int_t flag=-1;
02821
02822 if (nu.detector==Detector::kNear) {
02823 contained=z<zMax &&
02824 x>xMin && x<xMax &&
02825 y>yMin && y<yMax &&
02826 y>-x+xMin &&
02827 y<+x-xMin &&
02828 y<-x+uMax &&
02829 y>+x-uMax;
02830
02831 if (contained) {
02832
02833
02834
02835
02836
02837
02838 if (z<7) flag=1;
02839 else if (z>=7) flag=3;
02840 }
02841 else {
02842 if (z<7) flag=2;
02843 else if (z>=7) flag=4;
02844 }
02845 }
02846 else if (nu.detector==Detector::kFar) {
02847 contained=nu.planeTrkEnd>=planeMin &&
02848 nu.planeTrkEnd<=planeMax &&
02849 r<rMax;
02850 if (contained) flag=1;
02851 else flag=2;
02852 }
02853 else cout<<"Ahh, detector="<<nu.detector<<endl;
02854
02855 return flag;
02856 }
02857
02858
02859
02860 Int_t NuReco::GetContainmentFlagCC0093Std(const NuEvent& nu) const
02861 {
02862 const Float_t x=nu.xTrkEnd;
02863 const Float_t y=nu.yTrkEnd;
02864 const Float_t z=nu.zTrkEnd;
02865 const Float_t r=sqrt(x*x+y*y);
02866
02867 const Float_t xMin=-1.65;
02868 const Float_t xMax=+2.7;
02869 const Float_t yMin=-1.65;
02870 const Float_t yMax=+1.65;
02871 const Float_t zMax=+16;
02872 const Float_t rMin=+0.4;
02873 const Float_t rMinFD=+0.5;
02874
02875 Bool_t contained=false;
02876 Int_t flag=-1;
02877
02878
02879
02880
02881
02882
02883
02884
02885
02886
02887
02888
02889
02890
02891
02892
02893
02894
02895
02896
02897 if (nu.detector==Detector::kNear) {
02898
02899
02900 contained=z<zMax && r>rMin &&
02901 x>xMin && x<xMax &&
02902 y>yMin && y<yMax &&
02903 y>-x+xMin &&
02904 y<+x-xMin &&
02905 y<-x+3.55 &&
02906 y>+x-3.55;
02907
02908 if (contained) {
02909
02910
02911
02912
02913
02914
02915 if (z<7) flag=1;
02916 else if (z>=7) flag=3;
02917 }
02918 else {
02919 if (z<7) flag=2;
02920 else if (z>=7) flag=4;
02921 }
02922 }
02923 else if (nu.detector==Detector::kFar) {
02924 contained=!(nu.drTrkFidall<rMinFD || nu.dzTrkFidall<rMinFD);
02925 if (contained) flag=1;
02926 else flag=2;
02927 }
02928 else cout<<"Ahh, detector="<<nu.detector<<endl;
02929
02930 return flag;
02931 }
02932
02933
02934
02935 Int_t NuReco::GetContainmentFlagPitt(const NuEvent& nu) const
02936 {
02938
02944
02948
02949 Int_t containmentFlagPitt=0;
02950
02951 const Float_t x=nu.xTrkEnd;
02952 const Float_t y=nu.yTrkEnd;
02953 const Float_t z=nu.zTrkEnd;
02954 const Float_t u=nu.uTrkEnd;
02955 const Float_t v=nu.vTrkEnd;
02956 const Float_t rad=sqrt(x*x + y*y);
02957
02958
02959
02960
02961
02962
02963
02964
02965
02966
02967
02968
02969
02970
02971
02972
02973
02974
02975
02976
02977
02978
02979
02980
02981
02982
02983
02984
02985
02986
02987
02988
02989
02990
02991
02992
02993
02994
02995
02996
02997
02998
02999
03000
03001
03002
03003
03004
03005
03006
03007
03008
03009 if (nu.detector==Detector::kNear) {
03010 Float_t calZ=7*(Munits::m);
03011 Float_t specZ=15.6*(Munits::m);
03012 Float_t minR=0.8*(Munits::m);
03013
03014 Bool_t down_stop=false;
03015 Bool_t down_exit=false;
03016 Bool_t up_stop=false;
03017 Bool_t up_exit=false;
03018
03019
03020 if (z>calZ) {
03021 if (u>-0.85 && u<2.1 &&
03022 v>-2.1 && v<0.85 &&
03023 x>-1.5 && x<2.4 &&
03024 y>-1.4 && y<1.4 &&
03025 rad>minR &&
03026 z>calZ && z<=specZ) down_stop=true;
03027 else down_exit=true;
03028 }
03029
03030 else if (z<=calZ && z>=0) {
03031 if (u>0.3 && u<1.8 &&
03032 v>-1.8 && v<-0.3 &&
03033 x<2.4 &&
03034 rad>minR &&
03035 z<calZ && z>=0) up_stop=true;
03036 else up_exit=true;
03037 }
03038 else cout<<"Ahhh, z="<<z<<endl;
03039
03040
03041 if (down_stop) containmentFlagPitt = 3;
03042 if (up_stop ) containmentFlagPitt = 1;
03043 if (down_exit) containmentFlagPitt = 4;
03044 if (up_exit ) containmentFlagPitt = 2;
03045 }
03046 else if (nu.detector==Detector::kFar) {
03047
03048 const Float_t rMinFD=+0.5;
03049 Bool_t contained=!(nu.drTrkFidall<rMinFD || nu.dzTrkFidall<rMinFD);
03050 if (contained) containmentFlagPitt=1;
03051 else containmentFlagPitt=2;
03052 }
03053 else cout<<"Ahhh, detector="<<nu.detector<<endl;
03054
03055
03056 if (!(containmentFlagPitt==1 || containmentFlagPitt==2 ||
03057 containmentFlagPitt==3 || containmentFlagPitt==4)) {
03058 MSG("NuReco",Msg::kWarning)
03059 <<"containmentFlagPitt="<<containmentFlagPitt<<endl;
03060 }
03061
03062 return containmentFlagPitt;
03063 }
03064
03065
03066
03067
03068
03069
03070
03071
03072
03073
03074
03075
03076
03077
03078
03079
03080
03081
03082
03083
03084
03085
03086
03087
03088
03089
03090
03091
03092
03093
03094
03095
03096
03097
03098
03099
03100
03101
03102
03103
03104
03105
03106
03107
03108
03109
03110
03111
03112
03113
03114
03115
03116
03117
03118
03119
03120 Int_t NuReco::GetInitialState(const NtpStRecord& ntp,
03121 const NuEvent& nu) const
03122 {
03124
03125
03126
03127
03128 Int_t itr=nu.mc;
03129
03130 Int_t initial_state=0;
03131 TClonesArray* pointStdhepArray = NULL;
03132 pointStdhepArray=ntp.stdhep;
03133 TClonesArray& stdhepArray = *pointStdhepArray;
03134 Int_t nStdHep = stdhepArray.GetEntries();
03135
03136 int protneut = -1;
03137 int nubarnu = 0;
03138
03139 for(int i=0;i<nStdHep;i++){
03140
03141 const NtpMCStdHep* ntpStdHep=
03142 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
03143 if(ntpStdHep->mc==itr){
03144
03145 if(ntpStdHep->IstHEP==0){
03146 if(abs(ntpStdHep->IdHEP)==12 ||
03147 abs(ntpStdHep->IdHEP)==14 ||
03148 abs(ntpStdHep->IdHEP)==16){
03149 nubarnu = ntpStdHep->IdHEP/abs(ntpStdHep->IdHEP);
03150 }
03151 }
03152 if(ntpStdHep->IstHEP==11){
03153 if(ntpStdHep->IdHEP==2212) protneut = 1;
03154 else if(ntpStdHep->IdHEP==2112) protneut = 0;
03155 else if(abs(ntpStdHep->IdHEP)>1000000000) protneut = 2;
03156 else protneut = 3;
03157 }
03158 }
03159 }
03160
03161 if(protneut==1 && nubarnu==1) initial_state=1;
03162 if(protneut==0 && nubarnu==1) initial_state=2;
03163 if(protneut==1 && nubarnu==-1) initial_state=3;
03164 if(protneut==0 && nubarnu==-1) initial_state=4;
03165 if(protneut==2 && nubarnu==1) initial_state=5;
03166 if(protneut==3 && nubarnu==1) initial_state=6;
03167 if(protneut==2 && nubarnu==-1) initial_state=7;
03168 if(protneut==3 && nubarnu==-1) initial_state=8;
03169
03170 return initial_state;
03171 }
03172
03173
03174
03175 Int_t NuReco::GetNucleus(const NtpStRecord& ,
03176 const NuEvent& nu) const
03177 {
03179
03180
03181 Int_t z=nu.zMC;
03182 Int_t a=nu.aMC;
03183 Int_t nucleus=1;
03184
03185 switch (z) {
03186
03187
03188
03189 case 1:
03190 switch (a) {
03191 case 1:
03192 nucleus=256;
03193 break;
03194 case 2:
03195 nucleus=257;
03196 break;
03197 case 3:
03198 nucleus=258;
03199 break;
03200 default:
03201 nucleus=256;
03202 break;
03203 }
03204 break;
03205 case 6:
03206 switch (a) {
03207 case 11:
03208 nucleus=274;
03209 break;
03210 case 12:
03211 nucleus=275;
03212 break;
03213 case 13:
03214 nucleus=276;
03215 break;
03216 case 14:
03217 nucleus=277;
03218 break;
03219 default:
03220 nucleus=275;
03221 break;
03222 }
03223 break;
03224 case 7:
03225 switch (a) {
03226 case 13:
03227 nucleus=278;
03228 break;
03229 case 14:
03230 nucleus=279;
03231 break;
03232 case 15:
03233 nucleus=280;
03234 break;
03235 case 16:
03236 nucleus=281;
03237 break;
03238 case 17:
03239 nucleus=282;
03240 break;
03241 default:
03242 nucleus=279;
03243 break;
03244 }
03245 break;
03246 case 8:
03247 switch (a) {
03248 case 15:
03249 nucleus=283;
03250 break;
03251 case 16:
03252 nucleus=284;
03253 break;
03254 case 17:
03255 nucleus=285;
03256 break;
03257 case 18:
03258 nucleus=286;
03259 break;
03260 default:
03261 nucleus=284;
03262 break;
03263 }
03264 break;
03265 case 13:
03266 switch (a) {
03267 case 26:
03268 nucleus=303;
03269 break;
03270 case 27:
03271 nucleus=304;
03272 break;
03273 case 28:
03274 nucleus=305;
03275 break;
03276 case 29:
03277 nucleus=306;
03278 break;
03279 default:
03280 nucleus=304;
03281 break;
03282 }
03283 break;
03284 case 14:
03285 switch (a) {
03286 case 27:
03287 nucleus=307;
03288 break;
03289 case 28:
03290 nucleus=308;
03291 break;
03292 case 29:
03293 nucleus=309;
03294 break;
03295 case 30:
03296 nucleus=310;
03297 break;
03298 default:
03299 nucleus=308;
03300 break;
03301 }
03302 break;
03303 case 15:
03304 switch (a) {
03305 case 30:
03306 nucleus=311;
03307 break;
03308 case 31:
03309 nucleus=312;
03310 break;
03311 case 32:
03312 nucleus=313;
03313 break;
03314 case 33:
03315 nucleus=314;
03316 break;
03317 default:
03318 nucleus=312;
03319 break;
03320 }
03321 break;
03322 case 16:
03323 switch (a) {
03324 case 31:
03325 nucleus=315;
03326 break;
03327 case 32:
03328 nucleus=316;
03329 break;
03330 case 33:
03331 nucleus=317;
03332 break;
03333 case 34:
03334 nucleus=318;
03335 break;
03336 case 35:
03337 nucleus=319;
03338 break;
03339 case 36:
03340 nucleus=320;
03341 break;
03342 default:
03343 nucleus=316;
03344 break;
03345 }
03346 break;
03347 case 22:
03348 switch (a) {
03349 case 45:
03350 nucleus=347;
03351 break;
03352 case 46:
03353 nucleus=348;
03354 break;
03355 case 47:
03356 nucleus=349;
03357 break;
03358 case 48:
03359 nucleus=350;
03360 break;
03361 case 49:
03362 nucleus=351;
03363 break;
03364 case 50:
03365 nucleus=352;
03366 break;
03367 default:
03368 nucleus=350;
03369 break;
03370 }
03371 break;
03372 case 23:
03373 switch (a) {
03374 case 49:
03375 nucleus=353;
03376 break;
03377 case 50:
03378 nucleus=354;
03379 break;
03380 case 51:
03381 nucleus=355;
03382 break;
03383 case 52:
03384 nucleus=356;
03385 break;
03386 case 53:
03387 nucleus=357;
03388 break;
03389 default:
03390 nucleus=355;
03391 break;
03392 }
03393 break;
03394 case 24:
03395 switch (a) {
03396 case 49:
03397 nucleus=358;
03398 break;
03399 case 50:
03400 nucleus=359;
03401 break;
03402 case 51:
03403 nucleus=360;
03404 break;
03405 case 52:
03406 nucleus=361;
03407 break;
03408 case 53:
03409 nucleus=362;
03410 break;
03411 case 54:
03412 nucleus=363;
03413 break;
03414 default:
03415 nucleus=361;
03416 break;
03417 }
03418 break;
03419 case 25:
03420 switch (a) {
03421 case 53:
03422 nucleus=364;
03423 break;
03424 case 54:
03425 nucleus=365;
03426 break;
03427 case 55:
03428 nucleus=366;
03429 break;
03430 case 56:
03431 nucleus=367;
03432 break;
03433 case 57:
03434 nucleus=368;
03435 break;
03436 default:
03437 nucleus=366;
03438 break;
03439 }
03440 break;
03441 case 26:
03442 switch (a) {
03443 case 53:
03444 nucleus=369;
03445 break;
03446 case 54:
03447 nucleus=370;
03448 break;
03449 case 55:
03450 nucleus=371;
03451 break;
03452 case 56:
03453 nucleus=372;
03454 break;
03455 case 57:
03456 nucleus=373;
03457 break;
03458 case 58:
03459 nucleus=374;
03460 break;
03461 default:
03462 nucleus=372;
03463 break;
03464 }
03465 break;
03466 case 28:
03467 switch (a) {
03468 case 57:
03469 nucleus=382;
03470 break;
03471 case 58:
03472 nucleus=383;
03473 break;
03474 case 59:
03475 nucleus=384;
03476 break;
03477 case 60:
03478 nucleus=385;
03479 break;
03480 case 61:
03481 nucleus=386;
03482 break;
03483 case 62:
03484 nucleus=387;
03485 break;
03486 case 63:
03487 nucleus=388;
03488 break;
03489 case 64:
03490 nucleus=389;
03491 break;
03492 default:
03493 nucleus=383;
03494 break;
03495 }
03496 break;
03497 case 29:
03498 switch (a) {
03499 case 62:
03500 nucleus=390;
03501 break;
03502 case 63:
03503 nucleus=391;
03504 break;
03505 case 64:
03506 nucleus=392;
03507 break;
03508 case 65:
03509 nucleus=393;
03510 break;
03511 case 66:
03512 nucleus=394;
03513 break;
03514 case 67:
03515 nucleus=395;
03516 break;
03517 default:
03518 nucleus=392;
03519 break;
03520 }
03521 break;
03522 default:
03523 nucleus=1;
03524 break;
03525 }
03526
03527 return nucleus;
03528 }
03529
03530
03531
03532 Int_t NuReco::GetHadronicFinalState(const NtpStRecord& ntp,
03533 const NuEvent& nu) const
03534 {
03536
03537 Int_t hfs=0;
03538 Int_t proc=nu.iresonance;
03539 TClonesArray* pointStdhepArray = NULL;
03540 pointStdhepArray=ntp.stdhep;
03541 TClonesArray& stdhepArray = *pointStdhepArray;
03542 Int_t nStdHep = stdhepArray.GetEntries();
03543
03544 Int_t itr=nu.mc;
03545
03546 if(proc==1002){
03547 for(int i=0;i<nStdHep;i++){
03548
03549 const NtpMCStdHep* ntpStdHep=
03550 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
03551
03552 if(ntpStdHep->mc==itr && ntpStdHep->IstHEP==3 &&
03553 !(abs(ntpStdHep->IdHEP)==15)){
03554 hfs = ntpStdHep->IdHEP;
03555 break;
03556 }
03557 }
03558 hfs = hfs%1000;
03559 }
03560 else {
03561 for(int i=0;i<nStdHep;i++){
03562
03563 const NtpMCStdHep* ntpStdHep=
03564 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
03565 if(ntpStdHep->mc==itr && ntpStdHep->IstHEP==3){
03566 hfs = ntpStdHep->IdHEP;
03567 break;
03568 }
03569 }
03570 hfs = hfs%1000;
03571 }
03572 return hfs;
03573 }
03574
03575
03576
03577 Float_t NuReco::GetDirCosNu(const NuEvent& nu) const
03578 {
03580
03581 Float_t dirCosNu=-999;
03582 if (nu.detector==Detector::kFar) {
03583 const Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.);
03584 const Float_t bl_y = sqrt(1. - bl_z*bl_z);
03585 dirCosNu=nu.trkvtxdcosz*bl_z+nu.trkvtxdcosy*bl_y;
03586 }
03587 else if (nu.detector==Detector::kNear) {
03588 const Float_t nu_cos = -5.799E-2;
03589 const Float_t nu_sin = sqrt(1 -nu_cos*nu_cos);
03590 dirCosNu=nu.trkvtxdcosz*nu_sin + nu.trkvtxdcosy*nu_cos;
03591 }
03592 else cout<<"Ahhh, detector="<<nu.detector<<endl;
03593
03594
03595 return dirCosNu;
03596 }
03597
03598
03599
03600 void NuReco::GetDirCosNu(const NtpSRTrack& trk,NuEvent& nu) const
03601 {
03602
03603
03604 if (nu.detector==Detector::kFar) {
03605 const Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.);
03606 const Float_t bl_y = sqrt(1. - bl_z*bl_z);
03607 nu.dirCosNu=trk.vtx.dcosz*bl_z+trk.vtx.dcosy*bl_y;
03608 }
03609 else if (nu.detector==Detector::kNear) {
03610
03611
03612
03613
03614
03615
03616
03617
03618
03619
03620
03621
03622
03623 const Float_t nu_cos = -5.799E-2;
03624 const Float_t nu_sin = sqrt(1 -nu_cos*nu_cos);
03625 nu.dirCosNu=trk.vtx.dcosz*nu_sin + trk.vtx.dcosy*nu_cos;
03626 }
03627 else cout<<"Ahhh, detector="<<nu.detector<<endl;
03628 }
03629
03630
03631
03632 Float_t NuReco::GetRadius(Float_t x,Float_t y,const NuEvent& nu) const
03633 {
03634 Float_t zerox=0;
03635 Float_t zeroy=0;
03636 if (nu.detector==Detector::kNear) {
03637
03638 if (nu.anaVersion==NuCuts::kCC0093Std) {
03639 zerox=1.4885*(Munits::m);
03640 zeroy=0.1397*(Munits::m);
03641 }
03642 else {
03643 zerox=1.4828*(Munits::m);
03644 zeroy=0.2384*(Munits::m);
03645 }
03646 }
03647
03648 return sqrt(pow((x-zerox),2)+pow((y-zeroy),2));
03649 }
03650
03651
03652
03653 Float_t NuReco::GetSmallestDeepDistToEdge(const NuEvent& nu) const
03654 {
03655
03656
03657
03658 Float_t distUPI=0;
03659 Float_t distVPI=0;
03660 Float_t distUFI=0;
03661 Float_t distVFI=0;
03662 Float_t xedge=0;
03663 Float_t yedge=0;
03664
03665 PlaneOutline po;
03666
03667
03668 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
03669 PlaneCoverage::kNearPartial,
03670 distUPI,xedge,yedge);
03671 Bool_t isInsideUPI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
03672 PlaneCoverage::kNearPartial);
03673
03674
03675 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
03676 PlaneCoverage::kNearFull,
03677 distUFI,xedge,yedge);
03678 Bool_t isInsideUFI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
03679 PlaneCoverage::kNearFull);
03680
03681
03682 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
03683 PlaneCoverage::kNearPartial,
03684 distVPI,xedge,yedge);
03685 Bool_t isInsideVPI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
03686 PlaneCoverage::kNearPartial);
03687
03688
03689 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
03690 PlaneCoverage::kNearFull,
03691 distVFI,xedge,yedge);
03692 Bool_t isInsideVFI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
03693 PlaneCoverage::kNearFull);
03694
03695
03696
03697 Float_t smallestDist=999;
03698 if (isInsideUPI && isInsideUFI && isInsideVPI && isInsideVFI){
03699 if (smallestDist>distUPI) smallestDist=distUPI;
03700 if (smallestDist>distUFI) smallestDist=distUFI;
03701 if (smallestDist>distVPI) smallestDist=distVPI;
03702 if (smallestDist>distVFI) smallestDist=distVFI;
03703 }
03704
03705
03706
03707
03708 if (smallestDist==999) {
03709 smallestDist=0;
03710 if (smallestDist<distUPI && !isInsideUPI) smallestDist=distUPI;
03711 if (smallestDist<distUFI && !isInsideUFI) smallestDist=distUFI;
03712 if (smallestDist<distVPI && !isInsideVPI) smallestDist=distVPI;
03713 if (smallestDist<distVFI && !isInsideVFI) smallestDist=distVFI;
03714
03715
03716 smallestDist*=-1;
03717 }
03718
03719 if (smallestDist==999) cout<<"Ahhhhhhh"<<endl;
03720
03721 MAXMSG("NuReco",Msg::kDebug,500)
03722 <<"smallestDist="<<smallestDist
03723 <<", UPI="<<distUPI<<" (in="<<isInsideUPI<<")"
03724 <<", UFI="<<distUFI<<" (in="<<isInsideUFI<<")"
03725 <<", VPI="<<distVPI<<" (in="<<isInsideVPI<<")"
03726 <<", VFI="<<distVFI<<" (in="<<isInsideVFI<<")"<<endl;
03727
03728 return smallestDist;
03729 }
03730
03731
03732
03733 void NuReco::TestGetSmallestDeepDistToEdge() const
03734 {
03735 NuEvent nu;
03736
03737 TH2F* hYvsXDistToEdge=new TH2F
03738 ("hYvsXDistToEdge","hYvsXDistToEdge",
03739 600,-3,3, 600,-3,3);
03740 hYvsXDistToEdge->SetTitle("SmallestDeepDistToEdge");
03741 hYvsXDistToEdge->GetXaxis()->SetTitle("X (m)");
03742 hYvsXDistToEdge->GetXaxis()->CenterTitle();
03743 hYvsXDistToEdge->GetYaxis()->SetTitle("Y (m)");
03744 hYvsXDistToEdge->GetYaxis()->CenterTitle();
03745
03746 for (Float_t x=-3;x<3;x+=0.01) {
03747 for (Float_t y=-3;y<3;y+=0.01) {
03748
03749 nu.xEvtVtx=x;
03750 nu.yEvtVtx=y;
03751 Float_t dist=this->GetSmallestDeepDistToEdge(nu);
03752
03753 hYvsXDistToEdge->Fill(x,y,dist);
03754 }
03755 }
03756 }
03757
03758
03759
03760 TVector3 NuReco::xyz2uvz(Float_t x,Float_t y,Float_t z,
03761 VldContext vc) const
03762 {
03763
03764 UgliGeomHandle ugh(vc);
03765
03766
03767 TVector3 xyz(x,y,z);
03768 TVector3 uvz=ugh.xyz2uvz(xyz);
03769 return uvz;
03770 }
03771
03772
03773
03774 TVector3 NuReco::uvz2xyz(Float_t u,Float_t v,Float_t z,
03775 VldContext vc) const
03776 {
03777
03778 UgliGeomHandle ugh(vc);
03779
03780
03781 TVector3 uvz(u,v,z);
03782 TVector3 xyz=ugh.uvz2xyz(uvz);
03783 return xyz;
03784 }
03785
03786
