// $Id: hpidsingleleptoneff.cc,v 1.17 2009-10-15 14:42:36 halo Exp $ // // Author: Laura Fabbietti // Modified by Laura Fabbietti: 05/06/05 21:47:04 // Modified by W. Przygoda / T. Wojcik (Cracow) in March 2006 // Modified by Thomas.Eberl@ph.tum.de, last modified : 2006-11-08 10:42:49 // //_HADES_CLASS_DESCRIPTION //////////////////////////////////////////////////////////////////////////////// // // HPidSingleLeptonEff // // This task is used to generate ntuples of embedded leptons // with white distributions in 1/p, theta, phi. // There are ntuples with such leptons in the acceptance and // additionally identified and accepted leptons. // These ntuples are later divided by each other (external code) // to get (p,theta,phi) distributions of reconstruction efficiency // // Strategy: // // 1. Start from HGeantKine. Select primary leptons. Do not select // tracks which have a partner closer than 9.5 deg in the same sector. // ( "single tracks") // 2. Try to find the GEANT track in created Pairs. Definition for the // single track to be accepted is: // a. most common track == GEANT track // b. is primary track // c. has reconstructed lepton PID and correlation with a RICH ring // d. Geant PID == reconstructed PID // e. partner of the track can be real track, hadron track or secondary lepton. // Pairs from embedded lepton are rejected // f. PairCut flag != 1 (not direct cut) // g. common coorelation flag >= 76 (fully correlated) // 3. If the searched track is not found in a pair, it is // searched in HPidParticleSim (reconstructed as single track) // Same criteria as for the single track accepted in the pair search. // // ALL ANGLES ARE IN DEGREES: 0 < phi < 360 and 0 < theta < 90 as it is usual // all over the HADES analysis chain. // // ntuple description: // // members of the ACCEPTANCE and YIELD ntuples // p:theta:phi:pt calculated from HGeantKine // // EFFICIENCY ntuple: // // pRecon:thetaRecon:phiRecon: // reconstructed track properties // pGeant:thetaGeant:phiGeant:ptGeant: // GEANT track properties // // comDet: // common correlation flag (repaired for hidden leptons) // isCutNb: // cut nb of recursive cut for pairs // // ( pairs with isCutNb == 1 ==> direct cut are rejected), // // -1 for singles // angleHadronFitted: // environment of the reconstructed track // angleHadronNonFitted: // .. // angleLeptonFitted: // .. // angleLeptonNonFitted: // .. // isRing: // ring correlated with RK // innerChi2:outerChi2: // chi2 of inner and outer MDCs // rkChi2: // RK chi2 // metaMatchQality: // // metaMatchRK: // matching RK with META // vertX:vertY:vertZ: // event vertex // // isHiddenLepton: // 1 if a lepton track has been hidden by secondary in SHOWER // beta: // reconstructed beta of the track (start time corrected) // pairVertx:pairVerty:pairVertz: // vertex of the pair // opang: // reconstructed opening angle of the pair // pReconPartner // momentum of the partner in the pair // pidRecPartner // PID of the partner in the pair // //////////////////////////////////////////////////////////////////////////////// #include "hpidsingleleptoneff.h" #include "hades.h" #include "hcategory.h" #include "hevent.h" #include "hgeantkine.h" #include "hgeantmdc.h" #include "hgeantshower.h" #include "hgeanttof.h" #include "hiterator.h" #include "hmdcseg.h" #include "hpairdata.h" #include "hpairfl.h" #include "hpairgeantdata.h" #include "hpairsim.h" #include "hpidfl.h" #include "hpidparticlesim.h" #include "hreconstructor.h" #include "pairsdef.h" #include "piddef.h" #include "TFile.h" #include "TMath.h" #include "TNtuple.h" #include "TString.h" #include "TVector3.h" #include using namespace std; ClassImp(HPidSingleLeptonEff) // ----------------------------------------------------------------------------- HPidSingleLeptonEff::HPidSingleLeptonEff(const Char_t * file) : HReconstructor("PidPartFiller","Filler of HPidParticle category") { out = new TFile(file,"recreate"); bookNTuple(); maxNParticle = 100; } // ----------------------------------------------------------------------------- HPidSingleLeptonEff::~HPidSingleLeptonEff() { if( NULL != out ) { delete out; out = NULL; } if ( NULL != pitPartCat ) { delete pitPartCat; pitPartCat = NULL; } if ( NULL != pitPairCat ) { delete pitPairCat; pitPairCat = NULL; } if ( NULL != pitGeant ) { delete pitGeant; pitGeant = NULL; } if ( NULL != pitGeant1 ) { delete pitGeant1; pitGeant1 = NULL; } } // ----------------------------------------------------------------------------- HPidSingleLeptonEff::HPidSingleLeptonEff(const Text_t name[],const Text_t title[], const Char_t * file) : HReconstructor(name,title) { out = new TFile(file,"recreate"); bookNTuple(); maxNParticle = 100; } // ----------------------------------------------------------------------------- void HPidSingleLeptonEff::bookNTuple() { simYieldEle = new TNtuple("simYieldEle","simYieldEle","p:theta:phi:pt"); simYieldPos = new TNtuple("simYieldPos","simYieldPos","p:theta:phi:pt"); TString sVarList("pRecon:pGeant:thetaRecon:thetaGeant:" "phiRecon:phiGeant:ptGeant:isCutNb:" "commonTrack:bitFieldPassedCuts:" "angleHadronFitted:angleHadronNonFitted:" "angleLeptonFitted:angleLeptonNonFitted:" "isRing:innerChi2:rkChi2:metaMatchQality:" "metaMatchRK:vertX:vertY:vertZ:" "outerChi2:isHiddenLepton:beta:" "pairVertx:pairVerty:pairVertz:" "opang:pReconPartner:pidRecPartner" ); effEle = new TNtuple("effEle", "effEle", sVarList); effPos = new TNtuple("effPos", "effPos", sVarList); acceEle = new TNtuple("acceEle", "acceEle", "p:theta:phi:pt"); accePos = new TNtuple("accePos", "accePos", "p:theta:phi:pt"); } // ----------------------------------------------------------------------------- Bool_t HPidSingleLeptonEff::init(void) { if( NULL == (pPartCat = gHades->getCurrentEvent()->getCategory(catPidPart))) { Error("init", "No category catPidPart"); return kFALSE; } pitPartCat = (HIterator *)pPartCat->MakeIterator(); if( NULL == (pPairCat = gHades->getCurrentEvent()->getCategory(catPairs)) ) { Error("init", "No category catPairs"); return kFALSE; } pitPairCat = (HIterator *) pPairCat->MakeIterator(); if( NULL == (pGeantCat = gHades->getCurrentEvent()->getCategory(catGeantKine)) ) { Error("init", "No category catGeantKine "); return kFALSE; } pitGeant = (HIterator *)pGeantCat->MakeIterator(); pitGeant1 = (HIterator *)pGeantCat->MakeIterator(); // --------- MDC GEANT category if( NULL == gHades->getCurrentEvent()->getCategory(catMdcGeantRaw) ) { Error("init", "No category catMdcGeantRaw "); return kFALSE; } // --------- TOF GEANT category if( NULL == gHades->getCurrentEvent()->getCategory(catTofGeantRaw) ) { Error("init", "No category catTofGeantRaw "); return kFALSE; } // --------- SHOWER GEANT category if( NULL == gHades->getCurrentEvent()->getCategory(catShowerGeantRaw) ) { Error("init", "No category catTofGeantRaw "); return kFALSE; } // ---------- EVENT HEADER category if ( NULL == (pEvtHeader = (HEventHeader*) gHades->getCurrentEvent()->getHeader()) ) { Error("init", "Event header not found"); return kFALSE; } return kTRUE; } // ----------------------------------------------------------------------------- Int_t HPidSingleLeptonEff::execute(void) { // All the lepton pairs contained in the GeantKine container are formed // The yield is stored in a ntuple (simYieldN) as a function of the pair mass // rapidity and transverse momentum. Moreover the lepton pairs in the acceptance // are stored in a ntuple (accept). // The HPidDilepton class contains all the pairs that have been reconstructed. // For each of the GeantKine pair in the acceptance we checked if this pair // has been properly // recontructed and is contained in the HPidDilepton Container. // For the first correctly recontructed pair a second ntuple (eff) will // be filled with weight 1, for each additional reconstructed pair // that dont correpond to any Geant pair, a third ntuple (fake) will be filled // with weight 1. // All these informations will be stored in Ntuples to allow // later on the user to create of efficiency matrixes with different bining. HGeantKine *pKine = NULL; TVector3 pD; Float_t xMom,yMom,zMom; Int_t aTrackLepton, aIDLepton; Int_t aParLepton, aMedLepton, aMechLepton; Int_t* partIndex = new Int_t[maxNParticle]; for( Int_t i = 0; i < maxNParticle; ++i ) partIndex[i] = 0; //----------------------------------------------------------------- // do not allow tracks from pairs closer than 9 deg removecloseTracks(partIndex); //----------------------------------------------------------------- pitGeant->Reset(); while( NULL != (pKine = (HGeantKine*) pitGeant->Next()) ) { pKine->getParticle(aTrackLepton,aIDLepton); pKine->getCreator (aParLepton,aMedLepton,aMechLepton); if(aParLepton == 0 && (aIDLepton == 2 || aIDLepton == 3)) { if( aTrackLepton >= maxNParticle ) { Warning("execute","Track number of lepton exceeded maxNumber ! skipping!"); break; } pKine->getMomentum(xMom, yMom, zMom); pD.SetXYZ(xMom, yMom, zMom); //----------------------------------------------------------------- // fill ntuple for total yield, all angles are in degrees // in usual HADES range 0 < phi < 360 and 0 < theta < 90 if( aIDLepton == 3 ) simYieldEle->Fill(pD.Mag(),pD.Theta()*TMath::RadToDeg(),transformPhi(pD.Phi()),pD.Pt()); if( aIDLepton == 2 ) simYieldPos->Fill(pD.Mag(),pD.Theta()*TMath::RadToDeg(),transformPhi(pD.Phi()),pD.Pt()); //----------------------------------------------------------------- // check acceptance and if the track has been reconstructed. // fill corresponding ntuples for efficiency calculations. if (partIndex[aTrackLepton] == 0) checkEff(pKine); //----------------------------------------------------------------- } } delete [] partIndex; return 0; } // ----------------------------------------------------------------------------- void HPidSingleLeptonEff::removecloseTracks(Int_t* p) { // all the tracks coming fomr the white generator that // has a neighbour also coming from the white generator // that is closer than 9 deg are tagged if they // are inside the same sector HGeantKine *pKine1 = NULL; HGeantKine *pKine2 = NULL; Int_t aTrackLepton1, aIDLepton1; Int_t aTrackLepton2, aIDLepton2; Int_t aParLepton1, aMedLepton1, aMechLepton1; Int_t aParLepton2, aMedLepton2, aMechLepton2; pitGeant->Reset(); while( NULL != (pKine1 = (HGeantKine*) pitGeant->Next()) ) { pKine1->getParticle(aTrackLepton1, aIDLepton1); pKine1->getCreator (aParLepton1, aMedLepton1, aMechLepton1); if ( aParLepton1 == 0 && (aIDLepton1 == 2 || aIDLepton1 == 3) ) { if( aTrackLepton1 >= maxNParticle ) { Warning("removecloseTracks", "Track number of lepton exceeded maxNumber ! skipping!"); break; } pitGeant1->Reset(); while( NULL != (pKine2 = (HGeantKine*) pitGeant1->Next()) ) { if ( pKine1->getSector() != pKine2->getSector() ) continue; // only tracks inside the same sector can influence each other pKine2->getParticle(aTrackLepton2, aIDLepton2); pKine2->getCreator (aParLepton2, aMedLepton2, aMechLepton2); if(aParLepton2 == 0 && (aIDLepton2 == 2 || aIDLepton2 == 3) && aTrackLepton1 != aTrackLepton2 ) { if(calcOpeningAngleKine(pKine1, pKine2) < 9.0) { p[aTrackLepton1] = 1; } } } } } } // ----------------------------------------------------------------------------- Float_t HPidSingleLeptonEff::calcOpeningAngleKine(HGeantKine* kine1, HGeantKine* kine2) { //input 2 kine objects //output opening angle of trajectories TVector3 vec1; TVector3 vec2; Float_t xMom1,yMom1,zMom1; Float_t xMom2,yMom2,zMom2; if ( NULL != kine1 ) { kine1->getMomentum(xMom1, yMom1, zMom1); vec1.SetX(xMom1); vec1.SetY(yMom1); vec1.SetZ(zMom1); } else { Error("calcOpeningAngleKine", "First pointer is NULL"); return -1.; } if ( NULL != kine2 ) { kine2->getMomentum(xMom2, yMom2, zMom2); vec2.SetX(xMom2); vec2.SetY(yMom2); vec2.SetZ(zMom2); } else { Error("calcOpeningAngleKine", "Second pointer is NULL"); return -1.; } return TMath::RadToDeg() * vec1.Angle(vec2); } // ----------------------------------------------------------------------------- void HPidSingleLeptonEff::checkEff(HGeantKine* pKine) { // checks // 1. if GEANT track is in acceptance ==> fills the acceptance ntuple. // 2. loops over HPairSim category and counts number of reconstructed pairs // which found the fully correlated track in a pair containing hadron or leptons. // checks the direct cut flag for the pair. // 3. loops again over HPairSim and tries to find the searched track in a pair // 4. if track has not been found in a pair we do loop over HPidParticle category // and try to find the corresponding single track. // 5. lepton ntuple are filled in case 3 and 4. HPairSim* pair = NULL; HPidParticleSim* part = NULL; HPidParticleSim* part2 = NULL; const HPidHitData* pHitData = NULL; const HPidTrackData* pTrackData = NULL; TVector3 pD; TVector3 pairvertex; TVector3 pairdistance; Float_t xMom = 0.; Float_t yMom = 0.; Float_t zMom = 0.; Float_t thetaExp = 0.; Float_t phiExp = 0.; Float_t momExp = 0.; Float_t mmRK = -1.; Int_t aTrackLepton, aIDLepton; Int_t comDet = 0; Int_t whichParticle = -1; Int_t indexPart = -1; Bool_t isHiddenLepton = kFALSE; Bool_t flagPair = 0; Bool_t isGoodPair = 0; pKine->getMomentum(xMom, yMom, zMom); pKine->getParticle(aTrackLepton, aIDLepton); pD.SetXYZ(xMom, yMom, zMom); //----------------------------------------------------------------- // 1. if GEANT track is in acceptance and fills the acceptance ntuple if true. // All angles are in degrees // in usual HADES range 0 < phi < 360 and 0 < theta < 90 if( isGeantTrackInAcceptance(pKine) ) { if(aIDLepton == 2) accePos->Fill(pD.Mag(),pD.Theta() * TMath::RadToDeg(),transformPhi(pD.Phi()),pD.Pt()); if(aIDLepton == 3) acceEle->Fill(pD.Mag(),pD.Theta() * TMath::RadToDeg(),transformPhi(pD.Phi()),pD.Pt()); } else return; // nothing to do for this track //----------------------------------------------------------------- //----------------------------------------------------------------- // 3. loops again over HPairSim and tries to find the // searched track in a pair pitPairCat->Reset(); while ( NULL != ( pair = (HPairSim*) pitPairCat->Next()) ) { comDet = checkTrackId(pair, aTrackLepton, whichParticle, isHiddenLepton); if( comDet >= 76 ) // track == trackNb, lepton PID, primary, PID == rec PID { // not rejected flagPair = 1; if (pair->getIsCutFlag() != 1 && isGoodPair == 0 ) { // not direct cut, recursive cuting sets this variable to values > 1 // fully correlated track and first pair match found calculateTrackProperties(pair, thetaExp, phiExp, momExp, whichParticle); //----------------------------------------------------------------- // ----------- track & particle data retrieved here if (whichParticle == 1) indexPart = pair->getIndexParticle1(); else if (whichParticle == 2) indexPart = pair->getIndexParticle2(); else { Error("checkEff()","Wrong number for particle selection!"); continue; } part = (HPidParticleSim*) pPartCat->getObject(indexPart); pHitData = part->getHitData(); pTrackData = part->getTrackData(); mmRK = TMath::Sqrt(pTrackData->getRkMetadx() * pTrackData->getRkMetadx() + pTrackData->getRkMetady() * pTrackData->getRkMetady() + pTrackData->getRkMetadz() * pTrackData->getRkMetadz() ); HPairFL::calcVertex((HPidParticleSim*) pPartCat->getObject(pair->getIndexParticle1()), (HPidParticleSim*) pPartCat->getObject(pair->getIndexParticle2()), &pairvertex,&pairdistance); if (whichParticle == 1) part2 = (HPidParticleSim*) pPartCat->getObject(pair->getIndexParticle2()); else if (whichParticle == 2) part2 = (HPidParticleSim*) pPartCat->getObject(pair->getIndexParticle1()); //----------------------------------------------------------------- //----------------------------------------------------------------- // Fill ntuple for track found in a pair Float_t aVar[] = { momExp, pD.Mag(), thetaExp, pD.Theta() * TMath::RadToDeg(), phiExp, transformPhi(pD.Phi()), pD.Pt(), pair->getIsCutFlag(), // isCutNb comDet, // common detector pair->getBitFieldPassedCuts(), // for direct cut = info about passed cuts pTrackData->getAngleWithClosestCandidate(-1, 1), // hadron fitted, pTrackData->getAngleWithClosestCandidate(-1, -1), // hadron non fitted, pTrackData->getAngleWithClosestCandidate( 1, 1), // lepton fitted, pTrackData->getAngleWithClosestCandidate( 1, -1), // lepton non fitted, static_cast(pHitData->getRingCorrelation(4)), // has RK ring pHitData->getInnerMdcChiSquare(), // innerChi2, pTrackData->getRKChiSquare(), // rkChi2, pTrackData->getMetaMatchingQuality(), // metaMatchQuality, mmRK, // RK-META distance pEvtHeader->getVertex().getX(), // vertex X pEvtHeader->getVertex().getY(), // vertex Y pEvtHeader->getVertex().getZ(), // vertex Z pHitData->getOuterMdcChiSquare(), // outerChi, isHiddenLepton, // hidden lepton flag pTrackData->fTofRecBeta[4], // beta, pairvertex.x(), // pairVertx, pairvertex.y(), // pairVerty, pairvertex.z(), // pairVertz, pair->getOpeningAngleDeg(), // opening angle, part2->P(), // momentum of partner part2->getPid() // PID of partner }; if(aIDLepton == 3) effEle->Fill(aVar); if(aIDLepton == 2) effPos->Fill(aVar); //----------------------------------------------------------------- isGoodPair = 1; // we found the track, do not loop further } // end of iscutflag == 0 } // end of comDet > 0 } //----------------------------------------------------------------- //----------------------------------------------------------------- // 4. if track has not been found in a pair we do loop over // HPidParticle category and try to find the corresponding single track. if ( 0 == flagPair ) { pitPartCat->Reset(); while( NULL != (part = ( HPidParticleSim *)pitPartCat->Next()) ) { comDet = checkSingleTrack(part, aTrackLepton, isHiddenLepton); if(comDet >= 76) // track == trackNb, lepton PID, primary, PID == rec PID, fully correlated { //----------------------------------------------------------------- // ----------- track & particle data retrieved here pHitData = part->getHitData(); pTrackData = part->getTrackData(); mmRK = TMath::Sqrt(pTrackData->getRkMetadx() * pTrackData->getRkMetadx() + pTrackData->getRkMetady() * pTrackData->getRkMetady() + pTrackData->getRkMetadz() * pTrackData->getRkMetadz() ); //----------------------------------------------------------------- //----------------------------------------------------------------- // Fill ntuple for track not found in a pair Float_t aVar[] = { part->P(), pD.Mag(), part->thetaDeg(), pD.Theta() * TMath::RadToDeg(), part->phiDeg(), transformPhi(pD.Phi()), pD.Pt(), -1, // isCutNb comDet, -1, // bitFieldPassedCuts pTrackData->getAngleWithClosestCandidate(-1, 1), // hadron fitted, pTrackData->getAngleWithClosestCandidate(-1, -1), // hadron non fitted, pTrackData->getAngleWithClosestCandidate( 1, 1), // lepton fitted, pTrackData->getAngleWithClosestCandidate( 1, -1), // lepton non fitted, static_cast(pHitData->getRingCorrelation(4)), // has RK ring pHitData->getInnerMdcChiSquare(), // innerChi2, pTrackData->getRKChiSquare(), // rkChi2, pTrackData->getMetaMatchingQuality(), // metaMatchQuality, mmRK, // RK-META distance pEvtHeader->getVertex().getX(), // vertex X pEvtHeader->getVertex().getY(), // vertex Y pEvtHeader->getVertex().getZ(), // vertex Z pHitData->getOuterMdcChiSquare(), // outerChi, isHiddenLepton, // hidden lepton flag pTrackData->fTofRecBeta[4], // beta, -1000, // pairVertx -1000, // pairVerty -1000, // pairVertz -1, // opangle, -1, // pReconPartner -1 // PID_recPartner }; if(aIDLepton == 3) effEle ->Fill(aVar); if(aIDLepton == 2) effPos->Fill(aVar); //----------------------------------------------------------------- break; } // comDet >=76 } // end of HpidPartcle loop } // case there has been no pair found } // ----------------------------------------------------------------------------- Int_t HPidSingleLeptonEff::checkSingleTrack(HPidParticleSim* part, Int_t lepTrackNb, Bool_t& isHiddenLepton) { // Checks if the track euqal to the GEANT track, is primary, has lepton PID // and PID == GEANT PID. Return the most common detector correlation if successfull // or 0 if not. If the most common correlation misses the META ( < 76) and the // track was going to system 0 (SHOWER), we try to find if the lepton track // was hidden by a secondary. If so, the in commonCorrelation the META bit is set // and the isHiddenLepton flag becomes true. Int_t comDet = 0; HGeantKine* kine = NULL; isHiddenLepton = kFALSE; kine = (HGeantKine*) pGeantCat->getObject(lepTrackNb - 1); if(part->getGeantTrackSet()->getGeantTrackID() == lepTrackNb && // same track number part->getGeantTrackSet()->getGeantParentTrackID() == 0 && // primary part->getHitData()->getRingCorrelation(4) == kTRUE && // RK-RICH correlation (part->getPid() == 2 || part->getPid() == 3) && // lepton PID kine->getID() == part->getPid() // PID == GEANT PID ) { comDet = part->getGeantTrackSet()->getMostCommonCorrelation(); if( comDet < 76 && part->getSystem() == 0 ) { //---------------------------------------------------------------- // check if Track is sill registered in SHOWER // if so, change comDet and set isHiddenLepton UInt_t nShowerTracks = part->getGeantTrackSet()->getNShowerTracks(); if(nShowerTracks > 1){ for(UInt_t i = 0; i < nShowerTracks; i ++){ if(part->getGeantTrackSet()->getShowerTrack(i) == lepTrackNb){ comDet = comDet | kSeenInMETA; // set META bit isHiddenLepton = kTRUE; break; } } } //---------------------------------------------------------------- } else { isHiddenLepton = kFALSE; } return comDet; } else { return 0; } } // ----------------------------------------------------------------------------- Int_t HPidSingleLeptonEff::checkTrackId(HPairSim* p, Int_t lepTrackNb, Int_t& whichPart, Bool_t& isHiddenLepton) { // check if the searched track is contained in the pair // allow pairs with geant hadrons or geant secondary leptons // pairs with embedded leptons are rejected. Checks the // reconstructed PID for leptons and check reconstructed PID == GEANT PID // Int_t whichPart will be -1 if no track was accepted and 1 // or 2 if particle 1 or 2 of the pair was accepted. HPidParticleSim* part = NULL; HPairGeantData pg(p); Int_t comDet = 0; whichPart = -1; isHiddenLepton = kFALSE; part = (HPidParticleSim*) pPartCat->getObject(p->getIndexParticle1()); comDet = checkSingleTrack(part, lepTrackNb, isHiddenLepton); if ( comDet > 0 ) { // track == lepTrackNb, lepton PID, primary, PID == rec PID if ( (pg.getParentTrackNumber2() == 0 && (pg.getGPid2() != 2 && pg.getGPid2() != 3)) || // generator hadron or (pg.getParentTrackNumber2() != 0 )) // secondary lepton (decay in GEANT) or real tracks { whichPart = 1; return comDet; } else return 0; // reject pairs from embedded leptons } part = (HPidParticleSim*) pPartCat->getObject(p->getIndexParticle2()); comDet = checkSingleTrack(part,lepTrackNb,isHiddenLepton); if ( comDet > 0 ) { // track == lepTrackNb, lepton PID, primary, PID == rec PID if ( (pg.getParentTrackNumber1() == 0 && (pg.getGPid1() != 2 && pg.getGPid1() != 3)) || // generator hadron or (pg.getParentTrackNumber1() != 0 )) // secondary lepton (decay in GEANT) or real tracks { whichPart = 2; return comDet; } else return 0; // reject pairs from embedded leptons } return 0; } // ----------------------------------------------------------------------------- void HPidSingleLeptonEff::calculateTrackProperties(HPairSim* pair, Float_t& th, Float_t& ph, Float_t& p, Int_t whichParticle) { // calculates the pair theta, phi (in deg!!!) // Close pair candidates the angles to the closest // neighbours are stored separately HPairGeantData pg(pair); HPairData pd(pair); if ( 1 == whichParticle && 0 == pg.getParentTrackNumber1() ) { th = pd.getThetaDeg1(); ph = pd.getPhiDeg1(); p = pd.getMom1(); } else if ( 2 == whichParticle && 0 == pg.getParentTrackNumber2() ) { th = pd.getThetaDeg2(); ph = pd.getPhiDeg2(); p = pd.getMom2(); } else { Error("calculateTrackProperties()","Called with not identified Particle number!"); } } // ----------------------------------------------------------------------------- Bool_t HPidSingleLeptonEff::isGeantTrackInAcceptance(HGeantKine *pG) { // checks if thw Gent track is inside the Hades acceptance. // The acceptance is defined by a hit in the 4 MDCs and TOF // or SHOWER. Since Drift chambers can be taken out of the // HGeant geometry, for the outer drift chambers a hit in // MDC3 or MDC4 is required. Int_t nStatMDC1 = pG->getNMdcHits(0); Int_t nStatMDC2 = pG->getNMdcHits(1); Int_t nStatMDC3 = pG->getNMdcHits(2); Int_t nStatMDC4 = pG->getNMdcHits(3); Int_t nStatTof = pG->getNTofHits(); Int_t nStatShower = pG->getNShowerHits(); if(nStatMDC1 && nStatMDC2 && ( nStatMDC3 || nStatMDC4 ) && (nStatTof || nStatShower)) { return kTRUE; } return kFALSE; } // ----------------------------------------------------------------------------- Bool_t HPidSingleLeptonEff::finalize() { out->cd(); simYieldEle->Write(); simYieldPos->Write(); effEle ->Write(); effPos ->Write(); acceEle->Write(); accePos->Write(); out->Close(); return kTRUE; } // ----------------------------------------------------------------------------- Float_t HPidSingleLeptonEff::transformPhi(Float_t Phi) { Float_t dPhi; if( (dPhi = TMath::RadToDeg() * Phi) < 0.0 ) dPhi += 360.0; return dPhi; }