#include "hmdcfunc1.h" #include "hpiddilepton.h" #include "hpidparticle.h" #include "hpidparticlesim.h" #include "hpidtrackcand.h" #include "hpidcandidate.h" #include "hpiddilepton.h" #include "hmdcseg.h" #include "hmdchit.h" #include "hmdchitsim.h" #include "hmdcclusinf.h" #include "hmdcclusinfsim.h" #include "hrichhit.h" #include "hrichhitsim.h" #include "hkicktracksim.h" #include "hshowerhittof.h" #include "htofhit.h" #include "hgeantkine.h" #include "hgeantmdc.h" #include "hgeantshower.h" #include "hgeanttof.h" #include "hgeantrich.h" #include "hiterator.h" #include "hgeomvector.h" #include "TArrayI.h" #include "TVector3.h" #include "hpidcpparam.h" Int_t HMdcFunc1::getMdcClsSize(HPidParticle *part,Int_t module){ HMdcSeg * pMdcSeg=NULL; HPidTrackCand * pPidTrackCand=NULL; Int_t cls=0; if(part) { pPidTrackCand=getPidTrackCand(part); if(pPidTrackCand){ pMdcSeg=(HMdcSeg*)pPidTrackCand->getMdcSeg(); if(pMdcSeg) { cls=getMdcClsSize(pMdcSeg,module); } } } return cls; } Int_t HMdcFunc1::getMdcNWires(HPidParticle *part,Int_t module){ HMdcSeg * pMdcSeg=NULL; HPidTrackCand * pPidTrackCand=NULL; Int_t nwires=0; if(part) { pPidTrackCand=getPidTrackCand(part); if(pPidTrackCand){ pMdcSeg=(HMdcSeg*)pPidTrackCand->getMdcSeg(); if(pMdcSeg) { nwires=getMdcNWires(pMdcSeg,module); } } } return nwires; } Int_t HMdcFunc1::getMdcLevelCls(HPidParticle *part,Int_t module){ HMdcSeg * pMdcSeg=NULL; HPidTrackCand * pPidTrackCand=NULL; Int_t level=0; if(part) { pPidTrackCand=getPidTrackCand(part); if(pPidTrackCand){ pMdcSeg=(HMdcSeg*)pPidTrackCand->getMdcSeg(); if(pMdcSeg) { level=getMdcLevelCls(pMdcSeg,module); } } } return level; } Int_t HMdcFunc1::getIntCharge(HPidParticle *part){ HRichHit * pRichHit=NULL; HPidTrackCand * pPidTrackCand=NULL; Int_t int_charge=0; if(part) { pPidTrackCand=getPidTrackCand(part); if(pPidTrackCand){ pRichHit=(HRichHit*)pPidTrackCand->getRichHit(); if(pRichHit) { int_charge=pRichHit->getRingAmplitude(); } } } return int_charge; } Int_t HMdcFunc1::getMdcClsSize(HMdcSeg *seg,Int_t module){ Int_t cls=0; HMdcHit *hit; if(seg) { hit=getMdcHit(seg,module); if(hit) cls=getMdcClsSize(hit); } return cls; } Int_t HMdcFunc1::getMdcNWires(HMdcSeg *seg,Int_t module){ Int_t nwires=0; HMdcHit *hit; if(seg) { hit=getMdcHit(seg,module); if(hit) nwires=getMdcNWires(hit); } return nwires; } Int_t HMdcFunc1::getMdcLevelCls(HMdcSeg *seg,Int_t module){ Int_t level=0; HMdcHit *hit; if(seg) { hit=getMdcHit(seg,module); if(hit) level=getMdcLevelCls(hit); } return level; } HPidTrackCand* HMdcFunc1::getPidTrackCand(HPidParticle *part){ HPidCandidate* pidcand=NULL; HPidTrackCand* pidtrackcand=NULL; if(part) { pidcand=part->getPidCandidate(); if(pidcand) { pidtrackcand=pidcand->getCandidate(); } } return pidtrackcand; } HMdcHit* HMdcFunc1::getMdcHit(HMdcSeg *seg,Int_t module){ HMdcHit *hit=NULL; Int_t hit_inx=-1; if(seg) { hit_inx=seg->getHitInd(module); if(hit_inx>-1){ if(getMdcHitCat()) { // cout<getEntries()<getObject(hit_inx); // cout<getClusInfIndex(); pObj_clusinf=getClusInfObj(clsindex); if(pObj_clusinf){ clssize=pObj_clusinf->getClusSize(); } return clssize; } Int_t HMdcFunc1::getMdcNWires(HMdcHit *hit){ Int_t nwires=-1; Int_t clsindex=-1; HMdcClusInf *pObj_clusinf=NULL; clsindex=hit->getClusInfIndex(); pObj_clusinf=getClusInfObj(clsindex); if(pObj_clusinf){ nwires=pObj_clusinf->getNDrTimes(); } return nwires; } Int_t HMdcFunc1::getMdcLevelCls(HMdcHit *hit){ Int_t level=-1; Int_t clsindex=-1; HMdcClusInf *pObj_clusinf=NULL; clsindex=hit->getClusInfIndex(); pObj_clusinf=getClusInfObj(clsindex); if(pObj_clusinf){ level=pObj_clusinf->getLevelClFinding(); } return level; } HMdcClusInf* HMdcFunc1::getClusInfObj(Int_t clsindex){ //get ClusInfObj from category using index HMdcClusInf *pObj_clusinf=NULL; if(getClusInfCat()){ pObj_clusinf=(HMdcClusInf* )getClusInfCat()->getObject(clsindex); } else { printf("******MdcClusInfCat is NULL********\n"); printf("******use HMdcFunc1::setMdcClusInfCat()********\n"); } return pObj_clusinf; } Int_t HMdcFunc1::calculateLevelBin(Int_t level) { Int_t mdc_lev_bin =-1; if (level>=4&&level<=6) mdc_lev_bin=(Int_t )level-4; return mdc_lev_bin; } Int_t HMdcFunc1::calculateBin(Float_t mdc_angle,Int_t option) { // calculate the the bin of the angle (option=0 ... theta; =1 ... phi) Int_t mdc_theta_bin=-9; Int_t mdc_phi_bin=-9; Int_t theta_bin_size=10; Int_t phi_bin_size=10; Float_t theta_min=10.0; Float_t theta_max=90.0; Float_t phi_min=0.0; Float_t phi_max=360.0; switch (option) { case 0: // case theta if(mdc_angle>=theta_min&&mdc_angle<=theta_max) { mdc_theta_bin=((Int_t)(mdc_angle-theta_min))/theta_bin_size; // printf("mdc_angle=%f bin=%i \n",mdc_angle,mdc_theta_bin); return mdc_theta_bin; } else { printf ("theta out of range\n"); return -2; } case 1: if(mdc_angle>=phi_min&&mdc_angle<=phi_max){ mdc_phi_bin=(((Int_t) mdc_angle)%60)/phi_bin_size; if(mdc_phi_bin>2) mdc_phi_bin=5-mdc_phi_bin ;//e.g. 10-10 and 50-60 is bin 0 return mdc_phi_bin; }else{ printf ("phi out of range\n"); return -2; } default: printf("Third argument should be 0-1 for theta-phi angle\n"); return -1; } } Float_t HMdcFunc1::getNormalMdcPhi(Int_t iSector, Float_t fPhiMdc) { // Convert MDC's phi angle to the coordinate system used in other detectors // input mdc phi in rad and sector coordinate system => output in deg and hades coordinate system Float_t fPhi = TMath::RadToDeg() * fPhiMdc; switch(iSector) { case 0: break; case 1: case 2: case 3: case 4: fPhi += 60.0f * iSector; break; default: fPhi -= 60.0f; break; } return fPhi; } // ----------------------------------------------------------------------------- Float_t HMdcFunc1::getNormalMdcTheta(Float_t fThetaMdc) { // Convert MDC's theta angle to the coordinate system used in other detectors return TMath::RadToDeg() * fThetaMdc; } Int_t HMdcFunc1::getGeantId(HRichHitSim * pRichHitSim,HKickTrackSim * pKickTrack){ //this calculated GeantId in a way I consider it correct //return 2,3 or 102 and 103 if e+ e- // -2 ring-pi // -3 ring-p // -5 ring-another lepton // -4 ring-mdc corr ele-meta not // -1 ring-non sence in meta HGeantKine *kine=NULL; Int_t richTr1=-1; Int_t richTr2=-1; Int_t richTr3 =-1; Int_t iTrack=-1, iId=-1; Int_t flag_richmdcele=0; Int_t flag_missinmdc=0; if(pKickTrack!=NULL&&pRichHitSim!=NULL) { if(pRichHitSim->weigTrack1>0) richTr1=pRichHitSim->track1; if(pRichHitSim->weigTrack2>0) richTr2=pRichHitSim->track2; if(pRichHitSim->weigTrack3>0) richTr3=pRichHitSim->track3; //printf("THIS TRACK: RICH %i %i %i MDC (num %i) %i %i %i %i %i META (num %i) %i %i %i\n",richTr1,richTr2,richTr3,pKickTrack->getNMdcTracks(),pKickTrack->getMdcTrack(0),pKickTrack->getMdcTrack(1),pKickTrack->getMdcTrack(2),pKickTrack->getMdcTrack(3),pKickTrack->getMdcTrack(4),pKickTrack->getNMetaTracks(),pKickTrack->getMetaTrack(0),pKickTrack->getMetaTrack(1),pKickTrack->getMetaTrack(2)); flag_richmdcele=0; for (Int_t j=0;jgetNMdcTracks();j++){ //go over mdc iId=-1; if(pKickTrack->getMdcTrack(j)>0&&pKickTrack->getMdcTrack(j)<10000){ for (Int_t k=0;kgetNMetaTracks();k++){ //iId=-1; //completely new if(pKickTrack->getMdcTrack(j)==pKickTrack->getMetaTrack(k)){ //common MDC-META kine=HPidParticleSim::getGeantKine(pKickTrack->getMdcTrack(j)); kine->getParticle(iTrack, iId); if((pKickTrack->getMdcTrack(j)==richTr1)||(pKickTrack->getMdcTrack(j)==richTr2)||(pKickTrack->getMdcTrack(j)==richTr3)){ //common RING-MDC-META //iId=-1; //kine=HPidParticleSim::getGeantKine(pKickTrack->getMdcTrack(j)); //kine->getParticle(iTrack, iId); //printf("Rich-Meta-Mdc common track=%i part=%i\n",pKickTrack->getMdcTrack(j),iId); //printf("return %i\n",iId); return iId; } if(iId==9||iId==8) { //printf("Ring corr Meta-Mdc common track=%i part=%i\n",pKickTrack->getMdcTrack(j),iId); //printf("return -2 (pion)\n"); return -2; } if(iId==14){ // printf("Ring corr Meta-Mdc common track=%i part=%i\n",pKickTrack->getMdcTrack(j),iId); //printf("return -3 (proton)\n"); return -3; } if(iId==3||iId==2){ // printf("Ring corr Meta-Mdc common track=%i part=%i\n",pKickTrack->getMdcTrack(j),iId); //printf("return -5 (another lepton)\n"); return -5; } } if((pKickTrack->getMdcTrack(j)==richTr1)||(pKickTrack->getMdcTrack(j)==richTr2)||(pKickTrack->getMdcTrack(j)==richTr3)){ //common RING-MDC but not META because otherwise it is already return before kine=HPidParticleSim::getGeantKine(pKickTrack->getMdcTrack(j)); kine->getParticle(iTrack, iId); if(iId==2||iId==3){ //common RING-MDC electron (but not META) // printf("common RING-MDC lepton (but not META) tr=%i \n",pKickTrack->getMdcTrack(j)); //printf("set iId=-4\n"); flag_richmdcele=1; } else { printf("Strange RING-MDC is not lepton\n"); } } //one has to treat separately case that in MDC I do not see a second pair from conv but in META yes if((pKickTrack->getMetaTrack(k)==richTr1)||(pKickTrack->getMetaTrack(k)==richTr2)||(pKickTrack->getMetaTrack(k)==richTr3)){ for (Int_t l=0;lgetNMdcTracks();l++){ if(TMath::Abs(pKickTrack->getMdcTrack(l)-pKickTrack->getMetaTrack(k))==1){ // printf("common RING-MDC lepton1 ,common RING-META lepton2 from pair tr=%i in mdc lost \n",pKickTrack->getMetaTrack(l)); kine=HPidParticleSim::getGeantKine(pKickTrack->getMetaTrack(k)); kine->getParticle(iTrack, iId); //printf("common miss mdc track iTrack=%i set flag_missinmdc=%i\n ",pKickTrack->getMetaTrack(l),(100+iId)); //return (100+iId); //I do not return because it could still be double-double-single case flag_missinmdc=(100+iId); } } } } } } //printf("THIS TRACK: RICH %i %i %i MDC (num %i) %i %i %i %i %i META (num %i) %i %i %i\n",richTr1,richTr2,richTr3,pKickTrack->getNMdcTracks(),pKickTrack->getMdcTrack(0),pKickTrack->getMdcTrack(1),pKickTrack->getMdcTrack(2),pKickTrack->getMdcTrack(3),pKickTrack->getMdcTrack(4),pKickTrack->getNMetaTracks(),pKickTrack->getMetaTrack(0),pKickTrack->getMetaTrack(1),pKickTrack->getMetaTrack(2)); if(flag_missinmdc!=0){ return flag_missinmdc; } if(flag_richmdcele==1){ // printf("return %i\n",-4); return -4; } //printf("return %i\n",-1); return -1; } return 0; } Float_t HMdcFunc1::getKickAngle(HPidParticle * part){ //return META-MDC theta diff //HMdcSeg * pMdcSeg=NULL; //HPidTrackCand * pPidTrackCand=NULL; Float_t meta_theta=-100.0; Float_t mdc_the=100.0; // Int_t cls=0; // if(part) { // pPidTrackCand=getPidTrackCand(part); // if(pPidTrackCand){ // pMdcSeg=(HMdcSeg*)pPidTrackCand->getMdcSeg(); // mdc_the=getNormalMdcTheta(pMdcSeg->getTheta()); // if(pPidTrackCand->getSystem()==0){ // ((HShowerHitTof*) pPidTrackCand->getShowerHitTof())->getSphereCoord(dummy,dummy,meta_theta); // } else { // ((HTofHit*) pPidTrackCand->getTofHit())->getTheta(meta_theta); // } // } // } return (meta_theta- mdc_the); } Int_t HMdcFunc1::getGeantCommonTrack(HRichHitSim * pRichHitSim,HKickTrackSim * pKickTrack){ //this calculated CommonGeantTrack in a way I consider it is correct //return positive value for e+ e- common track and: // -2 ring-pi // -3 ring-p // -5 ring-another lepton // -4 ring-mdc corr ele-meta not // -1 ring-non sence in meta // -10 richhit is null pointer HGeantKine *kine=NULL; Int_t richTr1=-1; Int_t richTr2=-1; Int_t richTr3 =-1; Int_t iTrack=-1, iId=-1; Int_t flag_richmdcele=0; Int_t flag_missinmdc=0; if(pKickTrack!=NULL&&pRichHitSim!=NULL) { if(pRichHitSim->weigTrack1>0) richTr1=pRichHitSim->track1; if(pRichHitSim->weigTrack2>0) richTr2=pRichHitSim->track2; if(pRichHitSim->weigTrack3>0) richTr3=pRichHitSim->track3; flag_richmdcele=0; for (Int_t j=0;jgetNMdcTracks();j++){ //go over mdc iId=-1; if(pKickTrack->getMdcTrack(j)>0&&pKickTrack->getMdcTrack(j)<10000){ for (Int_t k=0;kgetNMetaTracks();k++){ //iId=-1; //completely new if(pKickTrack->getMdcTrack(j)==pKickTrack->getMetaTrack(k)){ //common MDC-META kine=HPidParticleSim::getGeantKine(pKickTrack->getMdcTrack(j)); kine->getParticle(iTrack, iId); if((pKickTrack->getMdcTrack(j)==richTr1)||(pKickTrack->getMdcTrack(j)==richTr2)||(pKickTrack->getMdcTrack(j)==richTr3)){ //common RING-MDC-META //iId=-1; //kine=HPidParticleSim::getGeantKine(pKickTrack->getMdcTrack(j)); //kine->getParticle(iTrack, iId); //printf("Rich-Meta-Mdc common track=%i part=%i\n",pKickTrack->getMdcTrack(j),iId); //printf("return %i\n",iId); return pKickTrack->getMdcTrack(j); } if(iId==9||iId==8) { //printf("Ring corr Meta-Mdc common track=%i part=%i\n",pKickTrack->getMdcTrack(j),iId); //printf("return -2 (pion)\n"); return -2; } if(iId==14){ // printf("Ring corr Meta-Mdc common track=%i part=%i\n",pKickTrack->getMdcTrack(j),iId); //printf("return -3 (proton)\n"); return -3; } if(iId==3||iId==2){ // printf("Ring corr Meta-Mdc common track=%i part=%i\n",pKickTrack->getMdcTrack(j),iId); //printf("return -5 (another lepton)\n"); return -5; } } if((pKickTrack->getMdcTrack(j)==richTr1)||(pKickTrack->getMdcTrack(j)==richTr2)||(pKickTrack->getMdcTrack(j)==richTr3)){ //common RING-MDC but not META because otherwise it is already return before kine=HPidParticleSim::getGeantKine(pKickTrack->getMdcTrack(j)); kine->getParticle(iTrack, iId); if(iId==2||iId==3){ //common RING-MDC electron (but not META) // printf("common RING-MDC lepton (but not META) tr=%i \n",pKickTrack->getMdcTrack(j)); //printf("set iId=-4\n"); flag_richmdcele=1; } else { printf("Strange RING-MDC is not lepton\n"); } } //one has to treat separately case that in MDC I do not see a second pair from conv but in META yes if((pKickTrack->getMetaTrack(k)==richTr1)||(pKickTrack->getMetaTrack(k)==richTr2)||(pKickTrack->getMetaTrack(k)==richTr3)){ for (Int_t l=0;lgetNMdcTracks();l++){ if(TMath::Abs(pKickTrack->getMdcTrack(l)-pKickTrack->getMetaTrack(k))==1){ // printf("common RING-MDC lepton1 ,common RING-META lepton2 from pair tr=%i in mdc lost \n",pKickTrack->getMetaTrack(l)); kine=HPidParticleSim::getGeantKine(pKickTrack->getMetaTrack(k)); kine->getParticle(iTrack, iId); //printf("common miss mdc track iTrack=%i set flag_missinmdc=%i\n ",pKickTrack->getMetaTrack(l),(100+iId)); //return (100+iId); //I do not return because it could still be double-double-single case // flag_missinmdc=(100+iId); flag_missinmdc=pKickTrack->getMetaTrack(k); } } } } } } // printf("THIS TRACK: RICH %i %i %i MDC (num %i) %i %i %i %i %i META (num %i) %i %i %i\n",richTr1,richTr2,richTr3,pKickTrack->getNMdcTracks(),pKickTrack->getMdcTrack(0),pKickTrack->getMdcTrack(1),pKickTrack->getMdcTrack(2),pKickTrack->getMdcTrack(3),pKickTrack->getMdcTrack(4),pKickTrack->getNMetaTracks(),pKickTrack->getMetaTrack(0),pKickTrack->getMetaTrack(1),pKickTrack->getMetaTrack(2)); if(flag_missinmdc!=0){ return flag_missinmdc; } if(flag_richmdcele==1){ // printf("return %i\n",-4); return -4; } //printf("return %i\n",-1); return -1; } return -10; } HGeantKine* HMdcFunc1::getKineObj(Int_t pTrack){ if(fGeantKineCat!=NULL) { return (HGeantKine*)(fGeantKineCat->getObject(pTrack-1)); } else { printf("HMdcFunc1::getKineObj there is no fGeantKineCat defined please use HMdcFunc1::setGeantKineCat before\n"); return NULL; } } Int_t HMdcFunc1::getGeantCommonTrack(HPidParticleSim *part){ HPidTrackCand* pPidTrackCand=NULL; HKickTrackSim *pKickTrack=NULL; HRichHitSim *pRichHitSim=NULL; if (part!=NULL) { pPidTrackCand=getPidTrackCand(part); if (pPidTrackCand!=NULL) { pKickTrack=(HKickTrackSim *) pPidTrackCand->getKickTrack(); pRichHitSim=(HRichHitSim *) pPidTrackCand->getRichHit(); return getGeantCommonTrack(pRichHitSim,pKickTrack); } } return -10; } Bool_t HMdcFunc1::isCompt(HGeantKine *kine){ Int_t pTrack_par=-1,pMed=-1,pMech=-1,pID_par=0; Int_t pTrack_particle=-1,pID_particle=0; if(isLepton(kine)) { kine->getParticle(pTrack_particle,pID_particle); kine->getCreator(pTrack_par,pMed,pMech); if(pTrack_par!=0){ getKineObj(pTrack_par)->getParticle(pTrack_par,pID_par); //pMech==7 means that it was compt gama->e- gama if(pID_par==1&&pMech==7){ return kTRUE ; } } } return kFALSE; } Bool_t HMdcFunc1::isPi0Conv(HGeantKine *kine){ Int_t pTrack_par=-1,pMed=-1,pMech=-1,pID_par=0,pTrack_par_par=-1,pMed_par=-1,pMech_par=-1,pID_par_par=-1; Int_t pTrack_particle=-1,pID_particle=0; if(isLepton(kine)) { kine->getParticle(pTrack_particle,pID_particle); kine->getCreator(pTrack_par,pMed,pMech); if(pTrack_par!=0){ getKineObj(pTrack_par)->getParticle(pTrack_par,pID_par); //pMech==6 means that it was pair creation gama->e+e- if(pID_par==1&&pMech==6){ getKineObj(pTrack_par)->getCreator(pTrack_par_par,pMed_par,pMech_par); if(pTrack_par_par!=0){ getKineObj(pTrack_par_par)->getParticle(pTrack_par_par,pID_par_par); // cout<<"in checkConv prarodic ="<getParticle(pTrack_particle,pID_particle); kine->getCreator(pTrack_par,pMed,pMech); if(pTrack_par!=0){ getKineObj(pTrack_par)->getParticle(pTrack_par,pID_par); //pMech==6 means that it was pair creation gama->e+e- if(pID_par==1&&pMech==6){ return kTRUE ; // getKineObj(pTrack_par)->getCreator(pTrack_par_par,pMed_par,pMech_par); // if(pTrack_par_par!=0){ // getKineObj(pTrack_par_par)->getParticle(pTrack_par_par,pID_par_par); // // cout<<"in checkConv prarodic ="<getParticle(pTrack_particle,pID_particle); kine->getCreator(pTrack_par,pMed,pMech); if(pTrack_par!=0){ getKineObj(pTrack_par)->getParticle(pTrack_par,pID_par); if(pID_par==7){ return kTRUE; } } } return kFALSE; } Bool_t HMdcFunc1::isEtaDalitz(HGeantKine *kine){ Int_t pTrack_par=-1,pMed=-1,pMech=-1,pID_par=0; Int_t pTrack_particle=-1,pID_particle=0; if(isLepton(kine)) { kine->getParticle(pTrack_particle,pID_particle); kine->getCreator(pTrack_par,pMed,pMech); if(pTrack_par!=0){ getKineObj(pTrack_par)->getParticle(pTrack_par,pID_par); if(pID_par==17){ return kTRUE; } } } return kFALSE; } Bool_t HMdcFunc1::isLepton(HGeantKine *kine){ Int_t pTrack_particle=-1,pID_particle=0; kine->getParticle(pTrack_particle,pID_particle); if (pID_particle==2||pID_particle==3) return kTRUE; return kFALSE; } Bool_t HMdcFunc1::isSingle(HMdcSegSim *seg){ //The method operates only on simulation data. //Using the geant track numbers of the tracks in mdcseg (5 tracks stored) decide whether //it is a single tracks from lepton or not //Check also whether it is inner segment. If not returns false. Int_t track_mul; //look only at inner segments :very important if(seg->getIOSeg()==0){ track_mul=seg->getNTracks(); //it gives number of tracks (max) stored for this segment if(track_mul==1){ //only segment that consist of one track can be a ideal single (MAYBE TO STRONG HAS TO BE RECONSIDER) if(isLepton(getKineObj(seg->getTrack(0)))) { //this is segment that consist only from one geant track that is an lepton => we call it SINGLE return kTRUE; } } } return kFALSE; } Bool_t HMdcFunc1::isSingle(HMdcHitSim *hit){ Int_t track_mul; Int_t clsindex=-1; HMdcClusInfSim *pObj_clusinf=NULL; clsindex=hit->getClusInfIndex(); pObj_clusinf=(HMdcClusInfSim*) getClusInfObj(clsindex); if(pObj_clusinf){ track_mul=pObj_clusinf->getNTracks(); if(track_mul==1){ if(isLepton(getKineObj(pObj_clusinf->getTrack(0)))) { return kTRUE; } } } return kFALSE; } Int_t HMdcFunc1::isSingleDouble(HMdcSegSim *seg,Int_t track,Int_t conv_mode){ //this method have following return values // 0 ... normal single // 1 ... normal double // 2 ... double => because possition of pair lepton to track in MdcGeant too close to MdcHit possition //-1 ... nothing of previous // return value 2 was introduced because the lack of complete propagation of all contributing track numbers // from contributing wires to segment for fitted one. This is in version v6_11 Int_t track_in_pair=-1; if(track>0){ if(isDouble(seg,conv_mode)) return 1; if(isSingle(seg)){ //here I have check position of pair lepton => if too close than I consider it as double, because //lost info after fitter remove some wires track_in_pair=getPairTrack(track); if(isMdcSegInGeant(seg,track_in_pair)) return 2; return 0; } } return -1; } Bool_t HMdcFunc1::isDouble(HMdcSegSim *seg,Int_t conv_mode){ //The method operates only on simulated data. //Using the geant track numbers of the tracks in mdcseg (5 tracks stored) decide whether //it is a double tracks from lepton or not //conv_more =0 seg is DOUBLE if there are 2 different leptons contributing to it //conv_more =1 seg is DOUBLE if there are 2 different leptons contributing to it and leptons are conversion pair Int_t track_mul=0; Int_t pID2=0,pID1=-1,pTrack_par1=-1,pTrack_par2=-2,dummy=-1; //look only at inner segments :very important if(seg->getIOSeg()==0){ track_mul=seg->getNTracks(); //it gives number of tracks (max) stored for this segment if(track_mul>1){ //only segment that consist more than one track can be a ideal double for(Int_t i=0;igetTrack(i))),isPi0Conv(getKineObj(seg->getTrack(j)))); // check i-th track in seg if(isLepton(getKineObj(seg->getTrack(i)))) { if(conv_mode==1){ //in this mode I check whether it is conversion lepton if(isPi0Conv(getKineObj(seg->getTrack(i)))){ getKineObj(seg->getTrack(i))->getParticle(dummy,pID1); getKineObj(seg->getTrack(i))->getCreator(pTrack_par1,dummy,dummy); } } else { getKineObj(seg->getTrack(i))->getParticle(dummy,pID1); getKineObj(seg->getTrack(i))->getCreator(pTrack_par1,dummy,dummy); } } // check j-th track in seg if(isLepton(getKineObj(seg->getTrack(j)))) { if(conv_mode==1){ //in this mode I check whether it is conversion lepton if(isPi0Conv(getKineObj(seg->getTrack(j)))){ getKineObj(seg->getTrack(j))->getParticle(dummy,pID2); getKineObj(seg->getTrack(j))->getCreator(pTrack_par2,dummy,dummy); } } else { getKineObj(seg->getTrack(j))->getParticle(dummy,pID2); getKineObj(seg->getTrack(j))->getCreator(pTrack_par2,dummy,dummy); } } // printf("track_mul=%i i=%i j=%i pid1=%i pid2=%i pTrack_par1=%i pTrack_par2=%i \n",track_mul,i,j,pID1,pID2,pTrack_par1,pTrack_par2); if((((pID1==2&&pID2==3)||(pID2==2&&pID1==3))&&(pTrack_par2==pTrack_par1))||(conv_mode!=1)){ //This is a DOUBLE mdcseg (either 2 leptons or conversion pair (depending on the conv_mode)) //if conv_mode!=1 then I do not check weather both leptons are e+e- from one mother particle //check for leptons or conv is done already above return kTRUE; } } } } } return kFALSE; } Bool_t HMdcFunc1::isDouble(HMdcHitSim *hit,Int_t conv_mode){ //The method operates only on simulated data. //Using the geant track numbers of the tracks in mdcseg (5 tracks stored) decide whether //it is a double tracks from lepton or not //conv_more =0 seg is DOUBLE if there are 2 different leptons contributing to it and they are pair //conv_more =1 seg is DOUBLE if there are 2 different leptons contributing to it and leptons are conversion pair //Update from v6_12 for complete list of track num related to hit one has to look at MdcClusInf Int_t debug=0; Int_t clsindex=-1; Int_t track_mul=0; Int_t pID2=0,pID1=-1,pTrack_par1=-1,pTrack_par2=-2,dummy=-1; HMdcClusInfSim *pObj_clusinf=NULL; clsindex=hit->getClusInfIndex(); pObj_clusinf=(HMdcClusInfSim*) getClusInfObj(clsindex); if(pObj_clusinf){ track_mul=pObj_clusinf->getNTracks(); //it gives number of tracks (max) stored for this segment if(track_mul>1&&pObj_clusinf->getNTimesInTrack(1)>2){ //only segment that consist more than one track can be a ideal double if (debug==1){ if (track_mul==2) printf("mul=%i t0=%i(%i) t1=%i(%i)\n",track_mul,pObj_clusinf->getTrack(0),pObj_clusinf->getNTimesInTrack(0),pObj_clusinf->getTrack(1),pObj_clusinf->getNTimesInTrack(1)); if (track_mul==3) printf("mul=%i t0=%i(%i) t1=%i(%i) t2=%i(%i) \n",track_mul,pObj_clusinf->getTrack(0),pObj_clusinf->getNTimesInTrack(0),pObj_clusinf->getTrack(1),pObj_clusinf->getNTimesInTrack(1),pObj_clusinf->getTrack(2),pObj_clusinf->getNTimesInTrack(2)); if (track_mul==4) printf("mul=%i t0=%i(%i) t1=%i(%i) t2=%i(%i) t3=%i(%i)\n",track_mul,pObj_clusinf->getTrack(0),pObj_clusinf->getNTimesInTrack(0),pObj_clusinf->getTrack(1),pObj_clusinf->getNTimesInTrack(1),pObj_clusinf->getTrack(2),pObj_clusinf->getNTimesInTrack(2),pObj_clusinf->getTrack(3),pObj_clusinf->getNTimesInTrack(3)); if (track_mul==5) printf("mul=%i t0=%i(%i) t1=%i(%i) t2=%i(%i) t3=%i(%i) t4=%i(%i)\n",track_mul,pObj_clusinf->getTrack(0),pObj_clusinf->getNTimesInTrack(0),pObj_clusinf->getTrack(1),pObj_clusinf->getNTimesInTrack(1),pObj_clusinf->getTrack(2),pObj_clusinf->getNTimesInTrack(2),pObj_clusinf->getTrack(3),pObj_clusinf->getNTimesInTrack(3),pObj_clusinf->getTrack(4),pObj_clusinf->getNTimesInTrack(4)); } for(Int_t i=0;igetTrack(i))),isPi0Conv(getKineObj(seg->getTrack(j)))); // check i-th track in seg if(isLepton(getKineObj(pObj_clusinf->getTrack(i)))) { if(conv_mode==1){ //in this mode I check whether it is conversion lepton if(isPi0Conv(getKineObj(pObj_clusinf->getTrack(i)))){ getKineObj(pObj_clusinf->getTrack(i))->getParticle(dummy,pID1); getKineObj(pObj_clusinf->getTrack(i))->getCreator(pTrack_par1,dummy,dummy); } } else { getKineObj(pObj_clusinf->getTrack(i))->getParticle(dummy,pID1); getKineObj(pObj_clusinf->getTrack(i))->getCreator(pTrack_par1,dummy,dummy); } } // check j-th track in seg if(isLepton(getKineObj(pObj_clusinf->getTrack(j)))) { if(conv_mode==1){ //in this mode I check whether it is conversion lepton if(isPi0Conv(getKineObj(pObj_clusinf->getTrack(j)))){ getKineObj(pObj_clusinf->getTrack(j))->getParticle(dummy,pID2); getKineObj(pObj_clusinf->getTrack(j))->getCreator(pTrack_par2,dummy,dummy); } } else { getKineObj(pObj_clusinf->getTrack(j))->getParticle(dummy,pID2); getKineObj(pObj_clusinf->getTrack(j))->getCreator(pTrack_par2,dummy,dummy); } } // printf("track_mul=%i i=%i j=%i pid1=%i pid2=%i pTrack_par1=%i pTrack_par2=%i \n",track_mul,i,j,pID1,pID2,pTrack_par1,pTrack_par2); if(((pID1==2&&pID2==3)||(pID2==2&&pID1==3))&&(pTrack_par2==pTrack_par1)){ //This is a DOUBLE mdcseg (either 2 leptons from the same mother or conversion pair //(depending on the conv_mode)) if(debug==1) printf(" Module=%i DOUBLE %i %i \n",hit->getModule(),pObj_clusinf->getTrack(i),pObj_clusinf->getTrack(j)); return kTRUE; } } } } } return kFALSE; } Bool_t HMdcFunc1::isMdcSegInGeant(HMdcSegSim *seg,Int_t track){ //this method look if mdcseg that corresponds to geant tracknumber track if found at the same place in GeantMdc Float_t xgeant=1000.0,ygeant=1000.0,dummy=1000.0,x_hit=-1000.0,y_hit=-1000.0; Float_t Xdiff=3*2.5,Ydiff=3*1.5; Int_t module_geant=-1; HGeantMdc *pObj_geantmdc=NULL; if(getGeantMdcCat()){ HIterator *iterator_mdcgeant= (HIterator*)getGeantMdcCat()->MakeIterator(); iterator_mdcgeant->Reset(); while ((pObj_geantmdc=(HGeantMdc*)iterator_mdcgeant->Next())!= 0) { if(pObj_geantmdc->getLayer()!=6||(pObj_geantmdc->getModule()==2||pObj_geantmdc->getModule()==3)) continue; if(pObj_geantmdc->getTrack()==track){ //this pObj_geantmdc object has the same track number as this segment pObj_geantmdc->getHit(xgeant,ygeant,dummy,dummy); module_geant=pObj_geantmdc->getModule(); if(getMdcHit(seg,module_geant)){ x_hit=getMdcHit(seg,module_geant)->getX(); y_hit=getMdcHit(seg,module_geant)->getY(); } // printf("module=%i xgeant=%f ygeant=%f x_hit=%f y_hit=%f\n",module_geant,xgeant,ygeant,x_hit,y_hit); if(TMath::Abs(x_hit-xgeant)=0){ getKineObj(track1)->getCreator(pTrack_par1,dummy,dummy); getKineObj(track1+1)->getCreator(pTrack_par2,dummy,dummy); getKineObj(track1+1)->getParticle(aTrack2,aID2); //printf("A tr1=%i cr1=%i cr2=%i ID2=%i tr2=%i \n",track1,pTrack_par1,pTrack_par2,aID2,aTrack2); if(pTrack_par2==pTrack_par1) return (track1+1); getKineObj(track1-1)->getCreator(pTrack_par2,dummy,dummy); // printf("B tr1=%i cr1=%i cr2=%i \n",track1,pTrack_par1,pTrack_par2); if(pTrack_par2==pTrack_par1) return (track1-1); /* getKineObj(track1+2)->getCreator(pTrack_par2,dummy,dummy); getKineObj(track1+2)->getParticle(aTrack2,aID2); printf("C tr1=%i cr1=%i cr2=%i ID2=%i tr2=%i \n",track1,pTrack_par1,pTrack_par2,aID2,aTrack2); if(pTrack_par2==pTrack_par1) return (track1+2); getKineObj(track1+3)->getCreator(pTrack_par2,dummy,dummy); getKineObj(track1+3)->getParticle(aTrack2,aID2); printf("D tr1=%i cr1=%i cr2=%i ID2=%i tr2=%i \n",track1,pTrack_par1,pTrack_par2,aID2,aTrack2); if(pTrack_par2==pTrack_par1) return (track1+3); Int_t tr=track1; Int_t i=0; Int_t pMech2=-2; while (tr>=0) { i=i+1; if(getKineObj(i)){ getKineObj(i)->getCreator(pTrack_par2,dummy,pMech2); getKineObj(i)->getParticle(aTrack2,aID2); printf("%i cr1=%i cr2=%i ID2=%i tr2=%i pMech2=%i \n",i,pTrack_par1,pTrack_par2,aID2,aTrack2,pMech2); if(pTrack_par2==pTrack_par1 && aTrack2!=track1) return (i); } else { tr=-1;printf("LAST track=%i\n",i-1); } } */ } //printf("return -1 for track1=%i \n",track1); return -1; } HGeantKine* HMdcFunc1::getPairTrack(HGeantKine* kine1){ //this method finds a pair lepton to lepton with track1 //I use the fact that it is either track1+1 or track1+1 Int_t pTrack_par1=-1,pTrack_par2=-2,dummy=0,track1=-1; if(kine1){ kine1->getParticle(track1,dummy); kine1->getCreator(pTrack_par1,dummy,dummy); if(getKineObj(track1+1)) getKineObj(track1+1)->getCreator(pTrack_par2,dummy,dummy); if(pTrack_par2==pTrack_par1) return getKineObj(track1+1); if(getKineObj(track1-1)) getKineObj(track1-1)->getCreator(pTrack_par2,dummy,dummy); if(pTrack_par2==pTrack_par1) return getKineObj(track1-1); } return NULL; } Float_t HMdcFunc1::getPairOpenAngle(HGeantKine* kine1,HGeantKine* kine2){ //input 2 kine objects //output opening angle of trajectories Double_t rad2deg = TMath::RadToDeg(); HGeomVector vec1; if (kine1){ Float_t xMom1,yMom1,zMom1; kine1->getMomentum(xMom1,yMom1,zMom1); vec1.setX(xMom1); vec1.setY(yMom1); vec1.setZ(zMom1); vec1/=vec1.length();//norm } HGeomVector vec2; if (kine2){ Float_t xMom2,yMom2,zMom2; kine2->getMomentum(xMom2,yMom2,zMom2); vec2.setX(xMom2); vec2.setY(yMom2); vec2.setZ(zMom2); vec2/=vec2.length();//norm } Float_t cosfOpeningAngle = vec1.scalarProduct(vec2); //cout<getMomentum(xMom1,yMom1,zMom1); kine2->getMomentum(xMom2,yMom2,zMom2); TVector3 vvec1(xMom1,yMom1,zMom1); TVector3 vvec2(xMom2,yMom2,zMom2); fOpeningAngle=vvec1.Angle(vvec2)* rad2deg; return fOpeningAngle; } Float_t HMdcFunc1::getPairInvMass(HGeantKine* kine1,HGeantKine* kine2){ //copy from HPidFL if(kine1&&kine2){ Float_t opang = TMath::DegToRad()*getPairOpenAngle(kine1,kine2); Float_t p1 = kine1->getTotalMomentum(); Float_t p2 = kine2->getTotalMomentum(); return 2.*sin(opang/2.)*sqrt(p1*p2); } return -1; } void HMdcFunc1::putInArray(Int_t cls,TArrayI * arr,Int_t index){ //filling the array arr if (indexGetSize()){ arr->AddAt(cls,index); } else { arr->Set(2*arr->GetSize()); arr->AddAt(cls,index); } } Bool_t HMdcFunc1::isInArray(Int_t cls,TArrayI * arr){ //checking if number cls is already in arr for (Int_t i=0;iGetSize();i++){ if(arr->At(i)==cls) return kTRUE; } return kFALSE; } void HMdcFunc1::resetArray(TArrayI *arr) { //put all entries to -1 for (Int_t i=0;iGetSize();i++) arr->AddAt(-1,i); } Bool_t HMdcFunc1::isInGeantMdc(Int_t tr){ //check whether track tr is in GeantMdc inner modules HGeantMdc *pObj_geantmdc=NULL; if(getGeantMdcCat()){ HIterator *iterator_mdcgeant= (HIterator*)getGeantMdcCat()->MakeIterator(); iterator_mdcgeant->Reset(); while ((pObj_geantmdc=(HGeantMdc*)iterator_mdcgeant->Next())!= 0) { if(pObj_geantmdc->getModule()==2||pObj_geantmdc->getModule()==3) continue; if(pObj_geantmdc->getTrack()==tr){ return kTRUE; } } } else { printf("To Call HMdcFunc1::isInGeantMdc(Int tr)one has to set HMdcFunc1::GeomMdc category\n"); } return kFALSE; } Bool_t HMdcFunc1::isInGeantRichMirror(Int_t tr){ //check whether track tr is in GeantMdc inner modules HGeantRichMirror *pObj_geantrichmirr=NULL; Int_t aTrk,aID; if(getGeantRichMirrorCat()){ HIterator *iterator_geant= (HIterator*)getGeantRichMirrorCat()->MakeIterator(); iterator_geant->Reset(); while ((pObj_geantrichmirr=(HGeantRichMirror*)iterator_geant->Next())!= 0) { pObj_geantrichmirr->getTrack(aTrk,aID); if(aTrk==tr){ return kTRUE; } } } else { printf("To Call HMdcFunc1::isInGeantRichMirr(Int tr)one has to set HMdcFunc1::GeomRichMirr category\n"); } return kFALSE; } Bool_t HMdcFunc1::isInGeantTof(Int_t tr){ //check whether track tr is in GeantTof HGeantTof *pObj_geant=NULL; if(getGeantTofCat()){ HIterator *iterator_geant= (HIterator*)getGeantTofCat()->MakeIterator(); iterator_geant->Reset(); while ((pObj_geant=(HGeantTof*)iterator_geant->Next())!= 0) { if(pObj_geant->getModule()>21) continue; //this would be TOFINO if(pObj_geant->getTrack()==tr){ return kTRUE; } } } else { printf("To Call HMdcFunc1::isInGeantTof(Int tr)one has to set HMdcFunc1::GeomTof category\n"); } return kFALSE; } Bool_t HMdcFunc1::isInGeantShower(Int_t tr){ //check whether track tr is in GeantShower HGeantShower *pObj_geant=NULL; if(getGeantShowCat()){ HIterator *iterator_geant= (HIterator*)getGeantShowCat()->MakeIterator(); iterator_geant->Reset(); while ((pObj_geant=(HGeantShower*)iterator_geant->Next())!= 0) { if(pObj_geant->getTrack()==tr){ return kTRUE; } } } else { printf("To Call HMdcFunc1::isInGeantShower(Int tr)one has to set HMdcFunc1::GeomShower category\n"); } return kFALSE; } Bool_t HMdcFunc1::isHadesAccepted(Int_t tr){ //look whether this track is seen in rich-mdc tof||show if(isInGeantRichMirror(tr)) { if(isInGeantMdc(tr)) { if(isInGeantTof(tr)) { return kTRUE; } if(isInGeantShower(tr)) { return kTRUE; } } } return kFALSE; } Bool_t HMdcFunc1::couldBeAccepted(Int_t tr){ if(tr>=0) return couldBeAccepted(getKineObj(tr)); return kFALSE; } Bool_t HMdcFunc1::couldBeAccepted(HGeantKine *kine){ //check whether this particle was created before rich mirror that it could be possible flying to HADES //acceptance and should be taken in an account for acceptance studies Float_t verX,verY,verZ; Int_t material,dum; if(kine){ kine->getVertex(verX,verY,verZ); //mirror max radius is 715 mm if(sqrt(verX*verX+verY*verY+verZ*verZ)>715) return kFALSE; kine->getCreator(dum,material,dum); if(material==3||material==4||material==8||material==9||material==10||material==17||material==19) return kTRUE; //materials Al,Steel,C4F10,CaF2,CKFshell,RichMirr,CarbonTarget could change in various gen of sim data // this is valid for gen5 Nov01 } return kFALSE; } Bool_t HMdcFunc1::passCprPdfCuts(HPidDilepton* pDil,Float_t prob,Char_t * name_of_pdf_file,Int_t mode_rs,Option_t *rej_condition){ //checking wheather both PidParticles fullfill CPR propabilities to be single if(passCprPdfCuts((HPidParticle *)pDil->getParticle1(),prob,name_of_pdf_file,mode_rs,rej_condition)==kFALSE) return kFALSE; if(passCprPdfCuts((HPidParticle *)pDil->getParticle2(),prob,name_of_pdf_file,mode_rs,rej_condition)==kFALSE) return kFALSE; return kTRUE; } Bool_t HMdcFunc1::passCprPdfCuts(HPidParticle* pPart,Float_t prob,Char_t * name_of_pdf_file,Int_t mode_rs,Option_t *rej_condition){ // checking wheather HPidParticle pass a probability to be a single // mode_rs =0 for simulation data // mode_rs =1 for real data TString condition = rej_condition; condition.ToUpper(); HPidCPParam par(name_of_pdf_file); Int_t sec =-1,nwires_part_0=-1,nwires_part_1=-1,level_0=-1,level_1=-1,the_b=-1,phi_b=-1,level_b0=-1,level_b1=-1,cls_part_0=-1,cls_part_1=-1,charge_part=-1; Float_t mdc_the=-1.0, mdc_phi=-1.0; sec=(getPidTrackCand(pPart))->getMdcSeg()->getSec(); //if(sec==0) {printf("sec==%i\n",sec);return kFALSE;} mdc_the=(getPidTrackCand(pPart))->getMdcSeg()->getTheta(); mdc_the=getNormalMdcTheta(mdc_the); the_b=calculateBin(mdc_the,0); mdc_phi=(getPidTrackCand(pPart))->getMdcSeg()->getPhi(); mdc_phi=getNormalMdcPhi(sec,mdc_phi); phi_b=calculateBin(mdc_phi,1); level_0=getMdcLevelCls(pPart,0); level_b0=calculateLevelBin(level_0); level_1=getMdcLevelCls(pPart,1); level_b1=calculateLevelBin(level_1); cls_part_0=getMdcClsSize(pPart,0); cls_part_1=getMdcClsSize(pPart,1); nwires_part_0=getMdcNWires(pPart,0); nwires_part_1=getMdcNWires(pPart,1); charge_part=getIntCharge(pPart); if(mode_rs==0){ //this is real data cls_part_0=corect_cls(cls_part_0,the_b,phi_b,level_b0); cls_part_1=corect_cls(cls_part_1,the_b,phi_b,level_b1); nwires_part_0=corect_nwires(nwires_part_0,the_b,phi_b,level_b0); nwires_part_1=corect_nwires(nwires_part_0,the_b,phi_b,level_b1); } // Float_t prop_single=par.getProbSingle(the_b,phi_b,level_b0,level_b1,cls_part_0,cls_part_1,nwires_part_0,nwires_part_1); Float_t prop_single0=par.getProbSingle(the_b,phi_b,level_b0,cls_part_0,nwires_part_0,0); Float_t prop_single1=par.getProbSingle(the_b,phi_b,level_b1,cls_part_1,nwires_part_1,1); //in a case having -1 made default if(strstr(condition.Data(),"OR")){ // cout<<"current condition "<prob||prop_single1>prob) return kTRUE; } if(strstr(condition.Data(),"AND")){ // cout<<"current condition "<prob&&prop_single1>prob) return kTRUE; } if(strstr(condition.Data(),"AND")==NULL && strstr(condition.Data(),"OR")==NULL ){ printf("THIS condition: %s is unknown possible: AND and OR",condition.Data()); } return kFALSE; } Int_t HMdcFunc1::corect_cls(Int_t cls,Int_t the_bin,Int_t phi,Int_t level){ //corect cls to keep sim/real correspondance //function the,phi,level //for the moment only the Float_t k[8]={1.18,1.17,1.10,1.10,1.13,1.12,1.15,1.10}; if( the_bin>=0&&the_bin<8){ cls=(Int_t) (cls*k[the_bin]); } return cls; } Int_t HMdcFunc1::corect_nwires(Int_t nwires,Int_t the_bin,Int_t phi,Int_t level){ Float_t k[8]={1.06,1.04,1.03,1.04,1.05,1.04,1.04,0.99}; if( the_bin>=0&&the_bin<8){ nwires=(Int_t) (nwires*k[the_bin]); } return nwires; } Int_t HMdcFunc1::whatCommonHits(HPidDilepton* dil){ //this check common hits of leptons in pairs //return values: //0 : no common hits //1 : R- common RICH hit //2 : M- common MDC hit //3 : RM- common RICH and MDC hits //4 : Me- common META hit //5 : RMe-common RICH and META hits //6 : MMe-common MDC and META hits //7 : should not happend //-1: default return HPidParticle *pPart1=NULL; HPidParticle *pPart2=NULL; Int_t sum=-1; if(dil){ pPart1=(HPidParticle *)dil->getParticle1(); pPart2=(HPidParticle *)dil->getParticle2(); // pPart1->print(); // pPart2->print(); //dil->print(); if(pPart1&& pPart2){ HPidTrackCand* pTrackCand1=pPart1->getTrackCand(); HPidTrackCand* pTrackCand2=pPart2->getTrackCand(); // pTrackCand1->print(); // pTrackCand2->print(); // printf("info dilepton hits:\n"); // printf("sectors:%i %i \n",pTrackCand1->getSector(),pTrackCand2->getSector()); // printf("RICH hits:%i %i \n",pTrackCand1->getRingId(),pTrackCand2->getRingId()); // printf("MDC hits:%i %i \n",pTrackCand1->getInnerMdcId(0),pTrackCand2->getInnerMdcId(0)); // printf("META hits:%i %i \n",pTrackCand1->getMetaId(),pTrackCand2->getMetaId()); Bool_t common_RICH=kFALSE; Bool_t common_MDC=kFALSE; Bool_t common_META=kFALSE; //check ring ID if(pTrackCand1&& pTrackCand2){ if(pTrackCand1->getRingId()==pTrackCand2->getRingId()){ //printf ("in dilepton common RICH ind1=%i ind2=%i\n",pTrackCand1->getRingId(),pTrackCand2->getRingId()); common_RICH=kTRUE; } //check MdcSeg Id if((pTrackCand1->getInnerMdcId(0)==pTrackCand2->getInnerMdcId(0))&&(pTrackCand1->getSector()==pTrackCand2->getSector())){ //printf ("in dilepton common MDC ind1=%i ind2=%i\n",pTrackCand1->getInnerMdcId(0),pTrackCand2->getInnerMdcId(0)); common_MDC=kTRUE; } //check META Id if(pTrackCand1->getMetaId()==pTrackCand2->getMetaId()){ // printf ("in dilepton common META ind1=%i ind2=%i\n",pTrackCand1->getMetaId(),pTrackCand2->getMetaId()); common_META=kTRUE; } sum=0; if(common_RICH) sum=sum+1; if(common_MDC) sum=sum+2; if(common_META) sum=sum+4; } } } return sum; } ClassImp(HMdcFunc1)