// File: hpseudotracking.cc
//
// Author: Jaro Bielcik <bielcik@gsi.de>
// Last update by Laura Fabbietti 06/02/01 19:30
// changed by: T. Eberl (Thomas.Eberl@ph.tum.de),
// Thu Jul 11 13:35:28 CEST 2001
using namespace std;
#include "hpseudotracking.h"
#include "hdebug.h"
#include "hades.h"
#include "tofdef.h"
#include "richdef.h"
#include "hrichdetector.h"
#include "hiterator.h"
#include "hspectrometer.h"
#include "htofdetector.h"
#include "hevent.h"
#include "hcategory.h"
#include "hlinearcategory.h"
#include <iostream>
#include <iomanip>
#include "htofhitsim.h"
#include "hgeantkine.h"
#include "TRandom.h"
#include "haddef.h"
#include "hmdccal1.h"
#include "hmdccal1sim.h"
#include "hrichhitsim.h"
#include "hgeantrich.h"
#include "hparticlesim.h"
#include "hgeantkine.h"
#include "hlocation.h"
#include "hshowerhittrack.h"
#include "hmdcdef.h"
#include "showerdef.h"
#include "TIterator.h"
#include <TArrayI.h>
#include "hkinesim.h"
#include "hmdcsegsim.h"
#include "hphysicsconstants.h"
#include <TLorentzVector.h>
#include "hevent.h"
HPseudoTracking::HPseudoTracking(void) {
rich_tracknF.Set(100);
rich_trackn.Set(100);
rich_Theta.Set(100);
rich_Phi.Set(100);
mdc_trackn=new TArrayI(100);
tof_trackn=new TArrayI(100);
shower_trackn=new TArrayI(100);
iter_rich = NULL;
iter_tof= NULL;
iter_mdc= NULL;
iter_shower= NULL;
iter_tofino= NULL;
iter= NULL;
evtCount= 0;
}
HPseudoTracking::HPseudoTracking(Text_t *name,Text_t *title,Int_t closeRej) :
HReconstructor(name,title) {
i_RJ = closeRej;
rich_tracknF.Set(100);
rich_trackn.Set(100);
rich_Theta.Set(100);
rich_Phi.Set(100);
mdc_trackn=new TArrayI(100);
tof_trackn=new TArrayI(100);
shower_trackn=new TArrayI(100);
iter_rich = NULL;
iter_tof= NULL;
iter_mdc= NULL;
iter_shower= NULL;
iter_tofino= NULL;
iter= NULL;
evtCount= 0;
}
HPseudoTracking::~HPseudoTracking(void) {
// why no deletion of heap objects ??
//if (rich_trackn) delete [] rich_trackn;
//if ( mdc_trackn) delete [] mdc_trackn;
//if (tof_trackn) delete [] tof_trackn;
//if (shower_trackn) delete [] shower_trackn;
}
Bool_t HPseudoTracking::init(void) {
HEvent *event=gHades->getCurrentEvent();
fGeantKineCat =(HLinearCategory*) gHades->getCurrentEvent()->getCategory(catGeantKine);
if (!fGeantKineCat) { cout<<"no GEANT Kine category availablen"<<endl; }
iter_kine = (HIterator *)fGeantKineCat ->MakeIterator("native");
// set HGeantKine category pointer
fHGeantKine=gHades->getCurrentEvent()->getCategory(catGeantKine);
if (!fHGeantKine){ fHGeantKine = new HLinearCategory("HGeantKine");
gHades->getCurrentEvent()->addCategory(catGeantKine,fHGeantKine,"HGeantKine");
}
iter = (HIterator *)fHGeantKine->MakeIterator("native");
cout<<" after hkine "<<endl;
// fParticle = gHades->getCurrentEvent()->getCategory(catParticle);
// if (!fParticle){
// fParticle = new HLinearCategory("HParticleSim");
// gHades->getCurrentEvent()->addCategory(catParticle,fParticle,"Particle");
// }
// set particle category pointer
// =========================================================
// new category creation for HParticle
fParticle=gHades->getCurrentEvent()->getCategory(catParticle);
if (!fParticle) {
fParticle = new HLinearCategory("HParticleSim",1000);
if (fParticle) gHades->getCurrentEvent()->addCategory(catParticle,fParticle,"PhyAna");
else {
Error("init","Unable to create HParticle container");
return kFALSE;
}
}
// ===========================================================
iter_part = (HIterator *)fParticle ->MakeIterator("native");
cout<<" after particle "<<endl;
//set tof category pointer and create iterator
fTofCat=gHades->getCurrentEvent()->getCategory(catTofHit);
if (!fTofCat) {
fTofCat=gHades->getSetup()->getDetector("Tof")->buildCategory(catTofHit);
if (!fTofCat) return kFALSE;
else gHades->getCurrentEvent()->
addCategory(catTofHit,fTofCat,"Tof");
}
iter_tof=(HIterator *)fTofCat->MakeIterator("native");
cout<<"after tof"<<endl;
//set mdc category pointer and create iterator
fMdcCat=event->getCategory(catMdcSeg);
if (!fMdcCat) {
Error("init","No MDC segment category defined");
return kFALSE;
}
else iter_mdc=(HIterator *)fMdcCat->MakeIterator();
// fMdcCat=gHades->getCurrentEvent()->getCategory(catMdcCal1);
// if (!fMdcCat) {
// fMdcCat=gHades->getSetup()->getDetector("Mdc")->buildCategory(catMdcCal1);
// if (!fMdcCat) return kFALSE;
// else gHades->getCurrentEvent()->
// addCategory(catMdcCal1,fMdcCat,"Mdc");
// }
// iter_mdc=(HIterator *)fMdcCat->MakeIterator();
cout<<" after mdc "<<endl;
// set shower category pointer and create iterator
fShowerCat=gHades->getCurrentEvent()->getCategory(catShowerHitTrack);
if (!fShowerCat) {
fShowerCat=gHades->getSetup()->getDetector("Shower")->buildCategory(catShowerHitTrack);
if (!fShowerCat) return kFALSE;
else gHades->getCurrentEvent()->
addCategory(catShowerHitTrack,fShowerCat,"Shower");
}
iter_shower=(HIterator *)fShowerCat->MakeIterator("native");
// set rich category pointer and create iterator
fRichCat=gHades->getCurrentEvent()->getCategory(catRichHit);
if (!fRichCat) {
fRichCat =gHades->getSetup()->getDetector("Rich")->buildCategory(catRichHit);
if (!fRichCat) return kFALSE;
else gHades->getCurrentEvent()->addCategory(catRichHit,fRichCat,"Rich");
}
iter_rich = (HIterator *)fRichCat->MakeIterator("native");
return kTRUE;
}
Bool_t HPseudoTracking::checkThetaPhi(Float_t theta,Float_t phi){
// this function checks that the recognized lepton is produced in
// the rich and that its polar azymthal angles agreee with the HGeant one
// since the Theta RIch is off of 2 degree a correciotn is necessary ( 07/05/2001)
Int_t nTrackKine=0;
Int_t nId=0;
HGeantKine * kine =0;
Float_t x,y,z;
x = y = z =0;
Float_t pX,pY,pZ,pT,pTot,pTheta,pPhi;
pX=pY=pZ=pT=pTot=pTheta=pPhi=0;
Float_t factor = 1.0F/57.2958F;
Int_t aTRUE =0;
iter_kine->Reset();
while((kine=(HGeantKine *)iter_kine->Next())!=0){
kine->getParticle(nTrackKine,nId);
kine->getVertex(x,y,z);
// is a lepton and is created inside the rich
if( (nId ==2 || nId ==3) && (z>-40.0F && z<400.0F && x*x+y*y+(z+500.0F)*(z+500.0F) < 810000.0F)){
kine->getMomentum(pX,pY,pZ);
pT=pX*pX+pY*pY;
pTot=sqrt(pT+pZ*pZ);
pT=sqrt(pT);
pTheta=asin(pT/pTot)/factor;
pPhi=asin(pY/pT)/factor;
if(pX<0.0F) pPhi=180.0F-pPhi;
else if(pY<0.0F) pPhi+=360.0F;
// cout<<" theta kine "<<pTheta<< " phi kine "<<pPhi<<endl;
if(TMath::Abs(theta+2.-pTheta)<5 && TMath::Abs(phi-pPhi)<5){
// cout<<" theta "<< theta<< " phi "<<phi<<endl;
// cout<<" theta Geant "<< pTheta<< " phi Geant "<<pPhi<<endl;
aTRUE = 1;
break;
} // endif
} // endif particle id
} // end while loop
return aTRUE;
}
void HPseudoTracking::fillRichTrackList(){
//filling RICH tracks in event to the list
HRichHitSim* pObj_rich = NULL;
Int_t cont = 0;
Int_t rich_track1,rich_track2;
rich_track1 = rich_track2 = 0;
iter_rich->Reset();
while ((pObj_rich=(HRichHitSim*)iter_rich->Next())!= 0) {
if ((cont+1)== rich_trackn.GetSize()) { // resize
rich_Theta.Set(2*cont);
rich_Phi.Set(2*cont);
rich_trackn.Set(2*cont);
}
rich_track1=pObj_rich->track1;
rich_track2=pObj_rich->track2;
// check that the hit is a lepton and that it is produced in the RICH volume
if(checkThetaPhi(pObj_rich->getTheta(),pObj_rich->getPhi()) ){
// check for track numbers that show up twice
// function adds particle values to arrays
if (rich_track1) scanRingTrack(pObj_rich,rich_track1);
if (rich_track2) scanRingTrack(pObj_rich,rich_track2);
}
}
if(rich_tracknF.GetSize()!=rich_trackn.GetSize()) rich_tracknF.Set(rich_trackn.GetSize());
rich_tracknF = rich_trackn; // copy one array to the other
Int_t k =0;
while((rich_tracknF.At(k))!=0){ // what is this doing ???
k++;
}
}
void HPseudoTracking::scanRingTrack(HRichHitSim *richHit, Int_t track){
// the track array is scanned to see whether there are track numbers
// which show up twice ( int this case 1 copy has to be removed)
// this can happen when 2 rings are very closed to each other
// so that they have some overlap,
// but distinguished by the ring finder; then a mixing of track numbers
// occours between the 2 rings and they are assigned identical track
// numbers: both tracks for each ring.
Bool_t atrue = 0;
Int_t i =0;
while((rich_trackn.At(i))!=0){
if (track==rich_trackn.At(i)) atrue = 1;
i++;
}
if(!atrue) {
rich_trackn.AddAt(track,i);
rich_Theta.AddAt(richHit->getTheta(),i);
rich_Phi.AddAt(richHit->getPhi(),i);
}
}
void HPseudoTracking::closeLepRejNew(){
// this method apply the ideal close pairs rejection
// for leptons identified as rings the opening angle
// is calculated on the base of the tracks.
Int_t i = 0;
Int_t rich_list_track = 0;
TList * momLep = new TList();
TArrayI parentList(10);
Float_t xMom,yMom,zMom;
xMom = yMom = zMom = 0;
Int_t parent, dummy = 0;
parent = dummy = 0;
rich_tracknF.Reset();
while ((rich_list_track=rich_trackn.At(i))!=0){
getObj(rich_list_track)->getMomentum(xMom,yMom,zMom);
getObj(rich_list_track)->getCreator(parent,dummy,dummy);
momLep->Add(new TVector3(xMom,yMom,zMom));
parentList.AddAt(parent,i);
i++;
}
for(Int_t i = 0;i<momLep->GetSize();i++){
for(Int_t j = 0;j<momLep->GetSize();j++){
Float_t angle = ( (TVector3*)momLep->At(i) )->Angle( *(TVector3*)momLep->At(j) )*57.3;
if (angle<15 && angle>0 && i!=j ){
rich_trackn.AddAt(0,i);
rich_trackn.AddAt(0,j);
break;
}
}
}
Int_t cont1 = 0;
for(Int_t i = 0;i<rich_trackn.GetSize();i++){
if (rich_trackn.At(i)) {
rich_tracknF.AddAt(rich_trackn.At(i),cont1);
cont1++;
}
}
}
void HPseudoTracking::closeLepRej(){
// same as closeLepRejNew(), but w/o use of TVector3
// this method apply the close pairs rejection
// for leptons identified as rings whose opening angle is closer then
// 15 degree, the 2 leptons have to come from the same parent particle.
rich_tracknF.Reset();
Float_t factor = 57.2958;
Float_t t1,t2,f1,f2 ;
t1 = t2 = f1 = f2 =0;
Int_t i = 0;
Int_t j = 0;
while ((rich_trackn.At(i))!=0){
while ((rich_trackn.At(j))!=0){
t1 = rich_Theta.At(i)/factor;
t2 = rich_Theta.At(j)/factor;
f1 = rich_Phi.At(i)/factor;
f2 = rich_Phi.At(j)/factor;
Float_t angle = acos(sin(t1)*sin(t2)*(cos(f1)*cos(f2)+sin(f1)*sin(f2))+cos(t1)*cos(t2))*factor;
if (angle<15 && angle>0 && i!=j ){
//cout<<"close pair found &&&&&&&&&&&&&&&&&&&&&&&&& "<<endl;
// only real pairs are rejected
rich_trackn.AddAt(0,i);
rich_trackn.AddAt(0,j);
break;
}
j++;
}
i++;
}
Int_t cont = 0;
for(Int_t i = 0;i<rich_trackn.GetSize();i++){
//printf("old array :%d n",rich_trackn.At(i));
if (rich_trackn.At(i)) {
rich_tracknF.AddAt(rich_trackn.At(i),cont);
//printf("new array :%d , at: %dn",rich_tracknF.At(cont),cont);
cont++;
}
}
}
void HPseudoTracking::fillMdcTrackList(){
//filling MDC tracks in event to the list
HMdcSegSim* pObj_mdc = NULL;
Int_t cont = 0;
Int_t mdc_track = 0;
iter_mdc->Reset();
Int_t numtracks = 0;
while ((pObj_mdc=(HMdcSegSim*)iter_mdc->Next())!= 0) {
numtracks = pObj_mdc->getNTracks();
for (Int_t i=0;i<numtracks;i++){
if ((cont +1)== mdc_trackn->GetSize() ) mdc_trackn->Set(2*cont);
mdc_track=pObj_mdc->getTrack(i);
if(mdc_track>0){
mdc_trackn->AddAt(mdc_track, cont);
cont++;
}
}
}// end of while
}
void HPseudoTracking::fillTofTrackList(){
//filling TOF tracks in event to the list
HTofHitSim* pObj_tof = NULL;
Int_t cont = 0;
Int_t tof_track= 0;
iter_tof->Reset();
while ((pObj_tof=(HTofHitSim*)iter_tof->Next())!= 0) {
if ((cont + 1) == tof_trackn->GetSize()) tof_trackn->Set(2*cont);
tof_track=pObj_tof->getNTrack1();
tof_trackn->AddAt(tof_track,cont);
cont++;
tof_track=pObj_tof->getNTrack2();
if(tof_track>0){
tof_trackn->AddAt(tof_track,cont);
cont++;
}
}
}
void HPseudoTracking::fillShowerTrackList(){
//filling Shower tracks in event to the list
HShowerHitTrack* pObj_shower = NULL;
Int_t cont= 0;
Int_t shower_track = 0;
iter_shower->Reset();
while ((pObj_shower=(HShowerHitTrack*)iter_shower->Next())!= 0) {
if ((cont+1) == shower_trackn->GetSize()) shower_trackn->Set(2*cont);
shower_track=pObj_shower->getTrack();
shower_trackn->AddAt(shower_track,cont);
cont++;
}
}
void HPseudoTracking::pseudoTracking() {
//define the lists to store the track numbers of Hits in various detectors
//create iterators over lists of tracks for detectors
Int_t shower_list_track; // track number in the list
Int_t rich_list_track;
Int_t tof_list_track;
Int_t mdc_list_track;
Int_t index=0;
//here I fill info which than go to TRACK class about tracks which has hit in all detectors
TArrayI* foundtrack = new TArrayI(100);
foundtrack->Reset();
index=0;
rich_list_track=0;
mdc_list_track=0;
tof_list_track=0;
shower_list_track=0;
// iterate over the lists with the track numbers and finding the track numbers which appear in
// RICH-MDC-(TOF||SHOWER)
// those one are pseudo-found
Int_t i,j,k,m;
i = j =k =m =0;
while ( (rich_list_track=rich_tracknF.At(i)) != 0 ){
mdc_list_track=0;
j = 0;
while ((mdc_list_track=mdc_trackn->At(j))!= 0){
if(rich_list_track==mdc_list_track){
tof_list_track=0;
k = 0;
while ((tof_list_track=tof_trackn->At(k))!= 0) {
if(rich_list_track==tof_list_track){
Int_t already_in=0;
for (Int_t ii=0;ii<foundtrack->GetSize();ii++){
if(mdc_list_track==foundtrack->At(ii)) already_in=1;
// cout<<"flag already_in : "<<already_in<<endl;
} //endof for
if (already_in==0){
if (index==foundtrack->GetSize()) foundtrack->Set(index+1);
foundtrack->AddAt(mdc_list_track,index);
// printf(" ARRAYrich-mdc-tof track=:%in",mdc_list_track);
count_leptons++;
index++;
} //end if
} //endif rich-tof
k++;
} // end tof loop
shower_list_track=0;
m = 0;
while ((shower_list_track=shower_trackn->At(m))!= 0) {
if(rich_list_track==shower_list_track){
// cout<<"match rich, mdc, shower "<< rich_list_track<<endl;
Int_t already_in=0;
for (Int_t ii=0;ii<foundtrack->GetSize();ii++){
if(mdc_list_track==foundtrack->At(ii)) already_in=1;
}
if (already_in==0){
if (index==foundtrack->GetSize()) foundtrack->Set(index+1);
foundtrack->AddAt(mdc_list_track,index);
count_leptons++;
index++;
}
} // end if rich-shower
m++;
} // end if shower loop
} // end if rich-mdc matching
j++;
} // loop over mdc tracks
i++;
}// loop over rich tracks
fillPart(foundtrack); // fill TArrayI with particles.
}
void HPseudoTracking::fillPart(TArrayI * pTrackArray){
// read all found track numbers in pTrackArray and look in HGeantKine for inforamtion
// about this particle and fill the HParticleSim container
// changes needed due to the migration to the new official HParticle class
HGeantKine* pKine = NULL;
HParticleSim *pParticle = NULL;
Int_t dummy;
dummy = 0;
Int_t trknb = 0;
Int_t paridx,parentidx,aPar,aMech,aMed;
paridx = parentidx = aPar = aMech= aMed= 0;
Float_t weight = 0;
Float_t gen = 0;
Float_t xx,yy,zz;
xx = yy = zz =0;
Int_t medium,mechanism;
medium = mechanism = 0;
Float_t px,py,pz,en;
TLorentzVector mom4;
px = py = pz = en = 0;
for (Int_t ii=0;ii<pTrackArray->GetSize();ii++){
if(pTrackArray->At(ii)>0){
px = py = pz = en = 0;
pKine = (HGeantKine*) fHGeantKine->getObject(pTrackArray->At(ii)-1);
pKine->getParticle(trknb,paridx);
pKine->getVertex(xx,yy,zz);
pKine->getGenerator(gen,weight);
pKine->getMomentum(px,py,pz);
pKine->getCreator(dummy,medium,mechanism);
en= calcEnergy(px,py,pz,paridx);
mom4.SetPxPyPzE(px,py,pz,en);
HLocation loc1;
pParticle = (HParticleSim*)((HLinearCategory*)getPartCat())->getNewSlot(loc1);
if(pParticle){
pParticle = new (pParticle) HParticleSim;
pParticle->setPid(paridx);
pParticle->setVect4(mom4);
pParticle->setWeight(weight);
pParticle->setMedium(medium);
pParticle->setMechanism(mechanism);
//pParticle->SetTrackNumber(trknb);
// later an HVertex object can be added
}
}
}
}
Float_t HPseudoTracking::calcEnergy(Float_t px, Float_t py,Float_t pz,int id){
//Float_t mass = PData::Mass(id)*1000;
Float_t mass =HPhysicsConstants::mass(id);
return sqrt(px*px+py*py+pz*pz+(mass*mass));
}
HGeantKine* HPseudoTracking::getObj(Int_t nTrack){
Int_t nTrackKine=0;
Int_t nId=0;
HGeantKine * kine =0;
HGeantKine * kine1 =0;
iter->Reset();
while((kine=(HGeantKine *)iter->Next())!=0){
kine->getParticle(nTrackKine,nId);
if (nTrackKine==nTrack){
kine1 = kine;
break;
}
}
return kine1;
}
Int_t HPseudoTracking::execute(void) {
evtCount ++;
fillRichTrackList();
if(i_RJ ) closeLepRej();
fillMdcTrackList();
fillTofTrackList();
fillShowerTrackList();
pseudoTracking();
//assignEvtNum(evtCount);
rich_trackn.Reset();
rich_Theta.Reset();
rich_Phi.Reset();
rich_tracknF.Reset();
mdc_trackn->Reset();
tof_trackn->Reset();
shower_trackn->Reset();
return 0;
}
// void HPseudoTracking::assignEvtNum(Int_t evt){
// HParticleSim * partSim =NULL;
// iter_part->Reset();
// while((partSim = (HParticleSim*)iter_part->Next())!=0){
// partSim->setEvtNum(evt);
// }
// }
Bool_t HPseudoTracking::finalize(void){
return kTRUE;
}
ClassImp(HPseudoTracking)
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.