00001
00002
00003
00004
00005
00007
00008 #include <set>
00009 #include <cmath>
00010
00011 #include "TClonesArray.h"
00012 #include "TFile.h"
00013 #include "TH2F.h"
00014 #include "TRandom3.h"
00015
00016 #include "BeamDataUtil/BDSpillAccessor.h"
00017 #include "Conventions/BeamType.h"
00018 #include "DataUtil/EnergyCorrections.h"
00019 #include "MessageService/MsgService.h"
00020 #include "MessageService/MsgFormat.h"
00021 #include "MCNtuple/NtpMCStdHep.h"
00022 #include "MCNtuple/NtpMCTruth.h"
00023 #include "CandNtupleSR/NtpSREvent.h"
00024 #include "CandNtupleSR/NtpSRShower.h"
00025 #include "CandNtupleSR/NtpSRStrip.h"
00026 #include "CandNtupleSR/NtpSRTrack.h"
00027 #include "StandardNtuple/NtpStRecord.h"
00028 #include "TruthHelperNtuple/NtpTHEvent.h"
00029 #include "TruthHelperNtuple/NtpTHShower.h"
00030 #include "TruthHelperNtuple/NtpTHTrack.h"
00031 #include "DataUtil/PlaneOutline.h"
00032 #include "Conventions/ReleaseType.h"
00033 #include "UgliGeometry/UgliGeomHandle.h"
00034 #include "DataUtil/infid.h"
00035 #include "DataUtil/EnergyResolution.h"
00036
00037 #include "NtupleUtils/LISieve.h"
00038 #include "NtupleUtils/NuEvent.h"
00039 #include "NtupleUtils/NuKDTree.h"
00040 #include "NtupleUtils/NuLibrary.h"
00041 #include "NtupleUtils/NuReco.h"
00042
00043 #include "NtupleUtils/NuIntranuke.h"
00044
00045 using std::endl;
00046 using std::cout;
00047 using std::map;
00048 using std::vector;
00049
00050 namespace EnergyCorrections {}
00051 using namespace EnergyCorrections;
00052
00053 CVSID("$Id: NuReco.cxx,v 1.135 2010/02/09 16:23:39 bckhouse Exp $");
00054
00055
00056
00057 NuReco::NuReco()
00058 {
00059 MSG("NuReco",Msg::kDebug)
00060 <<"Running NuReco Constructor..."<<endl;
00061
00062 MSG("NuReco",Msg::kDebug)
00063 <<"Finished NuReco Constructor"<<endl;
00064 }
00065
00066
00067 NuReco::~NuReco()
00068 {
00069 MSG("NuReco",Msg::kDebug)
00070 <<"Running NuReco Destructor..."<<endl;
00071
00072
00073 MSG("NuReco",Msg::kDebug)
00074 <<"Finished NuReco Destructor"<<endl;
00075 }
00076
00077
00078
00079 Double_t NuReco::GetEvtEnergy(NuEvent& nu, bool includekNN) const
00080 {
00081 Msg::LogLevel_t logLevel=Msg::kDebug;
00082
00083 MAXMSG("NuReco",logLevel,200)
00084 <<"GetEvtEnergy: evt="<<nu.evt<<", trks="<<nu.ntrk
00085 <<", shws="<<nu.nshw<<endl;
00086
00087
00088 this->SetBestTrk(nu);
00089
00090 this->SetBestShw(nu);
00091
00092 this->CalcEvtVariables(nu);
00093 this->CalcTrkVariables(nu);
00094 this->CalcTrkReclamation(nu);
00095
00096 this->CalcShwVariables(nu, includekNN);
00097
00098
00099 nu.energyCC=nu.trkEn+nu.shwEnCC;
00100 nu.energyNC=nu.shwEnNC;
00101 nu.energyRM=nu.trkEn;
00102 nu.energy=nu.energyCC;
00103
00104
00105 if(nu.anaVersion == NuCuts::kCC0720Std){
00106 nu.shwEn = nu.shwEnkNN;
00107 nu.energy = nu.trkEn + nu.shwEn;
00108 }
00109
00110
00111 nu.energyNoRw=nu.energy;
00112 this->CalcKinematicVariables(nu);
00113
00114
00115 this->CalcResolution(nu);
00116
00117 MAXMSG("NuReco",Msg::kInfo,10)
00118 <<endl
00119 <<"GetEvtEnergy:"
00120 <<" nshw="<<nu.nshw
00121 <<", ntrk="<<nu.ntrk
00122 <<", primshw="<<nu.primshw
00123 <<", primtrk="<<nu.primtrk
00124 <<endl
00125 <<"shwExists1,2,3="<<nu.shwExists1<<","
00126 <<nu.shwExists2<<","<<nu.shwExists3
00127 <<", trkExists1,2,3="<<nu.trkExists1<<","
00128 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
00129
00130 MAXMSG("NuReco",Msg::kInfo,10)
00131 <<"E="<<nu.energy<<", trkEn="<<nu.trkEn<<", shwEn="<<nu.shwEn
00132 <<endl
00133 <<"dircosnu="<<nu.dirCosNu
00134 <<", y="<<nu.y<<", q2="<<nu.q2
00135 <<", w2="<<nu.w2<<", x="<<nu.x<<endl;
00136
00137
00138 if (nu.ntrk>3 || nu.nshw>7) {
00139 MAXMSG("NuReco",Msg::kWarning,3)
00140 <<endl<<"Lots of trks/shws"
00141 <<", ntrk="<<nu.ntrk<<", nshw="<<nu.nshw
00142 <<", primtrk="<<nu.primtrk<<", primshw="<<nu.primshw
00143 <<", trkEn="<<nu.trkEn<<", shwEn="<<nu.shwEn
00144 <<endl
00145 <<"shwEn1="<<nu.shwEnCor1<<", shw2="<<nu.shwEnCor2
00146 <<", shw3="<<nu.shwEnCor3<<", shw4="<<nu.shwEnCor4
00147 <<", shw5="<<nu.shwEnCor5
00148
00149 <<endl
00150 <<"trkCv1="<<nu.trkEnCorCurv1<<", trkCv2="<<nu.trkEnCorCurv2
00151 <<", trkCv3="<<nu.trkEnCorCurv3<<endl
00152 <<"trkRg1="<<nu.trkEnCorRange1<<", trkRg2="<<nu.trkEnCorRange2
00153 <<", trkRg3="<<nu.trkEnCorRange3<<endl;
00154 if (nu.ntrk>3) {
00155 MAXMSG("NuReco",Msg::kError,300)
00156 <<endl<<"Lots of trks/shws"
00157 <<", ntrk="<<nu.ntrk<<", nshw="<<nu.nshw
00158 <<", primtrk="<<nu.primtrk<<", primshw="<<nu.primshw
00159 <<", trkEn="<<nu.trkEn<<", shwEn="<<nu.shwEn
00160 <<endl
00161 <<"shwEn1="<<nu.shwEnCor1<<", shw2="<<nu.shwEnCor2
00162 <<", shw3="<<nu.shwEnCor3<<", shw4="<<nu.shwEnCor4
00163 <<", shw5="<<nu.shwEnCor5
00164
00165 <<endl
00166 <<"trkCv1="<<nu.trkEnCorCurv1<<", trkCv2="<<nu.trkEnCorCurv2
00167 <<", trkCv3="<<nu.trkEnCorCurv3<<endl
00168 <<"trkRg1="<<nu.trkEnCorRange1<<", trkRg2="<<nu.trkEnCorRange2
00169 <<", trkRg3="<<nu.trkEnCorRange3<<endl;
00170 }
00171 }
00172
00173 return nu.energy;
00174 }
00175
00176
00177
00178 void NuReco::CalcKinematicVariables(NuEvent& nu) const
00179 {
00180 nu.dirCosNu=this->GetDirCosNu(nu);
00181
00183
00184
00185
00186
00187
00188
00189
00190 const Double_t M=(0.93827 + 0.93957)/2.0;
00191
00192
00193 if (nu.trkEn+nu.shwEn) nu.y=nu.shwEn/(nu.trkEn+nu.shwEn);
00194 nu.q2=2*nu.energy*nu.trkEn*(1.0-nu.dirCosNu);
00195 nu.w2=M*M-nu.q2+2*M*nu.shwEn;
00196 if (nu.shwEn) nu.x=nu.q2/(2*M*nu.shwEn);
00197
00198 }
00199
00200
00201
00202 void NuReco::CalcEvtVariables(NuEvent& nu) const
00203 {
00204
00205 nu.rEvtVtx=this->GetRadius(nu.xEvtVtx,nu.yEvtVtx,nu);
00206 nu.rEvtEnd=this->GetRadius(nu.xEvtEnd,nu.yEvtEnd,nu);
00207 nu.evtVtxUVDiffPl=nu.planeEvtBegu-nu.planeEvtBegv;
00208 nu.distToEdgeEvtVtx=this->GetSmallestDeepDistToEdge(nu);
00209 }
00210
00211
00212
00213 void NuReco::CalcTrkVariables(NuEvent& nu) const
00214 {
00215 if (nu.trkExists) {
00216 nu.rTrkVtx=this->GetRadius(nu.xTrkVtx,nu.yTrkVtx,nu);
00217 nu.rTrkEnd=this->GetRadius(nu.xTrkEnd,nu.yTrkEnd,nu);
00218
00219
00220 if (nu.qp) nu.sigqp_qp=nu.sigqp/nu.qp;
00221 if (nu.ndof) {
00222 nu.chi2PerNdof=nu.chi2/nu.ndof;
00223 nu.prob=TMath::Prob(nu.chi2,static_cast<Int_t>(nu.ndof));
00224 }
00225
00226
00227 nu.containmentFlagCC0093Std=this->GetContainmentFlagCC0093Std(nu);
00228 nu.containmentFlagCC0250Std=this->GetContainmentFlagCC0250Std(nu);
00229 nu.containmentFlagPitt=this->GetContainmentFlagPitt(nu);
00230 nu.containmentFlag=this->GetContainmentFlag(nu);
00231
00232
00233 nu.trkEnRange=this->GetTrackEnergyFromRange(nu);
00234 nu.trkEnCurv=this->GetTrackEnergyFromCurv(nu);
00235 nu.usedRange=this->UseRange(nu);
00236 nu.usedCurv=this->UseCurvature(nu);
00237 nu.trkEn=this->GetTrackEnergy(nu);
00238 nu.trkEnNoRw=nu.trkEn;
00239
00240
00241 nu.trkMomentumTransverse = this->GetTrackTransverseMomentum(nu);
00242 }
00243 else {
00244 MAXMSG("NuReco",Msg::kDebug,10)
00245 <<"No track, nu.primtrk="<<nu.primtrk<<endl;
00246 nu.rTrkVtx=-999;
00247 nu.rTrkEnd=-999;
00248
00249 nu.sigqp_qp=-999;
00250 nu.chi2PerNdof=-1;
00251 nu.prob=-1;
00252
00253 nu.containmentFlagCC0093Std=-1;
00254 nu.containmentFlagCC0250Std=-1;
00255 nu.containmentFlagPitt=-1;
00256 nu.containmentFlag=-1;
00257
00258 nu.trkEnRange=0;
00259 nu.trkEnCurv=0;
00260 nu.usedRange=0;
00261 nu.usedCurv=0;
00262 nu.trkEn=0;
00263 nu.trkEnNoRw=nu.trkEn;
00264
00265 nu.trkMomentumTransverse = -1;
00266 }
00267 }
00268
00269
00270
00271 void NuReco::CalcTrkReclamation(NuEvent& nu) const
00272 {
00273 if (nu.trkExists && nu.trkEnCurv==0) {
00274
00275 if (nu.detector==Detector::kFar) {
00276 MAXMSG("NuReco",Msg::kDebug,100)
00277 <<"Looks like track reclamation for the FD..."
00278 <<", usedCurv="<<nu.usedCurv<<", trkEnCurv="<<nu.trkEnCurv
00279 <<", trkEnRange="<<nu.trkEnRange
00280 <<", trkfitpass="<<nu.trkfitpass
00281 <<endl;
00282 }
00283
00284 MAXMSG("NuReco",Msg::kDebug,10)
00285 <<endl<<"Setting curv=range"
00286 <<", usedCurv="<<nu.usedCurv<<", trkEnCurv="<<nu.trkEnCurv
00287 <<", trkEnRange="<<nu.trkEnRange<<", trkfitpass="<<nu.trkfitpass
00288 <<endl;
00289 nu.trkEnCurv=nu.trkEnRange;
00290 nu.trkEn=nu.trkEnRange;
00291 nu.trkEnNoRw=nu.trkEn;
00292 }
00293 else if (!nu.trkExists && nu.trkEnCurv==0) {
00294 MAXMSG("NuReco",Msg::kDebug,10)
00295 <<"nu.trkExists="<<nu.trkExists<<", nu.trkEnCurv="<<nu.trkEnCurv
00296 <<endl;
00297 }
00298 }
00299
00300
00301
00302 void NuReco::CalcResolution(NuEvent& nu) const
00303 {
00304 if(nu.sigqp>=0){
00305 nu.trkResolution=EnergyResolution::MuonResolution(nu.trkEn,
00306 (nu.trkEn*nu.trkEn)*nu.sigqp,
00307 nu.usedRange,
00308 true,
00309
00310 static_cast<Detector::Detector_t>(nu.detector),
00311 nu.releaseType);
00312
00313 }
00314 nu.shwResolution=EnergyResolution::ShowerResolution(nu.shwEn,
00315 0.0,
00316 static_cast<Detector::Detector_t>(nu.detector),
00317 nu.releaseType);
00318
00319 if(nu.trkResolution==-1.) nu.trkResolution=0.0;
00320
00321 nu.resolution=sqrt( nu.trkResolution*nu.trkResolution + nu.shwResolution*nu.shwResolution );
00322
00323
00324 }
00325
00326
00327
00328 void NuReco::CalcShwVariables(NuEvent& nu, bool includekNN) const
00329 {
00330 if (nu.shwExists) {
00331 nu.shwEnCC=this->GetShowerEnergyCC(nu);
00332 nu.shwEnNC=this->GetShowerEnergyNC(nu);
00333 nu.shwEn=nu.shwEnCC;
00334 nu.shwEnNoRw=nu.shwEn;
00335 if(includekNN){
00336 nu.shwEnkNN = GetShowerEnergykNN(nu, &nu.shwEnReskNN);
00337 }
00338 else{
00339
00340 }
00341 }
00342 else {
00343 nu.shwEn=0;
00344 nu.shwEnNoRw=nu.shwEn;
00345 nu.shwEnCC=0;
00346 nu.shwEnNC=0;
00347 nu.shwEnkNN = 0;
00348 nu.shwEnReskNN = 0;
00349
00350 nu.shwEnCor=0;
00351 nu.shwEnNoCor=0;
00352 nu.shwEnLinCCNoCor=0;
00353 nu.shwEnLinCCCor=0;
00354 nu.shwEnWtCCNoCor=0;
00355 nu.shwEnWtCCCor=0;
00356 nu.shwEnLinNCNoCor=0;
00357 nu.shwEnLinNCCor=0;
00358 nu.shwEnWtNCNoCor=0;
00359 nu.shwEnWtNCCor=0;
00360 nu.shwEnMip=0;
00361 }
00362 }
00363
00364
00365
00366 void NuReco::CalcExtraTruthVariables(NuEvent& nu) const
00367 {
00368
00369 const VldContext& vc=this->GetVldContext(nu);
00370 UgliGeomHandle ugh(vc);
00371
00372
00373 const NuLibrary& lib=NuLibrary::Instance();
00374
00375
00376 TVector3 uvz=lib.reco.xyz2uvz(nu.vtxxMC,nu.vtxyMC,nu.vtxzMC,vc);
00377
00378
00379 nu.vtxuMC=uvz.X();
00380 nu.vtxvMC=uvz.Y();
00381
00382
00383 nu.planeTrkVtxMC=ugh.GetPlaneIdFromZ(nu.vtxzMC).GetPlane();
00384 nu.rTrkVtxMC=lib.reco.GetRadius(nu.vtxxMC,nu.vtxyMC,nu);
00385 }
00386
00387
00388
00389 void NuReco::ApplyReweights(NuEvent& nu) const
00390 {
00391 if (nu.simFlag==SimFlag::kData) {
00392 MAXMSG("NuReco",Msg::kInfo,1)
00393 <<"Not doing ApplyReweights for simFlag="<<nu.simFlag<<endl;
00394
00395 nu.trkEnRw=nu.trkEnNoRw;
00396 nu.shwEnRw=nu.shwEnNoRw;
00397 nu.energyRw=nu.trkEnRw+nu.shwEnRw;
00398 return;
00399 }
00400
00401
00402 nu.trkEnRw=nu.trkEnNoRw*nu.trkEnWeight;
00403 nu.shwEnRw=nu.shwEnNoRw*nu.shwEnWeight;
00404 nu.energyRw=nu.trkEnRw+nu.shwEnRw;
00405
00406
00407 if (nu.applyEnergyShifts) {
00408 MAXMSG("NuReco",Msg::kInfo,1)
00409 <<"Applying SKZP energy shifts"<<endl;
00410 nu.trkEn=nu.trkEnRw;
00411 nu.shwEn=nu.shwEnRw;
00412 nu.energy=nu.energyRw;
00413 }
00414 else {
00415 MAXMSG("NuReco",Msg::kInfo,1)
00416 <<"NOT applying SKZP energy shifts"<<endl;
00417 }
00418
00419
00420
00421 if (nu.anaVersion==NuCuts::kJJH1 ||
00422 nu.anaVersion==NuCuts::kJJE1 ||
00423 nu.anaVersion==NuCuts::kFullDST || NuCuts::IsNMBPQ(nu) ) {
00424 MAXMSG("NuReco",Msg::kInfo,1)
00425 <<"Setting the detector weight for NuMuBar selection"<<endl;
00426 nu.detectorWeight=nu.detectorWeightNMB;
00427 }
00428 else {
00429 MAXMSG("NuReco",Msg::kInfo,1)
00430 <<"Setting the detector weight for NuMu selection"<<endl;
00431 nu.detectorWeight=nu.detectorWeightNM;
00432 }
00433
00434
00435
00436 nu.rw=1;
00437
00438
00439 if (nu.applyBeamWeight) {
00440 MAXMSG("NuReco",Msg::kInfo,1)
00441 <<"Applying beam weight"<<endl;
00442 nu.rw*=nu.beamWeight;
00443 MAXMSG("NuReco",Msg::kInfo,1)
00444 <<"Beam weight is " << nu.rw <<endl;
00445 }
00446 else {
00447 MAXMSG("NuReco",Msg::kInfo,1)
00448 <<"NOT applying beam weight"<<endl;
00449 }
00450
00451
00452 if (nu.applyDetectorWeight) {
00453 MAXMSG("NuReco",Msg::kInfo,1)
00454 <<"Applying detector weight"<<endl;
00455 nu.rw*=nu.detectorWeight;
00456 }
00457 else {
00458 MAXMSG("NuReco",Msg::kInfo,1)
00459 <<"NOT applying detector weight"<<endl;
00460 }
00461
00462
00463 if (nu.applyGeneratorWeight) {
00464 MAXMSG("NuReco",Msg::kInfo,1)
00465 <<"Applying generator weight"<<endl;
00466 nu.rw*=nu.generatorWeight;
00467 }
00468 else {
00469 MAXMSG("NuReco",Msg::kInfo,1)
00470 <<"NOT applying generator weight"<<endl;
00471 }
00472
00473
00474 if (nu.apply1SigmaWeight) {
00475 MAXMSG("NuReco",Msg::kInfo,1)
00476 <<"Applying 1-sigma weight"<<endl;
00477 nu.rw*=(nu.fluxErrTotalErrorAfterTune-1.0)+1.0;
00478 }
00479 else {
00480 MAXMSG("NuReco",Msg::kInfo,1)
00481 <<"NOT applying 1-sigma weight"<<endl;
00482 }
00483
00484 nu.rwActual=nu.rw;
00485
00487
00488 nu.fluxErr=nu.fluxErrTotalErrorAfterTune;
00489
00490
00491 this->CalcResolution(nu);
00492
00493 }
00494
00495
00496
00497 Int_t NuReco::FDRCBoundary(const NuEvent& nu) const
00498 {
00500 Int_t litag=0;
00501
00502 if(nu.detector!=Detector::kFar) return litag;
00503
00504
00505
00506
00507
00508 Int_t numshower=nu.nshw;
00509 Float_t allph=nu.rawPhEvt;
00510
00511
00512 Int_t plbeg=nu.planeEvtHdrBeg;
00513 Int_t plend=nu.planeEvtHdrEnd;
00514
00515 if (numshower) {
00516 if (allph>1e6) litag+=10;
00517
00518 if (((plbeg==1 || plbeg==2) && (plend==63 || plend==64)) ||
00519 ((plbeg==65 || plbeg==66) && (plend==127 || plend==128)) ||
00520 ((plbeg==129 || plbeg==130) && (plend==191 || plend==192)) ||
00521 ((plbeg==193 || plbeg==194) && (plend==247 || plend==248))) litag++;
00522 if (((plbeg==250 || plbeg==251) && (plend==312 || plend==313)) ||
00523 ((plbeg==314 || plbeg==315) && (plend==376 || plend==377)) ||
00524 ((plbeg==378 || plbeg==379) && (plend==440 || plend==441)) ||
00525 ((plbeg==442 || plbeg==443) && (plend==484 || plend==485))) litag++;
00526 }
00527 return litag;
00528 }
00529
00530
00531
00532 Double_t NuReco::GetTrackEnergy(const NuEvent& nu) const
00533 {
00534
00535
00536 if (this->UseRange(nu)) {
00537 return this->GetTrackEnergyFromRange(nu);
00538 }
00539
00540 else if (this->UseCurvature(nu)){
00541 return this->GetTrackEnergyFromCurv(nu);
00542 }
00543 else {
00544 cout<<"Ahhhh, can't use range or curvature"<<endl;
00545 return -1;
00546 }
00547 }
00548
00549
00550
00551 Double_t NuReco::GetTrackEnergyFromRange(const NuEvent& nu) const
00552 {
00553
00554 if (true) {
00555 MAXMSG("NuReco",Msg::kInfo,1)
00556 <<"GetTrackEnergyFromRange: (re)doing EnergyCorrection"<<endl;
00557 Double_t trkEn=this->GetTrackEnergyFromRangeCor
00558 (nu.trkMomentumRange,nu);
00559 MAXMSG("NuReco",Msg::kDebug,3000)
00560 <<"Range: trkEn="<<trkEn<<", trkEnCorRange="<<nu.trkEnCorRange
00561 <<", trkMomentumRange="<<nu.trkMomentumRange<<endl;
00562 return trkEn;
00563 }
00564 else {
00565 MAXMSG("NuReco",Msg::kInfo,1)
00566 <<"GetTrackEnergyFromRange: NOT (re)doing EnergyCorrection"<<endl;
00567 return nu.trkEnRange;
00568 }
00569 }
00570
00571
00572
00573 void NuReco::CapTrackMomentum(Double_t& trkMom) const
00574 {
00575
00576 if (trkMom>(1000*Munits::GeV)) trkMom=1000*Munits::GeV;
00577 else if (trkMom<(-1000*Munits::GeV)) trkMom=-1000*Munits::GeV;
00578 }
00579
00580
00581
00582 Double_t NuReco::GetTrackEnergyFromRangeCor(Double_t trkMomentumRange,
00583 const NuEvent& nu) const
00584 {
00585
00586 if (trkMomentumRange<0) return -1;
00587
00588 this->CapTrackMomentum(trkMomentumRange);
00589
00590 Double_t mr=0;
00591 VldContext vc=this->GetVldContext(nu);
00592 mr=EnergyCorrections::FullyCorrectMomentumFromRange
00593 (trkMomentumRange,vc,nu.releaseType);
00594
00595 return sqrt(pow(mr,2)+pow(0.105658357*(Munits::GeV),2));
00596 }
00597
00598
00599
00600 Double_t NuReco::GetTrackEnergyFromCurv(const NuEvent& nu) const
00601 {
00602
00603 if (true) {
00604 MAXMSG("NuReco",Msg::kInfo,1)
00605 <<"GetTrackEnergyFromCurv: (re)doing EnergyCorrection"<<endl;
00606 Double_t trkEn=this->GetTrackEnergyFromCurvCor(nu.qp,nu);
00607 Double_t p=-1;
00608 if (nu.qp) p=1./nu.qp;
00609 MAXMSG("NuReco",Msg::kDebug,3000)
00610 <<"Curv: trkEn="<<trkEn<<", trkEnCorCurv="<<nu.trkEnCorCurv
00611 <<", 1/qp="<<p<<endl;
00612 return trkEn;
00613 }
00614 else {
00615 MAXMSG("NuReco",Msg::kInfo,1)
00616 <<"GetTrackEnergyFromCurv: NOT (re)doing EnergyCorrection"<<endl;
00617 return nu.trkEnCurv;
00618 }
00619 }
00620
00621
00622
00623 Double_t NuReco::GetTrackEnergyFromCurvCor(Double_t qp,
00624 const NuEvent& nu) const
00625 {
00629 if (qp!=0){
00630 Double_t p=1./qp;
00631 this->CapTrackMomentum(p);
00632 VldContext vc=this->GetVldContext(nu);
00633 p=EnergyCorrections::FullyCorrectSignedMomentumFromCurvature
00634 (p,vc,nu.releaseType);
00635 return sqrt(pow(p,2)+pow(0.105658357*(Munits::GeV),2));
00636 }
00637 else return 0.;
00638 }
00639
00640
00641
00642 Double_t NuReco::GetShowerEnergyCC(const NuEvent& nu) const
00643 {
00644
00645
00646
00647
00648
00649
00650 if (true) {
00651
00652 double shwEn = nu.shwEnLinCCNoCor;
00653 if (shwEn == -1) {
00654 MAXMSG("NuReco",Msg::kInfo,1)
00655 <<"GetShowerEnergy: This appears to be an old DST -- using shwEnNoCor instead of shwEnLinCCNoCor"<<endl;
00656 shwEn = nu.shwEnNoCor;
00657 }
00658 MAXMSG("NuReco",Msg::kInfo,1)
00659 <<"GetShowerEnergy: (re)doing EnergyCorrection"<<endl;
00660 Double_t shwEnCor=this->GetShowerEnergyCor(shwEn, CandShowerHandle::kCC,nu);
00661 MAXMSG("NuReco",Msg::kDebug,300)
00662 <<"Redo energy cor: shwEnNoCor="<<nu.shwEnNoCor
00663 <<", shwEnCor="<<shwEnCor<<endl;
00664 if (shwEnCor == -1) {
00665 MAXMSG("NuReco",Msg::kInfo,1)
00666 <<"GetShowerEnergy: Shower Energy Correction faulty -- using uncorrected"<<endl;
00667 return nu.shwEnNoCor;
00668 }
00669 return shwEnCor;
00670 }
00671 else {
00672 MAXMSG("NuReco",Msg::kInfo,1)
00673 <<"GetShowerEnergy: NOT (re)doing EnergyCorrection"<<endl;
00674 return nu.shwEnLinCCCor;
00675 }
00676 }
00677
00678
00679 Double_t NuReco::GetShowerEnergyNC(const NuEvent& nu) const
00680 {
00681
00682
00683
00684
00685
00686
00687 if (true) {
00688 MAXMSG("NuReco",Msg::kInfo,1)
00689 <<"GetShowerEnergy: (re)doing EnergyCorrection"<<endl;
00690 Double_t shwEnCor=this->GetShowerEnergyCor(nu.shwEnLinNCNoCor, CandShowerHandle::kNC,nu);
00691 MAXMSG("NuReco",Msg::kDebug,300)
00692 <<"Redo energy cor: shwEnNoCor="<<nu.shwEnNoCor
00693 <<", shwEnCor="<<shwEnCor<<endl;
00694 return shwEnCor;
00695 }
00696 else {
00697 MAXMSG("NuReco",Msg::kInfo,1)
00698 <<"GetShowerEnergy: NOT (re)doing EnergyCorrection"<<endl;
00699 return nu.shwEnLinNCCor;
00700 }
00701 }
00702
00703
00704 Double_t NuReco::GetShowerEnergyCor(Double_t shwEn,
00705 const CandShowerHandle::ShowerType_t& st,
00706 const NuEvent& nu) const
00707 {
00708
00709 if (shwEn<0) return -1;
00710
00711
00712 VldContext vc=this->GetVldContext(nu);
00713 return EnergyCorrections::FullyCorrectShowerEnergy(shwEn,st,vc,nu.releaseType);
00714 }
00715
00716
00717 bool NuReco::InitializeShowerEnergykNN(NuKDTree<double, 3>* kdTree,
00718 Detector::Detector_t det,
00719 vector<double>& C,
00720 double& cutoff_lo,
00721 double& cutoff_hi,
00722 int& neighbours) const
00723 {
00724 MSG("NuReco", Msg::kInfo) << "Initializing shower energy kNN..." << endl;
00725
00726
00727 TDirectory* restore = gDirectory;
00728
00729 const TString fname = (det == Detector::kFar) ? "kNNShwEnTrainFar.root"
00730 : "kNNShwEnTrainNear.root";
00731
00732 const NuGeneral& gen = NuLibrary::Instance().general;
00733 const TString releaseDir = gen.GetReleaseDirToUse("NtupleUtils");
00734
00735 const TString fullPath = releaseDir+"/NtupleUtils/data/"+fname;
00736
00737 TFile trainFile(fullPath, "READ");
00738
00739 restore->cd();
00740
00741 if(trainFile.IsZombie()){
00742 MSG("NuReco", Msg::kError) << "No kNN shower training file at "
00743 << fullPath
00744 << " variables will not be filled" << endl;
00745 return false;
00746 }
00747
00748 TTree* tr = (TTree*)trainFile.Get("train");
00749
00750 double trkShwEnNearDW, nplaneShw, sum2Shw, shwEnMC;
00751
00752 const int status = tr->SetBranchAddress("trkShwEnNearDW", &trkShwEnNearDW);
00753 if(status < 0) tr->SetBranchAddress("trkShwEnNear", &trkShwEnNearDW);
00754
00755 tr->SetBranchAddress("nplaneShw", &nplaneShw);
00756 tr->SetBranchAddress("sum2Shw", &sum2Shw);
00757 tr->SetBranchAddress("shwEnMC", &shwEnMC);
00758
00759 TRandom3 r;
00760 const int N = tr->GetEntries();
00761 for(int n = 0; n < N; ++n){
00762 tr->GetEntry(n);
00763 const NuKDPt<3> pt(trkShwEnNearDW,
00764 nplaneShw + r.Gaus(0, 1e-5),
00765 sum2Shw);
00766 kdTree->AddPt(shwEnMC, pt);
00767 }
00768
00769 kdTree->Commit();
00770
00771
00772 Registry* rEnCorr = (Registry*)trainFile.Get("energy_corrections");
00773
00774 C.resize(8);
00775 for(unsigned int n = 0; n < C.size(); ++n)
00776 C[n] = rEnCorr->GetDouble(TString::Format("C%d", n));
00777 cutoff_lo = rEnCorr->GetDouble("cutoff_lo");
00778 cutoff_hi = rEnCorr->GetDouble("cutoff_hi");
00779
00780 const bool success = rEnCorr->Get("num_neighbours", neighbours);
00781 if(!success){
00782 MSG("NuReco", Msg::kWarning)
00783 << "num_neighbours not found. Defaulting to 400" << endl;
00784 neighbours = 400;
00785 }
00786
00787 MSG("NuReco", Msg::kInfo) << "Done initializing shower energy kNN" << endl;
00788
00789 return true;
00790 }
00791
00792
00793 Double_t NuReco::GetShowerEnergykNN(const NuEvent& nu, Float_t* res) const
00794 {
00795
00796
00797 static bool treeOK = true;
00798 if(!treeOK){
00799 if(res) *res = -1;
00800 return -1;
00801 }
00802
00803
00804 static NuKDTree<double, 3> kdTree;
00805
00806 static vector<double> C;
00807 static double cutoff_lo;
00808 static double cutoff_hi;
00809
00810 static int num_neighbours;
00811
00812 if(!kdTree.Committed()){
00813 const Detector::Detector_t det = Detector::Detector_t(nu.detector);
00814 const bool success = InitializeShowerEnergykNN(&kdTree, det,
00815 C, cutoff_lo, cutoff_hi,
00816 num_neighbours);
00817
00818 if(!success){
00819 treeOK = false;
00820 if(res) *res = -1;
00821 return -1;
00822 }
00823 }
00824
00825 double sum2Shw = nu.shwEn;
00826 if(nu.nshw > 1) sum2Shw += nu.shwEnCor2;
00827
00828 const NuKDPt<3> pt(nu.trkShwEnNearDW, double(nu.nplaneShw), sum2Shw);
00829
00830
00831
00832
00833
00834 vector<NuKDPtWithPayload<double, 3> > nearest =
00835 kdTree.FindNearestPts(pt, num_neighbours + 10);
00836
00837
00838 double tot = 0;
00839 double sqrTot = 0;
00840 unsigned int N = num_neighbours;
00841 for(unsigned int n = 0; n < N; ++n){
00842 const double pl = nearest[n].payload;
00843
00844
00845
00846 if(nu.energyMC > 0 && fabs(pl-nu.energyMC) < 1e-5){
00847 ++N;
00848
00849
00850 assert(N <= nearest.size());
00851 continue;
00852 }
00853
00854 tot += pl;
00855 sqrTot += pl*pl;
00856 }
00857
00858 const double sd = sqrt(N*sqrTot-tot*tot)/num_neighbours;
00859 if(res) *res = sd;
00860
00861 double ret = tot/num_neighbours;
00862
00863
00864
00865 double clampedE = ret;
00866 if(clampedE < cutoff_lo) clampedE = cutoff_lo;
00867 if(clampedE > cutoff_hi) clampedE = cutoff_hi;
00868 const double y = TMath::Log(clampedE);
00869 const double y2 = y2*y; const double y3 = y3*y;
00870 const double y4 = y4*y; const double y5 = y5*y;
00871 const double y6 = y6*y; const double y7 = y7*y;
00872 const double corrFactor = C[0]+C[1]*y+C[2]*y2+C[3]*y3+
00873 C[4]*y4+C[5]*y5+C[6]*y6+C[7]*y7;
00874
00875 ret *= corrFactor;
00876
00877 return ret;
00878 }
00879
00880
00881 const NtpSRTrack* NuReco::GetLongestTrackTLikePlanes(const NtpStRecord& ntp,const NtpSREvent& evt) const
00882 {
00883 TClonesArray& trkTca=(*ntp.trk);
00884
00885
00886 MAXMSG("NuReco",Msg::kInfo,1)
00887 <<"Setting best track as longest track (track-like planes)"<<endl;
00888 Int_t trkToUse=0;
00889 Double_t trkLength=0;
00890 if (evt.ntrack>1){
00891 for(Int_t itrk=0;itrk<evt.ntrack;itrk++){
00892 const NtpSRTrack& trk=
00893 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
00894
00895 MAXMSG("NuReco",Msg::kDebug,200)
00896 <<"trk="<<itrk<<"/"<<evt.ntrack
00897 <<", pass="<<trk.fit.pass
00898 <<", ntrklike="<<trk.plane.ntrklike
00899 <<endl;
00900
00901
00902 if (trk.plane.ntrklike>trkLength) {
00903 trkLength=trk.plane.ntrklike;
00904 trkToUse=itrk;
00905 }
00906 }
00907 MAXMSG("NuReco",Msg::kDebug,200)
00908 <<"selected trkToUse="<<trkToUse<<endl;
00909 const NtpSRTrack& trk=
00910 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[trkToUse]]);
00911 return &trk;
00912 }
00913 else if (evt.ntrack==1){
00914 const NtpSRTrack& trk=
00915 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
00916 return &trk;
00917 }
00918 else return 0;
00919 }
00920
00921
00922
00923 Int_t NuReco::GetLongestTrack(const NuEvent& nu) const
00924 {
00925 if (nu.ntrk<=0) return -1;
00926 else if (nu.ntrk==1) return 1;
00927 else if (nu.ntrk==2) {
00928 Int_t trkIndex=1;
00929 Int_t longestTrack=nu.trkLength1;
00930 if (nu.trkLength2>longestTrack) {
00931 trkIndex=2;
00932 longestTrack=nu.trkLength2;
00933 }
00934 return trkIndex;
00935 }
00936 else {
00937 Int_t trkIndex=1;
00938 Int_t longestTrack=nu.trkLength1;
00939 if (nu.trkLength2>longestTrack) {
00940 trkIndex=2;
00941 longestTrack=nu.trkLength2;
00942 }
00943 if (nu.trkLength3>longestTrack) {
00944 trkIndex=3;
00945 longestTrack=nu.trkLength3;
00946 }
00947 return trkIndex;
00948 }
00949 }
00950
00951
00952
00953 const NtpSRTrack* NuReco::GetLongestTrack(const NtpStRecord& ntp,
00954 const NtpSREvent& evt) const
00955 {
00956 TClonesArray& trkTca=(*ntp.trk);
00957
00958
00959 MAXMSG("NuReco",Msg::kInfo,1)
00960 <<"Setting best track as longest track (end-beg+1)"<<endl;
00961 Int_t trkToUse=0;
00962 Double_t trkLength=0;
00963 if (evt.ntrack>1){
00964 for(Int_t itrk=0;itrk<evt.ntrack;itrk++){
00965 const NtpSRTrack& trk=
00966 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
00967
00968 Int_t absLength=abs(trk.plane.end-trk.plane.beg+1);
00969
00970
00971 if (absLength>trkLength) {
00972 trkLength=absLength;
00973 trkToUse=itrk;
00974 }
00975
00976 MAXMSG("NuReco",Msg::kDebug,200)
00977 <<"itrk="<<itrk<<" of "<<evt.ntrack
00978 <<", trkToUse="<<trkToUse
00979 <<", pass="<<trk.fit.pass
00980 <<", ntrklike="<<trk.plane.ntrklike
00981 <<", absLength="<<absLength
00982 <<", longest="<<trkLength
00983 <<endl;
00984 }
00985 MAXMSG("NuReco",Msg::kDebug,200)
00986 <<"Finished loop, selected trkToUse="<<trkToUse<<endl;
00987 const NtpSRTrack& trk=
00988 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[trkToUse]]);
00989 return &trk;
00990 }
00991 else if (evt.ntrack==1){
00992 const NtpSRTrack& trk=
00993 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
00994 return &trk;
00995 }
00996 else return 0;
00997 }
00998
00999
01000
01001 void NuReco::SetBestShw(NuEvent& nu) const
01002 {
01003
01004 if (nu.nshw<=0) {
01005 nu.shwExists=false;
01006 return;
01007 }
01008
01009
01010 Int_t bestShower=this->GetBestShower(nu);
01011
01012
01013 if (bestShower<1) {
01014 MAXMSG("NuReco",Msg::kDebug,3)
01015 <<"No best shower, run="<<nu.run<<", snarl="<<nu.snarl
01016 <<", primshw="<<nu.primshw
01017 <<", nshw="<<nu.nshw
01018 <<", shwExists1,2,3="<<nu.shwExists1<<","
01019 <<nu.shwExists2<<","<<nu.shwExists3<<endl;
01020
01021
01022 nu.shwExists=false;
01023 nu.ndigitShw=-1;
01024 nu.nstripShw=-1;
01025 nu.nplaneShw=-1;
01026 nu.shwEnCor=-1;
01027 nu.shwEnNoCor=-1;
01028 nu.shwEnCC=-1;
01029 nu.shwEnNC=-1;
01030 nu.shwEnLinCCNoCor=-1;
01031 nu.shwEnLinCCCor=-1;
01032 nu.shwEnLinNCNoCor=-1;
01033 nu.shwEnLinNCCor=-1;
01034 nu.shwEnWtCCNoCor=-1;
01035 nu.shwEnWtCCCor=-1;
01036 nu.shwEnWtNCNoCor=-1;
01037 nu.shwEnWtNCCor=-1;
01038
01039 nu.shwEnMip=-1;
01040 nu.planeShwBeg=-1;
01041 nu.planeShwEnd=-1;
01042 nu.planeShwMax=-1;
01043 nu.xShwVtx=-999;
01044 nu.yShwVtx=-999;
01045 nu.zShwVtx=-999;
01046 return;
01047 }
01048 else if (bestShower==1) {
01049 MAXMSG("NuReco",Msg::kDebug,100)
01050 <<"Best shower=first, primshw="<<nu.primshw
01051 <<", nshw="<<nu.nshw
01052 <<", shwExists1,2,3="<<nu.shwExists1<<","
01053 <<nu.shwExists2<<","<<nu.shwExists3<<endl;
01054 nu.shwExists=nu.shwExists1;
01055 nu.ndigitShw=nu.ndigitShw1;
01056 nu.nstripShw=nu.nstripShw1;
01057 nu.nplaneShw=nu.nplaneShw1;
01058 nu.shwEnCor=nu.shwEnCor1;
01059 nu.shwEnMip=nu.shwEnMip1;
01060 nu.shwEnNoCor=nu.shwEnNoCor1;
01061 nu.shwEnLinCCNoCor=nu.shwEnLinCCNoCor1;
01062 nu.shwEnLinCCCor=nu.shwEnLinCCCor1;
01063 nu.shwEnLinNCNoCor=nu.shwEnLinNCNoCor1;
01064 nu.shwEnLinNCCor=nu.shwEnLinNCCor1;
01065 nu.shwEnWtCCNoCor=nu.shwEnWtCCNoCor1;
01066 nu.shwEnWtCCCor=nu.shwEnWtCCCor1;
01067 nu.shwEnWtNCNoCor=nu.shwEnWtNCNoCor1;
01068 nu.shwEnWtNCCor=nu.shwEnWtNCCor1;
01069
01070 nu.planeShwBeg=nu.planeShwBeg1;
01071 nu.planeShwEnd=nu.planeShwEnd1;
01072 nu.planeShwMax=nu.planeShwMax1;
01073 nu.xShwVtx=nu.xShwVtx1;
01074 nu.yShwVtx=nu.yShwVtx1;
01075 nu.zShwVtx=nu.zShwVtx1;
01076 }
01077 else if (bestShower==2) {
01078 MAXMSG("NuReco",Msg::kWarning,3)
01079 <<"Best shower=second, primshw="<<nu.primshw
01080 <<", nshw="<<nu.nshw
01081 <<", shwExists1,2,3="<<nu.shwExists1<<","
01082 <<nu.shwExists2<<","<<nu.shwExists3<<endl;
01083 nu.shwExists=nu.shwExists2;
01084 nu.ndigitShw=nu.ndigitShw2;
01085 nu.nstripShw=nu.nstripShw2;
01086 nu.nplaneShw=nu.nplaneShw2;
01087 nu.shwEnCor=nu.shwEnCor2;
01088 nu.shwEnMip=nu.shwEnMip2;
01089 nu.shwEnNoCor=nu.shwEnNoCor2;
01090 nu.shwEnLinCCNoCor=nu.shwEnLinCCNoCor2;
01091 nu.shwEnLinCCCor=nu.shwEnLinCCCor2;
01092 nu.shwEnLinNCNoCor=nu.shwEnLinNCNoCor2;
01093 nu.shwEnLinNCCor=nu.shwEnLinNCCor2;
01094 nu.shwEnWtCCNoCor=nu.shwEnWtCCNoCor2;
01095 nu.shwEnWtCCCor=nu.shwEnWtCCCor2;
01096 nu.shwEnWtNCNoCor=nu.shwEnWtNCNoCor2;
01097 nu.shwEnWtNCCor=nu.shwEnWtNCCor2;
01098
01099 nu.planeShwBeg=nu.planeShwBeg2;
01100 nu.planeShwEnd=nu.planeShwEnd2;
01101 nu.planeShwMax=nu.planeShwMax2;
01102 nu.xShwVtx=nu.xShwVtx2;
01103 nu.yShwVtx=nu.yShwVtx2;
01104 nu.zShwVtx=nu.zShwVtx2;
01105 }
01106 else if (bestShower==3) {
01107 MAXMSG("NuReco",Msg::kWarning,3)
01108 <<"Best shower=third, primshw="<<nu.primshw
01109 <<", nshw="<<nu.nshw
01110 <<", shwExists1,2,3="<<nu.shwExists1<<","
01111 <<nu.shwExists2<<","<<nu.shwExists3<<endl;
01112 nu.shwExists=nu.shwExists3;
01113 nu.ndigitShw=nu.ndigitShw3;
01114 nu.nstripShw=nu.nstripShw3;
01115 nu.nplaneShw=nu.nplaneShw2;
01116 nu.shwEnCor=nu.shwEnCor3;
01117 nu.shwEnMip=nu.shwEnMip3;
01118 nu.shwEnNoCor=nu.shwEnNoCor3;
01119 nu.shwEnLinCCNoCor=nu.shwEnLinCCNoCor3;
01120 nu.shwEnLinCCCor=nu.shwEnLinCCCor3;
01121 nu.shwEnLinNCNoCor=nu.shwEnLinNCNoCor3;
01122 nu.shwEnLinNCCor=nu.shwEnLinNCCor3;
01123 nu.shwEnWtCCNoCor=nu.shwEnWtCCNoCor3;
01124 nu.shwEnWtCCCor=nu.shwEnWtCCCor3;
01125 nu.shwEnWtNCNoCor=nu.shwEnWtNCNoCor3;
01126 nu.shwEnWtNCCor=nu.shwEnWtNCCor3;
01127
01128 nu.planeShwBeg=nu.planeShwBeg3;
01129 nu.planeShwEnd=nu.planeShwEnd3;
01130 nu.planeShwMax=nu.planeShwMax3;
01131 nu.xShwVtx=nu.xShwVtx3;
01132 nu.yShwVtx=nu.yShwVtx3;
01133 nu.zShwVtx=nu.zShwVtx3;
01134 }
01135 else cout<<"Ahhhhhhhhhhhh"<<endl;
01136 }
01137
01138
01139
01140 void NuReco::SetBestTrk(NuEvent& nu) const
01141 {
01142
01143 if (nu.ntrk<=0) {
01144 nu.trkExists=false;
01145 return;
01146 }
01147
01148
01149 Int_t bestTrack=this->GetBestTrack(nu);
01150
01151
01152 if (bestTrack<1) {
01153 MAXMSG("NuReco",Msg::kInfo,3)
01154 <<"No best trk, primtrk="<<nu.primtrk
01155 <<", ntrk="<<nu.ntrk
01156 <<", trkExists1,2,3="<<nu.trkExists1<<","
01157 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01158
01159
01160 nu.trkExists=false;
01161 nu.trkIndex=-1;
01162 nu.ndigitTrk=-1;
01163 nu.nstripTrk=-1;
01164 nu.trkEnCorRange=-1;
01165 nu.trkEnCorCurv=-1;
01166 nu.trkShwEnNear=-1;
01167 nu.trkShwEnNearDW=-1;
01168 nu.trkMomentumRange=-1;
01169 nu.trkMomentumTransverse=-1;
01170 nu.containedTrk=0;
01171 nu.trkfitpass=-1;
01172 nu.trkvtxdcosz=-999;
01173 nu.trkvtxdcosy=-999;
01174 nu.trknplane=-999;
01175 nu.charge=0;
01176 nu.qp=-999;
01177 nu.qp_rangebiased=-999;
01178 nu.sigqp=-1;
01179 nu.qp_sigqp=-999;
01180 nu.chi2=-1;
01181 nu.ndof=0;
01182 nu.qpFraction=-1;
01183 nu.trkVtxUVDiffPl=-999;
01184 nu.trkLength=-1;
01185 nu.planeTrkNu=-1;
01186 nu.planeTrkNv=-1;
01187 nu.ntrklike=-1;
01188 nu.trkphsigcor=-1;
01189 nu.trkphsigmap=-1;
01190 nu.trkIdMC=0;
01191
01192 nu.jitter=-1;
01193 nu.jPID=-999;
01194 nu.majC=0;
01195
01196
01197
01198 nu.smoothMajC=-999;
01199
01200
01201
01202 nu.xTrkVtx=-999;
01203 nu.yTrkVtx=-999;
01204 nu.zTrkVtx=-999;
01205 nu.uTrkVtx=-999;
01206 nu.vTrkVtx=-999;
01207 nu.planeTrkVtx=-1;
01208 nu.planeTrkBeg=-1;
01209 nu.planeTrkBegu=-1;
01210 nu.planeTrkBegv=-1;
01211 nu.stripTrkBeg = -1;
01212 nu.stripTrkBegu = -1;
01213 nu.stripTrkBegv = -1;
01214 nu.stripTrkEnd = -1;
01215 nu.stripTrkEndu = -1;
01216 nu.stripTrkEndv = -1;
01217 nu.stripTrkBegIsu = false;
01218 nu.stripTrkEndIsu = false;
01219 nu.regionTrkVtx=-1;
01220 nu.edgeRegionTrkVtx=-1;
01221 nu.edgeRegionTrkEnd=-1;
01222 nu.phiTrkVtx=-999;
01223 nu.phiTrkEnd=-999;
01224 nu.parallelStripTrkVtx=-999;
01225 nu.parallelStripTrkEnd=-999;
01226 nu.stripTrkBegPerpFlag = false;
01227 nu.stripTrkEndPerpFlag = false;
01228 nu.stripHoveNumTrkVtx = -999;
01229 nu.stripHoveNumTrkEnd = -999;
01230
01231
01232 nu.xTrkEnd=-999;
01233 nu.yTrkEnd=-999;
01234 nu.zTrkEnd=-999;
01235 nu.uTrkEnd=-999;
01236 nu.vTrkEnd=-999;
01237 nu.planeTrkEnd=-1;
01238 nu.planeTrkEndu=-1;
01239 nu.planeTrkEndv=-1;
01240
01241 nu.drTrkFidall=-999;
01242 nu.dzTrkFidall=-999;
01243 nu.drTrkFidvtx=-999;
01244 nu.dzTrkFidvtx=-999;
01245 nu.drTrkFidend=-999;
01246 nu.dzTrkFidend=-999;
01247 nu.traceTrkFidall = -999;
01248 nu.traceTrkFidvtx = -999;
01249 nu.traceTrkFidend = -999;
01250
01251 nu.cosPrTrkVtx = -999;
01252 return;
01253 }
01254 else if (bestTrack==1) {
01255 MAXMSG("NuReco",Msg::kDebug,3)
01256 <<"Best trk=first, primtrk="<<nu.primtrk
01257 <<", ntrk="<<nu.ntrk
01258 <<", trkExists1,2,3="<<nu.trkExists1<<","
01259 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01260 nu.trkExists=nu.trkExists1;
01261 nu.trkIndex=nu.trkIndex1;
01262 nu.ndigitTrk=nu.ndigitTrk1;
01263 nu.nstripTrk=nu.nstripTrk1;
01264 nu.trkEnCorRange=nu.trkEnCorRange1;
01265 nu.trkEnCorCurv=nu.trkEnCorCurv1;
01266 nu.trkShwEnNear=nu.trkShwEnNear1;
01267 nu.trkShwEnNearDW=nu.trkShwEnNearDW1;
01268 nu.trkMomentumRange=nu.trkMomentumRange1;
01269 nu.trkMomentumTransverse=nu.trkMomentumTransverse1;
01270 nu.containedTrk=nu.containedTrk1;
01271 nu.trkfitpass=nu.trkfitpass1;
01272 nu.trkvtxdcosz=nu.trkvtxdcosz1;
01273 nu.trkvtxdcosy=nu.trkvtxdcosy1;
01274 nu.trknplane=nu.trknplane1;
01275 nu.charge=nu.charge1;
01276 nu.qp=nu.qp1;
01277 nu.qp_rangebiased=nu.qp_rangebiased1;
01278 nu.sigqp=nu.sigqp1;
01279 nu.qp_sigqp=nu.qp_sigqp1;
01280 nu.chi2=nu.chi21;
01281 nu.ndof=nu.ndof1;
01282 nu.qpFraction=nu.qpFraction1;
01283 nu.trkVtxUVDiffPl=nu.trkVtxUVDiffPl1;
01284 nu.trkLength=nu.trkLength1;
01285 nu.planeTrkNu=nu.planeTrkNu1;
01286 nu.planeTrkNv=nu.planeTrkNv1;
01287 nu.ntrklike=nu.ntrklike1;
01288 nu.trkphsigcor=nu.trkphsigcor1;
01289 nu.trkphsigmap=nu.trkphsigmap1;
01290 nu.trkIdMC=nu.trkIdMC1;
01291
01292
01293
01294 nu.xTrkVtx=nu.xTrkVtx1;
01295 nu.yTrkVtx=nu.yTrkVtx1;
01296 nu.zTrkVtx=nu.zTrkVtx1;
01297 nu.uTrkVtx=nu.uTrkVtx1;
01298 nu.vTrkVtx=nu.vTrkVtx1;
01299 nu.planeTrkVtx=nu.planeTrkVtx1;
01300 nu.planeTrkBeg=nu.planeTrkBeg1;
01301 nu.planeTrkBegu=nu.planeTrkBegu1;
01302 nu.planeTrkBegv=nu.planeTrkBegv1;
01303 nu.stripTrkBeg = nu.stripTrkBeg1;
01304 nu.stripTrkBegu = nu.stripTrkBegu1;
01305 nu.stripTrkBegv = nu.stripTrkBegv1;
01306 nu.stripTrkEnd = nu.stripTrkEnd1;
01307 nu.stripTrkEndu = nu.stripTrkEndu1;
01308 nu.stripTrkEndv = nu.stripTrkEndv1;
01309
01310 nu.stripTrkBegIsu = nu.stripTrkBegIsu1;
01311 nu.stripTrkEndIsu = nu.stripTrkEndIsu1;
01312 nu.regionTrkVtx=nu.regionTrkVtx1;
01313 nu.edgeRegionTrkVtx=nu.edgeRegionTrkVtx1;
01314 nu.edgeRegionTrkEnd=nu.edgeRegionTrkEnd1;
01315
01316 nu.phiTrkVtx=nu.phiTrkVtx1;
01317 nu.phiTrkEnd=nu.phiTrkEnd1;;
01318 nu.parallelStripTrkVtx= nu.parallelStripTrkVtx1;
01319 nu.parallelStripTrkEnd= nu.parallelStripTrkEnd1;
01320 nu.stripTrkBegPerpFlag = nu.stripTrkBegPerpFlag1;
01321 nu.stripTrkEndPerpFlag = nu.stripTrkEndPerpFlag1;
01322 nu.stripHoveNumTrkVtx = nu.stripHoveNumTrkVtx1;
01323 nu.stripHoveNumTrkEnd = nu.stripHoveNumTrkEnd1;
01324
01325 nu.xTrkEnd=nu.xTrkEnd1;
01326 nu.yTrkEnd=nu.yTrkEnd1;
01327 nu.zTrkEnd=nu.zTrkEnd1;
01328 nu.uTrkEnd=nu.uTrkEnd1;
01329 nu.vTrkEnd=nu.vTrkEnd1;
01330 nu.planeTrkEnd=nu.planeTrkEnd1;
01331 nu.planeTrkEndu=nu.planeTrkEndu1;
01332 nu.planeTrkEndv=nu.planeTrkEndv1;
01333
01334 nu.drTrkFidall = nu.drTrkFidall1;
01335 nu.dzTrkFidall = nu.dzTrkFidall1;
01336 nu.drTrkFidvtx = nu.drTrkFidvtx1;
01337 nu.dzTrkFidvtx = nu.dzTrkFidvtx1;
01338 nu.drTrkFidend = nu.drTrkFidend1;
01339 nu.dzTrkFidend = nu.dzTrkFidend1;
01340 nu.traceTrkFidall = nu.traceTrkFidall1;
01341 nu.traceTrkFidvtx = nu.traceTrkFidvtx1;
01342 nu.traceTrkFidend = nu.traceTrkFidend1;
01343
01344 nu.cosPrTrkVtx = nu.cosPrTrkVtx1;
01345 }
01346 else if (bestTrack==2) {
01347 MAXMSG("NuReco",Msg::kInfo,3)
01348 <<"Best trk=second, primtrk="<<nu.primtrk
01349 <<", ntrk="<<nu.ntrk
01350 <<", trkExists1,2,3="<<nu.trkExists1<<","
01351 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01352
01353 nu.trkExists=nu.trkExists2;
01354 nu.trkIndex=nu.trkIndex2;
01355 nu.ndigitTrk=nu.ndigitTrk2;
01356 nu.nstripTrk=nu.nstripTrk2;
01357 nu.trkEnCorRange=nu.trkEnCorRange2;
01358 nu.trkEnCorCurv=nu.trkEnCorCurv2;
01359 nu.trkShwEnNear=nu.trkShwEnNear2;
01360 nu.trkShwEnNearDW=nu.trkShwEnNearDW2;
01361 nu.trkMomentumRange=nu.trkMomentumRange2;
01362 nu.trkMomentumTransverse=nu.trkMomentumTransverse2;
01363 nu.containedTrk=nu.containedTrk2;
01364 nu.trkfitpass=nu.trkfitpass2;
01365 nu.trkvtxdcosz=nu.trkvtxdcosz2;
01366 nu.trkvtxdcosy=nu.trkvtxdcosy2;
01367 nu.trknplane=nu.trknplane2;
01368 nu.charge=nu.charge2;
01369 nu.qp=nu.qp2;
01370 nu.qp_rangebiased=nu.qp_rangebiased2;
01371 nu.sigqp=nu.sigqp2;
01372 nu.qp_sigqp=nu.qp_sigqp2;
01373 nu.chi2=nu.chi22;
01374 nu.ndof=nu.ndof2;
01375 nu.qpFraction=nu.qpFraction2;
01376 nu.trkVtxUVDiffPl=nu.trkVtxUVDiffPl2;
01377 nu.trkLength=nu.trkLength2;
01378 nu.planeTrkNu=nu.planeTrkNu2;
01379 nu.planeTrkNv=nu.planeTrkNv2;
01380 nu.ntrklike=nu.ntrklike2;
01381 nu.trkphsigcor=nu.trkphsigcor2;
01382 nu.trkphsigmap=nu.trkphsigmap2;
01383 nu.trkIdMC=nu.trkIdMC2;
01384
01385
01386
01387 nu.xTrkVtx=nu.xTrkVtx2;
01388 nu.yTrkVtx=nu.yTrkVtx2;
01389 nu.zTrkVtx=nu.zTrkVtx2;
01390 nu.uTrkVtx=nu.uTrkVtx2;
01391 nu.vTrkVtx=nu.vTrkVtx2;
01392 nu.planeTrkVtx=nu.planeTrkVtx2;
01393 nu.planeTrkBeg=nu.planeTrkBeg2;
01394 nu.planeTrkBegu=nu.planeTrkBegu2;
01395 nu.planeTrkBegv=nu.planeTrkBegv2;
01396 nu.stripTrkBeg = nu.stripTrkBeg2;
01397 nu.stripTrkBegu = nu.stripTrkBegu2;
01398 nu.stripTrkBegv = nu.stripTrkBegv2;
01399 nu.stripTrkEnd = nu.stripTrkEnd2;
01400 nu.stripTrkEndu = nu.stripTrkEndu2;
01401 nu.stripTrkEndv = nu.stripTrkEndv2;
01402
01403 nu.stripTrkBegIsu = nu.stripTrkBegIsu2;
01404 nu.stripTrkEndIsu = nu.stripTrkEndIsu2;
01405 nu.regionTrkVtx=nu.regionTrkVtx2;
01406 nu.edgeRegionTrkVtx=nu.edgeRegionTrkVtx2;
01407 nu.edgeRegionTrkEnd=nu.edgeRegionTrkEnd2;
01408 nu.phiTrkVtx=nu.phiTrkVtx2;
01409 nu.phiTrkEnd=nu.phiTrkEnd2;;
01410
01411 nu.parallelStripTrkVtx= nu.parallelStripTrkVtx2;
01412 nu.parallelStripTrkEnd= nu.parallelStripTrkEnd2;
01413 nu.stripTrkBegPerpFlag = nu.stripTrkBegPerpFlag2;
01414 nu.stripTrkEndPerpFlag = nu.stripTrkEndPerpFlag2;
01415 nu.stripHoveNumTrkVtx = nu.stripHoveNumTrkVtx2;
01416 nu.stripHoveNumTrkEnd = nu.stripHoveNumTrkEnd2;
01417
01418 nu.xTrkEnd=nu.xTrkEnd2;
01419 nu.yTrkEnd=nu.yTrkEnd2;
01420 nu.zTrkEnd=nu.zTrkEnd2;
01421 nu.uTrkEnd=nu.uTrkEnd2;
01422 nu.vTrkEnd=nu.vTrkEnd2;
01423 nu.planeTrkEnd=nu.planeTrkEnd2;
01424 nu.planeTrkEndu=nu.planeTrkEndu2;
01425 nu.planeTrkEndv=nu.planeTrkEndv2;
01426
01427 nu.drTrkFidall=nu.drTrkFidall2;
01428 nu.dzTrkFidall=nu.dzTrkFidall2;
01429 nu.drTrkFidvtx=nu.drTrkFidvtx2;
01430 nu.dzTrkFidvtx=nu.dzTrkFidvtx2;
01431 nu.drTrkFidend=nu.drTrkFidend2;
01432 nu.dzTrkFidend=nu.dzTrkFidend2;
01433 nu.traceTrkFidall = nu.traceTrkFidall2;
01434 nu.traceTrkFidvtx = nu.traceTrkFidvtx2;
01435 nu.traceTrkFidend = nu.traceTrkFidend2;
01436
01437 nu.cosPrTrkVtx = nu.cosPrTrkVtx2;
01438 }
01439 else if (bestTrack==3) {
01440 MAXMSG("NuReco",Msg::kInfo,3)
01441 <<"Best trk=third, primtrk="<<nu.primtrk
01442 <<", ntrk="<<nu.ntrk
01443 <<", trkExists1,2,3="<<nu.trkExists1<<","
01444 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01445
01446 nu.trkExists=nu.trkExists3;
01447 nu.trkIndex=nu.trkIndex3;
01448 nu.ndigitTrk=nu.ndigitTrk3;
01449 nu.nstripTrk=nu.nstripTrk3;
01450 nu.trkEnCorRange=nu.trkEnCorRange3;
01451 nu.trkEnCorCurv=nu.trkEnCorCurv3;
01452 nu.trkShwEnNear=nu.trkShwEnNear3;
01453 nu.trkShwEnNearDW=nu.trkShwEnNearDW3;
01454 nu.trkMomentumRange=nu.trkMomentumRange3;
01455 nu.trkMomentumTransverse=nu.trkMomentumTransverse3;
01456 nu.containedTrk=nu.containedTrk3;
01457 nu.trkfitpass=nu.trkfitpass3;
01458 nu.trkvtxdcosz=nu.trkvtxdcosz3;
01459 nu.trkvtxdcosy=nu.trkvtxdcosy3;
01460 nu.trknplane=nu.trknplane3;
01461 nu.charge=nu.charge3;
01462 nu.qp=nu.qp3;
01463 nu.qp_rangebiased=nu.qp_rangebiased3;
01464 nu.sigqp=nu.sigqp3;
01465 nu.qp_sigqp=nu.qp_sigqp3;
01466 nu.chi2=nu.chi23;
01467 nu.ndof=nu.ndof3;
01468 nu.qpFraction=nu.qpFraction3;
01469 nu.trkVtxUVDiffPl=nu.trkVtxUVDiffPl3;
01470 nu.trkLength=nu.trkLength3;
01471 nu.planeTrkNu=nu.planeTrkNu3;
01472 nu.planeTrkNv=nu.planeTrkNv3;
01473 nu.ntrklike=nu.ntrklike3;
01474 nu.trkphsigcor=nu.trkphsigcor3;
01475 nu.trkphsigmap=nu.trkphsigmap3;
01476 nu.trkIdMC=nu.trkIdMC3;
01477
01478
01479
01480 nu.xTrkVtx=nu.xTrkVtx3;
01481 nu.yTrkVtx=nu.yTrkVtx3;
01482 nu.zTrkVtx=nu.zTrkVtx3;
01483 nu.uTrkVtx=nu.uTrkVtx3;
01484 nu.vTrkVtx=nu.vTrkVtx3;
01485 nu.planeTrkVtx=nu.planeTrkVtx3;
01486 nu.planeTrkBeg=nu.planeTrkBeg3;
01487 nu.planeTrkBegu=nu.planeTrkBegu3;
01488 nu.planeTrkBegv=nu.planeTrkBegv3;
01489 nu.stripTrkBeg = nu.stripTrkBeg3;
01490 nu.stripTrkBegu = nu.stripTrkBegu3;
01491 nu.stripTrkBegv = nu.stripTrkBegv3;
01492 nu.stripTrkEnd = nu.stripTrkEnd3;
01493 nu.stripTrkEndu = nu.stripTrkEndu3;
01494 nu.stripTrkEndv = nu.stripTrkEndv3;
01495 nu.stripTrkBegIsu = nu.stripTrkBegIsu3;
01496 nu.stripTrkEndIsu = nu.stripTrkEndIsu3;
01497 nu.regionTrkVtx=nu.regionTrkVtx3;
01498 nu.edgeRegionTrkVtx=nu.edgeRegionTrkVtx3;
01499 nu.edgeRegionTrkEnd=nu.edgeRegionTrkEnd3;
01500 nu.phiTrkVtx=nu.phiTrkVtx3;
01501 nu.phiTrkEnd=nu.phiTrkEnd3;
01502
01503 nu.parallelStripTrkVtx= nu.parallelStripTrkVtx3;
01504 nu.parallelStripTrkEnd= nu.parallelStripTrkEnd3;
01505 nu.stripTrkBegPerpFlag = nu.stripTrkBegPerpFlag3;
01506 nu.stripTrkEndPerpFlag = nu.stripTrkEndPerpFlag3;
01507 nu.stripHoveNumTrkVtx = nu.stripHoveNumTrkVtx3;
01508 nu.stripHoveNumTrkEnd = nu.stripHoveNumTrkEnd3;
01509
01510 nu.xTrkEnd=nu.xTrkEnd3;
01511 nu.yTrkEnd=nu.yTrkEnd3;
01512 nu.zTrkEnd=nu.zTrkEnd3;
01513 nu.uTrkEnd=nu.uTrkEnd3;
01514 nu.vTrkEnd=nu.vTrkEnd3;
01515 nu.planeTrkEnd=nu.planeTrkEnd3;
01516 nu.planeTrkEndu=nu.planeTrkEndu3;
01517 nu.planeTrkEndv=nu.planeTrkEndv3;
01518
01519 nu.drTrkFidall=nu.drTrkFidall3;
01520 nu.dzTrkFidall=nu.dzTrkFidall3;
01521 nu.drTrkFidvtx=nu.drTrkFidvtx3;
01522 nu.dzTrkFidvtx=nu.dzTrkFidvtx3;
01523 nu.drTrkFidend=nu.drTrkFidend3;
01524 nu.dzTrkFidend=nu.dzTrkFidend3;
01525 nu.traceTrkFidall = nu.traceTrkFidall3;
01526 nu.traceTrkFidvtx = nu.traceTrkFidvtx3;
01527 nu.traceTrkFidend = nu.traceTrkFidend3;
01528
01529 nu.cosPrTrkVtx = nu.cosPrTrkVtx3;
01530 }
01531 else {
01532
01533
01534 MSG("NuReco", Msg::kWarning) << "Ahhhhhhhhhhhh"
01535 << ": Trying to set a best track higher than stored" << endl;
01536 }
01537
01538
01539
01540
01541
01542
01543
01544 this->SetBestTrkMajorityCurvature(nu);
01545
01546 this->SetBestTrkSAFit(nu);
01547 }
01548
01549
01550
01551 void NuReco::SetBestTrkMajorityCurvature(NuEvent& nu) const
01552 {
01553
01554 const NuLibrary& lib=NuLibrary::Instance();
01555
01556
01557 Int_t bestTrack=lib.reco.GetBestTrack(nu);
01558
01559
01560 if (bestTrack<1) {
01561 MAXMSG("NuReco",Msg::kInfo,3)
01562 <<"SetBestTrkMajorityCurvature: No best trk, primtrk="
01563 <<nu.primtrk
01564 <<", ntrk="<<nu.ntrk
01565 <<", trkExists1,2,3="<<nu.trkExists1<<","
01566 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01567 return;
01568 }
01569 else if (bestTrack==1) {
01570 MAXMSG("NuReco",Msg::kDebug,3)
01571 <<"SetBestTrkMajorityCurvature: Best trk=first, primtrk="
01572 <<nu.primtrk
01573 <<", ntrk="<<nu.ntrk
01574 <<", trkExists1,2,3="<<nu.trkExists1<<","
01575 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01576
01577 nu.jitter=nu.jitter1;
01578 nu.jPID=nu.jPID1;
01579 nu.majC=nu.majC1;
01580
01581
01582
01583 nu.smoothMajC=nu.smoothMajC1;
01584
01585
01586 }
01587 else if (bestTrack==2) {
01588 MAXMSG("NuReco",Msg::kDebug,3)
01589 <<"SetBestTrkMajorityCurvature: Best trk=second, primtrk="
01590 <<nu.primtrk
01591 <<", ntrk="<<nu.ntrk
01592 <<", trkExists1,2,3="<<nu.trkExists1<<","
01593 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01594
01595 nu.jitter=nu.jitter2;
01596 nu.jPID=nu.jPID2;
01597 nu.majC=nu.majC2;
01598
01599
01600
01601 nu.smoothMajC=nu.smoothMajC2;
01602
01603
01604 }
01605 else if (bestTrack==3) {
01606 MAXMSG("NuReco",Msg::kDebug,3)
01607 <<"SetBestTrkMajorityCurvature: Best trk=third, primtrk="
01608 <<nu.primtrk
01609 <<", ntrk="<<nu.ntrk
01610 <<", trkExists1,2,3="<<nu.trkExists1<<","
01611 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01612
01613 nu.jitter=nu.jitter3;
01614 nu.jPID=nu.jPID3;
01615 nu.majC=nu.majC3;
01616
01617
01618
01619 nu.smoothMajC=nu.smoothMajC3;
01620
01621
01622 }
01623 else cout<<"Ahhhhhhhhhhhh"<<endl;
01624 }
01625
01626
01627
01628 void NuReco::SetBestTrkSAFit(NuEvent& nu) const
01629 {
01630
01631 const NuLibrary& lib=NuLibrary::Instance();
01632
01633
01634 Int_t bestTrack=lib.reco.GetBestTrack(nu);
01635
01636
01637 if (bestTrack<1) {
01638 MAXMSG("NuReco",Msg::kInfo,3)
01639 <<"SetBestTrkSAFit: No best trk, primtrk="
01640 <<nu.primtrk
01641 <<", ntrk="<<nu.ntrk
01642 <<", trkExists1,2,3="<<nu.trkExists1<<","
01643 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01644 return;
01645 }
01646 else if (bestTrack==1) {
01647 MAXMSG("NuReco",Msg::kDebug,3)
01648 <<"SetBestTrkSAFit: Best trk=first, primtrk="
01649 <<nu.primtrk
01650 <<", ntrk="<<nu.ntrk
01651 <<", trkExists1,2,3="<<nu.trkExists1<<","
01652 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01653
01654 nu.trkfitpassSA=nu.trkfitpassSA1;
01655 nu.trkvtxdcoszSA=nu.trkvtxdcoszSA1;
01656 nu.chargeSA=nu.chargeSA1;
01657 nu.qpSA=nu.qpSA1;
01658 nu.sigqpSA=nu.sigqpSA1;
01659 nu.chi2SA=nu.chi2SA1;
01660 nu.ndofSA=nu.ndofSA1;
01661 nu.probSA=nu.probSA1;
01662 nu.xTrkVtxSA=nu.xTrkVtxSA1;
01663 nu.yTrkVtxSA=nu.yTrkVtxSA1;
01664 nu.zTrkVtxSA=nu.zTrkVtxSA1;
01665 nu.uTrkVtxSA=nu.uTrkVtxSA1;
01666 nu.vTrkVtxSA=nu.vTrkVtxSA1;
01667 }
01668 else if (bestTrack==2) {
01669 MAXMSG("NuReco",Msg::kDebug,3)
01670 <<"SetBestTrkSAFit: Best trk=second, primtrk="
01671 <<nu.primtrk
01672 <<", ntrk="<<nu.ntrk
01673 <<", trkExists1,2,3="<<nu.trkExists1<<","
01674 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01675
01676 nu.trkfitpassSA=nu.trkfitpassSA2;
01677 nu.trkvtxdcoszSA=nu.trkvtxdcoszSA2;
01678 nu.chargeSA=nu.chargeSA2;
01679 nu.qpSA=nu.qpSA2;
01680 nu.sigqpSA=nu.sigqpSA2;
01681 nu.chi2SA=nu.chi2SA2;
01682 nu.ndofSA=nu.ndofSA2;
01683 nu.probSA=nu.probSA2;
01684 nu.xTrkVtxSA=nu.xTrkVtxSA2;
01685 nu.yTrkVtxSA=nu.yTrkVtxSA2;
01686 nu.zTrkVtxSA=nu.zTrkVtxSA2;
01687 nu.uTrkVtxSA=nu.uTrkVtxSA2;
01688 nu.vTrkVtxSA=nu.vTrkVtxSA2;
01689 }
01690 else if (bestTrack==3) {
01691 MAXMSG("NuReco",Msg::kDebug,3)
01692 <<"SetBestTrkSAFit: Best trk=third, primtrk="
01693 <<nu.primtrk
01694 <<", ntrk="<<nu.ntrk
01695 <<", trkExists1,2,3="<<nu.trkExists1<<","
01696 <<nu.trkExists2<<","<<nu.trkExists3<<endl;
01697
01698 nu.trkfitpassSA=nu.trkfitpassSA3;
01699 nu.trkvtxdcoszSA=nu.trkvtxdcoszSA3;
01700 nu.chargeSA=nu.chargeSA3;
01701 nu.qpSA=nu.qpSA3;
01702 nu.sigqpSA=nu.sigqpSA3;
01703 nu.chi2SA=nu.chi2SA3;
01704 nu.ndofSA=nu.ndofSA3;
01705 nu.probSA=nu.probSA3;
01706 nu.xTrkVtxSA=nu.xTrkVtxSA3;
01707 nu.yTrkVtxSA=nu.yTrkVtxSA3;
01708 nu.zTrkVtxSA=nu.zTrkVtxSA3;
01709 nu.uTrkVtxSA=nu.uTrkVtxSA3;
01710 nu.vTrkVtxSA=nu.vTrkVtxSA3;
01711 }
01712 else cout<<"Ahhhhhhhhhhhh"<<endl;
01713 }
01714
01715
01716
01717 Int_t NuReco::GetTrackCharge(const NtpSRTrack& trk, bool isReverse) const
01718 {
01719 Int_t charge=this->GetTrackCharge(trk.momentum.qp, isReverse);
01720 if (trk.momentum.qp==0) {
01721 if (isReverse) {
01722 MAXMSG("NuReco",Msg::kInfo,3)
01723 <<"Note: assigning charge to be +1 since trk.momentum.qp="
01724 <<trk.momentum.qp<<", trkfitpass="<<trk.fit.pass
01725 <<" and this is RHC" << endl;
01726 }
01727 else {
01728 MAXMSG("NuReco",Msg::kInfo,3)
01729 <<"Note: assigning charge to be -1 since trk.momentum.qp="
01730 <<trk.momentum.qp<<", trkfitpass="<<trk.fit.pass
01731 <<" and this is FHC" << endl;
01732 }
01733 }
01734 return charge;
01735 }
01736
01737
01738
01739 Int_t NuReco::GetTrackCharge(Double_t qp, bool isReverse) const
01740 {
01741 Int_t charge=qp<0 ? -1 : 1;
01742 if (qp==0) {
01743 if (isReverse) charge=+1;
01744 else charge=-1;
01745 }
01746
01747 return charge;
01748 }
01749
01750
01751
01752 Int_t NuReco::GetPrimaryTrack(const NuEvent& nu) const
01753 {
01754 if (nu.trkExists1) return 1;
01755 else return 0;
01756 }
01757
01758
01759
01760 Int_t NuReco::GetPrimaryShowerCCStd(const NuEvent& nu) const
01761 {
01764 if (nu.nshw<=0) {
01765 MAXMSG("NuReco",Msg::kDebug,3)
01766 <<"GetPrimaryShower: No showers so return 0"<<endl;
01767 return 0;
01768 }
01769
01770 if (nu.primshw<0) {
01771
01772 return 0;
01773 }
01774
01775
01776 if (!nu.trkExists) return 1;
01777
01778 MAXMSG("NuReco",Msg::kInfo,1)
01779 <<"GetPrimaryShw: require trk vtx within 0.5 m || >2 GeV"<<endl;
01780
01781
01782
01783
01784 if (fabs(nu.zTrkVtx-nu.zShwVtx1)<(0.5*Munits::m) ||
01785 (fabs(nu.zTrkVtx-nu.zShwVtx1)>=(0.5*Munits::m) &&
01786 nu.shwEnCor1>(2*Munits::GeV))) {
01787 return 1;
01788 }
01789 else {
01790 return 0;
01791 }
01792 }
01793
01794
01795
01796 Int_t NuReco::GetPrimaryShowerStdReco(const NuEvent& nu) const
01797 {
01798 if (nu.nshw<=0) {
01799 MAXMSG("NuReco",Msg::kInfo,3)
01800 <<"GetPrimaryShower: No showers so return 0"<<endl;
01801 return 0;
01802 }
01803
01804 MAXMSG("NuReco",Msg::kInfo,1)
01805 <<"GetPrimaryShw: using evt.primshw"<<endl;
01806 if (nu.primshw<0) return 0;
01807 else {
01808 if (nu.primshw==0 && nu.shwExists1) return 1;
01809 else if (nu.primshw==1 && nu.shwExists2) {
01810 MAXMSG("NuReco",Msg::kWarning,300)
01811 <<"nu.primshw="<<nu.primshw<<endl;
01812 return 2;
01813 }
01814 else if (nu.primshw==2 && nu.shwExists3) {
01815 MAXMSG("NuReco",Msg::kWarning,300)
01816 <<"nu.primshw="<<nu.primshw<<endl;
01817 return 3;
01818 }
01819 else {
01820 cout<<"ahhhhhhhh nu.primshw="<<nu.primshw<<endl;
01821 return 0;
01822 }
01823 }
01824 }
01825
01826
01827
01828 Int_t NuReco::GetBestTrack(const NuEvent& nu) const
01829 {
01830
01831
01834
01835 if (nu.anaVersion==NuCuts::kJJH1 ||
01836 nu.anaVersion==NuCuts::kJJE1 || NuCuts::IsNMB(nu) ) {
01837 MAXMSG("NuReco",Msg::kInfo,1)
01838 <<"Setting best track as longest track"<<endl;
01839 return this->GetLongestTrack(nu);
01840 }
01841 else if (nu.anaVersion==NuCuts::kCC0093Std) {
01842 MAXMSG("NuReco",Msg::kInfo,1)
01843 <<"Setting best track as longest track"<<endl;
01844 return this->GetLongestTrack(nu);
01845 }
01846 else if (nu.anaVersion==NuCuts::kCC0250Std ||
01847 nu.anaVersion==NuCuts::kCC0325Std ||
01848 nu.anaVersion==NuCuts::kCC0720Std ||
01849 nu.anaVersion==NuCuts::kCC0720Test ) {
01850 MAXMSG("NuReco",Msg::kInfo,1)
01851 <<"Setting best track as longest track"<<endl;
01852 return this->GetLongestTrack(nu);
01853 }
01854 else {
01855 MAXMSG("NuReco",Msg::kInfo,1)
01856 <<"Setting best track as primary track (from reconstruction)"
01857 <<endl;
01858 return this->GetPrimaryTrack(nu);
01859 }
01860 }
01861
01862
01863
01864 const NtpSRTrack* NuReco::GetTrackWithIndexX(const NtpStRecord& ntp,
01865 const NtpSREvent& evt,
01866 Int_t index) const
01867 {
01868 TClonesArray& trkTca=(*ntp.trk);
01869
01870 if (evt.ntrack>=index+1 && index>=0){
01871 const NtpSRTrack* trk=
01872 dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[index]]);
01873 return trk;
01874 }
01875 else {
01876 return 0;
01877 }
01878 }
01879
01880
01881
01882 const NtpSRShower* NuReco::GetShowerWithIndexX(const NtpStRecord& ntp,
01883 const NtpSREvent& evt,
01884 Int_t index) const
01885 {
01886 TClonesArray& shwTca=(*ntp.shw);
01887
01888 if (evt.nshower>=index+1 && index>=0){
01889 const NtpSRShower* shw=
01890 dynamic_cast<NtpSRShower*>(shwTca[evt.shw[index]]);
01891 return shw;
01892 }
01893 else {
01894 return 0;
01895 }
01896 }
01897
01898
01899
01900 Int_t NuReco::GetBestShower(const NuEvent& nu) const
01901 {
01904
01905 if (nu.anaVersion==NuCuts::kJJH1 ||
01906 nu.anaVersion==NuCuts::kJJE1 || NuCuts::IsNMBPQ(nu) ) {
01907 MAXMSG("NuReco",Msg::kInfo,1)
01908 <<"Setting best shower as primary shower (from reconstruction)"
01909 <<endl;
01910 return this->GetPrimaryShowerStdReco(nu);
01911 }
01912 else if (nu.anaVersion==NuCuts::kCC0093Std) {
01913 MAXMSG("NuReco",Msg::kInfo,1)
01914 <<"Setting best shower according to CC std"<<endl;
01915 return this->GetPrimaryShowerCCStd(nu);
01916 }
01917 else if (nu.anaVersion==NuCuts::kCC0250Std ||
01918 nu.anaVersion==NuCuts::kCC0325Std ||
01919 nu.anaVersion==NuCuts::kCC0720Std ||
01920 nu.anaVersion==NuCuts::kCC0720Test ||
01921 NuCuts::IsNMBNQ(nu) ) {
01922 MAXMSG("NuReco",Msg::kInfo,1)
01923 <<"Setting best shower according to CC std"<<endl;
01924 return this->GetPrimaryShowerCCStd(nu);
01925 }
01926 else {
01927 MAXMSG("NuReco",Msg::kInfo,1)
01928 <<"Setting best shower as primary shower (from reconstruction)"
01929 <<endl;
01930 return this->GetPrimaryShowerStdReco(nu);
01931 }
01932 }
01933
01934
01935
01936 Bool_t NuReco::UseRange(const NuEvent& nu) const
01937 {
01938 if (nu.containmentFlag==1 ||
01939 nu.containmentFlag==3) return true;
01940 else return false;
01941 }
01942
01943
01944
01945 Bool_t NuReco::UseCurvature(const NuEvent& nu) const
01946 {
01947 if (nu.containmentFlag==2 ||
01948 nu.containmentFlag==4) return true;
01949 else return false;
01950 }
01951
01952
01953
01954 VldContext NuReco::GetVldContext(const NuEvent& nu) const
01955 {
01956
01957 VldTimeStamp time(nu.timeSec,nu.timeNanoSec);
01958
01959
01960 VldContext vc(static_cast<Detector::Detector_t>(nu.detector),
01961 static_cast<SimFlag::SimFlag_t>(nu.simFlag),time);
01962 return vc;
01963 }
01964
01965
01966
01967 void NuReco::GetEvtsPerSlc(const NtpStRecord& ntp,
01968 std::map<Int_t,Int_t>& evtsPerSlc) const
01969 {
01970 TClonesArray& evtTca=(*ntp.evt);
01971 const Int_t numEvts=evtTca.GetEntriesFast();
01972
01973
01974 for (Int_t ntpevt=0;ntpevt<numEvts;++ntpevt) {
01975 const NtpSREvent* pevt=
01976 dynamic_cast<NtpSREvent*>(evtTca[ntpevt]);
01977 const NtpSREvent& evt=(*pevt);
01978
01979
01980 evtsPerSlc[evt.slc]++;
01981 }
01982 }
01983
01984
01985
01986 void NuReco::GetEvtDeltaTs(const NtpStRecord& ntp,
01987 std::map<Int_t,Double_t>& deltaTs) const
01988 {
01989 Msg::LogLevel_t logLevel=Msg::kDebug;
01990
01991 const TClonesArray& evtTca=(*ntp.evt);
01992 const Int_t numEvts=evtTca.GetEntriesFast();
01993
01994 map<Int_t,Double_t> evtMedianTimes;
01995 multiset<Double_t> allMedianTimes;
01996 MsgFormat ffmt("%9.11f");
01997
01998
01999 for (Int_t i=0;i<numEvts;i++) {
02000 const NtpSREvent& evt=
02001 *dynamic_cast<NtpSREvent*>(evtTca[i]);
02002 Double_t medTime=this->GetEvtMedianTime(ntp,evt);
02003 evtMedianTimes[i]=medTime;
02004 allMedianTimes.insert(medTime);
02005 }
02006
02007
02008 typedef map<Int_t,Double_t>::iterator evtTimesIt;
02009 for (evtTimesIt it=evtMedianTimes.begin();
02010 it!=evtMedianTimes.end();++it){
02011
02012
02013 multiset<Double_t>::iterator tIt=allMedianTimes.find(it->second);
02014 if (tIt!=allMedianTimes.end()){
02015
02016 multiset<Double_t>::iterator tBeforeIt=tIt;
02017 --tBeforeIt;
02018 multiset<Double_t>::iterator tAfterIt=tIt;
02019 ++tAfterIt;
02020
02021 Double_t deltaAfter=-1;
02022 if (tAfterIt!=allMedianTimes.end()){
02023 MAXMSG("NuReco",logLevel,500)
02024 <<"t="<<ffmt(*tIt)
02025 <<", tAfterIt="<<ffmt(*tAfterIt)<<endl;
02026 deltaAfter=fabs((*tAfterIt)-(*tIt));
02027 }
02028
02029 Double_t deltaBefore=-1;
02030 if (tBeforeIt!=allMedianTimes.end()){
02031 MAXMSG("NuReco",logLevel,500)
02032 <<"t="<<ffmt(*tIt)
02033 <<", tBeforeIt="<<ffmt(*tBeforeIt)<<endl;
02034 deltaBefore=fabs((*tBeforeIt)-(*tIt));
02035 }
02036
02037 Double_t smallestDelta=deltaAfter;
02038 if (deltaAfter!=-1 && deltaBefore!=-1){
02039 if (deltaBefore<deltaAfter) smallestDelta=deltaBefore;
02040 }
02041 if (deltaBefore==-1){
02042 smallestDelta=deltaAfter;
02043 }
02044 if (deltaAfter==-1){
02045 smallestDelta=deltaBefore;
02046 }
02047
02048 MAXMSG("NuReco",logLevel,500)
02049 <<"smallestDelta="<<ffmt(smallestDelta)
02050 <<", deltaBefore="<<ffmt(deltaBefore)
02051 <<", deltaAfter="<<ffmt(deltaAfter)
02052 <<endl;
02053
02054
02055
02056 if (smallestDelta>=0) deltaTs[it->first]=smallestDelta;
02057 else cout<<"Ahhh bad delta T, dT="<<smallestDelta<<endl;
02058 }
02059 else cout<<"Ahh, not found time"<<endl;
02060 }
02061 }
02062
02063
02064
02065 Double_t NuReco::GetTrkQPFraction(const NtpSRTrack& trk) const
02066 {
02067 Int_t nPositive=0;
02068 Int_t nNegative=0;
02069 Int_t nZero=0;
02070
02071 for (Int_t i=0;i<trk.nstrip;i++){
02072
02073 Double_t stpfitqp=trk.stpfitqp[i];
02074
02075
02076 if (stpfitqp>0) nPositive++;
02077 else if (stpfitqp<0) nNegative++;
02078 else nZero++;
02079
02080 MAXMSG("NuReco",Msg::kDebug,1000)
02081 <<"stpfitqp="<<stpfitqp<<endl;
02082 }
02083
02084 Double_t qpFraction=-1;
02085 if (trk.nstrip>0) qpFraction=1.*nPositive/trk.nstrip;
02086
02087 if (trk.nstrip!=nPositive+nNegative) {
02088 MAXMSG("NuReco",Msg::kDebug,100)
02089 <<"??? stpfitqp: nPos="<<nPositive<<", nNeg="<<nNegative
02090 <<", nZero="<<nZero<<", sumN+P="<<nPositive+nNegative
02091 <<", trkn="<<trk.nstrip
02092 <<", qpFraction="<<qpFraction
02093 <<endl;
02094 }
02095 else {
02096 MAXMSG("NuReco",Msg::kDebug,100)
02097 <<"stpfitqp: nPos="<<nPositive<<", nNeg="<<nNegative
02098 <<", nZero="<<nZero<<", sumN+P="<<nPositive+nNegative
02099 <<", trkn="<<trk.nstrip
02100 <<", qpFraction="<<qpFraction
02101 <<endl;
02102 }
02103
02104 return qpFraction;
02105 }
02106
02107
02108
02109 Bool_t NuReco::IsTrkWithMuonFromNotNuNotPiNotKaonNotCharm
02110 (const NtpStRecord& ntp,const NtpSRTrack& trk) const
02111 {
02113 return false;
02114
02115 TClonesArray& thtrkTca=(*ntp.thtrk);
02116 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02117 if (numthtrks<=0) {
02118 MAXMSG("NuReco",Msg::kDebug,1)
02119 <<"No THTracks, so can't do IsMuonFromNotNuNotPiNotCharm..."
02120 <<endl;
02121 return false;
02122 }
02123 MAXMSG("NuReco",Msg::kDebug,1)
02124 <<"Found THTrack, IsMuonFromNotNuNotPiNotCharm..."<<endl;
02125 const NtpTHTrack& thtrk=
02126 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02127
02128
02129
02130 TClonesArray& mcTca=(*ntp.mc);
02131 const NtpMCTruth& mc=
02132 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
02133
02134 TClonesArray& hepTca=(*ntp.stdhep);
02135
02136
02137 Int_t muonEvent=-1;
02138 Bool_t isFromOther=false;
02139
02140 Int_t parent0=-999;
02141 Int_t parent1=-999;
02142
02143 Int_t parent0b=-999;
02144 Int_t parent1b=-999;
02145
02146 Int_t parent0c=-999;
02147 Int_t parent1c=-999;
02148
02149
02150 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02151 const NtpMCStdHep& stdhep=
02152 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02153
02154 if (abs(stdhep.IdHEP)==13){
02155 muonEvent=stdhep.mc;
02156 MAXMSG("NuReco",Msg::kDebug,3000)
02157 <<"Found muon, mc="<<muonEvent
02158 <<", ihep="<<ihep
02159 <<", id="<<stdhep.IdHEP
02160 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02161 <<", parent=["<<stdhep.parent[0]
02162 <<","<<stdhep.parent[1]<<"]"
02163 <<", E="<<stdhep.p4[3]
02164 <<endl;
02165
02166 if (stdhep.parent[0]-stdhep.parent[1]!=0) {
02167 MSG("NuReco",Msg::kWarning)
02168 <<"Ahhhhhhhh, multiple parents of a muon"
02169 <<", parent=["<<parent0<<","<<parent1<<"]"<<endl;
02170 this->PrintStdhep(ntp,mc);
02171 }
02172
02173 if (parent0==-999) {
02174 parent0=stdhep.parent[0];
02175 parent1=stdhep.parent[1];
02176 }
02177 else if (parent0!=-999 && parent0b==-999) {
02178 MAXMSG("NuReco",Msg::kDebug,3000)
02179 <<"Found second muon, 1st parent=["
02180 <<parent0<<","<<parent1<<"]"
02181 <<", 2nd parent=["<<stdhep.parent[0]
02182 <<","<<stdhep.parent[1]<<"]"
02183 <<endl;
02184 parent0b=stdhep.parent[0];
02185 parent1b=stdhep.parent[1];
02186 }
02187 else if (parent0b!=-999 && parent0c==-999) {
02188 MAXMSG("NuReco",Msg::kDebug,3000)
02189 <<endl
02190 <<"Found third muon, parentb=["
02191 <<parent0b<<","<<parent1b<<"]"
02192 <<", new parent=["<<stdhep.parent[0]
02193 <<","<<stdhep.parent[1]<<"]"
02194 <<endl;
02195 parent0c=stdhep.parent[0];
02196 parent1c=stdhep.parent[1];
02197
02198 }
02199 else if (parent0c!=-999) {
02200 MAXMSG("NuReco",Msg::kDebug,3000)
02201 <<endl
02202 <<"Ahhh, already found >=3 muons, parentb=["
02203 <<parent0b<<","<<parent1b<<"]"
02204 <<", new parent=["<<stdhep.parent[0]
02205 <<","<<stdhep.parent[1]<<"]"
02206 <<endl;
02207
02208 }
02209 }
02210 }
02211
02212
02213
02214
02215 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02216 const NtpMCStdHep& stdhep=
02217 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02218
02219 if (parent0 ==(Int_t)stdhep.index || parent1 ==(Int_t)stdhep.index ||
02220 parent0b==(Int_t)stdhep.index || parent1b==(Int_t)stdhep.index) {
02221
02222 if (abs(stdhep.IdHEP)!=14) {
02223 MAXMSG("NuReco",Msg::kDebug,3000)
02224 <<endl<<endl<<endl
02225 <<"**** Found non-nu muon parent..."
02226 <<", i="<<stdhep.index
02227 <<", id="<<stdhep.IdHEP
02228 <<", parent0="<<parent0
02229 <<", parent1="<<parent1
02230 <<", parent0b="<<parent0b
02231 <<", parent1b="<<parent1b
02232 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02233 <<endl<<endl<<endl;
02234 if (!((abs(stdhep.IdHEP)>200 && abs(stdhep.IdHEP)<400) ||
02235 (abs(stdhep.IdHEP)==411 || abs(stdhep.IdHEP)==421 ||
02236 abs(stdhep.IdHEP)==431 || abs(stdhep.IdHEP)==4122))){
02237 isFromOther=true;
02238 MAXMSG("NuReco",Msg::kInfo,3000)
02239 <<"From a ?? ?? !!!"
02240 <<" m="<<stdhep.mass
02241 <<", E="<<stdhep.p4[3]
02242 <<", id="<<stdhep.IdHEP
02243 <<endl<<endl;
02244
02245 this->PrintStdhep(ntp,mc);
02246 break;
02247 }
02248 }
02249 }
02250 }
02251
02252 if (isFromOther) return true;
02253 else return false;
02254 }
02255
02256
02257
02258 Bool_t NuReco::IsTrkWithMuonFromPionDecay(const NtpStRecord& ntp,
02259 const NtpSRTrack& trk) const
02260 {
02263
02264 TClonesArray& thtrkTca=(*ntp.thtrk);
02265 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02266 if (numthtrks<=0) {
02267 MAXMSG("NuReco",Msg::kDebug,1)
02268 <<"No THTracks, so can't do IsTrkWithMuonFromPionDecay..."<<endl;
02269 return false;
02270 }
02271 MAXMSG("NuReco",Msg::kDebug,1)
02272 <<"Found THTrack, IsTrkWithMuonFromPionDecay..."<<endl;
02273 const NtpTHTrack& thtrk=
02274 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02275
02276
02277
02278 TClonesArray& mcTca=(*ntp.mc);
02279 const NtpMCTruth& mc=
02280 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
02281
02282 TClonesArray& hepTca=(*ntp.stdhep);
02283
02284
02285 Int_t muonEvent=-1;
02286 Bool_t isFromPion=false;
02287
02288 Int_t parent0=-999;
02289 Int_t parent1=-999;
02290
02291 Int_t parent0b=-999;
02292 Int_t parent1b=-999;
02293
02294 Int_t parent0c=-999;
02295 Int_t parent1c=-999;
02296
02297
02298 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02299 const NtpMCStdHep& stdhep=
02300 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02301
02302 if (abs(stdhep.IdHEP)==13){
02303 muonEvent=stdhep.mc;
02304 MAXMSG("NuReco",Msg::kDebug,3000)
02305 <<"Found muon, mc="<<muonEvent
02306 <<", ihep="<<ihep
02307 <<", id="<<stdhep.IdHEP
02308 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02309 <<", parent=["<<stdhep.parent[0]
02310 <<","<<stdhep.parent[1]<<"]"
02311 <<", E="<<stdhep.p4[3]
02312 <<endl;
02313
02314 if (stdhep.parent[0]-stdhep.parent[1]!=0) {
02315 MSG("NuReco",Msg::kWarning)
02316 <<"Ahhhhhhhh, multiple parents of a muon"
02317 <<", parent=["<<parent0<<","<<parent1<<"]"<<endl;
02318 this->PrintStdhep(ntp,mc);
02319 }
02320
02321 if (parent0==-999) {
02322 parent0=stdhep.parent[0];
02323 parent1=stdhep.parent[1];
02324 }
02325 else if (parent0!=-999 && parent0b==-999) {
02326 MAXMSG("NuReco",Msg::kDebug,3000)
02327 <<"Found second muon, 1st parent=["
02328 <<parent0<<","<<parent1<<"]"
02329 <<", 2nd parent=["<<stdhep.parent[0]
02330 <<","<<stdhep.parent[1]<<"]"
02331 <<endl;
02332 parent0b=stdhep.parent[0];
02333 parent1b=stdhep.parent[1];
02334 }
02335 else if (parent0b!=-999 && parent0c==-999) {
02336 MAXMSG("NuReco",Msg::kDebug,3000)
02337 <<endl
02338 <<"Found third muon, parentb=["
02339 <<parent0b<<","<<parent1b<<"]"
02340 <<", new parent=["<<stdhep.parent[0]
02341 <<","<<stdhep.parent[1]<<"]"
02342 <<endl;
02343 parent0c=stdhep.parent[0];
02344 parent1c=stdhep.parent[1];
02345 }
02346 else if (parent0c!=-999) {
02347 MAXMSG("NuReco",Msg::kDebug,3000)
02348 <<endl
02349 <<"Ahhh, already found >=3 muons, parentb=["
02350 <<parent0b<<","<<parent1b<<"]"
02351 <<", new parent=["<<stdhep.parent[0]
02352 <<","<<stdhep.parent[1]<<"]"
02353 <<endl;
02354
02355 }
02356 }
02357 }
02358
02359
02360
02361
02362 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02363 const NtpMCStdHep& stdhep=
02364 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02365
02366 if (parent0 ==(Int_t)stdhep.index || parent1 ==(Int_t)stdhep.index ||
02367 parent0b==(Int_t)stdhep.index || parent1b==(Int_t)stdhep.index) {
02368
02369 if (abs(stdhep.IdHEP)!=14) {
02370 MAXMSG("NuReco",Msg::kDebug,3000)
02371 <<endl
02372 <<"**** Found non-nu muon parent..."
02373 <<", i="<<stdhep.index
02374 <<", id="<<stdhep.IdHEP
02375 <<", parent0="<<parent0
02376 <<", parent1="<<parent1
02377 <<", parent0b="<<parent0b
02378 <<", parent1b="<<parent1b
02379 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02380 <<endl;
02381 if (abs(stdhep.IdHEP)>200 && abs(stdhep.IdHEP)<300){
02382 isFromPion=true;
02383 MAXMSG("NuReco",Msg::kDebug,3000)
02384 <<"From a pion!!!"
02385 <<" m="<<stdhep.mass
02386 <<", E="<<stdhep.p4[3]<<endl<<endl;
02387 break;
02388 }
02389 }
02390 }
02391
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414
02415
02416
02417
02418
02419 }
02420
02421 if (isFromPion) return true;
02422 else return false;
02423 }
02424
02425
02426
02427 Bool_t NuReco::IsTrkWithMuonFromKaonDecay(const NtpStRecord& ntp,
02428 const NtpSRTrack& trk) const
02429 {
02432
02433 TClonesArray& thtrkTca=(*ntp.thtrk);
02434 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02435 if (numthtrks<=0) {
02436 MAXMSG("NuReco",Msg::kDebug,1)
02437 <<"No THTracks, so can't do IsTrkWithMuonFromKaonDecay..."<<endl;
02438 return false;
02439 }
02440 MAXMSG("NuReco",Msg::kDebug,1)
02441 <<"Found THTrack, IsTrkWithMuonFromKaonDecay..."<<endl;
02442 const NtpTHTrack& thtrk=
02443 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02444
02445
02446
02447 TClonesArray& mcTca=(*ntp.mc);
02448 const NtpMCTruth& mc=
02449 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
02450
02451 TClonesArray& hepTca=(*ntp.stdhep);
02452
02453
02454 Int_t muonEvent=-1;
02455 Bool_t isFromKaon=false;
02456
02457 Int_t parent0=-999;
02458 Int_t parent1=-999;
02459
02460 Int_t parent0b=-999;
02461 Int_t parent1b=-999;
02462
02463 Int_t parent0c=-999;
02464 Int_t parent1c=-999;
02465
02466
02467 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02468 const NtpMCStdHep& stdhep=
02469 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02470
02471 if (abs(stdhep.IdHEP)==13){
02472 muonEvent=stdhep.mc;
02473 MAXMSG("NuReco",Msg::kDebug,3000)
02474 <<"Found muon, mc="<<muonEvent
02475 <<", ihep="<<ihep
02476 <<", id="<<stdhep.IdHEP
02477 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02478 <<", parent=["<<stdhep.parent[0]
02479 <<","<<stdhep.parent[1]<<"]"
02480 <<", E="<<stdhep.p4[3]
02481 <<endl;
02482
02483 if (stdhep.parent[0]-stdhep.parent[1]!=0) {
02484 MSG("NuReco",Msg::kWarning)
02485 <<"Ahhhhhhhh, multiple parents of a muon"
02486 <<", parent=["<<parent0<<","<<parent1<<"]"<<endl;
02487 this->PrintStdhep(ntp,mc);
02488 }
02489
02490 if (parent0==-999) {
02491 parent0=stdhep.parent[0];
02492 parent1=stdhep.parent[1];
02493 }
02494 else if (parent0!=-999 && parent0b==-999) {
02495 MAXMSG("NuReco",Msg::kDebug,3000)
02496 <<"Found second muon, 1st parent=["
02497 <<parent0<<","<<parent1<<"]"
02498 <<", 2nd parent=["<<stdhep.parent[0]
02499 <<","<<stdhep.parent[1]<<"]"
02500 <<endl;
02501 parent0b=stdhep.parent[0];
02502 parent1b=stdhep.parent[1];
02503 }
02504 else if (parent0b!=-999 && parent0c==-999) {
02505 MAXMSG("NuReco",Msg::kDebug,3000)
02506 <<endl
02507 <<"Found third muon, parentb=["
02508 <<parent0b<<","<<parent1b<<"]"
02509 <<", new parent=["<<stdhep.parent[0]
02510 <<","<<stdhep.parent[1]<<"]"
02511 <<endl;
02512 parent0c=stdhep.parent[0];
02513 parent1c=stdhep.parent[1];
02514 }
02515 else if (parent0c!=-999) {
02516 MAXMSG("NuReco",Msg::kDebug,3000)
02517 <<endl
02518 <<"Ahhh, already found >=3 muons, parentb=["
02519 <<parent0b<<","<<parent1b<<"]"
02520 <<", new parent=["<<stdhep.parent[0]
02521 <<","<<stdhep.parent[1]<<"]"
02522 <<endl;
02523
02524 }
02525 }
02526 }
02527
02528
02529
02530
02531 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02532 const NtpMCStdHep& stdhep=
02533 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02534
02535 if (parent0 ==(Int_t)stdhep.index || parent1 ==(Int_t)stdhep.index ||
02536 parent0b==(Int_t)stdhep.index || parent1b==(Int_t)stdhep.index) {
02537
02538 if (abs(stdhep.IdHEP)!=14) {
02539 MAXMSG("NuReco",Msg::kDebug,3000)
02540 <<endl
02541 <<"**** Found non-nu muon parent..."
02542 <<", i="<<stdhep.index
02543 <<", id="<<stdhep.IdHEP
02544 <<", parent0="<<parent0
02545 <<", parent1="<<parent1
02546 <<", parent0b="<<parent0b
02547 <<", parent1b="<<parent1b
02548 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02549 <<endl;
02550 if (abs(stdhep.IdHEP)>300 && abs(stdhep.IdHEP)<400){
02551 isFromKaon=true;
02552 MAXMSG("NuReco",Msg::kDebug,3000)
02553 <<"From a kaon!!!"
02554 <<" m="<<stdhep.mass
02555 <<", E="<<stdhep.p4[3]<<endl<<endl;
02556 break;
02557 }
02558 }
02559 }
02560 }
02561
02562 if (isFromKaon) return true;
02563 else return false;
02564 }
02565
02566
02567
02568 Bool_t NuReco::IsTrkWithDimuon(const NtpStRecord& ntp,
02569 const NtpSRTrack& trk) const
02570 {
02571 TClonesArray& thtrkTca=(*ntp.thtrk);
02572 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02573 if (numthtrks<=0) {
02574 MAXMSG("NuReco",Msg::kDebug,1)
02575 <<"No THTracks, so can't do IsDimuon..."<<endl;
02576 return false;
02577 }
02578 MAXMSG("NuReco",Msg::kDebug,1)
02579 <<"Found THTrack, IsDimuon..."<<endl;
02580 const NtpTHTrack& thtrk=
02581 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02582
02583
02584
02585 TClonesArray& mcTca=(*ntp.mc);
02586 const NtpMCTruth& mc=
02587 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
02588
02589 TClonesArray& hepTca=(*ntp.stdhep);
02590
02591
02592 Int_t charmEvent=-1;
02593 Bool_t isDimuon=false;
02594
02595 Int_t child0=-999;
02596 Int_t child1=-999;
02597
02598
02599 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02600 const NtpMCStdHep& stdhep=
02601 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02602
02603 if (abs(stdhep.IdHEP)==411 || abs(stdhep.IdHEP)==421 ||
02604 abs(stdhep.IdHEP)==431 || abs(stdhep.IdHEP)==4122){
02605 charmEvent=stdhep.mc;
02606 MAXMSG("NuReco",Msg::kDebug,3000)
02607 <<"Found Charm particle, mc="<<charmEvent
02608 <<", ihep="<<ihep
02609 <<", id="<<stdhep.IdHEP
02610 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02611 <<endl;
02612
02613 child0=stdhep.child[0];
02614 child1=stdhep.child[1];
02615 }
02616 else if ((abs(stdhep.IdHEP)>4000 && abs(stdhep.IdHEP)<5000) ||
02617 (abs(stdhep.IdHEP)>400 && abs(stdhep.IdHEP)<500)) {
02618 MAXMSG("NuReco",Msg::kInfo,3000)
02619 <<endl<<endl
02620 <<"What is this particle??? Charm?, mc="<<charmEvent
02621 <<", ihep="<<ihep
02622 <<", id="<<stdhep.IdHEP
02623 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02624 <<endl<<endl<<endl;
02625 }
02626
02627 if ((Int_t) stdhep.index>=child0 && (Int_t) stdhep.index<=child1){
02628 MAXMSG("NuReco",Msg::kDebug,3000)
02629 <<"Found Charm child..."
02630 <<", i="<<stdhep.index
02631 <<", child0="<<child0
02632 <<", child1="<<child1
02633 <<", id="<<stdhep.IdHEP
02634 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02635 <<endl;
02636 if (abs(stdhep.IdHEP)==13){
02637 isDimuon=true;
02638 MAXMSG("NuReco",Msg::kDebug,3000)
02639 <<"It's a muon!!!"
02640 <<" m="<<stdhep.mass
02641 <<", E="<<stdhep.p4[3]<<endl<<endl;
02642 break;
02643 }
02644 }
02645 }
02646
02647 if (isDimuon) return true;
02648 else return false;
02649 }
02650
02651
02652
02653 Bool_t NuReco::PrintCharm(const NtpStRecord& ntp,
02654 const NtpSREvent& evt) const
02655 {
02656 TClonesArray& thevtTca=(*ntp.thevt);
02657 const Int_t numthevts=thevtTca.GetEntriesFast();
02658 if (numthevts<=0) {
02659 MAXMSG("NuReco",Msg::kDebug,1)
02660 <<"No THEvents, so can't GetTruthInfo..."<<endl;
02661 return false;
02662 }
02663 MAXMSG("NuReco",Msg::kDebug,1)
02664 <<"Found THEvent, GetTruthInfo..."<<endl;
02665 const NtpTHEvent& thevt=
02666 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
02667
02668
02669
02670 TClonesArray& mcTca=(*ntp.mc);
02671 const NtpMCTruth& mc=
02672 *(dynamic_cast<NtpMCTruth*>(mcTca[thevt.neumc]));
02673
02674 TClonesArray& hepTca=(*ntp.stdhep);
02675
02676
02677 Int_t charmEvent=-1;
02678
02679 Int_t child0=-999;
02680 Int_t child1=-999;
02681
02682
02683 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02684 const NtpMCStdHep& stdhep=
02685 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02686
02687 if (abs(stdhep.IdHEP)==411 || abs(stdhep.IdHEP)==421 ||
02688 abs(stdhep.IdHEP)==431 || abs(stdhep.IdHEP)==4122){
02689 charmEvent=stdhep.mc;
02690 MAXMSG("NuReco",Msg::kDebug,3000)
02691 <<"Found Charm particle, mc="<<charmEvent
02692 <<", ihep="<<ihep
02693 <<", id="<<stdhep.IdHEP
02694 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02695 <<endl;
02696
02697 child0=stdhep.child[0];
02698 child1=stdhep.child[1];
02699 }
02700
02701 if (child0==(Int_t)stdhep.index || child1==(Int_t)stdhep.index) {
02702 MAXMSG("NuReco",Msg::kDebug,3000)
02703 <<"Found Charm child..."
02704 <<", i="<<stdhep.index
02705 <<", child0="<<child0
02706 <<", child1="<<child1
02707 <<", id="<<stdhep.IdHEP
02708 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02709 <<endl;
02710 if (abs(stdhep.IdHEP)==13){
02711 MAXMSG("NuReco",Msg::kDebug,3000)
02712 <<"It's a muon!!!"
02713 <<" m="<<stdhep.mass
02714 <<", E="<<stdhep.p4[3]<<endl<<endl;
02715 }
02716 }
02717
02718 }
02719
02720
02721 if (charmEvent!=-1){
02722 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02723 const NtpMCStdHep& stdhep=
02724 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02725
02726 MAXMSG("NuReco",Msg::kInfo,3000)
02727 <<"ihep="<<ihep
02728 <<", mc="<<stdhep.mc
02729 <<", id="<<stdhep.IdHEP
02730 <<", Ist="<<stdhep.IstHEP
02731 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02732 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02733 <<", m="<<stdhep.mass
02734 <<", E="<<stdhep.p4[3]
02735 <<endl;
02736 }
02737 return true;
02738 }
02739 else return false;
02740 }
02741
02742
02743
02744 void NuReco::PrintStdhep(const NtpStRecord& ntp,
02745 const NtpMCTruth& mc) const
02746 {
02747 TClonesArray& hepTca=(*ntp.stdhep);
02748
02749
02750 MAXMSG("NuReco",Msg::kInfo,3000)
02751 <<"PrintStdhep: mc.stdhep=["<<mc.stdhep[0]
02752 <<","<<mc.stdhep[1]<<"]"<<endl;
02753
02754
02755 for (Int_t ihep=mc.stdhep[0];ihep<=mc.stdhep[1];++ihep) {
02756 const NtpMCStdHep& stdhep=
02757 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
02758
02759 MAXMSG("NuReco",Msg::kInfo,3000)
02760 <<"ihep="<<ihep
02761 <<", mc="<<stdhep.mc
02762 <<", id="<<stdhep.IdHEP
02763 <<", Ist="<<stdhep.IstHEP
02764 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
02765 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
02766 <<", m="<<stdhep.mass
02767 <<", E="<<stdhep.p4[3]
02768 <<endl;
02769 }
02770 }
02771
02772
02773
02774 Double_t NuReco::GetShwMedianTime(const NtpStRecord& ntp,
02775 const NtpSRShower& shw) const
02776 {
02777
02778
02779 multiset<Double_t> times;
02780 const TClonesArray& stpTca=(*ntp.stp);
02781
02782
02783 MAXMSG("NuReco",Msg::kDebug,200)
02784 <<"shw.nstrip="<<shw.nstrip<<endl;
02785 for (Int_t i=0;i<shw.nstrip;i++) {
02786 const NtpSRStrip& stp=
02787 *(dynamic_cast<NtpSRStrip*>(stpTca[shw.stp[i]]));
02788
02789 Double_t time=stp.time0;
02790 if (time<0 || time>1) time=stp.time1;
02791
02792 if (time>0 && time<=1) {
02793 times.insert(time);
02794 }
02795
02796
02797 MAXMSG("NuReco",Msg::kDebug,500)
02798 <<" Strip index="<<stp.index<<", stp="<<stp.strip
02799 <<", t="<<time
02800 <<", pl="<<stp.plane
02801 <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
02802 }
02803
02804
02805 multiset<Double_t>::const_iterator it=times.begin();
02806 advance(it,times.size()/2);
02807 Double_t medianTime=*it;
02808 MAXMSG("NuReco",Msg::kDebug,100)
02809 <<"Median time="<<medianTime<<endl;
02810
02811 return medianTime;
02812 }
02813
02814
02815
02816 Int_t NuReco::GetShwMaxPlane(const NtpStRecord& ntp,
02817 const NtpSRShower& shw) const
02818 {
02819 typedef map<UShort_t, Float_t> planePH_t;
02820 planePH_t planePH;
02821
02822 for(int i = 0; i < shw.nstrip; ++i){
02823 const NtpSRStrip* stp = (NtpSRStrip*)ntp.stp->At(shw.stp[i]);
02824
02825 planePH[stp->plane] += stp->ph0.pe+stp->ph1.pe;
02826 }
02827
02828 Float_t maxPH = 0;
02829 Int_t maxPlane = -1;
02830
02831 for(planePH_t::iterator it = planePH.begin(); it != planePH.end(); ++it){
02832 if(it->second > maxPH){
02833 maxPH = it->second;
02834 maxPlane = it->first;
02835 }
02836 }
02837
02838 MAXMSG("NuReco", Msg::kDebug, 100) << "Max plane=" << maxPlane << endl;
02839
02840 return maxPlane;
02841 }
02842
02843
02844
02845 Double_t NuReco::GetEvtMedianTime(const NtpStRecord& ntp,
02846 const NtpSREvent& evt) const
02847 {
02848
02849
02850 multiset<Double_t> times;
02851 const TClonesArray& stpTca=(*ntp.stp);
02852
02853
02854 MAXMSG("NuReco",Msg::kDebug,200)
02855 <<"evt.nstrip="<<evt.nstrip<<endl;
02856 for (Int_t i=0;i<evt.nstrip;i++) {
02857 const NtpSRStrip& stp=
02858 *(dynamic_cast<NtpSRStrip*>(stpTca[evt.stp[i]]));
02859
02860 Double_t time=stp.time0;
02861 if (time<0 || time>1) time=stp.time1;
02862
02863 if (time>0 && time<=1) {
02864 times.insert(time);
02865 }
02866
02867
02868 MAXMSG("NuReco",Msg::kDebug,500)
02869 <<" Strip index="<<stp.index<<", stp="<<stp.strip
02870 <<", t="<<time
02871 <<", pl="<<stp.plane
02872 <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
02873 }
02874
02875
02876 multiset<Double_t>::const_iterator it=times.begin();
02877 advance(it,times.size()/2);
02878 Double_t medianTime=*it;
02879 MAXMSG("NuReco",Msg::kDebug,100)
02880 <<"Median time="<<medianTime<<endl;
02881
02882 return medianTime;
02883 }
02884
02885
02886
02887 void NuReco::GetBestTruthIndex(const NtpStRecord& ntp,
02888 const NtpSREvent& evt,NuEvent& nu) const
02889 {
02890 TClonesArray& thevtTca=(*ntp.thevt);
02891 const Int_t numthevts=thevtTca.GetEntriesFast();
02892 if (numthevts<=0) {
02893 MAXMSG("NuReco",Msg::kDebug,1)
02894 <<"No THEvents, so can't GetBestTruthIndex..."<<endl;
02895 return;
02896 }
02897 MAXMSG("NuReco",Msg::kDebug,1)
02898 <<"Found THEvent, GetBestTruthIndex..."<<endl;
02899
02900 Int_t mcTrk=-1;
02901 Int_t mcShw=-1;
02902 Int_t mcEvt=-1;
02903
02904 const NtpTHEvent& thevt=
02905 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
02906 mcEvt=thevt.neumc;
02907
02908
02909
02910
02911
02912 TClonesArray& thshwTca=(*ntp.thshw);
02913 const Int_t numthshws=thshwTca.GetEntriesFast();
02914 if (numthshws>0) {
02915 Int_t bestShower=this->GetBestShower(nu);
02916 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,
02917 bestShower-1);
02918 if (pshw) {
02919 const NtpSRShower& shw=*pshw;
02920 const NtpTHShower& thshw=
02921 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
02922 mcShw=thshw.neumc;
02923 }
02924 }
02925
02926
02927 TClonesArray& thtrkTca=(*ntp.thtrk);
02928 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02929 if (numthtrks>0) {
02930 Int_t bestTrack=this->GetBestTrack(nu);
02931 const NtpSRTrack* ptrk=this->GetTrackWithIndexX(ntp,evt,
02932 bestTrack-1);
02933 if (ptrk) {
02934 const NtpSRTrack& trk=*ptrk;
02935 const NtpTHTrack& thtrk=
02936 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
02937 mcTrk=thtrk.neumc;
02938 }
02939 }
02940
02941
02942 nu.mcTrk=mcTrk;
02943 nu.mcShw=mcShw;
02944 nu.mcEvt=mcEvt;
02945 this->ChooseTruthIndexToUse(nu);
02946 }
02947
02948
02949
02950 void NuReco::GetPrimaryTruthIndex(const NtpStRecord& ntp,
02951 const NtpSREvent& evt,NuEvent& nu) const
02952 {
02953 TClonesArray& thevtTca=(*ntp.thevt);
02954 const Int_t numthevts=thevtTca.GetEntriesFast();
02955 if (numthevts<=0) {
02956 MAXMSG("NuReco",Msg::kDebug,1)
02957 <<"No THEvents, so can't GetPrimaryTruthIndex..."<<endl;
02958 return;
02959 }
02960 MAXMSG("NuReco",Msg::kDebug,1)
02961 <<"Found THEvent, GetPrimaryTruthIndex..."<<endl;
02962
02963 Int_t mcTrk=-1;
02964 Int_t mcShw=-1;
02965 Int_t mcEvt=-1;
02966
02967 const NtpTHEvent& thevt=
02968 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
02969 mcEvt=thevt.neumc;
02970
02971
02972
02973
02974
02975 TClonesArray& thshwTca=(*ntp.thshw);
02976 const Int_t numthshws=thshwTca.GetEntriesFast();
02977 if (numthshws>0) {
02978 Int_t bestShower=1;
02979 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,
02980 bestShower-1);
02981 if (pshw) {
02982 const NtpSRShower& shw=*pshw;
02983 const NtpTHShower& thshw=
02984 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
02985 mcShw=thshw.neumc;
02986 }
02987 }
02988
02989
02990 TClonesArray& thtrkTca=(*ntp.thtrk);
02991 const Int_t numthtrks=thtrkTca.GetEntriesFast();
02992 if (numthtrks>0) {
02993 Int_t bestTrack=1;
02994 const NtpSRTrack* ptrk=this->GetTrackWithIndexX(ntp,evt,
02995 bestTrack-1);
02996 if (ptrk) {
02997 const NtpSRTrack& trk=*ptrk;
02998 const NtpTHTrack& thtrk=
02999 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
03000 mcTrk=thtrk.neumc;
03001 }
03002 }
03003
03004
03005 nu.mcTrk=mcTrk;
03006 nu.mcShw=mcShw;
03007 nu.mcEvt=mcEvt;
03008 this->ChooseTruthIndexToUse(nu);
03009 }
03010
03011
03012 void NuReco::GetSecondaryTruthIndex(const NtpStRecord& ntp,
03013 const NtpSREvent& evt,
03014 NuEvent& nu) const
03015 {
03016 Int_t mcShw1 = -1;
03017 Int_t mcShw2 = -1;
03018 Int_t mcShw3 = -1;
03019 Int_t mcShw4 = -1;
03020 Int_t mcShw5 = -1;
03021
03022
03023 TClonesArray& thshwTca=(*ntp.thshw);
03024 const Int_t numthshws=thshwTca.GetEntriesFast();
03025 if (numthshws>0) {
03026 if (nu.shwExists1){
03027 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,0);
03028 if (pshw) {
03029 const NtpSRShower& shw=*pshw;
03030 const NtpTHShower& thshw=
03031 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
03032 mcShw1=thshw.neumc;
03033 }
03034 }
03035 if (nu.shwExists2){
03036 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,1);
03037 if (pshw) {
03038 const NtpSRShower& shw=*pshw;
03039 const NtpTHShower& thshw=
03040 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
03041 mcShw2=thshw.neumc;
03042 }
03043 }
03044 if (nu.shwExists3){
03045 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,2);
03046 if (pshw) {
03047 const NtpSRShower& shw=*pshw;
03048 const NtpTHShower& thshw=
03049 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
03050 mcShw3=thshw.neumc;
03051 }
03052 }
03053 if (nu.shwExists4){
03054 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,3);
03055 if (pshw) {
03056 const NtpSRShower& shw=*pshw;
03057 const NtpTHShower& thshw=
03058 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
03059 mcShw4=thshw.neumc;
03060 }
03061 }
03062 if (nu.shwExists5){
03063 const NtpSRShower* pshw=this->GetShowerWithIndexX(ntp,evt,4);
03064 if (pshw) {
03065 const NtpSRShower& shw=*pshw;
03066 const NtpTHShower& thshw=
03067 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
03068 mcShw5=thshw.neumc;
03069 }
03070 }
03071 }
03072
03073 nu.mcShw1 = mcShw1;
03074 nu.mcShw2 = mcShw2;
03075 nu.mcShw3 = mcShw3;
03076 nu.mcShw4 = mcShw4;
03077 nu.mcShw5 = mcShw5;
03078
03079
03080 Int_t mcTrk1 = -1;
03081 Int_t mcTrk2 = -1;
03082 Int_t mcTrk3 = -1;
03083
03084
03085 TClonesArray& thtrkTca=(*ntp.thtrk);
03086 const Int_t numthtrks=thtrkTca.GetEntriesFast();
03087 if (numthtrks>0) {
03088 if (nu.trkExists1){
03089 const NtpSRTrack* ptrk =
03090 this->GetTrackWithIndexX(ntp,evt,0);
03091 if (ptrk){
03092 const NtpSRTrack& trk=*ptrk;
03093 const NtpTHTrack& thtrk=
03094 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
03095 mcTrk1=thtrk.neumc;
03096 }
03097 }
03098 if (nu.trkExists2){
03099 const NtpSRTrack* ptrk =
03100 this->GetTrackWithIndexX(ntp,evt,1);
03101 if (ptrk){
03102 const NtpSRTrack& trk=*ptrk;
03103 const NtpTHTrack& thtrk=
03104 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
03105 mcTrk2=thtrk.neumc;
03106 }
03107 }
03108 if (nu.trkExists3){
03109 const NtpSRTrack* ptrk =
03110 this->GetTrackWithIndexX(ntp,evt,2);
03111 if (ptrk){
03112 const NtpSRTrack& trk=*ptrk;
03113 const NtpTHTrack& thtrk=
03114 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
03115 mcTrk3=thtrk.neumc;
03116 }
03117 }
03118 }
03119
03120 nu.mcTrk1 = mcTrk1;
03121 nu.mcTrk2 = mcTrk2;
03122 nu.mcTrk3 = mcTrk3;
03123 }
03124
03125
03126
03127 void NuReco::ChooseTruthIndexToUse(NuEvent& nu) const
03128 {
03129
03130
03131
03132
03133
03134
03135
03136
03137
03138
03139
03140 if (nu.mcTrk>=0) nu.mc=nu.mcTrk;
03141 else nu.mc=nu.mcShw;
03142 }
03143
03144
03145
03146 void NuReco::GetTruthInfo(const NtpStRecord& ntp,
03147 const NtpSREvent& evt,NuEvent& nu) const
03148 {
03149
03150
03151 this->GetPrimaryTruthIndex(ntp,evt,nu);
03152 this->GetSecondaryTruthIndex(ntp,evt,nu);
03153
03154
03155
03156
03157
03158 if (nu.mc<0) {
03159 MAXMSG("NuReco",Msg::kInfo,1)
03160 <<"Not Getting TruthInfo for simFlag="<<nu.simFlag<<endl;
03161 return;
03162 }
03163 MAXMSG("NuReco",Msg::kInfo,1)
03164 <<"Getting TruthInfo for simFlag="<<nu.simFlag<<endl;
03165
03166
03167
03168 TClonesArray& mcTca=(*ntp.mc);
03169 const NtpMCTruth& mc=
03170 *(dynamic_cast<NtpMCTruth*>(mcTca[nu.mc]));
03171
03172
03173 this->GetTruthInfo(ntp,mc,nu);
03174
03175 if (mc.iaction==1){
03176 MAXMSG("NuReco",Msg::kDebug,300)
03177 <<"y="<<mc.y<<endl
03178 <<"p4neu[3]="<<mc.p4neu[3]<<", p4tgt[3]="<<mc.p4tgt[3]
03179 <<", sum="<<mc.p4neu[3]+mc.p4tgt[3]
03180 <<endl
03181 <<"p4mu1[3]="<<mc.p4mu1[3]<<", mc.p4shw[3]="<<mc.p4shw[3]
03182 <<", sum="<<fabs(mc.p4mu1[3])+mc.p4shw[3]
03183 <<", labY="<<mc.p4shw[3]/(fabs(mc.p4mu1[3])+mc.p4shw[3])
03184 <<endl;
03185 }
03186 else if (mc.iaction==0){
03187 MAXMSG("NuReco",Msg::kDebug,300)
03188 <<"NC: y="<<mc.y<<endl
03189 <<"p4neu[3]="<<mc.p4neu[3]
03190 <<", p4tgt[3]="<<mc.p4tgt[3]
03191 <<endl
03192 <<"p4mu1[3]="<<mc.p4mu1[3]<<", p4shw[3]="<<mc.p4shw[3]
03193 <<endl;
03194 }
03195 }
03196
03197
03198
03199 void NuReco::GetTruthInfo(const NtpStRecord& ntp,
03200 const NtpMCTruth& mc,NuEvent& nu) const
03201 {
03202
03203 Double_t trueEn=mc.y*mc.p4neu[3];
03204 if (mc.iaction==1) trueEn=mc.p4neu[3];
03205
03206 Double_t muEn=(1-mc.y)*mc.p4neu[3];
03207 if (mc.iaction==0) muEn=0;
03208 Double_t hadEn=mc.y*mc.p4neu[3];
03209
03210 nu.energyMC=trueEn;
03211 nu.neuEnMC=mc.p4neu[3];
03212 nu.neuPxMC=mc.p4neu[0];
03213 nu.neuPyMC=mc.p4neu[1];
03214 nu.neuPzMC=mc.p4neu[2];
03215 nu.mu1EnMC=fabs(mc.p4mu1[3]);
03216 nu.mu1PxMC=mc.p4mu1[0];
03217 nu.mu1PyMC=mc.p4mu1[1];
03218 nu.mu1PzMC=mc.p4mu1[2];
03219 nu.tgtEnMC=mc.p4tgt[3];
03220 nu.tgtPxMC=mc.p4tgt[0];
03221 nu.tgtPyMC=mc.p4tgt[1];
03222 nu.tgtPzMC=mc.p4tgt[2];
03223 nu.zMC=static_cast<Int_t>(mc.z);
03224 nu.aMC=static_cast<Int_t>(mc.a);
03225 nu.yMC=mc.y;
03226 nu.y2MC=mc.p4shw[3]/(fabs(mc.p4mu1[3])+mc.p4shw[3]);
03227 nu.xMC=mc.x;
03228 nu.q2MC=mc.q2;
03229 nu.w2MC=mc.w2;
03230 nu.trkEnMC=mc.p4mu1[3];
03231 nu.trkEn2MC=muEn;
03232 nu.shwEnMC=mc.p4shw[3];
03233 nu.shwEn2MC=hadEn;
03234 nu.iaction=mc.iaction;
03235 nu.iresonance=mc.iresonance;
03236 nu.inu=mc.inu;
03237 nu.inunoosc=mc.inunoosc;
03238 nu.itg=mc.itg;
03239
03240
03241 nu.trkStartEnMC = this->GetMuonFirstHitEnergy(ntp, mc);
03242 nu.trkEndEnMC = this->GetMuonLastHitEnergy(ntp, mc);
03243
03244 nu.trkContainmentMC = this->GetMuonContainmentMC(ntp, mc);
03245
03246 nu.vtxxMC=mc.vtxx;
03247 nu.vtxyMC=mc.vtxy;
03248 nu.vtxzMC=mc.vtxz;
03249
03250 nu.Npz=mc.flux.npz;
03251 nu.NdxdzNea=mc.flux.ndxdznear;
03252 nu.NdydzNea=mc.flux.ndydznear;
03253 nu.NenergyN=mc.flux.nenergynear;
03254 nu.NWtNear=mc.flux.nwtnear;
03255 nu.NdxdzFar=mc.flux.ndxdzfar;
03256 nu.NdydzFar=mc.flux.ndydzfar;
03257 nu.NenergyF=mc.flux.nenergyfar;
03258 nu.NWtFar=mc.flux.nwtfar;
03259 nu.Ndecay=mc.flux.ndecay;
03260 nu.Vx=mc.flux.vx;
03261 nu.Vy=mc.flux.vy;
03262 nu.Vz=mc.flux.vz;
03263 nu.pdPx=mc.flux.pdpx;
03264 nu.pdPy=mc.flux.pdpy;
03265 nu.pdPz=mc.flux.pdpz;
03266 nu.ppdxdz=mc.flux.ppdxdz;
03267 nu.ppdydz=mc.flux.ppdydz;
03268 nu.pppz=mc.flux.pppz;
03269 nu.ppenergy=mc.flux.ppenergy;
03270 nu.ppmedium=mc.flux.ppmedium;
03271 nu.ppvx=mc.flux.ppvx;
03272 nu.ppvy=mc.flux.ppvy;
03273 nu.ppvz=mc.flux.ppvz;
03274 nu.ptype=mc.flux.ptype;
03275 nu.Necm=mc.flux.necm;
03276 nu.Nimpwt=mc.flux.nimpwt;
03277 nu.tvx=mc.flux.tvx;
03278 nu.tvy=mc.flux.tvy;
03279 nu.tvz=mc.flux.tvz;
03280 nu.tpx=mc.flux.tpx;
03281 nu.tpy=mc.flux.tpy;
03282 nu.tpz=mc.flux.tpz;
03283 nu.tptype=mc.flux.tptype;
03284 nu.tgen=mc.flux.tgen;
03285
03286
03287
03288
03289
03290 static NuIntranuke inuke;
03291 if(inuke.Reweight(&ntp,&mc))
03292 {
03293 nu.InukeNwts = inuke.GetNwts() ;
03294 nu.InukePiCExchgP = inuke.GetWeight(0) ;
03295 nu.InukePiCExchgN = inuke.GetWeight(1) ;
03296 nu.InukePiEScatP = inuke.GetWeight(2) ;
03297 nu.InukePiEScatN = inuke.GetWeight(3) ;
03298 nu.InukePiInEScatP = inuke.GetWeight(4) ;
03299 nu.InukePiInEScatN = inuke.GetWeight(5) ;
03300 nu.InukePiAbsorbP = inuke.GetWeight(6) ;
03301 nu.InukePiAbsorbN = inuke.GetWeight(7) ;
03302 nu.InukePi2PiP = inuke.GetWeight(8) ;
03303 nu.InukePi2PiN = inuke.GetWeight(9) ;
03304 nu.InukeNknockP = inuke.GetWeight(10) ;
03305 nu.InukeNknockN = inuke.GetWeight(11) ;
03306 nu.InukeNNPiP = inuke.GetWeight(12) ;
03307 nu.InukeNNPiN = inuke.GetWeight(13) ;
03308 nu.InukeFormTP = inuke.GetWeight(14) ;
03309 nu.InukeFormTN = inuke.GetWeight(15) ;
03310 nu.InukePiXsecP = inuke.GetWeight(16) ;
03311 nu.InukePiXsecN = inuke.GetWeight(17) ;
03312 nu.InukeNXsecP = inuke.GetWeight(18) ;
03313 nu.InukeNXsecN = inuke.GetWeight(19) ;
03314 nu.InukeNucrad = inuke.GetNucrad() ;
03315 nu.InukeWrad = inuke.GetWrad();
03316 }
03317
03318
03319
03320
03321 this->CalcExtraTruthVariables(nu);
03322
03323 nu.nucleusMC=this->GetNucleus(ntp,nu);
03324 nu.initialStateMC=this->GetInitialState(ntp,nu);
03325 nu.hadronicFinalStateMC=this->GetHadronicFinalState(ntp,nu);
03326 this->GetPreInukeFinalStateVars(ntp,nu);
03327 }
03328
03329
03330
03331 void NuReco::PrintTrueEnergy(const NtpStRecord& ntp,
03332 const NtpSREvent& evt,Double_t recoEn) const
03333 {
03334 TClonesArray& thevtTca=(*ntp.thevt);
03335 const Int_t numthevts=thevtTca.GetEntriesFast();
03336 if (numthevts<=0) {
03337 MAXMSG("NuReco",Msg::kInfo,1)
03338 <<"No THEvents, so can't print truth information..."<<endl;
03339 MAXMSG("NuReco",Msg::kInfo,20)
03340 <<"recoEn="<<recoEn<<endl;
03341 return;
03342 }
03343 MAXMSG("NuReco",Msg::kInfo,1)
03344 <<"Found THEvent, printing info..."<<endl;
03345 const NtpTHEvent& thevt=
03346 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
03347
03348
03349 TClonesArray& mcTca=(*ntp.mc);
03350 const NtpMCTruth& mc=
03351 *(dynamic_cast<NtpMCTruth*>(mcTca[thevt.neumc]));
03352 MAXMSG("NuReco",Msg::kDebug,100)
03353 <<"Found mc="<<mc.index<<", thevt.index="<<thevt.index<<endl;
03354
03355 Double_t nuEn=mc.p4neu[3];
03356
03357 Double_t trueEn=mc.y*mc.p4neu[3];
03358 if (mc.iaction==1) trueEn=mc.p4neu[3];
03359
03360 Double_t muEn=(1-mc.y)*mc.p4neu[3];
03361 Double_t hadEn=mc.y*mc.p4neu[3];
03362
03363 string siaction="NC";
03364 if (mc.iaction==1) siaction="CC";
03365
03366 if (mc.iaction==1) {
03367 MAXMSG("NuReco",Msg::kInfo,100)
03368 <<"RecoEn="<<recoEn
03369 <<", "<<siaction
03370 <<" truth: En="<<trueEn<<", mu="<<muEn<<", had="<<hadEn
03371 <<", (y="<<mc.y<<")"
03372 <<endl;
03373 }
03374 else {
03375 MAXMSG("NuReco",Msg::kInfo,20)
03376 <<"RecoEn="<<recoEn
03377 <<", "<<siaction
03378 <<" truth: had="<<hadEn
03379 <<", (Nu="<<nuEn<<", y="<<mc.y<<")"
03380 <<endl;
03381 }
03382 }
03383
03384
03385
03386
03387
03388
03389
03390
03391
03392
03393
03394
03395
03396
03397
03398
03399
03400
03401
03402
03403
03404
03405
03406
03407
03408
03409
03410
03411
03412
03413
03414
03415
03416
03417
03418
03419
03420
03421
03422
03423
03424
03425
03426
03427
03428
03429
03430
03431
03432
03433
03434
03435
03436
03437
03438
03439
03440
03441
03442
03443
03444 Int_t NuReco::GetContainmentFlag(const NuEvent& nu) const
03445 {
03446 Int_t flag=-1;
03447
03448
03449 if (nu.anaVersion==NuCuts::kCC0093Std) {
03450 MAXMSG("NuReco",Msg::kInfo,1)
03451 <<"Using CC0093Std containment criteria"<<endl;
03452 flag=nu.containmentFlagCC0093Std;
03453 }
03454 else if (nu.anaVersion==NuCuts::kCC0250Std) {
03455 MAXMSG("NuReco",Msg::kInfo,1)
03456 <<"Using CC0250Std containment criteria"<<endl;
03457 flag=nu.containmentFlagCC0250Std;
03458 }
03459 else if (nu.anaVersion==NuCuts::kCC0325Std ||
03460 nu.anaVersion==NuCuts::kCC0720Std ||
03461 nu.anaVersion==NuCuts::kCC0720Test||
03462 NuCuts::IsNMBNQ(nu)) {
03463 MAXMSG("NuReco",Msg::kInfo,1)
03464 <<"Using CC0325Std hybrid containment criteria"<<endl;
03465 flag=nu.containmentFlagCC0250Std;
03466
03467
03468 if (nu.detector==Detector::kNear) {
03469 if (nu.xTrkEnd>1.3*Munits::m) {
03470 flag=this->GetContainmentFlagStdReco(nu);
03471 }
03472 }
03473
03474
03475
03476
03477
03478
03479
03480
03481 }
03482 else {
03483 if (ReleaseType::IsBirch(nu.recoVersion)) {
03484 MAXMSG("NuReco",Msg::kInfo,1)
03485 <<"Using Pitt containment criteria"<<endl;
03486
03487 flag=nu.containmentFlagPitt;
03488 }
03489 else {
03490 MAXMSG("NuReco",Msg::kInfo,1)
03491 <<"Using std reco containment criteria"<<endl;
03492 flag=this->GetContainmentFlagStdReco(nu);
03493
03494
03495
03496 }
03497 }
03498
03499
03500 if (!(flag==1 || flag==2 || flag==3 || flag==4)) {
03501 MSG("NuReco",Msg::kWarning)<<"flag="<<flag<<endl;
03502 }
03503
03504
03505 return flag;
03506 }
03507
03508
03509
03510 Int_t NuReco::GetContainmentFlagStdReco(const NuEvent& nu) const
03511 {
03512 Int_t flag=-1;
03513
03514
03515
03516
03517
03518
03519 if (nu.detector==Detector::kNear) {
03520 if (nu.containedTrk) {
03521
03522 if (nu.zTrkEnd<7.0) flag=1;
03523 else flag=3;
03524 }
03525 else {
03526
03527 if (nu.zTrkEnd<7.0) flag=2;
03528 else flag=4;
03529 }
03530 }
03531 else if (nu.detector==Detector::kFar) {
03532
03533 if (nu.containedTrk) flag=1;
03534 else flag=2;
03535 }
03536 else cout<<"Ahh, detector="<<nu.detector<<endl;
03537
03538 return flag;
03539 }
03540
03541
03542
03543 Int_t NuReco::GetContainmentFlagCC0250Std(const NuEvent& nu) const
03544 {
03545 const Float_t x=nu.xTrkEnd;
03546 const Float_t y=nu.yTrkEnd;
03547 const Float_t z=nu.zTrkEnd;
03548 const Float_t r=sqrt(x*x+y*y);
03549
03550
03551
03552 const Float_t xMin=-1.65*(Munits::m);
03553 const Float_t xMax=+2.7*(Munits::m);
03554 const Float_t yMin=-1.65*(Munits::m);
03555 const Float_t yMax=+1.65*(Munits::m);
03556
03557 const Float_t uMax=+3.55*(Munits::m);
03558
03559
03560 const Float_t zMax=+15*(Munits::m);
03561
03562
03563 const Float_t planeMin=4;
03564 const Float_t planeMax=475;
03565 const Float_t rMax=TMath::Sqrt(14)*(Munits::m);
03566
03567 Bool_t contained=false;
03568 Int_t flag=-1;
03569
03570 if (nu.detector==Detector::kNear) {
03571 contained=z<zMax &&
03572 x>xMin && x<xMax &&
03573 y>yMin && y<yMax &&
03574 y>-x+xMin &&
03575 y<+x-xMin &&
03576 y<-x+uMax &&
03577 y>+x-uMax;
03578
03579 if (contained) {
03580
03581
03582
03583
03584
03585
03586 if (z<7) flag=1;
03587 else if (z>=7) flag=3;
03588 }
03589 else {
03590 if (z<7) flag=2;
03591 else if (z>=7) flag=4;
03592 }
03593 }
03594 else if (nu.detector==Detector::kFar) {
03595 contained=nu.planeTrkEnd>=planeMin &&
03596 nu.planeTrkEnd<=planeMax &&
03597 r<rMax;
03598 if (contained) flag=1;
03599 else flag=2;
03600 }
03601 else cout<<"Ahh, detector="<<nu.detector<<endl;
03602
03603 return flag;
03604 }
03605
03606
03607
03608 Int_t NuReco::GetContainmentFlagCC0093Std(const NuEvent& nu) const
03609 {
03610 const Float_t x=nu.xTrkEnd;
03611 const Float_t y=nu.yTrkEnd;
03612 const Float_t z=nu.zTrkEnd;
03613 const Float_t r=sqrt(x*x+y*y);
03614
03615 const Float_t xMin=-1.65;
03616 const Float_t xMax=+2.7;
03617 const Float_t yMin=-1.65;
03618 const Float_t yMax=+1.65;
03619 const Float_t zMax=+16;
03620 const Float_t rMin=+0.4;
03621 const Float_t rMinFD=+0.5;
03622
03623 Bool_t contained=false;
03624 Int_t flag=-1;
03625
03626
03627
03628
03629
03630
03631
03632
03633
03634
03635
03636
03637
03638
03639
03640
03641
03642
03643
03644
03645 if (nu.detector==Detector::kNear) {
03646
03647
03648 contained=z<zMax && r>rMin &&
03649 x>xMin && x<xMax &&
03650 y>yMin && y<yMax &&
03651 y>-x+xMin &&
03652 y<+x-xMin &&
03653 y<-x+3.55 &&
03654 y>+x-3.55;
03655
03656 if (contained) {
03657
03658
03659
03660
03661
03662
03663 if (z<7) flag=1;
03664 else if (z>=7) flag=3;
03665 }
03666 else {
03667 if (z<7) flag=2;
03668 else if (z>=7) flag=4;
03669 }
03670 }
03671 else if (nu.detector==Detector::kFar) {
03672 contained=!(nu.drTrkFidall<rMinFD || nu.dzTrkFidall<rMinFD);
03673 if (contained) flag=1;
03674 else flag=2;
03675 }
03676 else cout<<"Ahh, detector="<<nu.detector<<endl;
03677
03678 return flag;
03679 }
03680
03681
03682
03683 Int_t NuReco::GetContainmentFlagPitt(const NuEvent& nu) const
03684 {
03686
03692
03696
03697 Int_t containmentFlagPitt=0;
03698
03699 const Float_t x=nu.xTrkEnd;
03700 const Float_t y=nu.yTrkEnd;
03701 const Float_t z=nu.zTrkEnd;
03702 const Float_t u=nu.uTrkEnd;
03703 const Float_t v=nu.vTrkEnd;
03704 const Float_t rad=sqrt(x*x + y*y);
03705
03706
03707
03708
03709
03710
03711
03712
03713
03714
03715
03716
03717
03718
03719
03720
03721
03722
03723
03724
03725
03726
03727
03728
03729
03730
03731
03732
03733
03734
03735
03736
03737
03738
03739
03740
03741
03742
03743
03744
03745
03746
03747
03748
03749
03750
03751
03752
03753
03754
03755
03756
03757 if (nu.detector==Detector::kNear) {
03758 Float_t calZ=7*(Munits::m);
03759 Float_t specZ=15.6*(Munits::m);
03760 Float_t minR=0.8*(Munits::m);
03761
03762 Bool_t down_stop=false;
03763 Bool_t down_exit=false;
03764 Bool_t up_stop=false;
03765 Bool_t up_exit=false;
03766
03767
03768 if (z>calZ) {
03769 if (u>-0.85 && u<2.1 &&
03770 v>-2.1 && v<0.85 &&
03771 x>-1.5 && x<2.4 &&
03772 y>-1.4 && y<1.4 &&
03773 rad>minR &&
03774 z>calZ && z<=specZ) down_stop=true;
03775 else down_exit=true;
03776 }
03777
03778 else if (z<=calZ && z>=0) {
03779 if (u>0.3 && u<1.8 &&
03780 v>-1.8 && v<-0.3 &&
03781 x<2.4 &&
03782 rad>minR &&
03783 z<calZ && z>=0) up_stop=true;
03784 else up_exit=true;
03785 }
03786 else cout<<"Ahhh, z="<<z<<endl;
03787
03788
03789 if (down_stop) containmentFlagPitt = 3;
03790 if (up_stop ) containmentFlagPitt = 1;
03791 if (down_exit) containmentFlagPitt = 4;
03792 if (up_exit ) containmentFlagPitt = 2;
03793 }
03794 else if (nu.detector==Detector::kFar) {
03795
03796 const Float_t rMinFD=+0.5;
03797 Bool_t contained=!(nu.drTrkFidall<rMinFD || nu.dzTrkFidall<rMinFD);
03798 if (contained) containmentFlagPitt=1;
03799 else containmentFlagPitt=2;
03800 }
03801 else cout<<"Ahhh, detector="<<nu.detector<<endl;
03802
03803
03804 if (!(containmentFlagPitt==1 || containmentFlagPitt==2 ||
03805 containmentFlagPitt==3 || containmentFlagPitt==4)) {
03806 MSG("NuReco",Msg::kWarning)
03807 <<"containmentFlagPitt="<<containmentFlagPitt<<endl;
03808 }
03809
03810 return containmentFlagPitt;
03811 }
03812
03813
03814
03815
03816
03817
03818
03819
03820
03821
03822
03823
03824
03825
03826
03827
03828
03829
03830
03831
03832
03833
03834
03835
03836
03837
03838
03839
03840
03841
03842
03843
03844
03845
03846
03847
03848
03849
03850
03851
03852
03853
03854
03855
03856
03857
03858
03859
03860
03861
03862
03863
03864
03865
03866
03867 void NuReco::GetPreInukeFinalStateVars(const NtpStRecord& ntp,
03868 NuEvent& nu) const
03869 {
03870 Int_t itr=nu.mc;
03871
03872 TClonesArray* pointStdhepArray = NULL;
03873 pointStdhepArray=ntp.stdhep;
03874 TClonesArray& stdhepArray = *pointStdhepArray;
03875 Int_t nStdHep = stdhepArray.GetEntries();
03876
03877 vector<Int_t> vChildIndices;
03878
03879 for(int i=0;i<nStdHep;i++){
03880
03881 const NtpMCStdHep* ntpStdHep=
03882 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
03883 if(ntpStdHep->mc==itr){
03884
03885 if (ntpStdHep->IstHEP==3){
03886 Int_t childStart = ntpStdHep->child[0];
03887 Int_t childEnd = ntpStdHep->child[1];
03888 for (Int_t goodChild = childStart; goodChild<=childEnd; ++goodChild){
03889 vChildIndices.push_back(goodChild);
03890 }
03891 }
03892 }
03893 }
03894
03895 Int_t numPreInukeFSprot = 0;
03896 Int_t numPreInukeFSneut = 0;
03897 Float_t maxMomPreInukeFSprot = 0.;
03898 Float_t maxMomPreInukeFSneut = 0.;
03899
03900 for (UInt_t i=0;i<vChildIndices.size();i++){
03901
03902 const NtpMCStdHep* ntpStdHep=
03903 dynamic_cast<NtpMCStdHep*>(stdhepArray[vChildIndices[i]]);
03904 if(ntpStdHep->mc==itr){
03905 Float_t energy = ntpStdHep->p4[3];
03906 Float_t mass = ntpStdHep->mass;
03907 Float_t mom = 0.0;
03908 if (energy*energy - mass*mass > 0.0){
03909 mom = TMath::Sqrt(energy*energy - mass*mass);
03910 }
03911 if (ntpStdHep->IdHEP==2212){
03912
03913 numPreInukeFSprot++;
03914 if (mom>maxMomPreInukeFSprot){
03915 maxMomPreInukeFSprot = mom;
03916 }
03917 }
03918 else if (ntpStdHep->IdHEP==2112){
03919
03920 numPreInukeFSneut++;
03921 if (mom>maxMomPreInukeFSneut){
03922 maxMomPreInukeFSneut = mom;
03923 }
03924 }
03925 }
03926 }
03927
03928 nu.numPreInukeFSprotMC = numPreInukeFSprot;
03929 nu.numPreInukeFSneutMC = numPreInukeFSneut;
03930 nu.maxMomPreInukeFSprotMC = maxMomPreInukeFSprot;
03931 nu.maxMomPreInukeFSneutMC = maxMomPreInukeFSneut;
03932
03933 }
03934
03935
03936
03937 Int_t NuReco::GetInitialState(const NtpStRecord& ntp,
03938 const NuEvent& nu) const
03939 {
03941
03942
03943
03944
03945 Int_t itr=nu.mc;
03946
03947 Int_t initial_state=0;
03948 TClonesArray* pointStdhepArray = NULL;
03949 pointStdhepArray=ntp.stdhep;
03950 TClonesArray& stdhepArray = *pointStdhepArray;
03951 Int_t nStdHep = stdhepArray.GetEntries();
03952
03953 int protneut = -1;
03954 int nubarnu = 0;
03955
03956 for(int i=0;i<nStdHep;i++){
03957
03958 const NtpMCStdHep* ntpStdHep=
03959 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
03960 if(ntpStdHep->mc==itr){
03961
03962 if(ntpStdHep->IstHEP==0){
03963 if(abs(ntpStdHep->IdHEP)==12 ||
03964 abs(ntpStdHep->IdHEP)==14 ||
03965 abs(ntpStdHep->IdHEP)==16){
03966 nubarnu = ntpStdHep->IdHEP/abs(ntpStdHep->IdHEP);
03967 }
03968 }
03969 if(ntpStdHep->IstHEP==11){
03970 if(ntpStdHep->IdHEP==2212) protneut = 1;
03971 else if(ntpStdHep->IdHEP==2112) protneut = 0;
03972 else if(abs(ntpStdHep->IdHEP)>1000000000) protneut = 2;
03973 else protneut = 3;
03974 }
03975 }
03976 }
03977
03978 if(protneut==1 && nubarnu==1) initial_state=1;
03979 if(protneut==0 && nubarnu==1) initial_state=2;
03980 if(protneut==1 && nubarnu==-1) initial_state=3;
03981 if(protneut==0 && nubarnu==-1) initial_state=4;
03982 if(protneut==2 && nubarnu==1) initial_state=5;
03983 if(protneut==3 && nubarnu==1) initial_state=6;
03984 if(protneut==2 && nubarnu==-1) initial_state=7;
03985 if(protneut==3 && nubarnu==-1) initial_state=8;
03986
03987 return initial_state;
03988 }
03989
03990
03991
03992 Int_t NuReco::GetNucleus(const NtpStRecord& ,
03993 const NuEvent& nu) const
03994 {
03996
03997
03998 Int_t z=nu.zMC;
03999 Int_t a=nu.aMC;
04000 Int_t nucleus=1;
04001
04002 switch (z) {
04003
04004
04005
04006 case 1:
04007 switch (a) {
04008 case 1:
04009 nucleus=256;
04010 break;
04011 case 2:
04012 nucleus=257;
04013 break;
04014 case 3:
04015 nucleus=258;
04016 break;
04017 default:
04018 nucleus=256;
04019 break;
04020 }
04021 break;
04022 case 6:
04023 switch (a) {
04024 case 11:
04025 nucleus=274;
04026 break;
04027 case 12:
04028 nucleus=275;
04029 break;
04030 case 13:
04031 nucleus=276;
04032 break;
04033 case 14:
04034 nucleus=277;
04035 break;
04036 default:
04037 nucleus=275;
04038 break;
04039 }
04040 break;
04041 case 7:
04042 switch (a) {
04043 case 13:
04044 nucleus=278;
04045 break;
04046 case 14:
04047 nucleus=279;
04048 break;
04049 case 15:
04050 nucleus=280;
04051 break;
04052 case 16:
04053 nucleus=281;
04054 break;
04055 case 17:
04056 nucleus=282;
04057 break;
04058 default:
04059 nucleus=279;
04060 break;
04061 }
04062 break;
04063 case 8:
04064 switch (a) {
04065 case 15:
04066 nucleus=283;
04067 break;
04068 case 16:
04069 nucleus=284;
04070 break;
04071 case 17:
04072 nucleus=285;
04073 break;
04074 case 18:
04075 nucleus=286;
04076 break;
04077 default:
04078 nucleus=284;
04079 break;
04080 }
04081 break;
04082 case 13:
04083 switch (a) {
04084 case 26:
04085 nucleus=303;
04086 break;
04087 case 27:
04088 nucleus=304;
04089 break;
04090 case 28:
04091 nucleus=305;
04092 break;
04093 case 29:
04094 nucleus=306;
04095 break;
04096 default:
04097 nucleus=304;
04098 break;
04099 }
04100 break;
04101 case 14:
04102 switch (a) {
04103 case 27:
04104 nucleus=307;
04105 break;
04106 case 28:
04107 nucleus=308;
04108 break;
04109 case 29:
04110 nucleus=309;
04111 break;
04112 case 30:
04113 nucleus=310;
04114 break;
04115 default:
04116 nucleus=308;
04117 break;
04118 }
04119 break;
04120 case 15:
04121 switch (a) {
04122 case 30:
04123 nucleus=311;
04124 break;
04125 case 31:
04126 nucleus=312;
04127 break;
04128 case 32:
04129 nucleus=313;
04130 break;
04131 case 33:
04132 nucleus=314;
04133 break;
04134 default:
04135 nucleus=312;
04136 break;
04137 }
04138 break;
04139 case 16:
04140 switch (a) {
04141 case 31:
04142 nucleus=315;
04143 break;
04144 case 32:
04145 nucleus=316;
04146 break;
04147 case 33:
04148 nucleus=317;
04149 break;
04150 case 34:
04151 nucleus=318;
04152 break;
04153 case 35:
04154 nucleus=319;
04155 break;
04156 case 36:
04157 nucleus=320;
04158 break;
04159 default:
04160 nucleus=316;
04161 break;
04162 }
04163 break;
04164 case 22:
04165 switch (a) {
04166 case 45:
04167 nucleus=347;
04168 break;
04169 case 46:
04170 nucleus=348;
04171 break;
04172 case 47:
04173 nucleus=349;
04174 break;
04175 case 48:
04176 nucleus=350;
04177 break;
04178 case 49:
04179 nucleus=351;
04180 break;
04181 case 50:
04182 nucleus=352;
04183 break;
04184 default:
04185 nucleus=350;
04186 break;
04187 }
04188 break;
04189 case 23:
04190 switch (a) {
04191 case 49:
04192 nucleus=353;
04193 break;
04194 case 50:
04195 nucleus=354;
04196 break;
04197 case 51:
04198 nucleus=355;
04199 break;
04200 case 52:
04201 nucleus=356;
04202 break;
04203 case 53:
04204 nucleus=357;
04205 break;
04206 default:
04207 nucleus=355;
04208 break;
04209 }
04210 break;
04211 case 24:
04212 switch (a) {
04213 case 49:
04214 nucleus=358;
04215 break;
04216 case 50:
04217 nucleus=359;
04218 break;
04219 case 51:
04220 nucleus=360;
04221 break;
04222 case 52:
04223 nucleus=361;
04224 break;
04225 case 53:
04226 nucleus=362;
04227 break;
04228 case 54:
04229 nucleus=363;
04230 break;
04231 default:
04232 nucleus=361;
04233 break;
04234 }
04235 break;
04236 case 25:
04237 switch (a) {
04238 case 53:
04239 nucleus=364;
04240 break;
04241 case 54:
04242 nucleus=365;
04243 break;
04244 case 55:
04245 nucleus=366;
04246 break;
04247 case 56:
04248 nucleus=367;
04249 break;
04250 case 57:
04251 nucleus=368;
04252 break;
04253 default:
04254 nucleus=366;
04255 break;
04256 }
04257 break;
04258 case 26:
04259 switch (a) {
04260 case 53:
04261 nucleus=369;
04262 break;
04263 case 54:
04264 nucleus=370;
04265 break;
04266 case 55:
04267 nucleus=371;
04268 break;
04269 case 56:
04270 nucleus=372;
04271 break;
04272 case 57:
04273 nucleus=373;
04274 break;
04275 case 58:
04276 nucleus=374;
04277 break;
04278 default:
04279 nucleus=372;
04280 break;
04281 }
04282 break;
04283 case 28:
04284 switch (a) {
04285 case 57:
04286 nucleus=382;
04287 break;
04288 case 58:
04289 nucleus=383;
04290 break;
04291 case 59:
04292 nucleus=384;
04293 break;
04294 case 60:
04295 nucleus=385;
04296 break;
04297 case 61:
04298 nucleus=386;
04299 break;
04300 case 62:
04301 nucleus=387;
04302 break;
04303 case 63:
04304 nucleus=388;
04305 break;
04306 case 64:
04307 nucleus=389;
04308 break;
04309 default:
04310 nucleus=383;
04311 break;
04312 }
04313 break;
04314 case 29:
04315 switch (a) {
04316 case 62:
04317 nucleus=390;
04318 break;
04319 case 63:
04320 nucleus=391;
04321 break;
04322 case 64:
04323 nucleus=392;
04324 break;
04325 case 65:
04326 nucleus=393;
04327 break;
04328 case 66:
04329 nucleus=394;
04330 break;
04331 case 67:
04332 nucleus=395;
04333 break;
04334 default:
04335 nucleus=392;
04336 break;
04337 }
04338 break;
04339 default:
04340 nucleus=1;
04341 break;
04342 }
04343
04344 return nucleus;
04345 }
04346
04347
04348
04349 Int_t NuReco::GetHadronicFinalState(const NtpStRecord& ntp,
04350 const NuEvent& nu) const
04351 {
04353
04354 Int_t hfs=0;
04355 Int_t proc=nu.iresonance;
04356 TClonesArray* pointStdhepArray = NULL;
04357 pointStdhepArray=ntp.stdhep;
04358 TClonesArray& stdhepArray = *pointStdhepArray;
04359 Int_t nStdHep = stdhepArray.GetEntries();
04360
04361 Int_t itr=nu.mc;
04362
04363 if(proc==1002){
04364 for(int i=0;i<nStdHep;i++){
04365
04366 const NtpMCStdHep* ntpStdHep=
04367 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
04368
04369 if(ntpStdHep->mc==itr && ntpStdHep->IstHEP==3 &&
04370 !(abs(ntpStdHep->IdHEP)==15)){
04371 hfs = ntpStdHep->IdHEP;
04372 break;
04373 }
04374 }
04375 hfs = hfs%1000;
04376 }
04377 else {
04378 for(int i=0;i<nStdHep;i++){
04379
04380 const NtpMCStdHep* ntpStdHep=
04381 dynamic_cast<NtpMCStdHep*>(stdhepArray[i]);
04382 if(ntpStdHep->mc==itr && ntpStdHep->IstHEP==3){
04383 hfs = ntpStdHep->IdHEP;
04384 break;
04385 }
04386 }
04387 hfs = hfs%1000;
04388 }
04389 return hfs;
04390 }
04391
04392
04393
04394 Float_t NuReco::GetDirCosNu(const NuEvent& nu) const
04395 {
04397
04398 Float_t dirCosNu=-999;
04399 if (nu.detector==Detector::kFar) {
04400 const Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.);
04401 const Float_t bl_y = sqrt(1. - bl_z*bl_z);
04402 dirCosNu=nu.trkvtxdcosz*bl_z+nu.trkvtxdcosy*bl_y;
04403 }
04404 else if (nu.detector==Detector::kNear) {
04405 const Float_t nu_cos = -5.799E-2;
04406 const Float_t nu_sin = sqrt(1 -nu_cos*nu_cos);
04407 dirCosNu=nu.trkvtxdcosz*nu_sin + nu.trkvtxdcosy*nu_cos;
04408 }
04409 else cout<<"Ahhh, detector="<<nu.detector<<endl;
04410
04411
04412 return dirCosNu;
04413 }
04414
04415
04416
04417 void NuReco::GetDirCosNu(const NtpSRTrack& trk,NuEvent& nu) const
04418 {
04419
04420
04421
04422
04423
04424
04425
04426
04427
04428
04429
04430
04431
04432 if (nu.detector==Detector::kFar) {
04433 const Float_t bl_y = 0.057184957;
04434 const Float_t bl_z = sqrt(1 - bl_y*bl_y);
04435 nu.dirCosNu=trk.vtx.dcosz*bl_z+trk.vtx.dcosy*bl_y;
04436 }
04437 else if (nu.detector==Detector::kNear) {
04438
04439
04440
04441
04442
04443
04444
04445
04446
04447
04448
04449
04450
04451 const Float_t bl_y = -0.058297760;
04452 const Float_t bl_z = sqrt(1 - bl_y*bl_y);
04453 nu.dirCosNu=trk.vtx.dcosz*bl_z+trk.vtx.dcosy*bl_y;
04454 }
04455 else cout<<"Ahhh, detector="<<nu.detector<<endl;
04456 }
04457
04458
04459
04460 Float_t NuReco::GetRadius(Float_t x,Float_t y,const NuEvent& nu) const
04461 {
04462 Float_t zerox=0;
04463 Float_t zeroy=0;
04464 if (nu.detector==Detector::kNear) {
04465
04466 if (nu.anaVersion==NuCuts::kCC0093Std) {
04467 zerox=1.4885*(Munits::m);
04468 zeroy=0.1397*(Munits::m);
04469 }
04470 else {
04471 zerox=1.4828*(Munits::m);
04472 zeroy=0.2384*(Munits::m);
04473 }
04474 }
04475
04476 return sqrt(pow((x-zerox),2)+pow((y-zeroy),2));
04477 }
04478
04479
04480
04481 Float_t NuReco::GetSmallestDeepDistToEdge(const NuEvent& nu) const
04482 {
04483
04484
04485
04486 Float_t distUPI=0;
04487 Float_t distVPI=0;
04488 Float_t distUFI=0;
04489 Float_t distVFI=0;
04490 Float_t xedge=0;
04491 Float_t yedge=0;
04492
04493 PlaneOutline po;
04494
04495
04496 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
04497 PlaneCoverage::kNearPartial,
04498 distUPI,xedge,yedge);
04499 Bool_t isInsideUPI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
04500 PlaneCoverage::kNearPartial);
04501
04502
04503 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
04504 PlaneCoverage::kNearFull,
04505 distUFI,xedge,yedge);
04506 Bool_t isInsideUFI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kU,
04507 PlaneCoverage::kNearFull);
04508
04509
04510 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
04511 PlaneCoverage::kNearPartial,
04512 distVPI,xedge,yedge);
04513 Bool_t isInsideVPI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
04514 PlaneCoverage::kNearPartial);
04515
04516
04517 po.DistanceToEdge(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
04518 PlaneCoverage::kNearFull,
04519 distVFI,xedge,yedge);
04520 Bool_t isInsideVFI=po.IsInside(nu.xEvtVtx,nu.yEvtVtx,PlaneView::kV,
04521 PlaneCoverage::kNearFull);
04522
04523
04524
04525 Float_t smallestDist=999;
04526 if (isInsideUPI && isInsideUFI && isInsideVPI && isInsideVFI){
04527 if (smallestDist>distUPI) smallestDist=distUPI;
04528 if (smallestDist>distUFI) smallestDist=distUFI;
04529 if (smallestDist>distVPI) smallestDist=distVPI;
04530 if (smallestDist>distVFI) smallestDist=distVFI;
04531 }
04532
04533
04534
04535
04536 if (smallestDist==999) {
04537 smallestDist=0;
04538 if (smallestDist<distUPI && !isInsideUPI) smallestDist=distUPI;
04539 if (smallestDist<distUFI && !isInsideUFI) smallestDist=distUFI;
04540 if (smallestDist<distVPI && !isInsideVPI) smallestDist=distVPI;
04541 if (smallestDist<distVFI && !isInsideVFI) smallestDist=distVFI;
04542
04543
04544 smallestDist*=-1;
04545 }
04546
04547 if (smallestDist==999) cout<<"Ahhhhhhh"<<endl;
04548
04549 MAXMSG("NuReco",Msg::kDebug,500)
04550 <<"smallestDist="<<smallestDist
04551 <<", UPI="<<distUPI<<" (in="<<isInsideUPI<<")"
04552 <<", UFI="<<distUFI<<" (in="<<isInsideUFI<<")"
04553 <<", VPI="<<distVPI<<" (in="<<isInsideVPI<<")"
04554 <<", VFI="<<distVFI<<" (in="<<isInsideVFI<<")"<<endl;
04555
04556 return smallestDist;
04557 }
04558
04559
04560
04561 void NuReco::TestGetSmallestDeepDistToEdge() const
04562 {
04563 NuEvent nu;
04564
04565 TH2F* hYvsXDistToEdge=new TH2F
04566 ("hYvsXDistToEdge","hYvsXDistToEdge",
04567 600,-3,3, 600,-3,3);
04568 hYvsXDistToEdge->SetTitle("SmallestDeepDistToEdge");
04569 hYvsXDistToEdge->GetXaxis()->SetTitle("X (m)");
04570 hYvsXDistToEdge->GetXaxis()->CenterTitle();
04571 hYvsXDistToEdge->GetYaxis()->SetTitle("Y (m)");
04572 hYvsXDistToEdge->GetYaxis()->CenterTitle();
04573
04574 for (Float_t x=-3;x<3;x+=0.01) {
04575 for (Float_t y=-3;y<3;y+=0.01) {
04576
04577 nu.xEvtVtx=x;
04578 nu.yEvtVtx=y;
04579 Float_t dist=this->GetSmallestDeepDistToEdge(nu);
04580
04581 hYvsXDistToEdge->Fill(x,y,dist);
04582 }
04583 }
04584 }
04585
04586
04587
04588 TVector3 NuReco::xyz2uvz(Float_t x,Float_t y,Float_t z,
04589 VldContext vc) const
04590 {
04591
04592 UgliGeomHandle ugh(vc);
04593
04594
04595 TVector3 xyz(x,y,z);
04596 TVector3 uvz=ugh.xyz2uvz(xyz);
04597 return uvz;
04598 }
04599
04600
04601
04602 TVector3 NuReco::uvz2xyz(Float_t u,Float_t v,Float_t z,
04603 VldContext vc) const
04604 {
04605
04606 UgliGeomHandle ugh(vc);
04607
04608
04609 TVector3 uvz(u,v,z);
04610 TVector3 xyz=ugh.uvz2xyz(uvz);
04611 return xyz;
04612 }
04613
04614
04615
04616
04617 Double_t NuReco::GetShowerEnergyNearTrack(const NtpStRecord &ntp,
04618 const NtpSRTrack &trk,
04619 Float_t* nearshowereDW) const
04620 {
04621 Double_t nearshowere = 0;
04622
04623
04624 if(nearshowereDW) *nearshowereDW = 0;
04625
04626 bool hasgoodshower = false;
04627
04628 const double trkTime = trk.vtx.t;
04629
04630 const Detector::Detector_t det = ntp.GetHeader().GetVldContext().GetDetector();
04631
04632
04633 for(Int_t it = 0; it < ntp.shw->GetEntries(); it++){
04634
04635 const NtpSRShower *shw = dynamic_cast<NtpSRShower *>(ntp.shw->At(it));
04636
04637
04638 if(det == Detector::kNear){
04639 const double shwTime = GetShwMedianTime(ntp, *shw);
04640 const double dt = (shwTime-trkTime)*1e9;
04641 if(dt > +75 || dt < -25) continue;
04642 }
04643
04644
04645 Double_t vertex_sep_sqr = (trk.vtx.x-shw->vtx.x)*(trk.vtx.x-shw->vtx.x) +
04646 (trk.vtx.y-shw->vtx.y)*(trk.vtx.y-shw->vtx.y) +
04647 (trk.vtx.z-shw->vtx.z)*(trk.vtx.z-shw->vtx.z);
04648
04649
04650 if(vertex_sep_sqr < 1.0){
04651 nearshowere += shw->shwph.linCCgev;
04652 if(nearshowereDW) *nearshowereDW += shw->shwph.wtCCgev;
04653 hasgoodshower = true;
04654 }
04655 }
04656
04657
04658 if(!hasgoodshower && nearshowereDW) *nearshowereDW = -1;
04659
04660 return nearshowere;
04661 }
04662
04663
04664
04665
04666 Double_t NuReco::GetCosBetweenPr_Theta(const NtpSRTrack& trk) const
04667 {
04668 Double_t x = trk.vtx.x,
04669 y = trk.vtx.y,
04670
04671 px = trk.momentum.best * trk.vtx.dcosx,
04672 py = trk.momentum.best * trk.vtx.dcosy;
04673
04674
04675 Double_t magxmagp = sqrt((x*x + y*y)*(px*px + py*py));
04676
04677
04678 if(magxmagp != 0) return (x*px + y*py)/magxmagp;
04679 else return 0;
04680 }
04681
04682
04683
04684 Int_t NuReco::GetTrackFirstStrip(const NtpStRecord &ntp, const NtpSRTrack &trk) const
04685 {
04686 NtpSRStrip * firststrip =
04687 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[0]));
04688 return firststrip->strip;
04689 }
04690
04691
04692
04693 Int_t NuReco::GetTrackFirstUStrip(const NtpStRecord &ntp, const NtpSRTrack &trk) const
04694 {
04695
04696 for(int i = 0; i < trk.nstrip; i++) {
04697 NtpSRStrip * firststrip =
04698 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[i]));
04699 if(firststrip->planeview == PlaneView::kU) return firststrip->strip;
04700 }
04701
04702 return -1;
04703 }
04704
04705
04706
04707 Int_t NuReco::GetTrackFirstVStrip(const NtpStRecord &ntp, const NtpSRTrack &trk) const
04708 {
04709
04710 for(int i = 0; i < trk.nstrip; i++) {
04711 NtpSRStrip * firststrip =
04712 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[i]));
04713 if(firststrip->planeview == PlaneView::kV) return firststrip->strip;
04714 }
04715
04716 return -1;
04717 }
04718
04719
04720
04721 Bool_t NuReco::GetTrackFirstStripIsU(const NtpStRecord &ntp, const NtpSRTrack &trk) const
04722 {
04723 NtpSRStrip * firststrip =
04724 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[0]));
04725 if(firststrip->planeview == PlaneView::kU) return true;
04726 else return false;
04727 }
04728
04729
04730
04731 Int_t NuReco::GetTrackLastStrip(const NtpStRecord &ntp, const NtpSRTrack &trk) const
04732 {
04733 NtpSRStrip * laststrip =
04734 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[trk.nstrip-1]));
04735 return laststrip->strip;
04736 }
04737
04738
04739
04740 Int_t NuReco::GetTrackLastUStrip(const NtpStRecord &ntp, const NtpSRTrack &trk) const
04741 {
04742
04743 for(Int_t i = trk.nstrip-1; i >= 0; i--) {
04744 NtpSRStrip * laststrip =
04745 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[i]));
04746 if(laststrip->planeview == PlaneView::kU) return laststrip->strip;
04747 }
04748
04749 return -1;
04750 }
04751
04752
04753
04754 Int_t NuReco::GetTrackLastVStrip(const NtpStRecord &ntp, const NtpSRTrack &trk) const
04755 {
04756
04757 for(Int_t i = trk.nstrip-1; i >= 0; i--) {
04758 NtpSRStrip * laststrip =
04759 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[i]));
04760 if(laststrip->planeview == PlaneView::kV) return laststrip->strip;
04761 }
04762
04763 return -1;
04764 }
04765
04766
04767
04768 Bool_t NuReco::GetTrackLastStripIsU(const NtpStRecord &ntp, const NtpSRTrack &trk) const
04769 {
04770 NtpSRStrip * laststrip =
04771 dynamic_cast<NtpSRStrip *>(ntp.stp->At(trk.stp[trk.nstrip-1]));
04772 if(laststrip->planeview == PlaneView::kU) return true;
04773 else return false;
04774 }
04775
04776
04777
04778
04779
04780 Int_t NuReco::GetRegion(const NtpStRecord &rec, const NtpSRVertex& vtx) const
04781 {
04782
04783 Detector::Detector_t det = rec.GetHeader().GetVldContext().GetDetector();
04784 SimFlag::SimFlag_t simFlag = rec.GetHeader().GetVldContext().GetSimFlag();
04785
04786 return GetRegion(det, simFlag, vtx.x, vtx.y, vtx.z, vtx.plane);
04787 }
04788
04789
04790
04791
04792
04793 Int_t NuReco::GetRegion(Detector::Detector_t det, SimFlag::SimFlag_t simFlag, Float_t vtxX, Float_t vtxY, Float_t vtxZ, Int_t vtxPlane) const
04794 {
04795
04796
04797
04798
04799
04800
04801
04802
04803
04804
04805
04806
04807
04808
04809
04810
04811
04812
04813
04814
04815
04816
04817
04818
04819
04820
04821
04822
04823
04824
04825
04826
04827
04828
04829
04830
04831
04832
04833
04834
04835
04836
04837
04838
04839
04840
04841
04842
04843
04844
04845
04846
04847
04848
04849
04850
04851
04852
04853
04854
04855 enum Rock_DetectorRegion {
04856 kUnknown = -2,
04857 kNoTrack = -1,
04858 kFiducial = 0,
04859 kEdge = 1,
04860 kCoilHole = 2,
04861 kFront = 4,
04862 kGapSouth = 8,
04863 kGapNorth = 12,
04864 kBack = 16
04865 };
04866
04867
04868
04869 Bool_t nearFront = false, nearEdge = false;
04870 Bool_t nearGapSouth = false, nearGapNorth = false;
04871 Bool_t nearCoil = false, nearBack = false;
04872 Double_t radius2 = vtxX*vtxX + vtxY*vtxY;
04873 Double_t z = vtxZ - FidVol::getTrkVtxZOffset();
04874
04875
04876 if (det != Detector::kFar) return kUnknown;
04877
04878
04879
04880
04881
04882 Double_t getFarZ[4] = {};
04883 if (simFlag == SimFlag::kData) {
04884 getFarZ[0] = FidVol::getFarZData(0);
04885 getFarZ[1] = FidVol::getFarZData(1);
04886 getFarZ[2] = FidVol::getFarZData(2);
04887 getFarZ[3] = FidVol::getFarZData(3);
04888 } else {
04889 getFarZ[0] = FidVol::getFarZMC(0);
04890 getFarZ[1] = FidVol::getFarZMC(1);
04891 getFarZ[2] = FidVol::getFarZMC(2);
04892 getFarZ[3] = FidVol::getFarZMC(3);
04893 }
04894
04895
04896
04897 const int GAPPLANE = 249;
04898
04899 nearEdge = radius2 >= (FidVol::getFarRouter()*FidVol::getFarRouter());
04900 nearCoil = radius2 < FidVol::getFarRinner()*FidVol::getFarRinner();
04901 nearFront = z < getFarZ[0];
04902 nearGapSouth = z > getFarZ[1] && vtxPlane < GAPPLANE;
04903 nearGapNorth = vtxPlane > GAPPLANE && z < getFarZ[2];
04904 nearBack = vtxZ - FidVol::getEvtVtxZOffset() > getFarZ[3];
04905
04906 return kEdge * nearEdge +
04907 kCoilHole * nearCoil +
04908 kFront * nearFront +
04909 kGapSouth * nearGapSouth +
04910 kGapNorth * nearGapNorth +
04911 kBack * nearBack;
04912 }
04913
04914
04915
04916
04917 Int_t NuReco::GetPrimaryMuonStdHepIndex(const NtpStRecord& ntp, const NtpMCTruth& mc) const
04918 {
04919
04920 if (mc.iaction == 0) return -1;
04921
04922
04923 if (mc.inu != 14 && mc.inu != -14) return -1;
04924
04925
04926 TClonesArray& hepTca=(*ntp.stdhep);
04927
04928
04929 const NtpMCStdHep& neutrino = *dynamic_cast<NtpMCStdHep*>(hepTca[mc.stdhep[0]]);
04930
04931 const NtpMCStdHep& child = *dynamic_cast<NtpMCStdHep*>(hepTca[neutrino.child[0]]);
04932
04933
04934 Int_t childcount = neutrino.child[1] - neutrino.child[0];
04935 if (childcount > 0) {
04936
04937 this->PrintStdhep(ntp, mc);
04938 MSG("NuReco",Msg::kFatal) << "More than one primary child in interaction" << endl;
04939 return 0;
04940 }
04941
04942
04943
04944
04945
04946
04947
04948
04949
04950 if (child.IdHEP != 13 && child.IdHEP != -13) {
04951 MSG("NuReco", Msg::kFatal) << "Not a muon child" << endl;
04952 }
04953
04954
04955 return neutrino.child[0];
04956
04957
04958
04959
04960
04961
04962
04963
04964
04965
04966
04967
04968
04969
04970
04971
04972
04973
04974
04975
04976
04977
04978
04979
04980
04981
04982
04983
04984
04985
04986
04987
04988
04989 }
04990
04991
04992
04993 Double_t NuReco::GetMuonFirstHitEnergy(const NtpStRecord& ntp, const NtpMCTruth& mc) const
04994 {
04995 Int_t stdindex = this->GetPrimaryMuonStdHepIndex(ntp, mc);
04996 if (stdindex == -1) return -1;
04997
04998
04999 TClonesArray& hepTca=(*ntp.stdhep);
05000 const NtpMCStdHep& muon = *dynamic_cast<NtpMCStdHep*>(hepTca[stdindex]);
05001
05002
05003 return muon.dethit[0].mom[3];
05004 }
05005
05006
05007
05008 Double_t NuReco::GetMuonLastHitEnergy(const NtpStRecord& ntp, const NtpMCTruth& mc) const
05009 {
05010 Int_t stdindex = this->GetPrimaryMuonStdHepIndex(ntp, mc);
05011 if (stdindex == -1) return -1;
05012
05013
05014 TClonesArray& hepTca=(*ntp.stdhep);
05015 const NtpMCStdHep& muon = *dynamic_cast<NtpMCStdHep*>(hepTca[stdindex]);
05016
05017
05018 return muon.dethit[1].mom[3];
05019 }
05020
05021
05022
05023 Double_t NuReco::GetTrackTransverseMomentum(const NuEvent& nu) const
05024 {
05025 Double_t bestmomentum = -1;
05026
05027 return bestmomentum*TMath::Sin(TMath::ACos(nu.dirCosNu));
05028 }
05029
05030
05031
05032 Bool_t NuReco::GetMuonContainmentMC(const NtpStRecord& ntp, const NtpMCTruth& mc) const
05033 {
05034 Int_t muonindex = this->GetPrimaryMuonStdHepIndex(ntp, mc);
05035
05036 if (muonindex == -1) return false;
05037
05038
05039 TClonesArray& hepTca=(*ntp.stdhep);
05040
05041
05042 NtpMCStdHep* muon = dynamic_cast<NtpMCStdHep*>(hepTca[muonindex]);
05043
05044
05045
05046
05047
05048
05049
05050 Bool_t containment = false;
05051
05052
05053 while (true) {
05054
05055 NtpMCStdHep* child = 0;
05056
05057
05058 if (muon->child[0] != -1) {
05059
05060
05061 for (Int_t i = muon->child[0]; i <= muon->child[1]; ++i) {
05062
05063 NtpMCStdHep* childcand = dynamic_cast<NtpMCStdHep*>(hepTca[i]);
05064
05065
05066 if (childcand->IdHEP == 13 || childcand->IdHEP == -13) {
05067
05068 child = childcand;
05069 }
05070 }
05071
05072
05073
05074 if (child == 0) {
05075
05076 containment = true;
05077 break;
05078 } else {
05079
05080 muon = child;
05081 }
05082 } else {
05083
05084 containment = false;
05085 break;
05086 }
05087 }
05088
05089
05090
05091
05092
05093
05094
05095
05096
05097 return containment;
05098 }
05099
05100
05101
05102 Float_t NuReco::GetPhi(const NtpSRVertex &vtx) const
05103 {
05104 return TMath::ATan2(vtx.y, vtx.x);
05105 }
05106
05107
05108
05109 Int_t NuReco::GetEdgeRegion(const NtpStRecord &rec, const NtpSRVertex& vtx) const
05110 {
05111
05112 Detector::Detector_t det = rec.GetHeader().GetVldContext().GetDetector();
05113
05114
05115
05116
05117 return GetEdgeRegion(det, GetPhi(vtx));
05118 }
05119
05120
05121
05122 Int_t NuReco::GetEdgeRegion(Detector::Detector_t det, const Float_t phi) const
05123 {
05124
05125
05126
05127
05128
05129
05130
05131
05132
05133
05134
05135 if (det == Detector::kNear) return -1;
05136
05137 Float_t positiveangle = phi;
05138
05139 if(positiveangle < 0) positiveangle += 2*M_PI;
05140 const Float_t torad = M_PI/180;
05141
05142 Float_t theoneang = atan(1.625/3.975);
05143 if(positiveangle < theoneang) return 0;
05144 if(positiveangle < 90*torad - theoneang) return 1;
05145 if(positiveangle < 90*torad + theoneang) return 2;
05146 if(positiveangle < 180*torad - theoneang) return 3;
05147 if(positiveangle < 180*torad + theoneang) return 4;
05148 if(positiveangle < 270*torad - theoneang) return 5;
05149 if(positiveangle < 270*torad + theoneang) return 6;
05150 if(positiveangle < 360*torad - theoneang) return 7;
05151
05152
05153 return 0;
05154 }
05155
05156
05157
05158 void NuReco::CalculateEdgeRegion(NuEvent &nu, Int_t track) const
05159 {
05160 Detector::Detector_t det = static_cast<Detector::Detector_t>(nu.detector);
05161
05162 if (track == 1) {
05163 nu.edgeRegionTrkVtx1 = GetEdgeRegion(det, nu.phiTrkVtx1);
05164 nu.edgeRegionTrkEnd1 = GetEdgeRegion(det, nu.phiTrkEnd1);
05165 } else if (track == 2) {
05166 nu.edgeRegionTrkVtx2 = GetEdgeRegion(det, nu.phiTrkVtx2);
05167 nu.edgeRegionTrkEnd2 = GetEdgeRegion(det, nu.phiTrkEnd2);
05168 } else if (track == 3) {
05169 nu.edgeRegionTrkVtx3 = GetEdgeRegion(det, nu.phiTrkVtx3);
05170 nu.edgeRegionTrkEnd3 = GetEdgeRegion(det, nu.phiTrkEnd3);
05171 } else if (track == 0){
05172
05173 nu.edgeRegionTrkVtx = GetEdgeRegion(det, nu.phiTrkVtx);
05174 nu.edgeRegionTrkEnd = GetEdgeRegion(det, nu.phiTrkEnd);
05175 } else {
05176 MSG("NuReco",Msg::kError) << "Unrecognised track id " << track << endl;
05177 }
05178 }
05179
05180
05181
05182 Short_t NuReco::GetParallelStripInset(Detector::Detector_t det, Short_t edgeRegion, Short_t stripBegU, Short_t stripBegV) const
05183 {
05184 Short_t parallelStripVtx = -10;
05185
05186
05187 if (det == Detector::kNear) return parallelStripVtx;
05188
05189 if (edgeRegion == 1) parallelStripVtx = 191-stripBegU;
05190 else if(edgeRegion == 3) parallelStripVtx = 191-stripBegV;
05191 else if(edgeRegion == 5) parallelStripVtx = stripBegU;
05192 else if(edgeRegion == 7) parallelStripVtx = stripBegV;
05193 else parallelStripVtx = -10;
05194
05195 return parallelStripVtx;
05196 }
05197
05198
05199
05200 void NuReco::CalculateParallelStripInset(NuEvent &nu, Int_t track) const
05201 {
05202 Detector::Detector_t det = static_cast<Detector::Detector_t>(nu.detector);
05203
05204 if (track == 1) {
05205 nu.parallelStripTrkVtx1 = GetParallelStripInset(det, nu.edgeRegionTrkVtx1, nu.stripTrkBegu1, nu.stripTrkBegv1);
05206 nu.parallelStripTrkEnd1 = GetParallelStripInset(det, nu.edgeRegionTrkEnd1, nu.stripTrkEndu1, nu.stripTrkEndv1);
05207 } else if (track == 2) {
05208 nu.parallelStripTrkVtx2 = GetParallelStripInset(det, nu.edgeRegionTrkVtx2, nu.stripTrkBegu2, nu.stripTrkBegv2);
05209 nu.parallelStripTrkEnd2 = GetParallelStripInset(det, nu.edgeRegionTrkEnd2, nu.stripTrkEndu2, nu.stripTrkEndv2);
05210 } else if (track == 3) {
05211 nu.parallelStripTrkVtx3 = GetParallelStripInset(det, nu.edgeRegionTrkVtx3, nu.stripTrkBegu3, nu.stripTrkBegv3);
05212 nu.parallelStripTrkEnd3 = GetParallelStripInset(det, nu.edgeRegionTrkEnd3, nu.stripTrkEndu3, nu.stripTrkEndv3);
05213 } else if (track == 0){
05214
05215 nu.parallelStripTrkVtx = GetParallelStripInset(det, nu.edgeRegionTrkVtx, nu.stripTrkBegu, nu.stripTrkBegv);
05216 nu.parallelStripTrkEnd = GetParallelStripInset(det, nu.edgeRegionTrkEnd, nu.stripTrkEndu, nu.stripTrkEndv);
05217 } else {
05218 MSG("NuReco",Msg::kError) << "Unrecognised track id " << track << endl;
05219 }
05220 }
05221
05222
05223
05224 Char_t NuReco::GetPerpendicularStripFlag(Detector::Detector_t det, Short_t edgeRegion, Bool_t IsHitu) const
05225 {
05226
05227 const Bool_t Isu = IsHitu;
05228 const Bool_t Isv = !Isu;
05229
05230 Char_t flag = -10;
05231
05232
05233 if (det == Detector::kNear) return flag;
05234
05235 if(edgeRegion == 1 || edgeRegion == 5) flag = Isv;
05236 else if(edgeRegion == 3 || edgeRegion == 7) flag = Isu;
05237
05238 return flag;
05239 }
05240
05241
05242
05243 void NuReco::CalculatePerpendicularStripFlag(NuEvent &nu, Int_t track) const
05244 {
05245 Detector::Detector_t det = static_cast<Detector::Detector_t>(nu.detector);
05246
05247 if (track == 1) {
05248 nu.stripTrkBegPerpFlag1 = GetPerpendicularStripFlag(det, nu.edgeRegionTrkVtx1, nu.stripTrkBegIsu1);
05249 nu.stripTrkEndPerpFlag1 = GetPerpendicularStripFlag(det, nu.edgeRegionTrkEnd1, nu.stripTrkEndIsu1);
05250 } else if (track == 2) {
05251 nu.stripTrkBegPerpFlag2 = GetPerpendicularStripFlag(det, nu.edgeRegionTrkVtx2, nu.stripTrkBegIsu2);
05252 nu.stripTrkEndPerpFlag2 = GetPerpendicularStripFlag(det, nu.edgeRegionTrkEnd2, nu.stripTrkEndIsu2);
05253 } else if (track == 3) {
05254 nu.stripTrkBegPerpFlag3 = GetPerpendicularStripFlag(det, nu.edgeRegionTrkVtx3, nu.stripTrkBegIsu3);
05255 nu.stripTrkEndPerpFlag3 = GetPerpendicularStripFlag(det, nu.edgeRegionTrkEnd3, nu.stripTrkEndIsu3);
05256 } else if (track == 0){
05257
05258 nu.stripTrkBegPerpFlag = GetPerpendicularStripFlag(det, nu.edgeRegionTrkVtx, nu.stripTrkBegIsu);
05259 nu.stripTrkEndPerpFlag = GetPerpendicularStripFlag(det, nu.edgeRegionTrkEnd, nu.stripTrkEndIsu);
05260 } else {
05261 MSG("NuReco",Msg::kError) << "Unrecognised track id " << track << endl;
05262 }
05263 }
05264
05265
05266
05267 Short_t NuReco::GetHorizontalVerticalStripNumber(Detector::Detector_t det, Short_t edgeRegion, Short_t stripTrku, Short_t stripTrkv) const
05268 {
05269 Short_t hoveNum = -999;
05270
05271
05272 if (det == Detector::kNear) return hoveNum;
05273
05274 if (edgeRegion == 0) hoveNum = 192 - stripTrku + stripTrkv;
05275 else if(edgeRegion == 2) hoveNum = 383 - stripTrku - stripTrkv;
05276 else if(edgeRegion == 4) hoveNum = 192 + stripTrku - stripTrkv;
05277 else if(edgeRegion == 6) hoveNum = 1 + stripTrku + stripTrkv;
05278
05279 return hoveNum;
05280 }
05281
05282
05283
05284 void NuReco::CalculateHorizontalVerticalStripNumber(NuEvent &nu, Int_t track) const
05285 {
05286 Detector::Detector_t det = static_cast<Detector::Detector_t>(nu.detector);
05287
05288 if (track == 1) {
05289 nu.stripHoveNumTrkVtx1 = GetHorizontalVerticalStripNumber(det, nu.edgeRegionTrkVtx1, nu.stripTrkBegu1, nu.stripTrkBegv1);
05290 nu.stripHoveNumTrkEnd1 = GetHorizontalVerticalStripNumber(det, nu.edgeRegionTrkVtx1, nu.stripTrkBegu1, nu.stripTrkBegv1);
05291 } else if (track == 2) {
05292 nu.stripHoveNumTrkVtx2 = GetHorizontalVerticalStripNumber(det, nu.edgeRegionTrkVtx2, nu.stripTrkBegu2, nu.stripTrkBegv2);
05293 nu.stripHoveNumTrkEnd2 = GetHorizontalVerticalStripNumber(det, nu.edgeRegionTrkVtx2, nu.stripTrkBegu2, nu.stripTrkBegv2);
05294 } else if (track == 3) {
05295 nu.stripHoveNumTrkVtx3 = GetHorizontalVerticalStripNumber(det, nu.edgeRegionTrkVtx3, nu.stripTrkBegu3, nu.stripTrkBegv3);
05296 nu.stripHoveNumTrkEnd3 = GetHorizontalVerticalStripNumber(det, nu.edgeRegionTrkVtx3, nu.stripTrkBegu3, nu.stripTrkBegv3);
05297 } else if (track == 0){
05298
05299 nu.stripHoveNumTrkVtx = GetHorizontalVerticalStripNumber(det, nu.edgeRegionTrkVtx, nu.stripTrkBegu, nu.stripTrkBegv);
05300 nu.stripHoveNumTrkEnd = GetHorizontalVerticalStripNumber(det, nu.edgeRegionTrkVtx, nu.stripTrkBegu, nu.stripTrkBegv);
05301 } else {
05302 MSG("NuReco",Msg::kError) << "Unrecognised track id " << track << endl;
05303 }
05304 }
05305
05306
05307
05308 Int_t NuReco::GetTrueTrackId(const NtpStRecord &ntp, const NtpSRTrack &trk) const
05309 {
05310
05311
05312 TClonesArray& thtrkTca=(*ntp.thtrk);
05313
05314 if (thtrkTca.GetEntriesFast() <= 0) return 0;
05315
05316
05317 const NtpTHTrack* thtrk= dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
05318
05319 if(!thtrk) {
05320 MSG("NuReco",Msg::kError) << "Truth track information available but track "
05321 << trk.index << " not found!" << endl;
05322 return 0;
05323 }
05324
05325
05326 NtpMCStdHep * stdhep =
05327 dynamic_cast<NtpMCStdHep *>(ntp.stdhep->At(thtrk->trkstdhep));
05328
05329 if(!stdhep) {
05330 MSG("NuReco",Msg::kError) << "Truth Helper, but no stdhep!\n";
05331 return 0;
05332 }
05333
05334 return stdhep->IdHEP;
05335 }
05336
05337
05338
05339 void NuReco::RecalculateDerivativeTrackGeometryVariables(NuEvent &nu) const
05340 {
05341
05342 CalculateEdgeRegion(nu, 0);
05343 CalculateEdgeRegion(nu, 1);
05344 CalculateEdgeRegion(nu, 2);
05345 CalculateEdgeRegion(nu, 3);
05346
05347
05348 CalculateParallelStripInset(nu, 0);
05349 CalculateParallelStripInset(nu, 1);
05350 CalculateParallelStripInset(nu, 2);
05351 CalculateParallelStripInset(nu, 3);
05352
05353
05354 CalculatePerpendicularStripFlag(nu, 0);
05355 CalculatePerpendicularStripFlag(nu, 1);
05356 CalculatePerpendicularStripFlag(nu, 2);
05357 CalculatePerpendicularStripFlag(nu, 3);
05358
05359
05360 CalculateHorizontalVerticalStripNumber(nu, 0);
05361 CalculateHorizontalVerticalStripNumber(nu, 1);
05362 CalculateHorizontalVerticalStripNumber(nu, 2);
05363 CalculateHorizontalVerticalStripNumber(nu, 3);
05364 }
