//*-- Author : M. Wisniowski //*-- Modified : 2005-9-13 //////////////////////////////////////////////////////////////////////// // // HHypPEpEmProjector // // HHypPEpEmProjector projects any PP data. At the moment output contains // //////////////////////////////////////////////////////////////////////// /** Ntuple content tescription sector_p - sector of first proton sector_e1 - sector of first lepton sector_e2 - sector of secound lepton system_p - system of first proton (shower-tofino->0, tof->1) system_e1 - system of first lepton system_e2 - system of secound lepton P_p - momentum [MeV] of first proton Th_p - theta angle [rad] of first proton Ph_p - 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_p piont) (first lepton) z_e1 - z vertex of first lepton track r_e2 - the smallest distance between track and z-axis (z_p piont) (secound lepton) z_e2 - z vertex of secound lepton track z_p - z vertex of proton track (first proton) r_p - the smallest distance between track and z-axis (z_p piont) (first proton) 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_miss - squerd missing mass of two protons p_miss - missing momentum of two protons m2_miss_pe1e2 - four particle squerd missing mass (p p e e) p_miss_pe1e2 - 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_p - Runge Kutta chi squerd of first 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_p - inner MDC chi squerd of first 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_p - opening angle between first proton and other closest lepton > 0 if closest lepton fitted < 0 if closest lepton not fitted (multipl. by -1) OpAngleClosestHadr_p - opening angle between first proton and other closest hadron > 0 if closest hadron fitted < 0 if closest hadron not fitted (multipl. by -1) **/ using namespace std; #include "hhypPEpEmProjector.h" #include "hgeomvector.h" #include "hgeomvertexfit.h" #define DEBUG 0 // well known constants should not be defined #define c 0.299792 #define D2R 0.0174532925199432955 #define R2D 57.2957795130823229 #define P_mass 938.272309999999998 #define E_mass 0.51099906000000002 ClassImp(HHypPEpEmProjector) HHypPEpEmProjector::HHypPEpEmProjector(Char_t *name_i, Option_t par[]) :HHypBaseAlgorithm(name_i,par) { simuflag = 0; } HHypPEpEmProjector::~HHypPEpEmProjector() { } Bool_t HHypPEpEmProjector::execute() { Int_t icomb=0; Float_t deltaTof= 100000000; Int_t BestComb=-1; 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++; } 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 r_e1=0,z_e1=0,r_e2=0,z_e2=0,r_p=0,z_p=0; Float_t q_e1=0, q_e2=0; Float_t RKchiq_p=-10, RKchiq_e1=-10, RKchiq_e2=-10; // RK chi2 for track1 & track2 Short_t sector_p=-1, system_p=-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_p=-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_p=-10, OpAngleClosestHadr_p=-10; Int_t geant_grandparentID_p = -10, geant_parentID_p = -10, geantID_p = -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_p=-10, geninfo2_p=-10; Float_t geninfo1_e1=-10, geninfo2_e1=-10; Float_t geninfo1_e2=-10, geninfo2_e2=-10; Bool_t IsCleanE1=0; Bool_t IsCleanE2=0; // Resetting the list and start looping over the combinations // Loop is only done over the VALID combinations mylist->CombIteratorReset(); Bool_t IsBestComb=0; Int_t ncomb=mylist->getNcomb(); icomb=0; while (mylist->CombIterator()) { if(mylist->getProbAlg(icomb)<=0) continue; if( icomb == BestComb ) IsBestComb = 1; else IsBestComb = 0; TLorentzVector proton(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_p(0,0,0,0); TLorentzVector geant_e1(0,0,0,0); TLorentzVector geant_e2(0,0,0,0); TVector3 pp_vertex; TVector3 pp_distance; TVector3 ee_vertex; TVector3 ee_distance; Float_t dist_ee=100; if (mylist->getIterStatus() == kTRUE) { HPidTrackCand *PidTrackCand[4]; const HPidHitData *PidData = NULL; const HPidTrackData *pTrack = NULL; //-------------------- simulation ---------------------------------------------- if (simuflag == 1 ) { HPidTrackCandSim *my_p = (HPidTrackCandSim *) CatTrackCandSim-> getObject(mylist->getIdxPidTrackCand(icomb, 0)); HPidTrackCandSim *my_e1 = (HPidTrackCandSim *) CatTrackCandSim-> getObject(mylist->getIdxPidTrackCand(icomb, 2)); HPidTrackCandSim *my_e2 = (HPidTrackCandSim *) CatTrackCandSim-> getObject(mylist->getIdxPidTrackCand(icomb, 3)); const HPidGeantTrackSet *p1GeantSet = my_p->getGeantTrackSet(); const HPidGeantTrackSet *e1GeantSet = my_e1->getGeantTrackSet(); const HPidGeantTrackSet *e2GeantSet = my_e2->getGeantTrackSet(); geninfo1_p = p1GeantSet->getGenInfo1(); geninfo2_p = p1GeantSet->getGenInfo2(); geantID_p = p1GeantSet->getGeantPID(); geant_parentID_p = p1GeantSet->getGeantParentID(); geant_grandparentID_p = p1GeantSet->getGeantGrandParentID(); TVector3 v1(p1GeantSet->getGeantMomX(), p1GeantSet->getGeantMomY(), p1GeantSet->getGeantMomZ()); geant_p.SetVectM(v1,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 v2(e1GeantSet->getGeantMomX(), e1GeantSet->getGeantMomY(), e1GeantSet->getGeantMomZ()); geant_e1.SetVectM(v2,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 v3(e2GeantSet->getGeantMomX(), e2GeantSet->getGeantMomY(), e2GeantSet->getGeantMomZ()); geant_e2.SetVectM(v3,P_mass); } // simuflag //-------------------- simulation end ------------------------------------------- HCategory *pidTrackCandCat = gHades->getCurrentEvent()->getCategory(catPidTrackCand); if (pidTrackCandCat != NULL ) { for(Int_t ipar=0; ipar<3; 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 Int_t system = PidData->iSystem; Bool_t isClean= PidTrackCand[ipar]->isFlagBit(HPidTrackCand::kIsUsed); 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); Float_t RKchiq = pTrack->getRKChiSquare(); switch(ipar) // Energy loss correction { case 0: { P = enLossCorr.getCorrMom(14 ,P,Th*TMath::RadToDeg()); } case 1: { if(q==1) P = enLossCorr.getCorrMom(2 ,P,Th*TMath::RadToDeg()); else P = enLossCorr.getCorrMom(3 ,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; } default: cerr<<"ERROR : HHypPEpEmProjector::execute() "<Fill(fpp); } else { Float_t fpp[]={ vertex.getX() , vertex.getY() ,vertex.getZ() ,vertex.getChi2() , sector_p, sector_e1, sector_e2, system_p, system_e1, system_e2, proton.P(), proton.Theta(), proton.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_p, z_p, ee_vertex.X(), ee_vertex.Y(), ee_vertex.Z(), dist_ee, miss.M2(), miss.P(), (miss - electron - positron).M2(), (miss - electron - positron).P(), op_angle_e1e2, (electron+positron).M2(), miss_cm.CosTheta(), RKchiq_p, RKchiq_e1, RKchiq_e2, ncomb, dsf, triggerBit, triggerDecision, date, time, InnerMDCchiq_p, InnerMDCchiq_e1, InnerMDCchiq_e2, OpAngleClosestLept_e1, OpAngleClosestHadr_e1, OpAngleClosestLept_e2, OpAngleClosestHadr_e2, OpAngleClosestLept_p, OpAngleClosestHadr_p, IsCleanE1, IsCleanE2, IsBestComb }; pp->Fill(fpp); } } else cerr << algoName << " got no TLorentzVector " << endl; icomb++; } return kTRUE; } Bool_t HHypPEpEmProjector::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_p:sector_e1:sector_e2:system_p:system_e1:system_e2:P_p:Th_p:Ph_p: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_p:z_p:xvert_ee:yvert_ee:zvert_ee:dist_ee:gP_p:gTh_p:gPh_p:gP_e1:gTh_e1:gPh_e1:gP_e2:gTh_e2:gPh_e2:m2_miss:p_miss:m2_miss_pe1e2:p_miss_pe1e2:opAngle_e1e2:m2_inv_e1e2:cosThCM_miss:RKchiq_p:RKchiq_e1:RKchiq_e2:ncomb:dsf:triggerBit:triggerDecision:geantID_p:geant_parentID_p:geant_grandparentID_p:geantID_e1:geant_parentID_e1:geant_grandparentID_e1:geantID_e2:geant_parentID_e2:geant_grandparentID_e2:geninfo1_p:geninfo1_e1:geninfo1_e2:geninfo2_p:geninfo2_e1:geninfo2_e2:InnerMDCchiq_p:InnerMDCchiq_e1:InnerMDCchiq_e2:OpAngleClosestLept_e1:OpAngleClosestHadr_e1:OpAngleClosestLept_e2:OpAngleClosestHadr_e2:OpAngleClosestLept_p:OpAngleClosestHadr_p", "sector_p:sector_e1:sector_e2:system_p:system_e1:system_e2:P_p:Th_p:Ph_p: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_p:z_p:xvert_ee:yvert_ee:zvert_ee:dist_ee:gP_p:gTh_p:gPh_p:gP_e1:gTh_e1:gPh_e1:gP_e2:gTh_e2:gPh_e2:m2_miss:p_miss:m2_miss_pe1e2:p_miss_pe1e2:opAngle_e1e2:m2_inv_e1e2:cosThCM_miss:RKchiq_p:RKchiq_e1:RKchiq_e2:ncomb:dsf:triggerBit:triggerDecision:geantID_p:geant_parentID_p:geant_grandparentID_p:geantID_e1:geant_parentID_e1:geant_grandparentID_e1:geantID_e2:geant_parentID_e2:geant_grandparentID_e2:geninfo1_p:geninfo1_e1:geninfo1_e2:geninfo2_p:geninfo2_e1:geninfo2_e2:InnerMDCchiq_p:InnerMDCchiq_e1:InnerMDCchiq_e2:OpAngleClosestLept_e1:OpAngleClosestHadr_e1:OpAngleClosestLept_e2:OpAngleClosestHadr_e2:OpAngleClosestLept_p:OpAngleClosestHadr_p"); } else { pp = new TNtuple(input + TString("_pp"), "eVert_x:eVert_y:eVert_Z:eVert_chi2:sector_p:sector_e1:sector_e2:system_p:system_e1:system_e2:P_p:Th_p:Ph_p: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_p:z_p:xvert_ee:yvert_ee:zvert_ee:dist_ee:m2_missp:p_missp:m2_miss_pee:p_miss_pee:m2_pee:p_pee:opAngle_ee:m2_inv_ee:cosThCM_miss:RKchiq_p:RKchiq_e1:RKchiq_e2:ncomb:dsf:triggerBit:triggerDecision:date:time:InnerMDCchiq_p:InnerMDCchiq_e1:InnerMDCchiq_e2:OpAngleClosestLept_e1:OpAngleClosestHadr_e1:OpAngleClosestLept_e2:OpAngleClosestHadr_e2:OpAngleClosestLept_p:OpAngleClosestHadr_p:IsCleanE1:IsCleanE2:IsBestComb", "eVert_x:eVert_y:eVert_Z:eVert_chi2:sector_p:sector_e1:sector_e2:system_p:system_e1:system_e2:P_p:Th_p:Ph_p: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_p:z_p:xvert_ee:yvert_ee:zvert_ee:dist_ee:m2_missp:p_missp:m2_miss_pee:p_miss_pee:m2_pee:p_pee:opAngle_ee:m2_inv_ee:cosThCM_miss:RKchiq_p:RKchiq_e1:RKchiq_e2:ncomb:dsf:triggerBit:triggerDecision:date:time:InnerMDCchiq_p:InnerMDCchiq_e1:InnerMDCchiq_e2:OpAngleClosestLept_e1:OpAngleClosestHadr_e1:OpAngleClosestLept_e2:OpAngleClosestHadr_e2:OpAngleClosestLept_p:OpAngleClosestHadr_p:IsCleanE1:IsCleanE2:IsBestComb"); } //f->Close(); return kTRUE; } Bool_t HHypPEpEmProjector::reinit() { return kTRUE; } Bool_t HHypPEpEmProjector::finalize() { pp->Write(); return kTRUE; } Bool_t HHypPEpEmProjector::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 HHypPEpEmProjector::SetParamFile(TString pFile) { paramFile=pFile; return kTRUE; } Float_t HHypPEpEmProjector::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; }