//*-- Author : M. Wisniowski //*-- Modified : 2005-9-13 //////////////////////////////////////////////////////////////////////// // // HHypPPEpEmGammaProjector // // HHypPPEpEmGammaProjector projects any PP data. At the moment output contains // //////////////////////////////////////////////////////////////////////// /** Ntuple content tescription sector_1 - sector of first proton sector_2 - sector of secound proton sector_e1 - sector of first lepton sector_e2 - sector of secound lepton system_p1 - system of first proton (shower-tofino->0, tof->1) system_p2 - system of secound proton system_e1 - system of first lepton system_e2 - system of secound lepton P_1 - momentum [MeV] of first proton Th_1 - theta angle [rad] of first proton Ph_1 - phi angle [rad] of secound proton P_2 - momentum [MeV] of secound proton Th_2 - theta angle [rad] of secound proton Ph_2 - phi angle [rad] of secound proton q_e1 - charge of first lepton P_e1 - momentum [MeV] of first lepton Th_e1 - theta angle [rad] of first lepton Ph_e1 - phi angle [rad] of first lepton q_e2 - charge of secound lepton P_e2 - momentum [MeV] of secound lepton Th_e2 - theta angle [rad] of secound lepton Ph_e2 - phi angle [rad] of secound lepton r_e1 - the smallest distance between track and z-axis (z_p1 piont) (first lepton) z_e1 - z vertex of first lepton track r_e2 - the smallest distance between track and z-axis (z_p1 piont) (secound lepton) z_e2 - z vertex of secound lepton track z_p1 - z vertex of proton track (first proton) r_p1 - the smallest distance between track and z-axis (z_p1 piont) (first proton) z_p2 - z vertex of proton track (secound proton) r_p2 - the smallest distance between track and z-axis (z_p1 piont) (secound proton) xvert_pp - x component of the vertex calculated from two protons yvert_pp - y component of the vertex calculated from two protons zvert_pp - z component of the vertex calculated from two protons dist_pp - closet distance between proton-proton tracks xvert_ee - x component of the vertex calculated from e-e pair yvert_ee - y component of the vertex calculated from e-e pair zvert_ee - z component of the vertex calculated from e-e pair dist_ee - closet distance between e-e tracks m2_miss2 - squerd missing mass of two protons p_miss2 - missing momentum of two protons dphi_pp - | Ph_p1 - Ph_p2 | -> for el. scattering m2_miss2e1e2 - four particle squerd missing mass (p p e e) p_miss2e1e2 - four particle squerd missing momentum (p p e e) opAngle_e1e2 - e-e pair openining angle [deg] m2_inv_e1e2 - invariant mass of e-e pair [MeV/c2] tanT1tanT2 - tan()*tan() -> for el. scattering cosThCM_miss - cos(Theta) of missing particle RKchiq_p1 - Runge Kutta chi squerd of first proton RKchiq_p2 - Runge Kutta chi squerd of secound proton RKchiq_e1 - Runge Kutta chi squerd of first lepton RKchiq_e2 - Runge Kutta chi squerd of secound lepton ncomb - number of track-track combination in the event dsf - down-scale factor triggerBit - tells what kind of trigger was used (trigger box setting) triggerDecision - secound level trigger decision date - date when file was taken time - time when file was taken InnerMDCchiq_p1 - inner MDC chi squerd of first proton InnerMDCchiq_p2 - inner MDC chi squerd of secound proton InnerMDCchiq_e1 - inner MDC chi squerd of first lepton InnerMDCchiq_e2 - inner MDC chi squerd of secound lepton OpAngleClosestLept_e1 - opening angle between first lepton and other closest lepton > 0 if closest lepton fitted < 0 if closest lepton not fitted (multipl. by -1) OpAngleClosestHadr_e1 - opening angle between first lepton and other closest hadron > 0 if closest hadron fitted < 0 if closest hadron not fitted (multipl. by -1) OpAngleClosestLept_e2 - opening angle between secound lepton and other closest lepton > 0 if closest lepton fitted < 0 if closest lepton not fitted OpAngleClosestHadr_e2 - opening angle between secound lepton and other closest hadron > 0 if closest hadron fitted < 0 if closest hadron not fitted OpAngleClosestLept_p1 - opening angle between first proton and other closest lepton > 0 if closest lepton fitted < 0 if closest lepton not fitted (multipl. by -1) OpAngleClosestHadr_p1 - opening angle between first proton and other closest hadron > 0 if closest hadron fitted < 0 if closest hadron not fitted (multipl. by -1) OpAngleClosestLept_p2 - opening angle between secound proton and other closest lepton > 0 if closest lepton fitted < 0 if closest lepton not fitted OpAngleClosestHadr_p2 - opening angle between secound proton and other closest lepton > 0 if closest hadron fitted < 0 if closest hadron not fitted **/ using namespace std; #include "hhypPPEpEmGammaProjector.h" #include "hgeomvector.h" #include "hgeomvertexfit.h" #define DEBUG 0 #define c 0.299792 #define D2R 0.0174532925199432955 #define R2D 57.2957795130823229 #define P_mass 938.272309999999998 #define E_mass 0.51099906000000002 #define Pip_mass 139.57018 ClassImp(HHypPPEpEmGammaProjector) HHypPPEpEmGammaProjector::HHypPPEpEmGammaProjector(Char_t *name_i, Option_t par[]) :HHypBaseAlgorithm(name_i,par) { simuflag = 0; } HHypPPEpEmGammaProjector::~HHypPPEpEmGammaProjector() { } Bool_t HHypPPEpEmGammaProjector::execute() { Short_t triggerBit = gHades->getCurrentEvent()->getHeader()->getTBit(); HEventHeader *evHeader = gHades->getCurrentEvent()->getHeader(); HVertex vertex = gHades->getCurrentEvent()->getHeader()->getVertex(); UInt_t date = evHeader->getDate(); UInt_t time = evHeader->getTime(); UInt_t dsf = evHeader->getDownscaling(); UInt_t triggerDecision = evHeader->getTriggerDecision(); if (!beam) { cerr << algoName << " needs beam particle! " << endl; return kFALSE; } Float_t P_1=0 , P_2=0, Th_1=0, Th_2=0, Ph_1=0, Ph_2=0; Float_t r_e1=0,z_e1=0,r_e2=0,z_e2=0,r_p1=0,z_p1=0,r_p2=0,z_p2=0; Float_t q_e1=0, q_e2=0; Float_t RKchiq_p1=-10, RKchiq_p2=-10, RKchiq_e1=-10, RKchiq_e2=-10; // RK chi2 for track1 & track2 Short_t sector_1=-1, sector_2=-1, system_1=-1, system_2=-1; // number of sector for first and secound particle Short_t sector_e1=-1, sector_e2=-1, system_e1=-1, system_e2=-1; // number of sector for e+ e- Float_t InnerMDCchiq_p1=-10, InnerMDCchiq_p2=-10; Float_t InnerMDCchiq_e1=-10, InnerMDCchiq_e2=-10; Float_t OpAngleClosestLept_e1=-10, OpAngleClosestHadr_e1=-10; Float_t OpAngleClosestLept_e2=-10, OpAngleClosestHadr_e2=-10; Float_t OpAngleClosestLept_p1=-10, OpAngleClosestHadr_p1=-10; Float_t OpAngleClosestLept_p2=-10, OpAngleClosestHadr_p2=-10; Int_t geant_grandparentID_p1 = -10, geant_parentID_p1 = -10, geantID_p1 = -100; Int_t geant_grandparentID_p2 = -10, geant_parentID_p2 = -10, geantID_p2 = -100; Int_t geant_grandparentID_e1 = -10, geant_parentID_e1 = -10, geantID_e1 = -100; Int_t geant_grandparentID_e2 = -10, geant_parentID_e2 = -10, geantID_e2 = -100; Float_t geninfo1_p1=-10, geninfo2_p1=-10; Float_t geninfo1_p2=-10, geninfo2_p2=-10; Float_t geninfo1_e1=-10, geninfo2_e1=-10; Float_t geninfo1_e2=-10, geninfo2_e2=-10; TVector3 v1(0,0,0),v2(0,0,0); TLorentzVector miss2(0,0,0,0), R1(0,0,0,0), R2(0,0,0,0), part_1(0,0,0,0), part_2(0,0,0,0); Float_t dEdx_p1 = 0; Float_t dEdx_p2 = 0; Float_t dEdx_e1 = 0; Float_t dEdx_e2 = 0; Int_t BestComb = 0; Int_t IsBestComb = 0; Int_t icomb=0; Float_t deltaTof= 100000000; mylist->CombIteratorReset(); while (mylist->CombIterator()) { Float_t deltaTof_tmp; mylist->getUserValue(DELTATOF_CHI2, deltaTof_tmp); if( deltaTof_tmp < deltaTof ) { deltaTof = deltaTof_tmp; BestComb = icomb; } icomb++; } // Resetting the list and start looping over the combinations // Loop is only done over the VALID combinations mylist->CombIteratorReset(); Int_t ncomb=mylist->getNcomb(); icomb=0; while (mylist->CombIterator()) { if( icomb == BestComb ) IsBestComb = 1; else IsBestComb = 0; if(mylist->getProbAlg(icomb)<=0) {icomb++; continue;} TLorentzVector proton1(0,0,0,0); TLorentzVector proton2(0,0,0,0); TLorentzVector pip1(0,0,0,0); TLorentzVector pip2(0,0,0,0); TLorentzVector positron(0,0,0,0); // in case of ppe-e- / ppe+e+ plays role of e- TLorentzVector electron(0,0,0,0); TLorentzVector positronRICH(0,0,0,0); // keep Theta and Phi taken from RICH TLorentzVector electronRICH(0,0,0,0); //TLorentzVector gamma(0,0,0,0); TLorentzVector geant_p1(0,0,0,0); TLorentzVector geant_p2(0,0,0,0); TLorentzVector geant_e1(0,0,0,0); TLorentzVector geant_e2(0,0,0,0); TVector3 pp_vertex; TVector3 pp_distance; Float_t dist_pp=100; TVector3 ee_vertex; TVector3 ee_distance; Float_t dist_ee=100; Int_t clean_e1=0, clean_e2=0, clean_p1=0, clean_p2=0; if (mylist->getIterStatus() == kTRUE) { HPidTrackCand *PidTrackCand[4]; const HPidHitData *PidData = NULL; const HPidTrackData *pTrack = NULL; //-------------------- simulation ---------------------------------------------- if (simuflag == 1 ) { HPidTrackCandSim *my_p1 = (HPidTrackCandSim *) CatTrackCandSim-> getObject(mylist->getIdxPidTrackCand(icomb, 0)); HPidTrackCandSim *my_p2 = (HPidTrackCandSim *) CatTrackCandSim-> getObject(mylist->getIdxPidTrackCand(icomb, 1)); HPidTrackCandSim *my_e1 = (HPidTrackCandSim *) CatTrackCandSim-> getObject(mylist->getIdxPidTrackCand(icomb, 2)); HPidTrackCandSim *my_e2 = (HPidTrackCandSim *) CatTrackCandSim-> getObject(mylist->getIdxPidTrackCand(icomb, 3)); const HPidGeantTrackSet *p1GeantSet = my_p1->getGeantTrackSet(); const HPidGeantTrackSet *p2GeantSet = my_p2->getGeantTrackSet(); const HPidGeantTrackSet *e1GeantSet = my_e1->getGeantTrackSet(); const HPidGeantTrackSet *e2GeantSet = my_e2->getGeantTrackSet(); geninfo1_p1 = p1GeantSet->getGenInfo1(); geninfo2_p1 = p1GeantSet->getGenInfo2(); geantID_p1 = p1GeantSet->getGeantPID(); geant_parentID_p1 = p1GeantSet->getGeantParentID(); geant_grandparentID_p1 = p1GeantSet->getGeantGrandParentID(); TVector3 v1(p1GeantSet->getGeantMomX(), p1GeantSet->getGeantMomY(), p1GeantSet->getGeantMomZ()); geant_p1.SetVectM(v1,P_mass); geninfo1_p2 = p2GeantSet->getGenInfo1(); geninfo2_p2 = p2GeantSet->getGenInfo2(); geantID_p2 = p2GeantSet->getGeantPID(); geant_parentID_p2 = p2GeantSet->getGeantParentID(); geant_grandparentID_p2 = p2GeantSet->getGeantGrandParentID(); TVector3 v2(p2GeantSet->getGeantMomX(), p2GeantSet->getGeantMomY(), p2GeantSet->getGeantMomZ()); geant_p2.SetVectM(v2,P_mass); geninfo1_e1 = e1GeantSet->getGenInfo1(); geninfo2_e1 = e1GeantSet->getGenInfo2(); geantID_e1 = e1GeantSet->getGeantPID(); geant_parentID_e1 = e1GeantSet->getGeantParentID(); geant_grandparentID_e1 = e1GeantSet->getGeantGrandParentID(); TVector3 v3(e1GeantSet->getGeantMomX(), e1GeantSet->getGeantMomY(), e1GeantSet->getGeantMomZ()); geant_e1.SetVectM(v3,P_mass); geninfo1_e2 = e2GeantSet->getGenInfo1(); geninfo2_e2 = e2GeantSet->getGenInfo2(); geantID_e2 = e2GeantSet->getGeantPID(); geant_parentID_e2 = e2GeantSet->getGeantParentID(); geant_grandparentID_e2 = e2GeantSet->getGeantGrandParentID(); TVector3 v4(e2GeantSet->getGeantMomX(), e2GeantSet->getGeantMomY(), e2GeantSet->getGeantMomZ()); geant_e2.SetVectM(v4,P_mass); } // simuflag //-------------------- simulation end ------------------------------------------- HCategory *pidTrackCandCat = gHades->getCurrentEvent()->getCategory(catPidTrackCand); if (pidTrackCandCat != NULL ) { for(Int_t ipar=0; ipar<4; ipar++) { PidTrackCand[ipar] = (HPidTrackCand *) pidTrackCandCat->getObject(mylist->getIdxPidTrackCand(icomb, ipar)); if (PidTrackCand[ipar] != NULL) { PidData = PidTrackCand[ipar]->getHitData(); pTrack = PidTrackCand[ipar]->getTrackData(); Float_t closestLepton = pTrack->getAngleWithClosestLeptonCandidate() *(int(pTrack->closestLeptonCandidateIsFitted())-0.5)*2; // Fitted tracks positive , not fitted negative value of opening angle Float_t closestHadron = pTrack->getAngleWithClosestHadronCandidate() *(int(pTrack->closestHadronCandidateIsFitted())-0.5)*2; // Fitted tracks positive , not fitted negative value of opening angle Float_t dEdx = PidData -> getInnerMdcdEdx(); Int_t system = PidData->iSystem; Float_t InnerMDCchiq = PidData->fInnerMdcChiSquare; Int_t sector = PidData->getSector();; Int_t q = pTrack->getPolarity(4); Float_t P = pTrack->fMomenta[4]; Float_t Th = D2R*pTrack->getRKTheta(); Float_t Ph = D2R*pTrack->getRKPhi(); Float_t Th_rich= D2R*PidData->getRichTheta(); Float_t Ph_rich= D2R*PidData->getRichPhi(); Float_t R = pTrack->getTrackR(4); Float_t Z = pTrack->getTrackZ(4); TVector3 p; TVector3 p_rich; Float_t RKchiq=0; switch(ipar) // Energy loss correction { case 0: { P = enLossCorr.getCorrMom(14 ,P,Th*TMath::RadToDeg()); break; } case 1: { P = enLossCorr.getCorrMom(14 ,P,Th*TMath::RadToDeg()); break; } case 2: { if(q==1) P = enLossCorr.getCorrMom(2 ,P,Th*TMath::RadToDeg()); else P = enLossCorr.getCorrMom(3 ,P,Th*TMath::RadToDeg()); break; } case 3: { if(q==1) P = enLossCorr.getCorrMom(2 ,P,Th*TMath::RadToDeg()); else P = enLossCorr.getCorrMom(3 ,P,Th*TMath::RadToDeg()); break; } default: cerr<<"ERROR : HHypPPEpEmGammaProjector::execute() "<isFlagBit(HPidTrackCand::kIsUsed); break; } case 1: { P_2=P; Th_2=Th; Ph_2=Ph; proton2.SetVectM(p,P_mass); pip2.SetVectM(p,Pip_mass); system_2=system; sector_2=sector; RKchiq_p2=RKchiq; r_p2=R; z_p2=Z; InnerMDCchiq_p2=InnerMDCchiq; v2.SetXYZ(P*sin(Th)*cos(Ph),P*sin(Th)*sin(Ph),P*cos(Th)); dEdx_p2=dEdx; OpAngleClosestLept_p2=closestLepton; OpAngleClosestHadr_p2=closestHadron; clean_p2 = PidTrackCand[ipar]->isFlagBit(HPidTrackCand::kIsUsed); break; } case 2: { positron.SetVectM(p,E_mass); positronRICH.SetVectM(p_rich,E_mass); system_e1=system; sector_e1=sector; dEdx_e1=dEdx; RKchiq_e1=RKchiq; q_e1=q; r_e1=R; z_e1=Z; InnerMDCchiq_e1=InnerMDCchiq; OpAngleClosestLept_e1=closestLepton; OpAngleClosestHadr_e1=closestHadron; clean_e1 = PidTrackCand[ipar]->isFlagBit(HPidTrackCand::kIsUsed); break; } case 3: { electron.SetVectM(p,E_mass); electronRICH.SetVectM(p_rich,E_mass); system_e2=system; sector_e2=sector; dEdx_e2=dEdx; RKchiq_e2=RKchiq; q_e2=q; r_e2=R; z_e2=Z; InnerMDCchiq_e2=InnerMDCchiq; OpAngleClosestLept_e2=closestLepton; OpAngleClosestHadr_e2=closestHadron; clean_e2 = PidTrackCand[ipar]->isFlagBit(HPidTrackCand::kIsUsed); break; } default: cerr<<"ERROR : HHypPPEpEmGammaProjector::execute() "<Fill(fpp); } else { Float_t fpp[]={ vertex.getX() , vertex.getY() ,vertex.getZ() ,vertex.getChi2() , sector_1, sector_2, sector_e1, sector_e2, system_1, system_2, system_e1, system_e2, proton1.P(), proton1.Theta(), proton1.Phi(), proton2.P(), proton2.Theta(), proton2.Phi(), q_e1, positron.P(), positron.Theta(), positron.Phi(), positronRICH.Theta(), positronRICH.Phi(), q_e2, electron.P(), electron.Theta(), electron.Phi(), electronRICH.Theta(), electronRICH.Phi(), r_e1, z_e1, r_e2, z_e2, r_p1, z_p1, r_p2, z_p2, pp_vertex.X(), pp_vertex.Y(), pp_vertex.Z(), dist_pp, ee_vertex.X(), ee_vertex.Y(), ee_vertex.Z(), dist_ee, miss2.M2(), miss2.P(), fabs(proton1.Phi() - proton2.Phi()), (miss2 - electron - positron).M2(), (miss2 - electron - positron).P(), op_angle_e1e2, (electron+positron).M2(), tanTh1*tanTh2, miss2_cm.CosTheta(), RKchiq_p1, RKchiq_p2, RKchiq_e1, RKchiq_e2, ncomb, dsf, triggerBit, triggerDecision, date, time, InnerMDCchiq_p1, InnerMDCchiq_p2, InnerMDCchiq_e1, InnerMDCchiq_e2, OpAngleClosestLept_e1, OpAngleClosestHadr_e1, OpAngleClosestLept_e2, OpAngleClosestHadr_e2, OpAngleClosestLept_p1, OpAngleClosestHadr_p1, OpAngleClosestLept_p2, OpAngleClosestHadr_p2, dEdx_p1, dEdx_p2, dEdx_e1, dEdx_e2, IsBestComb, clean_e1, clean_e2, clean_p1, clean_p2 }; pp->Fill(fpp); } } else cerr << algoName << " got no TLorentzVector " << endl; icomb++; } return kTRUE; } Bool_t HHypPPEpEmGammaProjector::init() { enLossCorr.setDefaultPar("jan04"); simCat = gHades->getCurrentEvent()->getCategory(catGeantKine); if (!simCat) { simuflag = 0; } else { simuflag = 1; //cout << "Projector uses SIMULATION" << endl; CatTrackCandSim = NULL; // Category if ((CatTrackCandSim = gHades->getCurrentEvent()->getCategory(catPidTrackCand)) == NULL) { Error("init", "Cannot get catPidTrackCandSim cat"); return kFALSE; } } // need to get name from channel TString input(channel->Get(initList)); TFile *f=GetHFile(); f->cd(); if(simuflag==1) { pp = new TNtuple(input + TString("_pp"), "sector_1:sector_2:sector_e1:sector_e2:system_1:system_2:system_e1:system_e2:P_1:Th_1:Ph_1:P_2:Th_2:Ph_2:q_e1:P_e1:Th_e1:Ph_e1:q_e2:P_e2:Th_e2:Ph_e2:r_e1:z_e1:r_e2:z_e2:r_p1:z_p1:r_p2:z_p2:xvert_pp:yvert_pp:zvert_pp:dist_pp:xvert_ee:yvert_ee:zvert_ee:dist_ee:gP_1:gTh_1:gPh_1:gP_2:gTh_2:gPh_2:gP_e1:gTh_e1:gPh_e1:gP_e2:gTh_e2:gPh_e2:m2_miss2:p_miss2:dphi_pp:m2_miss2e1e2:p_miss2e1e2:opAngle_e1e2:m2_inv_e1e2:tanT1tanT2:cosThCM_miss:RKchiq_p1:RKchiq_p2:RKchiq_e1:RKchiq_e2:ncomb:dsf:triggerBit:triggerDecision:geantID_p1:geant_parentID_p1:geant_grandparentID_p1:geant_ID_p2:geant_parentID_p2:geant_grandparentID_p2:geantID_e1:geant_parentID_e1:geant_grandparentID_e1:geantID_e2:geant_parentID_e2:geant_grandparentID_e2:geninfo1_p1:geninfo1_p2:geninfo1_e1:geninfo1_e2:geninfo2_p1:geninfo2_p2:geninfo2_e1:geninfo2_e2:InnerMDCchiq_p1:InnerMDCchiq_p2:InnerMDCchiq_e1:InnerMDCchiq_e2:OpAngleClosestLept_e1:OpAngleClosestHadr_e1:OpAngleClosestLept_e2:OpAngleClosestHadr_e2:OpAngleClosestLept_p1:OpAngleClosestHadr_p1:OpAngleClosestLept_p2:OpAngleClosestHadr_p2", "sector_1:sector_2:sector_e1:sector_e2:system_1:system_2:system_e1:system_e2:P_1:Th_1:Ph_1:P_2:Th_2:Ph_2:q_e1:P_e1:Th_e1:Ph_e1:q_e2:P_e2:Th_e2:Ph_e2:r_e1:z_e1:r_e2:z_e2:r_p1:z_p1:r_p2:z_p2:xvert_pp:yvert_pp:zvert_pp:dist_pp:xvert_ee:yvert_ee:zvert_ee:dist_ee:gP_1:gTh_1:gPh_1:gP_2:gTh_2:gPh_2:gP_e1:gTh_e1:gPh_e1:gP_e2:gTh_e2:gPh_e2:m2_miss2:p_miss2:dphi_pp:m2_miss2e1e2:p_miss2e1e2:opAngle_e1e2:m2_inv_e1e2:tanT1tanT2:cosThCM_miss:RKchiq_p1:RKchiq_p2:RKchiq_e1:RKchiq_e2:ncomb:dsf:triggerBit:triggerDecision:geantID_p1:geant_parentID_p1:geant_grandparentID_p1:geantID_p2:geant_parentID_p2:geant_grandparentID_p2:geantID_e1:geant_parentID_e1:geant_grandparentID_e1:geantID_e2:geant_parentID_e2:geant_grandparentID_e2:geninfo1_p1:geninfo1_p2:geninfo1_e1:geninfo1_e2:geninfo2_p1:geninfo2_p2:geninfo2_e1:geninfo2_e2:InnerMDCchiq_p1:InnerMDCchiq_p2:InnerMDCchiq_e1:InnerMDCchiq_e2:OpAngleClosestLept_e1:OpAngleClosestHadr_e1:OpAngleClosestLept_e2:OpAngleClosestHadr_e2:OpAngleClosestLept_p1:OpAngleClosestHadr_p1:OpAngleClosestLept_p2:OpAngleClosestHadr_p2"); } else { pp = new TNtuple(input + TString("_pp"), "eVert_x:eVert_y:eVert_Z:eVert_chi2:sector_1:sector_2:sector_e1:sector_e2:system_1:system_2:system_e1:system_e2:P_1:Th_1:Ph_1:P_2:Th_2:Ph_2:q_e1:P_e1:Th_e1:Ph_e1:Th_rich_e1:Ph_rich_e1:q_e2:P_e2:Th_e2:Ph_e2:Th_rich_e2:Ph_rich_e2:r_e1:z_e1:r_e2:z_e2:r_p1:z_p1:r_p2:z_p2:xvert_pp:yvert_pp:zvert_pp:dist_pp:xvert_ee:yvert_ee:zvert_ee:dist_ee:m2_miss2:p_miss2:dphi_pp:m2_miss2e1e2:p_miss2e1e2:opAngle_e1e2:m2_inv_e1e2:tanT1tanT2:cosThCM_miss:RKchiq_p1:RKchiq_p2:RKchiq_e1:RKchiq_e2:ncomb:dsf:triggerBit:triggerDecision:date:time:InnerMDCchiq_p1:InnerMDCchiq_p2:InnerMDCchiq_e1:InnerMDCchiq_e2:OpAngleClosestLept_e1:OpAngleClosestHadr_e1:OpAngleClosestLept_e2:OpAngleClosestHadr_e2:OpAngleClosestLept_p1:OpAngleClosestHadr_p1:OpAngleClosestLept_p2:OpAngleClosestHadr_p2:dEdx_p1:dEdx_p2:dEdx_e1:dEdx_e2:BestComb:clean_e1:clean_e2:clean_p1:clean_p2", "eVert_x:eVert_y:eVert_Z:eVert_chi2:sector_1:sector_2:sector_e1:sector_e2:system_1:system_2:system_e1:system_e2:P_1:Th_1:Ph_1:P_2:Th_2:Ph_2:q_e1:P_e1:Th_e1:Ph_e1:Th_rich_e1:Ph_rich_e1:q_e2:P_e2:Th_e2:Ph_e2:Th_rich_e2:Ph_rich_e2:r_e1:z_e1:r_e2:z_e2:r_p1:z_p1:r_p2:z_p2:xvert_pp:yvert_pp:zvert_pp:dist_pp:xvert_ee:yvert_ee:zvert_ee:dist_ee:m2_miss2:p_miss2:dphi_pp:m2_miss2e1e2:p_miss2e1e2:opAngle_e1e2:m2_inv_e1e2:tanT1tanT2:cosThCM_miss:RKchiq_p1:RKchiq_p2:RKchiq_e1:RKchiq_e2:ncomb:dsf:triggerBit:triggerDecision:date:time:InnerMDCchiq_p1:InnerMDCchiq_p2:InnerMDCchiq_e1:InnerMDCchiq_e2:OpAngleClosestLept_e1:OpAngleClosestHadr_e1:OpAngleClosestLept_e2:OpAngleClosestHadr_e2:OpAngleClosestLept_p1:OpAngleClosestHadr_p1:OpAngleClosestLept_p2:OpAngleClosestHadr_p2:dEdx_p1:dEdx_p2:dEdx_e1:dEdx_e2:BestComb:clean_e1:clean_e2:clean_p1:clean_p2"); } //f->Close(); return kTRUE; } Bool_t HHypPPEpEmGammaProjector::reinit() { return kTRUE; } Bool_t HHypPPEpEmGammaProjector::finalize() { pp->Write(); return kTRUE; } Bool_t HHypPPEpEmGammaProjector::IsOpposit(Short_t sec1, Short_t sec2) { if(sec1==0 && sec2==3) return 1; else if(sec1==1 && sec2==4) return 1; else if(sec1==2 && sec2==5) return 1; else if(sec1==3 && sec2==0) return 1; else if(sec1==4 && sec2==1) return 1; else if(sec1==5 && sec2==2) return 1; else return 0; } Bool_t HHypPPEpEmGammaProjector::SetParamFile(TString pFile) { paramFile=pFile; return kTRUE; } Float_t HHypPPEpEmGammaProjector::calcVertex(HPidTrackCand *p1, HPidTrackCand *p2, TVector3 *vertex, TVector3 *distance) { // calcVertex should return // 1. the vertex of two tracks (no weights included, so it returns the // center of closest approach vector) // 2. a vector with direction and magnitude of the distance // (using stefanos algebra to calculate the magnitude, root cross product // to give the direction) HGeomVector hoff[2]; HGeomVector hdir[2]; HGeomVector hvertex; HGeomVertexFit hfitter; TVector3 dir[2]; Float_t dist; Float_t det1, det2; // extract coordinates from p1, p2, fill them into HGeomVector // to use Manuels fitter hoff[0].setXYZ( p1->getTrackData()->getTrackR(4)*TMath::Cos(p1->getTrackData()->getRKPhi()*D2R + TMath::PiOver2()), p1->getTrackData()->getTrackR(4)*TMath::Sin(p1->getTrackData()->getRKPhi()*D2R + TMath::PiOver2()), p1->getTrackData()->getTrackZ(4)); hoff[1].setXYZ( p2->getTrackData()->getTrackR(4)*TMath::Cos(p1->getTrackData()->getRKPhi()*D2R + TMath::PiOver2()), p1->getTrackData()->getTrackR(4)*TMath::Sin(p1->getTrackData()->getRKPhi()*D2R + TMath::PiOver2()), p2->getTrackData()->getTrackZ(4)); dir[0].SetMagThetaPhi(p1->getTrackData()->getMomenta(4),p1->getTrackData()->getRKTheta()*D2R, p1->getTrackData()->getRKPhi()*D2R);// = p1->Vect(); dir[1].SetMagThetaPhi(p2->getTrackData()->getMomenta(4),p2->getTrackData()->getRKTheta()*D2R, p2->getTrackData()->getRKPhi()*D2R);// = p2->Vect(); hdir[0].setXYZ(dir[0].X(),dir[0].Y(),dir[0].Z()); hdir[1].setXYZ(dir[1].X(),dir[1].Y(),dir[1].Z()); hfitter.reset(); for (Int_t i = 0; i < 2; i++) { hfitter.addLine(hoff[i],hdir[i]); } hfitter.getVertex(hvertex); vertex->SetXYZ(hvertex.getX(),hvertex.getY(),hvertex.getZ()); // Function to calculate the distance between two lines in the space // c.f. Stefano det1 = ( (hoff[0].getX()-hoff[1].getX()) * (hdir[0].getY()*hdir[1].getZ()-hdir[0].getZ()*hdir[1].getY()) - (hoff[0].getY()-hoff[1].getY()) * (hdir[0].getX()*hdir[1].getZ()-hdir[0].getZ()*hdir[1].getX()) + (hoff[0].getZ()-hoff[1].getZ()) * (hdir[0].getX()*hdir[1].getY()-hdir[0].getY()*hdir[1].getX()) ); det2 = TMath::Sqrt( (hdir[0].getX()*hdir[1].getY() - hdir[0].getY()*hdir[1].getX()) * (hdir[0].getX()*hdir[1].getY() - hdir[0].getY()*hdir[1].getX()) + (hdir[0].getX()*hdir[1].getZ() - hdir[0].getZ()*hdir[1].getX()) * (hdir[0].getX()*hdir[1].getZ() - hdir[0].getZ()*hdir[1].getX()) + (hdir[0].getY()*hdir[1].getZ() - hdir[0].getZ()*hdir[1].getY()) * (hdir[0].getY()*hdir[1].getZ() - hdir[0].getZ()*hdir[1].getY()) ); // Create a distance vector and scale it with dist *distance = dir[0].Cross(dir[1]); if (det2==0) { dist = -1.; } else { dist = TMath::Abs(det1/det2); distance->SetMag(dist); } return dist; }