#include "hmdcfunc1.h" #include "hpiddilepton.h" #include "hpidparticle.h" #include "hpidparticlesim.h" #include "hpidtrackcand.h" #include "hpidcandidate.h" #include "hpiddilepton.h" #include "hpidfl.h" #include "hmdcseg.h" #include "hmdcsegsim.h" #include "hmdchit.h" #include "hmdchitsim.h" #include "hmdcclusinf.h" #include "hmdcclusinfsim.h" #include "hmdcdef.h" #include "hmdctrackddef.h" #include "richdef.h" #include "hgeantdef.h" #include "hrichhit.h" #include "hrichhitsim.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" Bool_t HMdcFunc1::initCategoryPointers() { // Init the the category pointers for the given event in Hades. // Needs to called before one can retrieve objects fMdcSegCat=gHades->getCurrentEvent()->getCategory(catMdcSeg); if (!fMdcSegCat) { ::Error("initCategories()","There is no catMdcSeg !"); return kFALSE; } fRichHitCat=gHades->getCurrentEvent()->getCategory(catRichHit); if (!fRichHitCat) { ::Error("initCategories()","There is no catRichHit !"); return kFALSE; } fMdcHitCat=gHades->getCurrentEvent()->getCategory(catMdcHit); if (!fMdcHitCat) { ::Error("initCategories()","There is no catMdcHit !"); return kFALSE; } fClusInf=gHades->getCurrentEvent()->getCategory(catMdcClusInf); if(!fClusInf) { ::Error("initCategories()","There is no catMdcClusInf !"); return kFALSE; } fGeantKineCat=gHades->getCurrentEvent()->getCategory(catGeantKine); if (!fGeantKineCat) { ::Error("initCategories()","There is no catGeantKine !"); } if(fGeantKineCat) { // simulation or embedding ! fGeantMdcCat=gHades->getCurrentEvent()->getCategory(catMdcGeantRaw); if (!fGeantMdcCat) { ::Error("initCategories()","There is no catMdcGeantRaw !"); return kFALSE; } fGeantTofCat=gHades->getCurrentEvent()->getCategory(catTofGeantRaw); if (!fGeantTofCat) { ::Error("initCategories()","There is no catTofGeantRaw !"); return kFALSE; } fGeantShowCat=gHades->getCurrentEvent()->getCategory(catShowerGeantRaw); if (!fGeantShowCat) { ::Error("initCategories()","There is no catShowerGeantRaw !"); return kFALSE; } fGeantRichMirrorCat=gHades->getCurrentEvent()->getCategory(catRichGeantRaw+2); if (!fGeantRichMirrorCat) { ::Error("initCategories()","There is no catRichGeantRaw+2 !"); return kFALSE; } } return kTRUE; } HMdcSeg* HMdcFunc1::getMdcSegFromPidTrackCand(HPidTrackCand* cand) { // Retrieve the Pointer of the inner segment corresponding // to the HPidTrackCand object. Return NULL if something wrong if(cand) { HPidHitData* hitdata = cand->getHitData(); if(hitdata) { HMdcSeg* seg=(HMdcSeg*)fMdcSegCat->getObject(hitdata->getIndInnerSeg()); if(seg){ return seg; } else { ::Error("getMdcSegFromPidTrackCand()","Zero pointer for MdcSeg retrieved!"); return 0; } } else { ::Error("getMdcSegFromPidTrackCand()","Zero pointer for PidHitData retrieved!"); return 0; } } else { ::Error("getMdcSegFromPidTrackCand()","Zero pointer for PidTrackCand retrieved!"); return 0; } return 0; } HRichHit* HMdcFunc1::getRichHitFromPidTrackCand(HPidTrackCand* cand) { // Retrieve the Pointer of the RichHit corresponding // to the HPidTrackCand object. Return NULL if something wrong if(cand) { HPidHitData* hitdata = cand->getHitData(); if(hitdata) { HRichHit* rich=(HRichHit*)fMdcSegCat->getObject(hitdata->getIndRICH()); if(rich){ return rich; } else { ::Error("getRichHitFromPidTrackCand()","Zero pointer for RichHit retrieved!"); return 0; } } else { ::Error("getRichHitFromPidTrackCand()","Zero pointer for PidHitData retrieved!"); return 0; } } else { ::Error("getRichHitFromPidTrackCand()","Zero pointer for PidTrackCand retrieved!"); return 0; } return 0; } Int_t HMdcFunc1::getMdcClsSize(HPidParticle* part,Int_t module) { // return the clustersize from HMdcClusInf for the // inner segment of a HPidParticle. // module == 0 -> first module in segment // module == 1 -> second module in segment // return 0 if something wrong HMdcSeg* pMdcSeg =NULL; HPidTrackCand* pPidTrackCand=NULL; Int_t cls =0; if(part) { pPidTrackCand=getPidTrackCand(part); if(pPidTrackCand) { pMdcSeg=getMdcSegFromPidTrackCand(pPidTrackCand); if(pMdcSeg) { cls=getMdcClsSize(pMdcSeg,module); } } } return cls; } Int_t HMdcFunc1::getMdcNWires(HPidParticle* part,Int_t module) { // return the number of wires from HMdcClusInf for the // inner segment of a HPidParticle. // module == 0 -> first module in segment // module == 1 -> second module in segment // return 0 if something wrong HMdcSeg* pMdcSeg =NULL; HPidTrackCand* pPidTrackCand=NULL; Int_t nwires =0; if(part) { pPidTrackCand=getPidTrackCand(part); if(pPidTrackCand) { pMdcSeg=getMdcSegFromPidTrackCand(pPidTrackCand); if(pMdcSeg) { nwires=getMdcNWires(pMdcSeg,module); } } } return nwires; } Int_t HMdcFunc1::getMdcLevelCls(HPidParticle* part,Int_t module) { // return the level of cluster finding from HMdcClusInf for the // inner segment of a HPidParticle. // module == 0 -> first module in segment // module == 1 -> second module in segment // return 0 if something wrong HMdcSeg* pMdcSeg =NULL; HPidTrackCand* pPidTrackCand=NULL; Int_t level =0; if(part) { pPidTrackCand=getPidTrackCand(part); if(pPidTrackCand) { pMdcSeg=getMdcSegFromPidTrackCand(pPidTrackCand); if(pMdcSeg) { level=getMdcLevelCls(pMdcSeg,module); } } } return level; } Int_t HMdcFunc1::getIntCharge(HPidParticle* part) { // return the integrated charge from HRichHit for the // of a HPidParticle. // return 0 if something wrong HRichHit* pRichHit =NULL; HPidTrackCand* pPidTrackCand=NULL; Int_t int_charge =0; if(part) { pPidTrackCand=getPidTrackCand(part); if(pPidTrackCand) { pRichHit=getRichHitFromPidTrackCand(pPidTrackCand); if(pRichHit) { int_charge=pRichHit->getRingAmplitude(); } } } return int_charge; } Int_t HMdcFunc1::getMdcClsSize(HMdcSeg* seg,Int_t module) { // return the clustersize from HMdcClusInf for the // segment. // module == 0 -> first module in segment // module == 1 -> second module in segment // return 0 if something wrong Int_t cls=0; HMdcHit* hit=0; if(seg) { hit =getMdcHit(seg,module); if(hit) cls=getMdcClsSize(hit); } return cls; } Int_t HMdcFunc1::getMdcNWires(HMdcSeg* seg,Int_t module) { // return the number of wires from HMdcClusInf for the // segment. // module == 0 -> first module in segment // module == 1 -> second module in segment // return 0 if something wrong Int_t nwires=0; HMdcHit* hit=0; if(seg) { hit =getMdcHit(seg,module); if(hit) nwires=getMdcNWires(hit); } return nwires; } Int_t HMdcFunc1::getMdcLevelCls(HMdcSeg* seg,Int_t module) { // return the level of cluster finding from HMdcClusInf for the // segment. // module == 0 -> first module in segment // module == 1 -> second module in segment // return 0 if something wrong Int_t level=0; HMdcHit* hit=0; if(seg) { hit =getMdcHit(seg,module); if(hit) level=getMdcLevelCls(hit); } return level; } HPidTrackCand* HMdcFunc1::getPidTrackCand(HPidParticle* part) { // return HPidTrackCand form an HPidParticle // using the indices HPidParticel->HPidCandidate->HPidTrackCand // return 0 if something wrong HPidCandidate* pidcand =NULL; HPidTrackCand* pidtrackcand=NULL; if(part) { pidcand=part->getPidCandidate(); if(pidcand) { pidtrackcand=pidcand->getTrackCandidate(); } } return pidtrackcand; } HMdcHit* HMdcFunc1::getMdcHit(HMdcSeg* seg,Int_t module) { // return the HMfdcHit for the // segment. // module == 0 -> first module in segment // module == 1 -> second module in segment // return 0 if something wrong HMdcHit* hit=NULL; Int_t hit_inx=-1; if(seg) { hit_inx=seg->getHitInd(module); if(hit_inx>-1) { if(getMdcHitCat()) { hit=(HMdcHit *)getMdcHitCat()->getObject(hit_inx); return hit; } else { ::Error("getMdcHit()","******MdcHitCat is NULL********use HMdcFunc1::setMdcHitCat()"); } } } return hit; } Int_t HMdcFunc1::getMdcClsSize(HMdcHit* hit) { // return the clustersize from HMdcClusInf for the // HMdcHit. // return -1 if something wrong Int_t clssize =-1; Int_t clsindex =-1; HMdcClusInf* pObj_clusinf=NULL; clsindex =hit->getClusInfIndex(); pObj_clusinf=getClusInfObj(clsindex); if(pObj_clusinf) { clssize=pObj_clusinf->getClusSize(); } return clssize; } Int_t HMdcFunc1::getMdcNWires(HMdcHit* hit) { // return the number of wires from HMdcClusInf for the // HMdcHit. // return -1 if something wrong 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) { // return the level of cluster finding from HMdcClusInf for the // HMdcHit. // return -1 if something wrong 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 { ::Error("getClusInfObj()","******MdcClusInfCat is NULL********use HMdcFunc1::setMdcClusInfCat()"); } return pObj_clusinf; } Int_t HMdcFunc1::calculateLevelBin(Int_t level) { // calculates the bin of the the level parameter // access the corrections 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; return mdc_theta_bin; } else { ::Error("calculateBin()","theta out of range"); 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 { ::Error("calculateBin()","phi out of range"); return -2; } default: ::Error("calculateBin()","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=HPidFL::getGeantKineObjFromTrkNbFast(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=HPidFL::getGeantKineObjFromTrkNbFast(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=HPidFL::getGeantKineObjFromTrkNbFast(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=HPidFL::getGeantKineObjFromTrkNbFast(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) { ::Warning("getKickAngle()","This function return hard wired numbers!"); //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); } */ HGeantKine* HMdcFunc1::getKineObj(Int_t pTrack) { // return the HGenatKine object pointer // for a given track number. returns NULL // if something wrong if(fGeantKineCat!=NULL||pTrack<1) { return (HGeantKine*)(fGeantKineCat->getObject(pTrack-1)); } else { ::Error("getKineObj()","HMdcFunc1::getKineObj there is no fGeantKineCat defined please use HMdcFunc1::setGeantKineCat before\n"); return NULL; } } Int_t HMdcFunc1::getGeantCommonTrack(HPidParticleSim* part) { //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 HPidTrackCand* pPidTrackCand=NULL; HRichHitSim* pRichHitSim =NULL; if(part!=NULL) { const HPidGeantTrackSet* pidset=part->getGeantTrackSet(); pPidTrackCand=getPidTrackCand(part); if(pPidTrackCand!=NULL) { HGeantKine *kine =NULL; Int_t richTr1 = -1; Int_t richTr2 = -1; Int_t richTr3 = -1; Int_t iTrack = -1; Int_t iId = -1; Int_t flag_richmdcele= 0; Int_t flag_missinmdc = 0; pRichHitSim=(HRichHitSim*) getRichHitFromPidTrackCand(pPidTrackCand); if(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; Int_t system =(Int_t)part->getSystem(); for(UInt_t j=0;jgetNInnerMdcTracks();j++) { //go over mdc iId=-1; if(pidset->getInnerMdcTrack(j)>0) { Int_t nmetatracks=0; if(system==0){ nmetatracks=pidset->getNShowerTracks(); } else if (system==1){ nmetatracks=pidset->getNTOFTracks(); } else { ::Error("getGeantCommonTrack()","System is not 0 or 1!"); return -10; } for(Int_t k=0;kgetShowerTrack(k) : metatrack=pidset->getTOFTrack(k); if(pidset->getInnerMdcTrack(j)==metatrack) { //common MDC-META kine=HPidFL::getGeantKineObjFromTrkNbFast(pidset->getInnerMdcTrack(j)); kine->getParticle(iTrack, iId); if((pidset->getInnerMdcTrack(j)==richTr1)|| (pidset->getInnerMdcTrack(j)==richTr2)|| (pidset->getInnerMdcTrack(j)==richTr3)) { //common RING-MDC-META return pidset->getInnerMdcTrack(j); } if(iId==9||iId==8) { return -2; } if(iId==14) { return -3; } if(iId==3||iId==2) { return -5; } } if((pidset->getInnerMdcTrack(j)==richTr1)|| (pidset->getInnerMdcTrack(j)==richTr2)|| (pidset->getInnerMdcTrack(j)==richTr3)) { //common RING-MDC but not META because otherwise it is already return before kine=HPidFL::getGeantKineObjFromTrkNbFast(pidset->getInnerMdcTrack(j)); kine->getParticle(iTrack, iId); if(iId==2||iId==3){ //common RING-MDC electron (but not META) flag_richmdcele=1; } else { ::Error("getGeantCommonTrack()","Strange RING-MDC is not lepton"); } } //one has to treat separately case that in MDC I do not see a second pair from conv but in META yes if((metatrack==richTr1)|| (metatrack==richTr2)|| (metatrack==richTr3)) { for(UInt_t l=0;lgetNInnerMdcTracks();l++) { if(TMath::Abs(pidset->getInnerMdcTrack(l)-metatrack)==1) { kine=HPidFL::getGeantKineObjFromTrkNbFast(metatrack); kine->getParticle(iTrack, iId); //I do not return because it could still be double-double-single case // flag_missinmdc=(100+iId); flag_missinmdc=metatrack; } } } } } } if(flag_missinmdc!=0) { return flag_missinmdc; } if(flag_richmdcele==1){ return -4; } return -1; } } } return -10; } Bool_t HMdcFunc1::isCompt(HGeantKine* kine) { // returns kTRUE if kine is a lepton // originating from gama->e- gama Int_t pTrack_par =-1; Int_t pMed =-1; Int_t pMech =-1; Int_t pID_par = 0; Int_t pTrack_particle=-1; Int_t 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) { // returns kTRUE if kine is a lepton // originating from pi0 ->gama->e+e- Int_t pTrack_par =-1; Int_t pMed =-1; Int_t pMech =-1; Int_t pID_par = 0; Int_t pTrack_par_par =-1; Int_t pMed_par =-1; Int_t pMech_par =-1; Int_t pID_par_par =-1; Int_t pTrack_particle=-1; Int_t 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); if(pID_par_par==7){ return kTRUE ; } } } } } return kFALSE; } Bool_t HMdcFunc1::isConv(HGeantKine* kine) { // returns kTRUE if kine is a lepton // originating from gama->e+e- Int_t pTrack_par =-1; Int_t pMed =-1; Int_t pMech =-1; Int_t pID_par = 0; Int_t pTrack_particle=-1; Int_t 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){ return kTRUE ; } } } return kFALSE; } Bool_t HMdcFunc1::isPi0Dalitz(HGeantKine* kine) { // returns kTRUE if kine is a lepton // originating from pi0 Int_t pTrack_par =-1; Int_t pMed =-1; Int_t pMech =-1; Int_t pID_par = 0; Int_t pTrack_particle=-1; Int_t 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==7){ return kTRUE; } } } return kFALSE; } Bool_t HMdcFunc1::isEtaDalitz(HGeantKine* kine) { // returns kTRUE if kine is a lepton // originating from eta dalitz Int_t pTrack_par =-1; Int_t pMed =-1; Int_t pMech =-1; Int_t pID_par = 0; Int_t pTrack_particle=-1; Int_t 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) { // returns kTRUE if kine is a lepton Int_t pTrack_particle=-1; Int_t 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. // look only at inner segments :very important if(seg->getIOSeg()==0) { Int_t 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) { // returns kTRUE if hit has only one track // conribution and is a lepton Int_t track_mul = 0; 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; Int_t pID1 =-1; Int_t pTrack_par1=-1; Int_t pTrack_par2=-2; Int_t 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)))) { 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); } } 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; Int_t pID1 =-1; Int_t pTrack_par1=-1; Int_t pTrack_par2=-2; Int_t 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)))) { 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); } } 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; Float_t ygeant = 1000.0; Float_t dummy = 1000.0; Float_t x_hit =-1000.0; Float_t y_hit =-1000.0; Float_t Xdiff = 3*2.5; Float_t 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(); } 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); if(pTrack_par2==pTrack_par1) return (track1+1); getKineObj(track1-1)->getCreator(pTrack_par2,dummy,dummy); if(pTrack_par2==pTrack_par1) return (track1-1); } 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; Int_t pTrack_par2=-2; Int_t dummy = 0; Int_t 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); Float_t fOpeningAngle=0.; if (TMath::Abs(cosfOpeningAngle) <= 1) { fOpeningAngle=TMath::ACos(cosfOpeningAngle) * rad2deg; } Float_t xMom1,yMom1,zMom1; Float_t xMom2,yMom2,zMom2; kine1->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){ delete iterator_mdcgeant; return kTRUE; } } delete iterator_mdcgeant; } else { ::Error("isInGeantMdc()","To Call HMdcFunc1::isInGeantMdc(Int tr)one has to set HMdcFunc1::GeomMdc category"); } 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){ delete iterator_geant; return kTRUE; } } delete iterator_geant; } else { ::Error("isInGeantRichMirror()","To Call HMdcFunc1::isInGeantRichMirr(Int tr)one has to set HMdcFunc1::GeomRichMirr category"); } 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){ delete iterator_geant; return kTRUE; } } delete iterator_geant; } else { ::Error("isInGeantTof()","To Call HMdcFunc1::isInGeantTof(Int tr)one has to set HMdcFunc1::GeomTof category"); } 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){ delete iterator_geant; return kTRUE; } } delete iterator_geant; } else { ::Error("isInGeantShower()","To Call HMdcFunc1::isInGeantShower(Int tr)one has to set HMdcFunc1::GeomShower category"); } 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 ::Warning("couldBeAccepted()","problem: materials could change! check the use of the function!"); 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; Int_t nwires_part_0=-1; Int_t nwires_part_1=-1; Int_t level_0 =-1; Int_t level_1 =-1; Int_t the_b =-1; Int_t phi_b =-1; Int_t level_b0 =-1; Int_t level_b1 =-1; Int_t cls_part_0 =-1; Int_t cls_part_1 =-1; Int_t charge_part =-1; Float_t mdc_the =-1.0; Float_t mdc_phi =-1.0; HMdcSeg* seg=getMdcSegFromPidTrackCand(getPidTrackCand(pPart)); sec =seg->getSec(); mdc_the =seg->getTheta(); mdc_phi =seg->getPhi(); mdc_the =getNormalMdcTheta(mdc_the); mdc_phi =getNormalMdcPhi(sec,mdc_phi); the_b =calculateBin(mdc_the,0); phi_b =calculateBin(mdc_phi,1); level_0 =getMdcLevelCls(pPart,0); level_1 =getMdcLevelCls(pPart,1); level_b0 =calculateLevelBin(level_0); 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_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")) { if (prop_single0>prob||prop_single1>prob) return kTRUE; } if(strstr(condition.Data(),"AND")) { if (prop_single0>prob&&prop_single1>prob) return kTRUE; } if(strstr(condition.Data(),"AND")==NULL && strstr(condition.Data(),"OR")==NULL ) { ::Error("passCprPdfCuts()","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 ::Warning("corect_cls()","Hardwired corrections!!!!check!"); 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) { //corect nwires to keep sim/real correspondance //function the,phi,level //for the moment only the ::Warning("corect_nwires()","Hardwired corrections!!!!check!"); 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(); if(pPart1&& pPart2) { HPidTrackCand* pTrackCand1=pPart1->getTrackCand(); HPidTrackCand* pTrackCand2=pPart2->getTrackCand(); Bool_t common_RICH=kFALSE; Bool_t common_MDC=kFALSE; Bool_t common_META=kFALSE; //check ring ID if(pTrackCand1&& pTrackCand2) { HPidHitData* hitdata1=pTrackCand1->getHitData(); HPidHitData* hitdata2=pTrackCand2->getHitData(); if(hitdata1->getIndRICH()==hitdata1->getIndRICH()){ common_RICH=kTRUE; } //check MdcSeg Id if((hitdata1->getIndInnerSeg()==hitdata2->getIndInnerSeg())&& (hitdata1->getSector()==hitdata2->getSector())){ common_MDC=kTRUE; } //check META Id Int_t system1=hitdata1->getSystem(); Int_t system2=hitdata2->getSystem(); Int_t ind1=-1; Int_t ind2=-1; if(system1==0){ ind1=hitdata1->getIndShower(); }else{ ind1=hitdata1->getIndTOF(); } if(system2==0){ ind2=hitdata2->getIndShower(); }else{ ind2=hitdata2->getIndTOF(); } if(ind1==ind2){ 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)