// $Id: hrichcorrelatorsim.cc,v 1.20 2006/08/12 12:53:41 halo Exp $
// Last update by Thomas Eberl: 04/08/02 14:22:26
//
#include "hrichcorrelatorsim.h"
#include "hrichcorrelatorpar.h"
#include "hrichgeometrypar.h"
#include "hrichutilfunc.h"
#include "hgeantkine.h"
#include "hgeantrich.h"
#include "hgeantmdc.h"
#include "hruntimedb.h"
#include "hevent.h"
#include "hspectrometer.h"
#include "hdetector.h"
#include "hrichdetector.h"
#include "hmdcdetector.h"
#include "htofdetector.h"
#include "hshowerdetector.h"
#include "hshowertrack.h"
#include "hcategory.h"
#include "hiterator.h"
#include "hmatrixcatiter.h"
#include "hrichhit.h"
#include "hrichhitsim.h"
#include "hmdcseg.h"
#include "hmdcsegsim.h"
#include "htofhit.h"
#include "htofhitsim.h"
#include "hshowerhittof.h"
#include "tofdef.h"
#include "hmdcdef.h"
#include "showertofinodef.h"
#include "showerdef.h"
#include "hdebug.h"
#include "hades.h"
#include "richdef.h"
#include "hhitmatchsim.h"
#include "hhitmatchheadersim.h"
#include "hlinearcategory.h"
#include "kickdef.h"
#include <fstream>
#include "hgeantheader.h"
#include "hrecevent.h"
#include "hpartialevent.h"
//_HADES_CLASS_DESCRIPTION
///////////////////////////////////////////////////////////
// HRichCorrelatorSim
//
// executes HRichCorrelator and uses additional info from GEANT
// fills HHitMatchSim
///////////////////////////////////////////////////////////
ClassImp(HRichCorrelatorSim)
HRichCorrelatorSim::HRichCorrelatorSim(Text_t *name,Text_t *title) :
HRichCorrelator(name,title)
{
}
HRichCorrelatorSim::HRichCorrelatorSim(Text_t *name,Text_t *title, char* filename, Bool_t style) :
HRichCorrelator(name,title,filename,style)
{
pFileName = filename;
isComplex = style;
}
HRichCorrelatorSim::HRichCorrelatorSim()
{
}
HRichCorrelatorSim::~HRichCorrelatorSim(void)
{
}
Bool_t HRichCorrelatorSim::init() {
if (gHades) {
HEvent *event=gHades->getCurrentEvent();
HRuntimeDb *rtdb=gHades->getRuntimeDb();
HSpectrometer *spec=gHades->getSetup();
if (event && rtdb) {
HTofDetector *tof=(HTofDetector*)spec->getDetector("Tof");
HShowerDetector *shower = (HShowerDetector*)spec->getDetector("Shower");
HMdcDetector *mdc =(HMdcDetector*)(spec->getDetector("Mdc"));
HRichDetector *rich = (HRichDetector*)spec->getDetector("Rich");
// init paramters for this task
HRichCorrelatorPar *pCorrPar =
(HRichCorrelatorPar*)rtdb->getContainer("RichCorrelatorParameters");
if (!pCorrPar)
{
pCorrPar = new HRichCorrelatorPar;
rtdb->addContainer(pCorrPar);
}
setCorrelationPar(pCorrPar);
if (!pCorrPar) return kFALSE;
// HRichGeometryPar *pGeoPar =
// (HRichGeometryPar*)rtdb->getContainer("RichGeometryParameters");
// if (!pGeoPar)
// {
// pGeoPar = new HRichGeometryPar;
// rtdb->addContainer(pGeoPar);
// }
// setGeometryPar(pGeoPar);
// if (!pGeoPar) return kFALSE;
// HGEANT KINE INFO
fGeantKineCat =(HLinearCategory*) event->getCategory(catGeantKine);
if (!fGeantKineCat)
{
Error("HRichCorrelatorSim::init",
"no GEANT Kine category available");
}
iter_kine = (HIterator *)fGeantKineCat ->MakeIterator("native");
// HGEANT RICH MIRROR INFO
fGeantRichMirrorCat = (HMatrixCategory*)(gHades->getCurrentEvent()
->getCategory(catRichGeantRaw+2));
if (!fGeantRichMirrorCat)
{
Warning("HRichCorrelatorSim::init",
"no GEANT RICH MIRROR category available");
}
if (fGeantRichMirrorCat) iter_mirror = (HIterator *)fGeantRichMirrorCat->MakeIterator("native");
// HGEANT MDC RAW INFO
fGeantMdcCat = (HMatrixCategory*)(gHades->getCurrentEvent()
->getCategory(catMdcGeantRaw));
if (!fGeantMdcCat)
{
Warning("HRichCorrelatorSim::init",
"no GEANT MDC RAW category available");
}
if (fGeantMdcCat) iter_mdcgeant = (HIterator *)fGeantMdcCat->MakeIterator("native");
// HGEANT RICH RAW INFO
// fGeantRichCat =(HLinearCategory*) event
// ->getCategory(catRichGeantRaw);
// if (!fGeantRichCat)
// {
// Error("HRichCorrelatorSim::init",
// "no GEANT RICH RAW category available");
// }
// iter_geantrichraw = (HIterator *) fGeantRichCat
// ->MakeIterator("native");
if (rich) { //Has the user set up a RICH?
// RICH HIT container
fRichPID=gHades->getCurrentEvent()->getCategory(catRichHit);
if (!fRichPID) {
fRichPID=rich->buildMatrixCat("HRichHitSim",1);
//cout<<" i have built hrichhitsim "<<endl;
if (!fRichPID) return kFALSE;
else gHades->getCurrentEvent()
->addCategory(catRichHit, fRichPID, "Rich");
}
fRichIter = (HIterator*) fRichPID->MakeIterator();
if(!fRichIter) Error("init","no RichHitSim iter defined");
pRichHitFitCat=gHades->getCurrentEvent()->getCategory(catRichHitFit);
if (!pRichHitFitCat){
Warning("init()","no Rich hit fit category in input");
}
if (pRichHitFitCat) pRichHitFitIter =
(HIterator*) pRichHitFitCat->MakeIterator();
// //Setup output
pHitMatchCat=event->getCategory(catMatchHit);
if (!pHitMatchCat) {
pHitMatchCat=rich->buildLinearCat("HHitMatchSim");
// cout<<"i have built hhitmatchsim"<<endl;
if (!pHitMatchCat) {
Error("init","No HIT MATCH category defined");
return kFALSE;
}
else event->addCategory(catMatchHit, pHitMatchCat, "Rich");
}
pIterMatchHit = (HIterator*)getHitMatchCat()->MakeIterator("native");
pHitMatchHeaderCat=event->getCategory(catMatchHitHeader);
if (!pHitMatchHeaderCat) {
pHitMatchHeaderCat=rich->buildLinearCat("HHitMatchHeaderSim");
//cout<<"i have built hhitmatchheadersim"<<endl;
if (!pHitMatchHeaderCat) {
Error("init","No HIT MATCH HEADER SIM category defined");
return kFALSE;
}
else event->addCategory(catMatchHitHeader, pHitMatchHeaderCat
, "Rich");
}
pIterMatchHitHeader = (HIterator*)getHitMatchHeaderCat()
->MakeIterator("native");
}
if (mdc) { //Has the user set up a MDC?
fMdcSeg = gHades->getCurrentEvent()->getCategory(catMdcSeg);
if (!fMdcSeg){
fMdcSeg = mdc->buildMatrixCategory("HMdcSegSim",0.5F);
//cout<<"i have built a hmdcsegsim"<<endl;
if (!fMdcSeg) {
Error("init","No MDC segment sim category defined");
return kFALSE;
}
}
else {
gHades->getCurrentEvent()
->addCategory(catMdcSeg,fMdcSeg,"Mdc");
fMdcSegIter=(HIterator *)fMdcSeg->MakeIterator();
}
}
if (tof) { //Has the user set up a TOF?
fTofHits=event->getCategory(catTofHit);
if (!fTofHits) {
fTofHits=tof->buildMatrixCategory("HTofHitSim",0.5F);
//cout<<"i have built ftofhits"<<endl;
if (!fTofHits) {
Error("init","No TOF input");
return kFALSE;
}else{
gHades->getCurrentEvent()
->addCategory(catTofHit,fTofHits,"Tof");
}
}
fTofIter=(HIterator *)fTofHits->MakeIterator();
}
if (shower) { //Has the user set up a Shower?
// fShowerHits=event->getCategory(catShowerHit);
// if (!fShowerHits) {
// fShowerHits = shower->buildCategory(catShowerHitTof);
// if (!fShowerHits) {
// Error("init","No Shower input available");
// return kFALSE;
// }
// }
// fShowerIter = (HIterator *)fShowerHits->MakeIterator();
// fShowerTofHits=event->getCategory(catShowerHitTof);
// if (!fShowerTofHits) {
// fShowerTofHits = shower->buildCategory(catShowerHitTof);
// if (!fShowerTofHits) {
// Error("init","No Shower input available");
// return kFALSE;
// }
// }
// fShowerTofIter = (HIterator *)fShowerTofHits->MakeIterator();
fShowerTrack=gHades->getCurrentEvent()->getCategory(catShowerTrack);
if (!fShowerTrack) {
fShowerTrack=shower->buildCategory(catShowerTrack);
//cout<<"i have built hshowertrack "<<endl;
if (!fShowerTrack) return kFALSE;
else gHades->getCurrentEvent()
->addCategory(catShowerTrack, fShowerTrack, "Shower");
}
iter_showertrack = (HIterator*)fShowerTrack->MakeIterator();
fShowerTofHits=event->getCategory(catShowerHitTofTrack);
if (!fShowerTofHits)
{
fShowerTofHits = new HLinearCategory("HShowerHitTofTrack", 1000);
//cout<<"i have built hshowerhittoftrack"<<endl;
if (!fShowerTofHits) return kFALSE;
else gHades->getCurrentEvent()
->addCategory(catShowerHitTofTrack,fShowerTofHits, "Tofino");
}
fShowerTofIter = (HIterator *)fShowerTofHits->MakeIterator();
}
// --------- Tracks from Kickplane --------------
fKickTrackCat = gHades->getCurrentEvent()->getCategory(catKickTrack);
if (!fKickTrackCat) {
Error("init","KickTrack not found in input file");
return kFALSE;
}
iterTracks = (HIterator *)fKickTrackCat->MakeIterator("native");
//cout<<"itertracks :"<<iterTracks<<endl;
} else {
Error ("init","! (event && rtdb)");
return kFALSE; //! (event && rtdb)
}
iniCounters();
pFileOut = new TFile(pFileName,"RECREATE");
//cout<<"end of init()"<<endl;
return kTRUE;
} else {
Error ("init","! (gHades)");
return kFALSE; //! (gHades)
}
}
void HRichCorrelatorSim::iniCounters(void)
{
HRichCorrelator::iniCounters();
// cout<<"in iniCountersSim"<<endl;
ng_Part=ngGNb_RM=ngGNb_RT_S=ngGNb_MT_S=ngGNb_RMT_S=0;
ngGNbRMonly=ngGNbRTSonly=ngGNbMTSonly=ngGNbRMTSonly=0;
ngNb_RW=ngNb_MW=0;
fg_RW=fg_MW=0.;
ngNb_Fakes=0;
ngNb_ConfTracklet=0;
for(Int_t k=0;k<MAXPARIDS;k++) ng_ParId[k]=0;
}
Bool_t HRichCorrelatorSim::reinit()
{
// cout<<"in reinit function"<<endl;
HRichCorrelatorPar* pCorrPar = (HRichCorrelatorPar*)getCorrelationPar();
// cout<<pCorrPar<<" is address"<<endl;
// cout<<pCorrPar->hasChanged()<<" has changed"<<endl;
// HRichGeometryPar* pGeo = (HRichGeometryPar*)getGeometryPar();
// if (pCorrPar->hasChanged())
//{
iniCuts();
iniSwitches();
//iniControlHistos();
// histArray = new TObjArray(5);
//histArray->Add(new TH1D("corr_obj_multi","corr_obj_multi",20,0,20));
pCorrPar -> printParam();
//}
return kTRUE;
}
void HRichCorrelatorSim::diag(HHitMatchSim * pHM)
{
// function that prints the found GEANT track numbers
// for different detector hits in a HHitMatchSim obj
cout<<"*** HRichCorrelatorSim::diag *********************** "<<endl;
HMdcSegSim *m=0;
HRichHitSim *r=0;
HShowerHitTof *st=0;
HTofHitSim *t=0;
// get objs from categories according to their index
// which is stored in HHitMatchSim
if(pHM->getMdcInd()>-1) m = ((HMdcSegSim*)getMdcSegCat()
->getObject(pHM->getMdcInd()));
if(pHM->getRichInd()>-1) r = ((HRichHitSim*)getRichHitCat()
->getObject(pHM->getRichInd()));
if(pHM->getShowInd()>-1) st = ((HShowerHitTof*)
getShowerTofHitCat()
->getObject(pHM->getShowInd()));
if(pHM->getTofInd()>-1) t = ((HTofHitSim*)getTofHitCat()
->getObject(pHM->getTofInd()));
if(r)
{
cout<<"RICH track/weight"<<endl;
cout<<r->track1<<"/"<<r->weigTrack1<<" ";
cout<<r->track2<<"/"<<r->weigTrack2<<" ";
cout<<r->track3<<"/"<<r->weigTrack3<<endl;
}
if(m)
{
Int_t nNbMdcSegTr=m->getNTracks();
cout<<"MDC track/weight"<<endl;
for (Int_t i=0;i<nNbMdcSegTr;i++)
{
cout<<m->getTrack(i)<<"/"<<(Int_t)m->getNTimes(i)<<" ";
}
cout<<endl;
}
if(st)
{
cout<<"SHOWER track/"<<endl;
for (Int_t i=0;i<SHOWERMAXTRACKS;i++)
{
cout<<showerTracks[i]<<"/"<<" ";
}
cout<<endl;
}
if(t)
{
cout<<"TOF track/"<<endl;
cout<<t->getNTrack1()<<"/"<<" ";
cout<<t->getNTrack2()<<"/"<<" ";
cout<<endl;
}
}
Int_t HRichCorrelatorSim::execute()
{
// cout<<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "
// <<nCounterProcessedNbEvents<<endl;
if (!HRichCorrelator::execute()) //treat sim data like exp data
//execute correlation search w/o
//GEANT info first
{// now add info from GEANT to the found corr objs
HHitMatchSim *pHitMatch=0;
pIterMatchHit->Reset();
// loop over corr objs in category
while((pHitMatch = (HHitMatchSim *)pIterMatchHit->Next()))
{
//get pointer to track info obj inside HHitMatchSim
//this obj stores all additional info from GEANT
(pHitMatch->getTrackInfoObj())->reset();
cleanTracks();//clean track helper arrays
// store the GEANT track numbers of the detector hits
// in temp helper arrays
if(pHitMatch->getRichInd()>=0) assignRichTracks(pHitMatch->getRichInd());
if(pHitMatch->getMdcInd()>=0) assignMdcTracks(pHitMatch->getMdcInd());
if(pHitMatch->getShowInd()>=0) assignShowerTrack(pHitMatch->getShowInd());
if(pHitMatch->getTofInd()>=0) assignTofTracks(pHitMatch->getTofInd());
// try to match the GEANT track numbers of the different detector hits
matchTracks(pHitMatch);
if ((richTracks[0] == -5 || richTracks[0] == -1) && richTracks[1] >0)
{
//cout<<"ALERT, ALERT ******************************************"
// <<nCounterProcessedNbEvents<<endl;
//diag(pHitMatch);
//pHitMatch->dumpToStdoutSim();
//cout<<"////////////////////////////////////////////////////////////////////////"<<endl;
}
}
// now fill one header object per event
// check for fakes
HHitMatchHeaderSim *pHMH = 0;
pIterMatchHitHeader->Reset();
//Int_t nnn=0;
while((pHMH = (HHitMatchHeaderSim *)pIterMatchHitHeader->Next()))
{
pHMH->resetSim();
//cout<<"evt nb "<<nCounterProcessedNbEvents<<" pHMH # "<<nnn<<endl;
//nnn++;
HHitMatchSim *pHitMatch=0;
pIterMatchHit->Reset();
Int_t nMaxParticles=pHMH->getNbCorrObjs()*MAXPARTICLES;
Int_t *nDParticles=new Int_t[nMaxParticles];
for(Int_t i=0;i<nMaxParticles;i++) nDParticles[i]=-2;
// local (for this evt) counter variables
Int_t nl_Part=0;
Int_t nl_RM=0;
Int_t nl_RT_S=0;
Int_t nl_MT_S=0;
Int_t nl_RMT_S=0;
Int_t nGNbRMonly=0;
Int_t nGNbRTSonly=0;
Int_t nGNbMTSonly=0;
Int_t nGNbRMTSonly=0;
Float_t fl_RW=0.;
Float_t fl_MW=0.;
Int_t nlNb_Fakes=0;
Int_t nlNb_ConfTracklet=0;
Int_t nl_ParId[MAXPARIDS];
for(Int_t k=0;k<MAXPARIDS;k++) nl_ParId[k]=0;
// loop over corr objs in category
while((pHitMatch = (HHitMatchSim *)pIterMatchHit->Next()))
{
HTrackInfo *t = pHitMatch->getTrackInfoObj();
pHitMatch->resetSim();
// check if ring is made by lepton
pHitMatch->setGLepRing(isGRing(pHitMatch));
// check if ring is made by lepton which is seen
// by other subdetectors
pHitMatch->setGCLepRing(isGCorrRing(t));
// number of different GEANT particles stored in
// tracklet that were seen by at least by 2 subdet
pHitMatch->setNbGPart(t->getPartNr());
//check if tracklet contains
//a ring that stems from a lepton that can be further
//correlated and went through the mirror
// Int_t nMatchedMirrorRingTrkNb = isLepOnMirror(pHitMatch);
Int_t nMatchedMirrorRingTrkNb = -3;
//check if particle was seen by GEANT in first layer of MDC
//meaning that it didn't perish somewhere in the Rich
if (nMatchedMirrorRingTrkNb>-1) isLepOnMDC(pHitMatch,nMatchedMirrorRingTrkNb);
Int_t isFirstGeantNbInTracklet=1;
for(Int_t i=0;i<MAXPARTICLES;i++)
{
Int_t nTrkNb=t->getTrkNb(i);
if (nTrkNb>-1)
{
if(isFirstGeantNbInTracklet) nlNb_ConfTracklet++;
// two tracklets can contain the same GEANT particle !
isFirstGeantNbInTracklet=0;
pHitMatch->setFakeFlag(0);
Int_t nParId=t->getParId(i);
if(!pHitMatch->isFake() && (nParId==2 || nParId==3) )
{
pHitMatch->setLeptonFlag(1);
}
Int_t index=fillUniqueIndex(nTrkNb,
nMaxParticles,
nDParticles);
if(index>0)
{//it is a new GEANT particle in this evt!!
// index was incremented by 1 in fillUnique....!
// only tracklets containing different GEANT
// particles are counted here !!!
//Int_t nParId=t->getParId(i);
// cout<<"it's a new particle: trk: "<<nTrkNb
// <<" part id: "<<nParId<<endl;
Float_t fRW=t->getRichWeight(i);
Float_t fMW=t->getMdcWeight(i);
Int_t nMatchRM=t->getMatchedRM(i);
Int_t nMatchRT_S=t->getMatchedRT_S(i);
Int_t nMatchMT_S=t->getMatchedMT_S(i);
Int_t nMatchRMT_S=t->getMatchedRMT_S(i);
nl_Part++;// count new part in this evt
if(nParId<MAXPARIDS) nl_ParId[nParId]+=1;
if(fRW>-1) fl_RW+=fRW;
if(fMW>-1) fl_MW+=fMW;
if(nMatchRM != -1) nl_RM++;
if(nMatchRT_S != -1) nl_RT_S++;
if(nMatchMT_S != -1) nl_MT_S++;
if(nMatchRMT_S != -1) nl_RMT_S++;
if(nMatchRM!=-1 && nMatchRT_S==-1 &&
nMatchMT_S==-1 && nMatchRMT_S==-1 ) nGNbRMonly++;
if(nMatchRM==-1 && nMatchRT_S!=-1 &&
nMatchMT_S==-1 && nMatchRMT_S==-1 ) nGNbRTSonly++;
if(nMatchRM==-1 && nMatchRT_S==-1 &&
nMatchMT_S!=-1 && nMatchRMT_S==-1 ) nGNbMTSonly++;
if(nMatchRM!=-1 && nMatchRT_S!=-1 &&
nMatchMT_S!=-1 && nMatchRMT_S!=-1 ) nGNbRMTSonly++;
}// endif new particle
}// endif valid trk nb
}//endfor particle in track info object
if(pHitMatch->isFake()) nlNb_Fakes++;
// DIAGNOSTICS
// diagOut(pHitMatch,t);
// -------------------
}// end while hhitmatchsim
//Int_t nParttmp=countDiscreteIndexes(nDParticles,nMaxParticles);
delete nDParticles;
if(nl_Part>0)
{
fl_RW/=nl_Part;//avrg Rich weight
fl_MW/=nl_Part;//avrg Mdc weight
}
pHMH->setNbPart(nl_Part);
ng_Part+=nl_Part;
ngNb_ConfTracklet+=nlNb_ConfTracklet;
pHMH->setNbConfTracklets(nlNb_ConfTracklet);
ngNb_Fakes+=nlNb_Fakes;
pHMH->setNbFakes(nlNb_Fakes);
pHMH->setAvrgRichWeight(fl_RW);
if(fl_RW>0.)
{
fg_RW+=fl_RW;//add avg Rich weight for this evt
ngNb_RW++;//count that there was a particle in Rich
}
pHMH->setAvrgMdcWeight(fl_MW);
if(fl_MW>0.)
{
fg_MW+=fl_MW;
ngNb_MW++;
}
for(Int_t n=0;n<MAXPARIDS;n++)
{
ng_ParId[n]+=nl_ParId[n];
if(nl_ParId[n]>0) pHMH->setParId(n,nl_ParId[n]);
}
// pHMH->setGeantNb_RM(nl_RM);
// ngGNb_RM+=nl_RM;
// pHMH->setGeantNb_RT_S(nl_RT_S);
// ngGNb_RT_S+=nl_RT_S;
// pHMH->setGeantNb_MT_S(nl_MT_S);
// ngGNb_MT_S+=nl_MT_S;
// pHMH->setGeantNb_RMT_S(nl_RMT_S);
// ngGNb_RMT_S+=nl_RMT_S;
pHMH->setGeantNb_RM(nGNbRMonly);
ngGNbRMonly+=nGNbRMonly;
pHMH->setGeantNb_RT_S(nGNbRTSonly);
ngGNbRTSonly+=nGNbRTSonly;
pHMH->setGeantNb_MT_S(nGNbMTSonly);
ngGNbMTSonly+=nGNbMTSonly;
pHMH->setGeantNb_RMT_S(nGNbRMTSonly);
ngGNbRMTSonly+=nGNbRMTSonly;
HRecEvent* pGeantEvent = (HRecEvent*)(gHades->getCurrentEvent());
HPartialEvent* pPartialEvent = pGeantEvent->getPartialEvent(catSimul);
HGeantHeader* pGeantHeader = (HGeantHeader*)(pPartialEvent->getSubHeader());
pHMH->setImpactParam(pGeantHeader->getImpactParameter());
//pHMH->dumpToStdoutSim();
}//this hitmatch header, there is only one !
}//end of base execute
// cout<<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "<<endl;
return 0;
}
void HRichCorrelatorSim::diagOut(HHitMatchSim* pHitMatch,HTrackInfo* t)
{
cout<<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "
<<nCounterProcessedNbEvents<<endl;
pHitMatch->dumpToStdoutSim();
diag(pHitMatch);
// t->dumpToStdout();
for(Int_t i=0;i<MAXPARTICLES;i++)
{
if (t->getTrkNb(i)>-1)
{
HRichUtilFunc::dumpKine(HRichUtilFunc::
getKineObj(t->getTrkNb(i),
(HLinearCategory*)
getGeantKineCat()));
}
}
cout<<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "<<endl;
}
void HRichCorrelatorSim::cleanTracks()
{
// temp helper arrays to hold GEANT track numbers and their weights
// need to be reset as they are data members
// GEANT starts track numbers at 1, therefore init them to 0
for(Int_t i = 0;i<RICHMAXTRACKS;i++){
richTracks[i] = 0;
richWeights[i] = 0;
}
for(Int_t i = 0;i<MDCMAXTRACKS;i++){
mdcTracks[i] = 0;
mdcWeights[i] = 0;
}
for(Int_t i = 0;i<TOFMAXTRACKS;i++) tofTracks[i]=0;
for(Int_t i = 0;i<SHOWERMAXTRACKS;i++) showerTracks[i] = 0;
}
void HRichCorrelatorSim::assignRichTracks(Int_t richInd)
{
// copy GEANT track info from the hit
HRichHitSim *pRichHit = ((HRichHitSim*)getRichHitCat()
->getObject(richInd));
richTracks[0] = pRichHit->track1;
richWeights[0] = pRichHit->weigTrack1;
richTracks[1] = pRichHit->track2;
richWeights[1] = pRichHit->weigTrack2;
richTracks[2] = pRichHit->track3;
richWeights[2] = pRichHit->weigTrack3;
}
void HRichCorrelatorSim::assignMdcTracks(Int_t mdcInd)
{
// according to info from Vladimir Pechenov
// copy GEANT track info from the hit
HMdcSegSim *pMdcSeg = ((HMdcSegSim*)getMdcSegCat()
->getObject(mdcInd));
for(Int_t i =0;i<pMdcSeg->getNTracks();i++)
{
mdcTracks[i] = pMdcSeg->getTrack(i);
mdcWeights[i] = (Int_t)pMdcSeg->getNTimes(i);
}
}
void HRichCorrelatorSim::assignShowerTrack(Int_t showerInd)
{
// look in HShowerTrack and HShowerHitTof
// and match the info
// basically what is done in HShowerHitTofMatcher
// but there for each GEANT track number one clone is created
// which is not useful here --- > we do it ourselves
HShowerTrack *pTrack=0;
iter_showertrack->Reset();
HShowerHitTof * pHit = ((HShowerHitTof*) getShowerTofHitCat()
->getObject(showerInd));
pTrack = (HShowerTrack *)iter_showertrack->Next();
do {
if (pTrack && pHit->getAddress()==pTrack->getAddress()) {
*showerTracks = pTrack->getTrack();
break;
}
}
while((pTrack = (HShowerTrack *)iter_showertrack->Next())) ;
//searching other tracks with the same address number
Int_t i=0;
while((pTrack = (HShowerTrack *)iter_showertrack->Next())) {
if (i<SHOWERMAXTRACKS-1 && pHit->getAddress()==pTrack->getAddress()) {
showerTracks[i++] = pTrack->getTrack();
} else break;
}
}
void HRichCorrelatorSim::assignTofTracks(Int_t tofInd)
{
HTofHitSim *pTofHit = ((HTofHitSim*)getTofHitCat()
->getObject(tofInd));
tofTracks[0] = pTofHit ->getNTrack1();
tofTracks[1] = pTofHit ->getNTrack2();
}
void HRichCorrelatorSim::matchTracks(HHitMatchSim *pHitMatch)
{
Bool_t kRICH = kFALSE;
Bool_t kMDC = kFALSE;
Bool_t kTOF = kFALSE;
Bool_t kSHOWER = kFALSE;
for(Int_t i = 0;i<RICHMAXTRACKS;i++)
{
if (richTracks[i] > 0) {kRICH=kTRUE; break;}
}
for(Int_t i = 0;i<MDCMAXTRACKS;i++)
{
if (mdcTracks[i] > 0) {kMDC=kTRUE; break;}
}
for(Int_t i = 0;i<TOFMAXTRACKS;i++)
{
if (tofTracks[i] > 0) {kTOF=kTRUE; break;}
}
for(Int_t i = 0;i<SHOWERMAXTRACKS;i++)
{
if (showerTracks[i] > 0) {kSHOWER=kTRUE; break;}
}
// case RICH-MDC only
if(kRICH && kMDC && !kSHOWER && !kTOF)
{
corrTracksRM(pHitMatch);
}
// case MDC and TOF or SHOWER
if(!kRICH && kMDC && (kSHOWER || kTOF) )
{
corrTracksMT_S(pHitMatch);
}
// case RICH and TOF or SHOWER
if(kRICH && !kMDC && (kSHOWER || kTOF) )
{
corrTracksRT_S(pHitMatch);
}
// case RICH - MDC and TOF or SHOWER
if(kRICH && kMDC && (kSHOWER || kTOF) )
{
corrTracksRM(pHitMatch);
corrTracksRT_S(pHitMatch);
corrTracksMT_S(pHitMatch);
evalBiCorrs(pHitMatch);
}
}
void HRichCorrelatorSim::evalBiCorrs(HHitMatchSim *pHitMatch)
{
// all possible RICH MDC META bi correlations are done
// now check whether there is a tri correlation of GEANT track numbers
HTrackInfo* t = pHitMatch->getTrackInfoObj();
for (Int_t i=0;i<MAXPARTICLES;i++)
{
Int_t nRM = t->getMatchedRM(i);
Int_t nRT_S = t->getMatchedRT_S(i);
Int_t nMT_S = t->getMatchedMT_S(i);
if (nRM==1 && nMT_S==1 && nRT_S==1) t->setMatchedRMT_S(i);
}
}
void HRichCorrelatorSim::corrTracksRT_S(HHitMatchSim *pHitMatch)
{
// look if we find the same GEANT track number for
// the rich hit and the hit in shower or tof
HTrackInfo* t = pHitMatch->getTrackInfoObj();
for (Int_t i=0;i<RICHMAXTRACKS;i++)
{
for (Int_t j=0;j<SHOWERMAXTRACKS;j++)
{
// store same GEANT nb found in RICH Hit and SHOWER Hit
if(richTracks[i]>0 && richTracks[i] == showerTracks[j])
{
setGeantTrackInfoRT_S(i,t);
break;
}
}
for (Int_t j=0;j<TOFMAXTRACKS;j++)
{
if(richTracks[i]>0 && richTracks[i] == tofTracks[j])
{
setGeantTrackInfoRT_S(i,t);
break;
}
}
}
}
void HRichCorrelatorSim::corrTracksMT_S(HHitMatchSim *pHitMatch){
HTrackInfo* t = pHitMatch->getTrackInfoObj();
for (Int_t i=0;i<MDCMAXTRACKS;i++)
{
for (Int_t j=0;j<SHOWERMAXTRACKS;j++)
{
if(mdcTracks[i]>0 && mdcTracks[i] == showerTracks[j])
{
setGeantTrackInfoMT_S(i,t);
break;
}
}
for (Int_t j=0;j<TOFMAXTRACKS;j++)
{
if(mdcTracks[i]>0 && mdcTracks[i] == tofTracks[j])
{
setGeantTrackInfoMT_S(i,t);
break;
}
}
}
}
void HRichCorrelatorSim::corrTracksRM(HHitMatchSim* pHitMatch)
{
HTrackInfo* t = pHitMatch->getTrackInfoObj();
for (Int_t i=0;i<RICHMAXTRACKS;i++)
{
for (Int_t j=0;j<MDCMAXTRACKS;j++)
{
if(richTracks[i]>0 && richTracks[i] == mdcTracks[j])
{
setGeantTrackInfoRM(i,j,t);
break;
}
}
}
}
Float_t HRichCorrelatorSim::calcRichWeight(Int_t ind)
{
// use the number of fired pads in a ring per tracknumber as
// a confidence value for this track number
Float_t sumWeightsRich=0.;
for(Int_t i=0;i<RICHMAXTRACKS;i++) sumWeightsRich += richWeights[i];
Float_t f = (Float_t) richWeights[ind];
return f/sumWeightsRich;
}
Float_t HRichCorrelatorSim::calcMdcWeight(Int_t ind)
{
// use number of drift times that have contributed to this segment
// and that come from one GEANT number as confidence value
Float_t sumWeightsMdc=0.;
for(Int_t j =0;j<MDCMAXTRACKS;j++) sumWeightsMdc += mdcWeights[j];
Float_t f = (Float_t) mdcWeights[ind];
return f/sumWeightsMdc;
}
void HRichCorrelatorSim::setGeantTrackInfoRT_S(Int_t index,HTrackInfo *t)
{
//check if geant nb is already stored
//we store only a maximum of MAXPARTICLES particles,
//the rest is discarded
Int_t trknb = richTracks[index];
if(trknb<=0) return;
Int_t newSlot=-1;
for (Int_t i=0;i<MAXPARTICLES;i++)
{
if(t->getTrkNb(i)<0) //assumes init of array elements to -1;
{// array element i is still -1, we can use it for a new particle
newSlot=i;
break;
}
if(t->getTrkNb(i) == trknb)
{// GEANT nb is already stored in Trackinfo
t->setMatchedRT_S(i);
if(t->getRichWeight(i)==-1) t->setRichWeight(i,calcRichWeight(index));
newSlot=-1;
break;
}
}
// newSlot is either -1 (nothing to be done)
// or i<MAXPARTICLES indicating a new as yet unused slot in the particle array
if (newSlot>=0) //insert new Geant Nb at newSlot
{
HGeantKine* kine=(HGeantKine*)HRichUtilFunc::getKineObj(trknb,(HLinearCategory*)getGeantKineCat());
if (kine==0)
{
Warning("setGeantTrackInfoRT_S","track number stored in rich hit not found in GEANT kine container");
cout<<"track number that creates problem: "<<trknb<<endl;
return;
}
Int_t aTrack, aID;
Int_t aPar, aMed, aMech;
Float_t ptot;
kine->getParticle(aTrack,aID);
kine->getCreator(aPar,aMed,aMech);
ptot=kine->getTotalMomentum();
if (aTrack!=trknb) Error("HRichCorrelatorSim::setGeantTrackInfoRT_S","track numbers not equal");
t->setTrkNb(newSlot,trknb);
t->setParId(newSlot,aID);
t->setMech(newSlot,aMech);
t->setCreatorId(newSlot,HRichUtilFunc::getParentParID(kine,(HLinearCategory*)getGeantKineCat()));
t->setCreatorTrkNb(newSlot,aPar);
t->setMed(newSlot,aMed);
Float_t xvert,yvert,zvert;
kine->getVertex(xvert,yvert,zvert);
t->setVertx(newSlot,xvert);
t->setVerty(newSlot,yvert);
t->setVertz(newSlot,zvert);
if (aID==2 || aID==3)
{
t->setPi0Dalitz(newSlot,isPi0Dalitz(trknb));
t->setConvLep(newSlot,isConvLep(trknb));
}
t->setTotMom(newSlot,ptot);
t->setRichWeight(newSlot,calcRichWeight(index));
t->setMatchedRT_S(newSlot);
}
}
void HRichCorrelatorSim::setGeantTrackInfoMT_S(Int_t index,HTrackInfo *t)
{
//check if geant nb is already stored
Int_t trknb = mdcTracks[index];
Int_t newSlot=-1;
for (Int_t i=0;i<MAXPARTICLES;i++)
{
if(t->getTrkNb(i)<0) //assumes init of array elements to -1;
{
newSlot=i;
break;
}
if(t->getTrkNb(i) == trknb)
{
t->setMatchedMT_S(i);
if(t->getMdcWeight(i)==-1) t->setMdcWeight(i,calcMdcWeight(index));
newSlot=-1;
break;
}
}
if (newSlot>=0) //insert new Geant Nb
{
HGeantKine* kine=(HGeantKine*)HRichUtilFunc::getKineObj(trknb,(HLinearCategory*)getGeantKineCat());
Int_t aTrack, aID;
Int_t aPar, aMed, aMech;
Float_t ptot;
kine->getParticle(aTrack,aID);
kine->getCreator(aPar,aMed,aMech);
ptot=kine->getTotalMomentum();
if (aTrack!=trknb) Error("HRichCorrelatorSim::setGeantTrackInfoRT_S","track numbers not equal");
t->setTrkNb(newSlot,trknb);
t->setParId(newSlot,aID);
t->setMech(newSlot,aMech);
t->setCreatorId(newSlot,HRichUtilFunc::getParentParID(kine,(HLinearCategory*)getGeantKineCat()));
t->setCreatorTrkNb(newSlot,aPar);
Float_t xvert,yvert,zvert;
kine->getVertex(xvert,yvert,zvert);
t->setVertx(newSlot,xvert);
t->setVerty(newSlot,yvert);
t->setVertz(newSlot,zvert);
t->setMed(newSlot,aMed);
if (aID==2 || aID==3)
{
t->setPi0Dalitz(newSlot,isPi0Dalitz(trknb));
t->setConvLep(newSlot,isConvLep(trknb));
}
t->setTotMom(newSlot,ptot);
t->setMdcWeight(newSlot,calcMdcWeight(index));
t->setMatchedMT_S(newSlot);
}
}
void HRichCorrelatorSim::setGeantTrackInfoRM(Int_t richindex,Int_t mdcindex,
HTrackInfo *t)
{
//this func needs two indexes as we need to find the weight for both rich and mdc
//check if geant nb is already stored
Int_t trknb = richTracks[richindex];
if (trknb<0) return;
Int_t newSlot=-1;
for (Int_t i=0;i<MAXPARTICLES;i++)
{
if(t->getTrkNb(i)<0) //assumes init of array elements to -1;
{
newSlot=i;
break;
}
if(t->getTrkNb(i) == trknb)
{
t->setMatchedRM(i);
if(t->getMdcWeight(i)==-1) t->setMdcWeight(i,calcMdcWeight(mdcindex));
if(t->getRichWeight(i)==-1) t->setRichWeight(i,calcRichWeight(richindex));
newSlot=-1;
break;
}
}
if (newSlot>=0) //insert new Geant Nb
{
HGeantKine* kine=(HGeantKine*)HRichUtilFunc::getKineObj(trknb,(HLinearCategory*)getGeantKineCat());
if (kine==0)
{
Warning("setGeantTrackInfoRT_S","track number stored in rich hit not found in GEANT kine container");
cout<<"track number that creates problem: "<<trknb<<endl;
return;
}
Int_t aTrack, aID;
Int_t aPar, aMed, aMech;
Float_t ptot;
kine->getParticle(aTrack,aID);
kine->getCreator(aPar,aMed,aMech);
ptot=kine->getTotalMomentum();
if (aTrack!=trknb) Error("HRichCorrelatorSim::setGeantTrackInfoRT_S","track numbers not equal");
t->setTrkNb(newSlot,trknb);
t->setParId(newSlot,aID);
t->setMech(newSlot,aMech);
t->setCreatorId(newSlot,HRichUtilFunc::getParentParID(kine,(HLinearCategory*)getGeantKineCat()));
t->setCreatorTrkNb(newSlot,aPar);
Float_t xvert,yvert,zvert;
kine->getVertex(xvert,yvert,zvert);
t->setVertx(newSlot,xvert);
t->setVerty(newSlot,yvert);
t->setVertz(newSlot,zvert);
t->setMed(newSlot,aMed);
if (aID==2 || aID==3)
{
t->setPi0Dalitz(newSlot,isPi0Dalitz(trknb));
t->setConvLep(newSlot,isConvLep(trknb));
}
t->setTotMom(newSlot,ptot);
t->setRichWeight(newSlot,calcRichWeight(richindex));
t->setMdcWeight(newSlot,calcMdcWeight(mdcindex));
t->setMatchedRM(newSlot);
}
}
Int_t HRichCorrelatorSim::isPi0Dalitz(Int_t trk)
{
// returns 1 if track nb corresponds to a lepton from pi0 Dalitz decay
// else returns zero
HGeantKine * kine =0;
// loop over kine container
iter_kine->Reset();
while((kine=(HGeantKine *)iter_kine->Next())!=0)
{
Int_t aTrack, aID;
Int_t aPar, aMed, aMech;
Float_t ptot;
kine->getParticle(aTrack,aID);
kine->getCreator(aPar,aMed,aMech);
ptot=kine->getTotalMomentum();
if (aTrack==trk)
{
if ( aMech==5 && (aID == 2 || aID==3) && (aMed>=8 && aMed<=19)) // direct particle decay
{
HGeantKine *kine_parent1 = HRichUtilFunc::findParent(kine,(HLinearCategory*) getGeantKineCat());
if(kine_parent1){//parent should be pion
// Float_t theta1,phi1;
// theta1=phi1=0.;
// HRichUtilFunc::calcParticleAnglesKine(kine_parent1,theta1,phi1);
Int_t aTrackp1, aIDp1;
kine_parent1->getParticle(aTrackp1,aIDp1);
//if(aIDp1==7 && (theta1>=10. && theta1<=90.))//pion
if(aIDp1==7)//pion, but no emission angle cut
{
return 1;
}
else return 0;
///////////////////////////////////////////////////////////////////////////
// get Second Pion Dalitz Lepton
///////////////////////////////////////////////////////////////////////////
// { //in acceptance
// HGeantKine *secondlepton =
// HRichUtilFunc::getSecondPionDalitzLepton(kine,(HLinearCategory*) getGeantKineCat());
// if (secondlepton)
// {
// Float_t deltaptot =
// TMath::Abs(ptot-secondlepton->getTotalMomentum());
// Float_t opangle =
// HRichUtilFunc::calcOpeningAngleKine(kine,secondlepton);
// deltaP_opangle->Fill(deltaptot,opangle);
// }
// med_dalitz->Fill(aMed);
// mech_dalitz->Fill(aMech);
// par_dalitz->Fill(aIDp1);
// par_med_dalitz->Fill(aID,aMed);
// par_mech_dalitz->Fill(aID,aMech);
// mech_mom_dalitz->Fill(aMech,ptot);
// }//pi0
///////////////////////////////////////////////////////////////////////////
}//kine_parent1
}//Mec==5 and lepton
// if (kDumpIt) dumpKine(kine);
// if (kDumpIt) dumpKine(findParent(kine));
}
} // end while kine loop
return 0;
}
Int_t HRichCorrelatorSim::isConvLep(Int_t trk)
{
// returns 1 if track nb corresponds to a lepton from gamma conv
// else returns zero
HGeantKine * kine =0;
// loop over kine container
iter_kine->Reset();
while((kine=(HGeantKine *)iter_kine->Next())!=0)
{
Int_t aTrack, aID;
Int_t aPar, aMed, aMech;
Float_t ptot;
kine->getParticle(aTrack,aID);
kine->getCreator(aPar,aMed,aMech);
ptot=kine->getTotalMomentum();
if (aTrack==trk)
{
if ( aMech==6 && (aID == 2 || aID==3) && (aMed>=8 && aMed<=19)) // photon pair production in target or Rich
{
HGeantKine *kine_parent1 = HRichUtilFunc::findParent(kine,(HLinearCategory*) getGeantKineCat());
if(kine_parent1){
Int_t aTrackp1, aIDp1;
Int_t aParp1,aMedp1,aMechp1;
kine_parent1->getParticle(aTrackp1,aIDp1);
kine_parent1->getCreator(aParp1,aMedp1,aMechp1);
if(aIDp1==1 && aMechp1==5)//gamma
{
HGeantKine *kine_parent2 = HRichUtilFunc::findParent(kine_parent1,(HLinearCategory*) getGeantKineCat());
if(kine_parent2){
Int_t aTrackp2, aIDp2;
kine_parent2->getParticle(aTrackp2,aIDp2);
// Float_t theta2,phi2;
// theta2=phi2=0.;
// HRichUtilFunc::calcParticleAnglesKine(kine_parent2,theta2,phi2);
//if(aIDp2==7 && (theta2>=10. && theta2<=90.))
if(aIDp2==7)
{//pi0 in acceptance
return 1;
}
else return 0;
///////////////////////////////////////////////////////////////////////////
// get Second Pion Decay Gamma
///////////////////////////////////////////////////////////////////////////
// HGeantKine *secondgamma = HRichUtilFunc::getSecondPionDecayGamma(kine_parent1,(HLinearCategory*) getGeantKineCat());
// if (secondgamma) g_conv_opangle->Fill(HRichUtilFunc::calcOpeningAngleKine(kine_parent1,secondgamma));
// med_conv->Fill(aMed);
// mech_conv->Fill(aMech);
// par_conv->Fill(aIDp2);
// par_med_conv->Fill(aID,aMed);
// par_mech_conv->Fill(aID,aMech);
// mech_mom_conv->Fill(aMech,ptot);
// }//pi0
////////////////////////////////////////////////////////////////////////////
}//kine_parent2
}//gamma
}//kine_parent1
}//Mec==6 and lepton
// if (kDumpIt) dumpKine(kine);
// if (kDumpIt) dumpKine(findParent(kine));
} // end if trk nb comparison
} // end while kine loop
return 0;
}
Int_t HRichCorrelatorSim::isLepOnMirror(HHitMatchSim *h)
{
Int_t RingCanBeInCorrelation=-1;//return value
HGeantRichMirror *mirr=0;
iter_mirror->Reset();
HRichGeometryPar* pGeo = (HRichGeometryPar*) getGeometryPar();
Float_t fYShift = pGeo->getSectorShift(); // to account for shifted volumes in HGEANT ?!
fYShift = fYShift/TMath::Cos(20./57.3);
HRichPadTab* pPadsTab = ((HRichGeometryPar*)pGeo)->getPadsPar();
HRichPad *pPad = 0;
// x,y in pad units of ring recognized on padplane
HRichHitSim* r=0;
Int_t lep1Count = 0;
Int_t lep2Count = 0;
if(h->getRichInd()>-1) r = ((HRichHitSim*)getRichHitCat()
->getObject(h->getRichInd()));
else return -1;//no rich hit in tracklet
if (r)
{
Int_t xRing=r->getRingCenterX();
Int_t yRing=r->getRingCenterY();
// cout<<" ++++++++++++++++++++++++++++++++++++++++++"<<endl;
//cout<<" ring X "<<r->getRingCenterX()<<" ring Y "<<r->getRingCenterY()<<endl;
Int_t secRing=r->getSector();
Float_t weight1 =0.;
Float_t weight2 =0.;
Float_t weightTemp = 0.;
// iterate over hits on mirror and compare hit ring x,y with
// center of reflected photons, calculate pad hit by photon ctr
while((mirr=(HGeantRichMirror *)iter_mirror->Next())!=0)
{
Int_t aTrk,aID;
mirr->getTrack(aTrk,aID);//tracknumber and pid from HGeant
// cout<<" mirror track "<<aTrk<<endl;
weightTemp = 0;
weightTemp = matchRingMirrtracks(r,aTrk,aID);
if(weightTemp>0){
Float_t x=mirr->getYRing()/10.;//in mm of the padplane
Float_t y=mirr->getXRing()/10. + fYShift;
pPad=pPadsTab->getPad(x, // they are swapped, sic !
y );// divide by 10 to account for mm->cm
if (!pPad) return -1;
Int_t sec=mirr->getSector();
Int_t xMirrPhotCM=pPad->getPadX();
Int_t yMirrPhotCM=pPad->getPadY();
// cout<<" mirror X "<<xMirrPhotCM<<" mirror Y "<<yMirrPhotCM<<endl;
Int_t phot=mirr->getNumPhot();
//cout<<" number of photons "<<phot<<endl;
if ((TMath::Abs(xMirrPhotCM-xRing)<2 &&
TMath::Abs(yMirrPhotCM-yRing)<2 &&
sec==secRing) && lep1Count==0){
RingCanBeInCorrelation=aTrk;
h->setLeptonOnMirror(1);
h->setNumPhot1(phot);
weight1 = weightTemp;
lep1Count++;
HTrackInfo* t=(HTrackInfo*) h->getTrackInfoObj();
for (Int_t i=0;i<MAXPARTICLES;i++)
{
if( t->getTrkNb(i)>-1 && t->getTrkNb(i)==aTrk)
{
if (t->getMatchedRM(i)==1 || t->getMatchedRT_S(i)==1 ||
t->getMatchedRMT_S(i)==1)
{
if (t->getParId(i)==2 || t->getParId(i)==3)
{
// ring in tracklet was made by lepton
// in correlation with other detector
// store the number of photons
// before reflection on mirror
t->setNumPhot(i,phot);
//cout<<"photon nb set: "<<phot<<endl;
} else Error("HRichCorrelatorSim::isLepOnMirror(HHitMatchSim *h)","it's not a lepton !!");
}
} else break;
}
}
else{
if ((TMath::Abs(xMirrPhotCM-xRing)<4 &&
TMath::Abs(yMirrPhotCM-yRing)<4 &&
sec==secRing) && lep2Count==0){
h->setLeptonOnMirror(2);
h->setNumPhot2(phot);
weight2 = weightTemp;
lep2Count++;
}
}
}
}
// cout<<" nr of lep1 "<<lep1Count<< " nr of lep2 "<<lep2Count<<endl;
if (lep1Count==0){
h->setLeptonOnMirror(0);
h->setNumPhot2(0);
}
else{
if(lep2Count>0){
h->setWeightRatio(weight1>weight2?weight1/weight2:weight2/weight1);
// cout<<" weight ratio "<< h->getWeightRatio()<<endl;
}
}
}else Error("HRichCorrelatorSim::isLepOnMirror(HHitMatchSim *h)","no valid Rich hit found");
// cout<<" RingCanBeInCorrelation "<<RingCanBeInCorrelation<<endl;
return RingCanBeInCorrelation;
}
Int_t HRichCorrelatorSim::matchRingMirrtracks(HRichHitSim *pH,Int_t trackMirr,Int_t idMirr){
Int_t t1,t2,t3;
t1=t2=t3=-1;
t1=pH->track1;t2=pH->track2;t3=pH->track3;
// cout<<" rich track 1 "<<t1<< " rich track2 "<<t2<<" rich track3 "<<t3<<endl;
//cout<<" rich weight 1 "<<pH->weigTrack1<<" rich weight 2 "<<pH->weigTrack2<<" rich weight 3 "<<pH->weigTrack3<<endl;
Int_t p1,p2,p3;p1=p2=p3=-1;
if (t1>-1) p1=HRichUtilFunc::getParID(t1,(HLinearCategory*)
getGeantKineCat());
if (t2>-1) p2=HRichUtilFunc::getParID(t2,(HLinearCategory*)
getGeantKineCat());
if (t3>-1) p3=HRichUtilFunc::getParID(t3,(HLinearCategory*)
getGeantKineCat());
//cout<<" rich ID 1 "<<p1<<"rich ID 2 "<<p2<<" rich ID 3 "<<p3<<endl;
//cout<<" mirror track "<<trackMirr<<" mirror id "<<idMirr<<endl;
if (t1==trackMirr && idMirr==p1) return pH->weigTrack1;
if (t2==trackMirr && idMirr==p2) return pH->weigTrack2;
if (t3==trackMirr && idMirr==p3) return pH->weigTrack3;
return 0;
}
Int_t HRichCorrelatorSim::isLepOnMDC(HHitMatchSim* h,Int_t trknb)
{
//check whether the particle with GEANT track number trk
//has reached the MDC module 0, layer 0
//this means that it didn't hit a spoke or edge in the RICH
//and can therefore be correlated and detected, finally going
//in efficiency.
//returns 1 if particle was seen in MDC, or 0 else
//trknb is return value of isLepOnMirror, the track number of
//the first found lepton passing through the mirror
//matching a track number in the ring
Int_t ParticleSeenInMDC=-1;// -1: not seen in MDC
if (h->getNumPhot1()>0) //has Rich ring made of lepton photons on mirr
{
HGeantMdc *gmdc = 0;
iter_mdcgeant->Reset();
while((gmdc=(HGeantMdc *)iter_mdcgeant->Next())!=0)
{
//cout<<gmdc->getTrack()<<" "<<(Int_t)gmdc->getSector()<<" ";
//cout<<(Int_t)gmdc->getModule()<<" "<<(Int_t)gmdc->getLayer()<<endl;
Int_t mod = (Int_t) gmdc->getModule();
Int_t lay = (Int_t) gmdc->getLayer();
if (gmdc->getTrack()==trknb && mod==0 && lay==0)
{//same track number of particle that went through
//the mirror
ParticleSeenInMDC=1;//enough for return value
h->setGLepInMDC(ParticleSeenInMDC);
//cout<<"flag indicating that lep reached MDC "<<h->getGLepInMDC()<<endl;
//check if this trk nb is in different subdetector hits of this tracklet
HTrackInfo* t=(HTrackInfo*) h->getTrackInfoObj();
for (Int_t i=0;i<MAXPARTICLES;i++)
{
if( t->getTrkNb(i)>-1 && t->getTrkNb(i)==trknb) //trk nb occurred in a corr hit
{
if (t->getMatchedRM(i)==1 || t->getMatchedRT_S(i)==1 ||
t->getMatchedRMT_S(i)==1)//not a fake
{
t->setGCLepInMDC(i,ParticleSeenInMDC);
break;
}
}
}
}
}//end while
}
return ParticleSeenInMDC;
}
Bool_t HRichCorrelatorSim::finalize() {
HRichCorrelator::finalize();
TFile *f = (TFile*) gHades->getOutputFile();
TString filename("testfile.stat");
if (f){
cout<<((TFile*)f)->GetName()<<endl;
filename = ((TFile*)gHades->getOutputFile())->GetName();
}
filename.Remove(filename.Length()-5, filename.Length());
filename.Append(".stat");
std::ofstream statusOut(filename.Data(),ios::app
// dump to stdout
cout<<"--- SIM INFO ---"<<endl;
cout<<"Nb of GEANT particles :"<<ng_Part<<endl;
cout<<"Particle ID/# ... ";
for(Int_t n=0;n<MAXPARIDS;n++)
{
if(ng_ParId[n]>0) cout<<n<<"/"<<ng_ParId[n]<<" ";
}
cout<<endl;
if(ngNb_RW>0.) cout<<"Avrg Rich Weight :"<<fg_RW/ngNb_RW<<endl;
if(ngNb_MW>0.) cout<<"Avrg Mdc Weight :"<<fg_MW/ngNb_MW<<endl;
cout<<"Nb of GEANT confirmed tracklets: "<<ngNb_ConfTracklet<<endl;
cout<<"Nb of fake tracklets: "<<ngNb_Fakes<<endl;
// cout<<"Nb of Geant confirmed RM tracklets: "<<ngGNb_RM<<endl;
// cout<<"Nb of Geant confirmed R (T||S) tracklets: "<<ngGNb_RT_S<<endl;
// cout<<"Nb of Geant confirmed M (T||S) tracklets: "<<ngGNb_MT_S<<endl;
// cout<<"Nb of Geant confirmed RM(T||S) tracklets: "<<ngGNb_RMT_S<<endl;
cout<<"Nb of Geant confirmed RM tracklets: "<<ngGNbRMonly<<endl;
cout<<"Nb of Geant confirmed R (T||S) tracklets: "<<ngGNbRTSonly<<endl;
cout<<"Nb of Geant confirmed M (T||S) tracklets: "<<ngGNbMTSonly<<endl;
cout<<"Nb of Geant confirmed RM(T||S) tracklets: "<<ngGNbRMTSonly<<endl;
// append to status file
statusOut<<"--- SIM INFO ---"<<endl;
statusOut<<"Nb of Geant particles :"<<ng_Part<<endl;
statusOut<<"Particle ID/# ... ";
for(Int_t n=0;n<MAXPARIDS;n++)
{
if(ng_ParId[n]>0) statusOut<<n<<"/"<<ng_ParId[n]<<" ";
}
statusOut<<endl;
if(ngNb_RW>0.) statusOut<<"Avrg Rich Weight :"<<fg_RW/ngNb_RW<<endl;
if(ngNb_MW>0.) statusOut<<"Avrg Mdc Weight :"<<fg_MW/ngNb_MW<<endl;
statusOut<<"Nb of GEANT confirmed tracklets: "<<ngNb_ConfTracklet<<endl;
statusOut<<"Nb of fake tracklets: "<<ngNb_Fakes<<endl;
// statusOut<<"Nb of Geant confirmed RM tracklets: "<<ngGNb_RM<<endl;
// statusOut<<"Nb of Geant confirmed R (T||S) tracklets: "<<ngGNb_RT_S<<endl;
// statusOut<<"Nb of Geant confirmed M (T||S) tracklets: "<<ngGNb_MT_S<<endl;
// statusOut<<"Nb of Geant confirmed RM(T||S) tracklets: "<<ngGNb_RMT_S<<endl;
statusOut<<"Nb of Geant confirmed RM tracklets: "<<ngGNbRMonly<<endl;
statusOut<<"Nb of Geant confirmed R (T||S) tracklets: "<<ngGNbRTSonly<<endl;
statusOut<<"Nb of Geant confirmed M (T||S) tracklets: "<<ngGNbMTSonly<<endl;
statusOut<<"Nb of Geant confirmed RM(T||S) tracklets: "<<ngGNbRMTSonly<<endl;
statusOut.close();
return kTRUE;
}
Int_t HRichCorrelatorSim::isGCorrRing(HTrackInfo *t)
{
// this function checks whether a GEANT particle correlation
// occurs for min 1 ring generating lepton in this tracklet
Int_t n=0;//counter for confirmed and min bi-correlated lepton rings
Bool_t isCorrRing=kFALSE;
for (Int_t i=0;i<MAXPARTICLES;i++)
{
if( t->getTrkNb(i)>-1 )
{
if (t->getMatchedRM(i)==1 || t->getMatchedRT_S(i)==1 ||
t->getMatchedRMT_S(i)==1)
{
if (t->getParId(i)==2 || t->getParId(i)==3) n++;
isCorrRing=kTRUE;
}
} else break;
}
if (isCorrRing) return n;// can be 0 or >0; 0 => not a lepton ring !
else return -1; //doesn't contain a rich hit in Geant correlation
}
Int_t HRichCorrelatorSim::isGRing(HHitMatchSim *h)
{
// this function checks whether the ring that went in the
// tracklet contains pads fired by photons that come from a lepton
// one should also check whether the lepton reached the mirror
// and the first MDC module
Int_t n=0;//counter for confirmed lepton rings
// n==0: this ring is a fake, no lepton contributed
// n==-1: no Rich hit in tracklet
// n>0: n>0 leptons really contributed photons to the hit pattern
HRichHitSim* r=0;
if(h->getRichInd()>-1) r = ((HRichHitSim*)getRichHitCat()
->getObject(h->getRichInd()));
else return -1;
Int_t t1,t2,t3,w1,w2,w3;
t1=t2=t3=w1=w2=w3=-1;
if(r)
{
t1=r->track1;t2=r->track2;t3=r->track3;
w1=r->weigTrack1;w2=r->weigTrack2;w3=r->weigTrack3;
Int_t p1,p2,p3;p1=p2=p3=-1;
if (t1>-1) p1=HRichUtilFunc::getParID(t1,(HLinearCategory*)
getGeantKineCat());
if (t2>-1) p2=HRichUtilFunc::getParID(t2,(HLinearCategory*)
getGeantKineCat());
if (t3>-1) p3=HRichUtilFunc::getParID(t3,(HLinearCategory*)
getGeantKineCat());
if (p1==2 || p1==3) n++;
if (p2==2 || p2==3) n++;
if (p3==2 || p3==3) n++;
}
else Error("HRichCorrelatorSim::isGRing","no rich hit found");
return n;
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.