/*
Origin: marian.ivanov@cern.ch
-
-Generate complex MC information - used for Comparison later on
-How to use it?
-
-gSystem->Load("libPWG1.so")
-AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",0,0)
-t->Exec();
+Container classes with MC infomation
*/
ClassImp(AliMCInfo);
ClassImp(AliGenV0Info)
ClassImp(AliGenKinkInfo)
-ClassImp(AliGenInfoMaker)
-
-
-
-AliTPCParam * GetTPCParam(){
- AliTPCParamSR * par = new AliTPCParamSR;
- par->Update();
- return par;
-}
-
-
-//_____________________________________________________________________________
-Float_t TPCBetheBloch(Float_t bg)
-{
- //
- // Bethe-Bloch energy loss formula
- //
- const Double_t kp1=0.76176e-1;
- const Double_t kp2=10.632;
- const Double_t kp3=0.13279e-4;
- const Double_t kp4=1.8631;
- const Double_t kp5=1.9479;
-
- Double_t dbg = (Double_t) bg;
-
- Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
-
- Double_t aa = TMath::Power(beta,kp4);
- Double_t bb = TMath::Power(1./dbg,kp5);
-
- bb=TMath::Log(kp3+bb);
-
- return ((Float_t)((kp2-aa-bb)*kp1/aa));
-}
-
-
TObject(),
fTrackRef(info.fTrackRef),
fTrackRefOut(info.fTrackRefOut),
- fTRdecay(info.fTRdecay),
+ fTRdecay(info.GetTRdecay()),
fPrimPart(info.fPrimPart),
fParticle(info.fParticle),
- fMass(info.fMass),
+ fMass(info.GetMass()),
fCharge(info.fCharge),
fLabel(info.fLabel),
fEventNr(info.fEventNr),
// Update information - calculates derived variables
//
- fMCPd[0] = fMCd.fParticle.Px();
- fMCPd[1] = fMCd.fParticle.Py();
- fMCPd[2] = fMCd.fParticle.Pz();
- fMCPd[3] = fMCd.fParticle.P();
- //
- fMCX[0] = fMCd.fParticle.Vx();
- fMCX[1] = fMCd.fParticle.Vy();
- fMCX[2] = fMCd.fParticle.Vz();
+ fMCPd[0] = fMCd.GetParticle().Px();
+ fMCPd[1] = fMCd.GetParticle().Py();
+ fMCPd[2] = fMCd.GetParticle().Pz();
+ fMCPd[3] = fMCd.GetParticle().P();
+ //
+ fMCX[0] = fMCd.GetParticle().Vx();
+ fMCX[1] = fMCd.GetParticle().Vy();
+ fMCX[2] = fMCd.GetParticle().Vz();
fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
//
- fPdg[0] = fMCd.fParticle.GetPdgCode();
- fPdg[1] = fMCm.fParticle.GetPdgCode();
+ fPdg[0] = fMCd.GetParticle().GetPdgCode();
+ fPdg[1] = fMCm.GetParticle().GetPdgCode();
//
- fLab[0] = fMCd.fParticle.GetUniqueID();
- fLab[1] = fMCm.fParticle.GetUniqueID();
+ fLab[0] = fMCd.GetParticle().GetUniqueID();
+ fLab[1] = fMCm.GetParticle().GetUniqueID();
//
//
//
Double_t x1[3],p1[3];
- TParticle & pdaughter = fMCd.fParticle;
+ TParticle& pdaughter = fMCd.GetParticle();
x1[0] = pdaughter.Vx();
x1[1] = pdaughter.Vy();
x1[2] = pdaughter.Vz();
//
Double_t x2[3],p2[3];
//
- TParticle & pmother = fMCm.fParticle;
+ TParticle& pmother = fMCm.GetParticle();
x2[0] = pmother.Vx();
x2[1] = pmother.Vy();
x2[2] = pmother.Vz();
fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
- Double_t mass1 = fMCd.fMass;
- Double_t mass2 = fMCm.fMass;
+ Double_t mass1 = fMCd.GetMass();
+ Double_t mass2 = fMCm.GetMass();
Double_t e1 = TMath::Sqrt(mass1*mass1+
// Update information
// some redundancy - faster acces to the values in analysis code
//
- fMCPd[0] = fMCd.fParticle.Px();
- fMCPd[1] = fMCd.fParticle.Py();
- fMCPd[2] = fMCd.fParticle.Pz();
- fMCPd[3] = fMCd.fParticle.P();
+ fMCPd[0] = fMCd.GetParticle().Px();
+ fMCPd[1] = fMCd.GetParticle().Py();
+ fMCPd[2] = fMCd.GetParticle().Pz();
+ fMCPd[3] = fMCd.GetParticle().P();
//
- fMCX[0] = fMCd.fParticle.Vx();
- fMCX[1] = fMCd.fParticle.Vy();
- fMCX[2] = fMCd.fParticle.Vz();
+ fMCX[0] = fMCd.GetParticle().Vx();
+ fMCX[1] = fMCd.GetParticle().Vy();
+ fMCX[2] = fMCd.GetParticle().Vz();
fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
//
- fPdg[0] = fMCd.fParticle.GetPdgCode();
- fPdg[1] = fMCm.fParticle.GetPdgCode();
+ fPdg[0] = fMCd.GetParticle().GetPdgCode();
+ fPdg[1] = fMCm.GetParticle().GetPdgCode();
//
- fLab[0] = fMCd.fParticle.GetUniqueID();
- fLab[1] = fMCm.fParticle.GetUniqueID();
+ fLab[0] = fMCd.GetParticle().GetUniqueID();
+ fLab[1] = fMCm.GetParticle().GetUniqueID();
//
//
//
Double_t x1[3],p1[3];
- TParticle & pdaughter = fMCd.fParticle;
+ TParticle& pdaughter = fMCd.GetParticle();
x1[0] = pdaughter.Vx();
x1[1] = pdaughter.Vy();
x1[2] = pdaughter.Vz();
//
Double_t x2[3],p2[3];
//
- TParticle & pmother = fMCm.fParticle;
+ TParticle& pmother = fMCm.GetParticle();
x2[0] = pmother.Vx();
x2[1] = pmother.Vy();
x2[2] = pmother.Vz();
p2[1] = pmother.Py();
p2[2] = pmother.Pz();
//
- AliTrackReference & pdecay = fMCm.fTRdecay;
+ const AliTrackReference & pdecay = fMCm.GetTRdecay();
x2[0] = pdecay.X();
x2[1] = pdecay.Y();
x2[2] = pdecay.Z();
}
-
-////////////////////////////////////////////////////////////////////////
-AliGenInfoMaker::AliGenInfoMaker():
- fDebug(0), //! debug flag
- fEventNr(0), //! current event number
- fLabel(0), //! track label
- fNEvents(0), //! number of events to process
- fFirstEventNr(0), //! first event to process
- fNParticles(0), //! number of particles in TreeK
- fTreeGenTracks(0), //! output tree with generated tracks
- fTreeKinks(0), //! output tree with Kinks
- fTreeV0(0), //! output tree with V0
- fTreeHitLines(0), //! tree with hit lines
- fFileGenTracks(0), //! output file with stored fTreeGenTracks
- fLoader(0), //! pointer to the run loader
- fTreeD(0), //! current tree with digits
- fTreeTR(0), //! current tree with TR
- fStack(0), //! current stack
- fGenInfo(0), //! array with pointers to gen info
- fNInfos(0), //! number of tracks with infos
- fParamTPC(0), //! AliTPCParam
- fTPCPtCut(0.03),
- fITSPtCut(0.1),
- fTRDPtCut(0.1),
- fTOFPtCut(0.1)
-{
-}
-
-////////////////////////////////////////////////////////////////////////
-AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
- Int_t nEvents, Int_t firstEvent):
- fDebug(0), //! debug flag
- fEventNr(0), //! current event number
- fLabel(0), //! track label
- fNEvents(0), //! number of events to process
- fFirstEventNr(0), //! first event to process
- fNParticles(0), //! number of particles in TreeK
- fTreeGenTracks(0), //! output tree with generated tracks
- fTreeKinks(0), //! output tree with Kinks
- fTreeV0(0), //! output tree with V0
- fTreeHitLines(0), //! tree with hit lines
- fFileGenTracks(0), //! output file with stored fTreeGenTracks
- fLoader(0), //! pointer to the run loader
- fTreeD(0), //! current tree with digits
- fTreeTR(0), //! current tree with TR
- fStack(0), //! current stack
- fGenInfo(0), //! array with pointers to gen info
- fNInfos(0), //! number of tracks with infos
- fParamTPC(0), //! AliTPCParam
- fTPCPtCut(0.03),
- fITSPtCut(0.1),
- fTRDPtCut(0.1),
- fTOFPtCut(0.1)
-
-{
- //
- //
- //
- fFirstEventNr = firstEvent;
- fEventNr = firstEvent;
- fNEvents = nEvents;
- sprintf(fFnRes,"%s",fnRes);
- //
- fLoader = AliRunLoader::Open(fnGalice);
- if (gAlice){
- delete gAlice->GetRunLoader();
- delete gAlice;
- gAlice = 0x0;
- }
- if (fLoader->LoadgAlice()){
- cerr<<"Error occured while l"<<endl;
- }
- Int_t nall = fLoader->GetNumberOfEvents();
- if (nEvents==0) {
- nEvents =nall;
- fNEvents=nall;
- fFirstEventNr=0;
- }
-
- if (nall<=0){
- cerr<<"no events available"<<endl;
- fEventNr = 0;
- return;
- }
- if (firstEvent+nEvents>nall) {
- fEventNr = nall-firstEvent;
- cerr<<"restricted number of events availaible"<<endl;
- }
- AliMagF * magf = gAlice->Field();
- AliTracker::SetFieldMap(magf,0);
-}
-
-
-AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
-{
- //
- if (i<fNParticles) {
- if (fGenInfo[i]) return fGenInfo[i];
- fGenInfo[i] = new AliMCInfo;
- fNInfos++;
- return fGenInfo[i];
- }
- else
- return 0;
-}
-
-////////////////////////////////////////////////////////////////////////
-AliGenInfoMaker::~AliGenInfoMaker()
-{
-
- if (fLoader){
- fLoader->UnloadgAlice();
- gAlice = 0;
- delete fLoader;
- }
-}
-
-Int_t AliGenInfoMaker::SetIO()
-{
- //
- //
- CreateTreeGenTracks();
- if (!fTreeGenTracks) return 1;
- // AliTracker::SetFieldFactor();
-
- fParamTPC = GetTPCParam();
- //
- return 0;
-}
-
-////////////////////////////////////////////////////////////////////////
-Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
-{
- //
- //
- // SET INPUT
- fLoader->SetEventNumber(eventNr);
- //
- fLoader->LoadHeader();
- fLoader->LoadKinematics();
- fStack = fLoader->Stack();
- //
- fLoader->LoadTrackRefs();
- fLoader->LoadHits();
- fTreeTR = fLoader->TreeTR();
- //
- AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
- tpcl->LoadDigits();
- fTreeD = tpcl->TreeD();
- return 0;
-}
-
-Int_t AliGenInfoMaker::CloseIOEvent()
-{
- fLoader->UnloadHeader();
- fLoader->UnloadKinematics();
- fLoader->UnloadTrackRefs();
- AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
- tpcl->UnloadDigits();
- return 0;
-}
-
-Int_t AliGenInfoMaker::CloseIO()
-{
- fLoader->UnloadgAlice();
- return 0;
-}
-
-
-
-////////////////////////////////////////////////////////////////////////
-Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
-{
- fNEvents = nEvents;
- fFirstEventNr = firstEventNr;
- return Exec();
-}
-
-////////////////////////////////////////////////////////////////////////
-Int_t AliGenInfoMaker::Exec()
-{
- TStopwatch timer;
- timer.Start();
- Int_t status =SetIO();
- if (status>0) return status;
- //
-
- for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
- fEventNr++) {
- SetIO(fEventNr);
- fNParticles = fStack->GetNtrack();
- //
- fGenInfo = new AliMCInfo*[fNParticles];
- for (UInt_t i = 0; i<fNParticles; i++) {
- fGenInfo[i]=0;
- }
- //
- cout<<"Start to process event "<<fEventNr<<endl;
- cout<<"\tfNParticles = "<<fNParticles<<endl;
- if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
- if (TreeTRLoop()>0) return 1;
- //
- if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
- if (TreeDLoop()>0) return 1;
- //
- if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
- if (TreeKLoop()>0) return 1;
- if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
- //
- if (BuildKinkInfo()>0) return 1;
- if (BuildV0Info()>0) return 1;
- //if (BuildHitLines()>0) return 1;
- if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
- //
- for (UInt_t i = 0; i<fNParticles; i++) {
- if (fGenInfo[i]) delete fGenInfo[i];
- }
- delete []fGenInfo;
- CloseIOEvent();
- }
- //
- CloseIO();
- CloseOutputFile();
-
- cerr<<"Exec finished"<<endl;
-
- timer.Stop();
- timer.Print();
- return 0;
-}
-
-////////////////////////////////////////////////////////////////////////
-void AliGenInfoMaker::CreateTreeGenTracks()
-{
- fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
- if (!fFileGenTracks) {
- cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
- return;
- }
- fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
- AliMCInfo * info = new AliMCInfo;
- fTreeGenTracks->Branch("MC","AliMCInfo",&info);
- delete info;
- //
- AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
- fTreeKinks = new TTree("genKinksTree","genKinksTree");
- fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
- delete kinkinfo;
- //
- AliGenV0Info *v0info = new AliGenV0Info;
- fTreeV0 = new TTree("genV0Tree","genV0Tree");
- fTreeV0->Branch("MC","AliGenV0Info",&v0info);
- delete v0info;
- //
- //
- AliTrackPointArray * points0 = new AliTrackPointArray;
- AliTrackPointArray * points1 = new AliTrackPointArray;
- AliTrackPointArray * points2 = new AliTrackPointArray;
- fTreeHitLines = new TTree("HitLines","HitLines");
- fTreeHitLines->Branch("TPC.","AliTrackPointArray",&points0,32000,10);
- fTreeHitLines->Branch("TRD.","AliTrackPointArray",&points1,32000,10);
- fTreeHitLines->Branch("ITS.","AliTrackPointArray",&points2,32000,10);
- Int_t label=0;
- fTreeHitLines->Branch("Label",&label,"label/I");
- //
- fTreeGenTracks->AutoSave();
- fTreeKinks->AutoSave();
- fTreeV0->AutoSave();
- fTreeHitLines->AutoSave();
-}
-////////////////////////////////////////////////////////////////////////
-void AliGenInfoMaker::CloseOutputFile()
-{
- if (!fFileGenTracks) {
- cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
- return;
- }
- fFileGenTracks->cd();
- fTreeGenTracks->Write();
- delete fTreeGenTracks;
- fTreeKinks->Write();
- delete fTreeKinks;
- fTreeV0->Write();
- delete fTreeV0;
-
- fFileGenTracks->Close();
- delete fFileGenTracks;
- return;
-}
-
-////////////////////////////////////////////////////////////////////////
-Int_t AliGenInfoMaker::TreeKLoop()
-{
-//
-// open the file with treeK
-// loop over all entries there and save information about some tracks
-//
-
- AliStack * stack = fStack;
- if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
-
- if (fDebug > 0) {
- cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
- <<fEventNr<<endl;
- }
- Int_t ipdg = 0;
- TParticlePDG *ppdg = 0;
- // not all generators give primary vertex position. Take the vertex
- // of the particle 0 as primary vertex.
- TDatabasePDG pdg; //get pdg table
- //thank you very much root for this
- TBranch * br = fTreeGenTracks->GetBranch("MC");
- TParticle *particle = stack->ParticleFromTreeK(0);
- fVPrim[0] = particle->Vx();
- fVPrim[1] = particle->Vy();
- fVPrim[2] = particle->Vz();
- for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
- // load only particles with TR
- AliMCInfo * info = GetInfo(iParticle);
- if (!info) continue;
- //////////////////////////////////////////////////////////////////////
- info->fLabel = iParticle;
- //
- info->fParticle = *(stack->Particle(iParticle));
- info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
- info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
- info->fVDist[2] = info->fParticle.Vz()-fVPrim[2];
- info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
- info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
- //
- //
- ipdg = info->fParticle.GetPdgCode();
- info->fPdg = ipdg;
- ppdg = pdg.GetParticle(ipdg);
- info->fEventNr = fEventNr;
- info->Update();
- //////////////////////////////////////////////////////////////////////
- br->SetAddress(&info);
- fTreeGenTracks->Fill();
- }
- fTreeGenTracks->AutoSave();
- if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
- return 0;
-}
-
-
-////////////////////////////////////////////////////////////////////////////////////
-Int_t AliGenInfoMaker::BuildKinkInfo()
-{
- //
- // Build tree with Kink Information
- //
- AliStack * stack = fStack;
- if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
-
- if (fDebug > 0) {
- cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
- <<fEventNr<<endl;
- }
- // Int_t ipdg = 0;
- //TParticlePDG *ppdg = 0;
- // not all generators give primary vertex position. Take the vertex
- // of the particle 0 as primary vertex.
- TDatabasePDG pdg; //get pdg table
- //thank you very much root for this
- TBranch * br = fTreeKinks->GetBranch("MC");
- //
- AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
- //
- for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
- // load only particles with TR
- AliMCInfo * info = GetInfo(iParticle);
- if (!info) continue;
- if (info->fCharge==0) continue;
- if (info->fTRdecay.P()<0.13) continue; //momenta cut
- if (info->fTRdecay.R()>500) continue; //R cut - decay outside barrel
- TParticle & particle = info->fParticle;
- Int_t first = particle.GetDaughter(0);
- if (first ==0) continue;
-
- Int_t last = particle.GetDaughter(1);
- if (last ==0) last=first;
- AliMCInfo * info2 = 0;
- AliMCInfo * dinfo = 0;
-
- for (Int_t id2=first;id2<=last;id2++){
- info2 = GetInfo(id2);
- if (!info2) continue;
- TParticle &p2 = info2->fParticle;
- Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
- if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
- if (!(p2.GetPDG())) continue;
- dinfo =info2;
-
- kinkinfo->fMCm = (*info);
- kinkinfo->fMCd = (*dinfo);
- if (kinkinfo->fMCm.fParticle.GetPDG()==0 || kinkinfo->fMCd.fParticle.GetPDG()==0) continue;
- kinkinfo->Update();
- br->SetAddress(&kinkinfo);
- fTreeKinks->Fill();
- }
- /*
- if (dinfo){
- kinkinfo->fMCm = (*info);
- kinkinfo->fMCd = (*dinfo);
- kinkinfo->Update();
- br->SetAddress(&kinkinfo);
- fTreeKinks->Fill();
- }
- */
- }
- fTreeGenTracks->AutoSave();
- if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
- return 0;
-}
-
-
-////////////////////////////////////////////////////////////////////////////////////
-Int_t AliGenInfoMaker::BuildV0Info()
+//_____________________________________________________________________________
+Float_t AliMCInfo::TPCBetheBloch(Float_t bg)
{
//
- // Build tree with V0 Information
- //
- AliStack * stack = fStack;
- if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
-
- if (fDebug > 0) {
- cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
- <<fEventNr<<endl;
- }
- // Int_t ipdg = 0;
- //TParticlePDG *ppdg = 0;
- // not all generators give primary vertex position. Take the vertex
- // of the particle 0 as primary vertex.
- TDatabasePDG pdg; //get pdg table
- //thank you very much root for this
- TBranch * br = fTreeV0->GetBranch("MC");
- //
- AliGenV0Info * v0info = new AliGenV0Info;
- //
+ // Bethe-Bloch energy loss formula
//
- for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
- // load only particles with TR
- AliMCInfo * info = GetInfo(iParticle);
- if (!info) continue;
- if (info->fCharge==0) continue;
- //
- //
- TParticle & particle = info->fParticle;
- Int_t mother = particle.GetMother(0);
- if (mother <=0) continue;
- //
- TParticle * motherparticle = stack->Particle(mother);
- if (!motherparticle) continue;
- //
- Int_t last = motherparticle->GetDaughter(1);
- if (last==(int)iParticle) continue;
- AliMCInfo * info2 = info;
- AliMCInfo * dinfo = GetInfo(last);
- if (!dinfo) continue;
- if (!dinfo->fParticle.GetPDG()) continue;
- if (!info2->fParticle.GetPDG()) continue;
- if (dinfo){
- v0info->fMCm = (*info);
- v0info->fMCd = (*dinfo);
- v0info->fMotherP = (*motherparticle);
- v0info->Update(fVPrim);
- br->SetAddress(&v0info);
- fTreeV0->Fill();
- }
- }
- fTreeV0->AutoSave();
- if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
- return 0;
-}
-
-
-
-
-////////////////////////////////////////////////////////////////////////
-Int_t AliGenInfoMaker::BuildHitLines()
-{
+ const Double_t kp1=0.76176e-1;
+ const Double_t kp2=10.632;
+ const Double_t kp3=0.13279e-4;
+ const Double_t kp4=1.8631;
+ const Double_t kp5=1.9479;
-//
-// open the file with treeK
-// loop over all entries there and save information about some tracks
-//
+ Double_t dbg = (Double_t) bg;
- AliStack * stack = fStack;
- if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
-
- if (fDebug > 0) {
- cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
- <<fEventNr<<endl;
- }
-// Int_t ipdg = 0;
-// // TParticlePDG *ppdg = 0;
-// // not all generators give primary vertex position. Take the vertex
-// // of the particle 0 as primary vertex.
-// TDatabasePDG pdg; //get pdg table
-// //thank you very much root for this
-// AliTrackPointArray *tpcp = new AliTrackPointArray;
-// AliTrackPointArray *trdp = new AliTrackPointArray;
-// AliTrackPointArray *itsp = new AliTrackPointArray;
-// Int_t label =0;
-// //
-// TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
-// TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");
-// TBranch * brits = fTreeHitLines->GetBranch("ITS.");
-// TBranch * brlabel = fTreeHitLines->GetBranch("Label");
-// brlabel->SetAddress(&label);
-// brtpc->SetAddress(&tpcp);
-// brtrd->SetAddress(&trdp);
-// brits->SetAddress(&itsp);
-// //
-// AliDetector *dtpc = gAlice->GetDetector("TPC");
-// AliDetector *dtrd = gAlice->GetDetector("TRD");
-// AliDetector *dits = gAlice->GetDetector("ITS");
-
-// for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
-// // load only particles with TR
-// AliMCInfo * info = GetInfo(iParticle);
-// if (!info) continue;
-// Int_t primpart = info->fPrimPart;
-// ipdg = info->fParticle.GetPdgCode();
-// label = iParticle;
-// //
-// gAlice->ResetHits();
-// tpcp->Reset();
-// itsp->Reset();
-// trdp->Reset();
-// tpcp->fLabel1 = ipdg;
-// trdp->fLabel1 = ipdg;
-// itsp->fLabel1 = ipdg;
-// if (dtpc->TreeH()->GetEvent(primpart)){
-// dtpc->LoadPoints(primpart);
-// tpcp->Reset(dtpc,iParticle);
-// }
-// if (dtrd->TreeH()->GetEvent(primpart)){
-// dtrd->LoadPoints(primpart);
-// trdp->Reset(dtrd,iParticle);
-// }
-// if (dits->TreeH()->GetEvent(primpart)){
-// dits->LoadPoints(primpart);
-// itsp->Reset(dits,iParticle);
-// }
-// //
-// fTreeHitLines->Fill();
-// dtpc->ResetPoints();
-// dtrd->ResetPoints();
-// dits->ResetPoints();
-// }
-// delete tpcp;
-// delete trdp;
-// delete itsp;
-// fTreeHitLines->AutoSave();
-// if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
- return 0;
-}
+ Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
+ Double_t aa = TMath::Power(beta,kp4);
+ Double_t bb = TMath::Power(1./dbg,kp5);
-////////////////////////////////////////////////////////////////////////
-Int_t AliGenInfoMaker::TreeDLoop()
-{
- //
- // open the file with treeD
- // loop over all entries there and save information about some tracks
- //
-
- Int_t nInnerSector = fParamTPC->GetNInnerSector();
- Int_t rowShift = 0;
- Int_t zero=fParamTPC->GetZeroSup()+6;
- //
- //
- AliSimDigits digitsAddress, *digits=&digitsAddress;
- fTreeD->GetBranch("Segment")->SetAddress(&digits);
-
- Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
- if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
- for (Int_t i=0; i<sectorsByRows; i++) {
- if (!fTreeD->GetEvent(i)) continue;
- Int_t sec,row;
- fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
- if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
- // here I expect that upper sectors follow lower sectors
- if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
- //
- digits->ExpandTrackBuffer();
- digits->First();
- do {
- Int_t iRow=digits->CurrentRow();
- Int_t iColumn=digits->CurrentColumn();
- Short_t digitValue = digits->CurrentDigit();
- if (digitValue >= zero) {
- Int_t label;
- for (Int_t j = 0; j<3; j++) {
- // label = digits->GetTrackID(iRow,iColumn,j);
- label = digits->GetTrackIDFast(iRow,iColumn,j)-2;
- if (label >= (Int_t)fNParticles) { //don't label from bakground event
- continue;
- }
- if (label >= 0 && label <= (Int_t)fNParticles) {
- if (fDebug > 6 ) {
- cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
- <<sec<<" "
- <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
- <<" "<<row<<endl;
- }
- AliMCInfo * info = GetInfo(label);
- if (info){
- info->fTPCRow.SetRow(row+rowShift);
- }
- }
- }
- }
- } while (digits->Next());
- }
+ bb=TMath::Log(kp3+bb);
- if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
- return 0;
-}
-
-
-////////////////////////////////////////////////////////////////////////
-Int_t AliGenInfoMaker::TreeTRLoop()
-{
- //
- // loop over TrackReferences and store the first one for each track
- //
- TTree * treeTR = fTreeTR;
- Int_t nPrimaries = (Int_t) treeTR->GetEntries();
- if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
- //
- //
- //track references for TPC
- TClonesArray* tpcArrayTR = new TClonesArray("AliTrackReference");
- TClonesArray* itsArrayTR = new TClonesArray("AliTrackReference");
- TClonesArray* trdArrayTR = new TClonesArray("AliTrackReference");
- TClonesArray* tofArrayTR = new TClonesArray("AliTrackReference");
- TClonesArray* runArrayTR = new TClonesArray("AliTrackReference");
- //
- if (treeTR->GetBranch("TPC")) treeTR->GetBranch("TPC")->SetAddress(&tpcArrayTR);
- if (treeTR->GetBranch("ITS")) treeTR->GetBranch("ITS")->SetAddress(&itsArrayTR);
- if (treeTR->GetBranch("TRD")) treeTR->GetBranch("TRD")->SetAddress(&trdArrayTR);
- if (treeTR->GetBranch("TOF")) treeTR->GetBranch("TOF")->SetAddress(&tofArrayTR);
- if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&runArrayTR);
- //
- //
- //
- for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
- treeTR->GetEntry(iPrimPart);
- //
- // Loop over TPC references
- //
- for (Int_t iTrackRef = 0; iTrackRef < tpcArrayTR->GetEntriesFast(); iTrackRef++) {
- AliTrackReference *trackRef = (AliTrackReference*)tpcArrayTR->At(iTrackRef);
- //
- if (trackRef->TestBit(BIT(2))){
- //if decay
- if (trackRef->P()<fTPCPtCut) continue;
- Int_t label = trackRef->GetTrack();
- AliMCInfo * info = GetInfo(label);
- if (!info) info = MakeInfo(label);
- info->fTRdecay = *trackRef;
- }
- //
- if (trackRef->P()<fTPCPtCut) continue;
- Int_t label = trackRef->GetTrack();
- AliMCInfo * info = GetInfo(label);
- if (!info) info = MakeInfo(label);
- if (!info) continue;
- info->fPrimPart = iPrimPart;
- TClonesArray & arr = *(info->fTPCReferences);
- new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
- }
- //
- // Loop over ITS references
- //
- for (Int_t iTrackRef = 0; iTrackRef < itsArrayTR->GetEntriesFast(); iTrackRef++) {
- AliTrackReference *trackRef = (AliTrackReference*)itsArrayTR->At(iTrackRef);
- //
- //
- if (trackRef->P()<fTPCPtCut) continue;
- Int_t label = trackRef->GetTrack();
- AliMCInfo * info = GetInfo(label);
- if ( (!info) && trackRef->Pt()<fITSPtCut) continue;
- if (!info) info = MakeInfo(label);
- if (!info) continue;
- info->fPrimPart = iPrimPart;
- TClonesArray & arr = *(info->fITSReferences);
- new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
- }
- //
- // Loop over TRD references
- //
- for (Int_t iTrackRef = 0; iTrackRef < trdArrayTR->GetEntriesFast(); iTrackRef++) {
- AliTrackReference *trackRef = (AliTrackReference*)trdArrayTR->At(iTrackRef);
- //
- if (trackRef->P()<fTPCPtCut) continue;
- Int_t label = trackRef->GetTrack();
- AliMCInfo * info = GetInfo(label);
- if ( (!info) && trackRef->Pt()<fTRDPtCut) continue;
- if (!info) info = MakeInfo(label);
- if (!info) continue;
- info->fPrimPart = iPrimPart;
- TClonesArray & arr = *(info->fTRDReferences);
- new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
- }
- //
- // Loop over TOF references
- //
- for (Int_t iTrackRef = 0; iTrackRef < tofArrayTR->GetEntriesFast(); iTrackRef++) {
- AliTrackReference *trackRef = (AliTrackReference*)tofArrayTR->At(iTrackRef);
- Int_t label = trackRef->GetTrack();
- AliMCInfo * info = GetInfo(label);
- if (!info){
- if (trackRef->Pt()<fTPCPtCut) continue;
- if ( (!info) && trackRef->Pt()<fTOFPtCut) continue;
- }
- if (!info) info = MakeInfo(label);
- if (!info) continue;
- info->fPrimPart = iPrimPart;
- TClonesArray & arr = *(info->fTOFReferences);
- new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
- }
- //
- // get dacay position
- for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
- AliTrackReference *trackRef = (AliTrackReference*)runArrayTR->At(iTrackRef);
- //
- Int_t label = trackRef->GetTrack();
- AliMCInfo * info = GetInfo(label);
- if (!info) continue;
- if (!trackRef->TestBit(BIT(2))) continue; //if not decay
- // if (TMath::Abs(trackRef.X());
- info->fTRdecay = *trackRef;
- }
- }
- //
- tpcArrayTR->Delete();
- delete tpcArrayTR;
- trdArrayTR->Delete();
- delete trdArrayTR;
- tofArrayTR->Delete();
- delete tofArrayTR;
- itsArrayTR->Delete();
- delete itsArrayTR;
- runArrayTR->Delete();
- delete runArrayTR;
- //
- return 0;
-}
-
-////////////////////////////////////////////////////////////////////////
-Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
- AliTPCParam *paramTPC) const {
-
- Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
- Int_t index[4];
- paramTPC->Transform0to1(x,index);
- paramTPC->Transform1to2(x,index);
- return x[0];
+ return ((Float_t)((kp2-aa-bb)*kp1/aa));
}
-////////////////////////////////////////////////////////////////////////
-
-
Int_t First() const ;
void Reset();
-//private:
+private:
UChar_t fDig[kgRowBytes];
ClassDef(AliTPCdigitRow,1) // container for digit pattern
};
////////////////////////////////////////////////////////////////////////
class AliMCInfo: public TObject {
-
+ friend class AliGenInfoMaker;
+ friend class AliRecInfoMaker;
+ friend class AliESDRecInfo;
public:
AliMCInfo();
~AliMCInfo();
AliMCInfo(const AliMCInfo& info);
void Update();
-
-
+ Int_t GetEventNr() const {return fEventNr;}
+ const AliTrackReference& GetTrackRef() const {return fTrackRef;}
+ const AliTrackReference& GetTrackRefOut() const {return fTrackRefOut;}
+ const AliTrackReference& GetTRdecay() const {return fTRdecay;}
+ TParticle& GetParticle() {return fParticle;}
+ Float_t TPCBetheBloch(Float_t bg);
+ //
+ Int_t GetPrimPart() const {return fPrimPart;}
+ Float_t GetMass() const {return fMass;}
+ Float_t GetCharge() const {return fCharge;}
+ Int_t GetLabel() const {return fLabel;}
+
+ Int_t GetMCtracks() const {return fMCtracks;}
+ Int_t GetPdg() const {return fPdg;}
+ const Float_t* GetDecayCoord() const {return fDecayCoord;}
+ const Double_t* GetVDist() const {return fVDist;}
+
+ Bool_t IsTPCdecay() const {return fTPCdecay;}
+
+ Int_t GetRowsWithDigitsInn() const {return fRowsWithDigitsInn;}
+ Int_t GetRowsWithDigits() const {return fRowsWithDigits;}
+ Int_t GetRowsTrackLength() const {return fRowsTrackLength;}
+ Float_t GetPrim() const { return fPrim;}
+
+ AliTPCdigitRow & GetTPCRow() {return fTPCRow;}
+ Int_t GetNTPCRef() const {return fNTPCRef;}
+ Int_t GetNITSRef() const {return fNITSRef;}
+ Int_t GetNTRDRef() const {return fNTRDRef;}
+ Int_t GetNTOFRef() const {return fNTOFRef;}
+ const TClonesArray *GetTPCReferences() const { return fTPCReferences;}
+ const TClonesArray * GetTRDReferences() const { return fTRDReferences;}
+ const TClonesArray * GetITSReferences() const { return fITSReferences;}
+ const TClonesArray * GetTOFReferences() const { return fTOFReferences;}
+private:
AliTrackReference fTrackRef; // track reference saved in the output tree
AliTrackReference fTrackRefOut; // decay track reference saved in the output tree
AliTrackReference fTRdecay; // track reference at decay point
TClonesArray * fTRDReferences; //container with TRD references
TClonesArray * fTOFReferences; //container with TRD references
//
- ClassDef(AliMCInfo,1) // container for
+ ClassDef(AliMCInfo,1); // container for
};
class AliGenV0Info: public TObject {
public:
AliGenV0Info(); //
- AliMCInfo fMCd; //info about daughter particle - second particle for V0
- AliMCInfo fMCm; //info about mother particle - first particle for V0
- TParticle fMotherP; //particle info about mother particle
- void Update(Float_t vertex[3]); // put some derived info to special field
+ void Update(Float_t vertex[3]);
+ AliMCInfo & GetPlus() {return fMCd;}
+ AliMCInfo & GetMinus() {return fMCm;}
+ TParticle & GetMopther() {return fMotherP;}
+ Double_t GetMCDist1() const { return fMCDist1;}
+ Double_t GetMCDist2() const {return fMCDist2;}
+ const Double_t* GetMCPdr() const {return fMCPdr;}
+ const Double_t* GetMCPd() const {return fMCPd;}
+ const Double_t* GetMCX() const {return fMCX;}
+ // const Double_t fMCXr;
+ //
+// Double_t fMCPm[3];
+// Double_t fMCAngle[3];
+// Double_t fMCRr;
+// Double_t fMCR;
+// Int_t fPdg[2];
+// Int_t fLab[2];
+// //
+// Double_t fInvMass;
+// Float_t fPointAngleFi;
+// Float_t fPointAngleTh;
+// Float_t fPointAngle;
+
+ void SetInfoP(AliMCInfo &plus) {fMCd=plus;}
+ void SetInfoM(AliMCInfo &minus){fMCm=minus;}
+ void SetMother(TParticle&mother){fMotherP=mother;}
+private:
+ AliMCInfo fMCd; //info about daughter particle - second particle for V0
+ AliMCInfo fMCm; //info about mother particle - first particle for V0
+ TParticle fMotherP; //particle info about mother particle
Double_t fMCDist1; //info about closest distance according closest MC - linear DCA
Double_t fMCDist2; //info about closest distance parabolic DCA
//
- Double_t fMCPdr[3]; //momentum at vertex daughter - according approx at DCA
- Double_t fMCPd[4]; //exact momentum from MC info
- Double_t fMCX[3]; //exact position of the vertex
- Double_t fMCXr[3]; //rec. position according helix
+ Double_t fMCPdr[3]; //momentum at vertex daughter - according approx at DCA
+ Double_t fMCPd[4]; //exact momentum from MC info
+ Double_t fMCX[3]; //exact position of the vertex
+ Double_t fMCXr[3]; //rec. position according helix
//
Double_t fMCPm[3]; //momentum at the vertex mother
Double_t fMCAngle[3]; //three angels
AliGenKinkInfo(); //default cosntructor
void Update(); // put some derived info to special field
Float_t GetQt(); //
-
+ AliMCInfo & GetPlus() {return fMCd;}
+ AliMCInfo & GetMinus() {return fMCm;}
+ void SetInfoDaughter(AliMCInfo &daughter) {fMCd=daughter;}
+ void SetInfoMother(AliMCInfo &mother){fMCm=mother;}
+private:
AliMCInfo fMCd; //info about daughter particle - second particle for V0
AliMCInfo fMCm; //info about mother particle - first particle for V0
Double_t fMCDist1; //info about closest distance according closest MC - linear DCA
ClassDef(AliGenKinkInfo,1) // container for
};
-
-
-////////////////////////////////////////////////////////////////////////
-//
-// Start of implementation of the class AliGenInfoMaker
-//
-////////////////////////////////////////////////////////////////////////
-
-class AliGenInfoMaker {
-
-public:
- AliGenInfoMaker();
- AliGenInfoMaker(const char * fnGalice, const char* fnRes ="genTracks.root",
- Int_t nEvents=1, Int_t firstEvent=0);
- virtual ~AliGenInfoMaker();
- Int_t Exec();
- Int_t Exec(Int_t nEvents, Int_t firstEventNr);
- void CreateTreeGenTracks();
- void CloseOutputFile();
- Int_t TreeKLoop();
- Int_t TreeTRLoop();
- Int_t TreeDLoop();
- Int_t BuildKinkInfo(); // build information about MC kinks
- Int_t BuildV0Info(); // build information about MC kinks
- Int_t BuildHitLines(); // build information about MC kinks
- void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
- void SetNEvents(Int_t i) {fNEvents = i;}
- void SetDebug(Int_t level) {fDebug = level;}
- Int_t SetIO(Int_t eventNr);
- Int_t CloseIOEvent();
- Int_t CloseIO();
- Int_t SetIO();
- Float_t TR2LocalX(AliTrackReference *trackRef,
- AliTPCParam *paramTPC) const;
- AliMCInfo * GetInfo(UInt_t i) const {return (i<fNParticles)? fGenInfo[i]:0;}
- AliMCInfo * MakeInfo(UInt_t i);
-
-private:
- Int_t fDebug; //! debug flag
- Int_t fEventNr; //! current event number
- Int_t fLabel; //! track label
- Int_t fNEvents; //! number of events to process
- Int_t fFirstEventNr; //! first event to process
- UInt_t fNParticles; //! number of particles in TreeK
- TTree *fTreeGenTracks; //! output tree with generated tracks
- TTree *fTreeKinks; //! output tree with Kinks
- TTree *fTreeV0; //! output tree with V0
- TTree *fTreeHitLines; //! tree with hit lines
- char fFnRes[1000]; //! output file name with stored tracks
- TFile *fFileGenTracks; //! output file with stored fTreeGenTracks
- //
- AliRunLoader * fLoader; //! pointer to the run loader
- TTree * fTreeD; //! current tree with digits
- TTree * fTreeTR; //! current tree with TR
- AliStack *fStack; //! current stack
- //
- AliMCInfo ** fGenInfo; //! array with pointers to gen info
- Int_t fNInfos; //! number of tracks with infos
- //
- AliTPCParam* fParamTPC; //! AliTPCParam
- Float_t fVPrim[3]; //! primary vertex position
- // the fVDist[3] contains size of the 3-vector
- // cuts
- //
- Double_t fTPCPtCut; // do not store particles with generated pT less than this
- Double_t fITSPtCut; // do not store particles with generated pT less than this
- Double_t fTRDPtCut; // do not store particles with generated pT less than this
- Double_t fTOFPtCut; // do not store particles with generated pT less than this
-
- ClassDef(AliGenInfoMaker,0) // class which creates and fills tree with TPCGenTrack objects
-};
-
-
-
-
-
-AliTPCParam * GetTPCParam();
-Float_t TPCBetheBloch(Float_t bg);
-
#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+///////////////////////////////////////////////////////////////////////////
+/*
+
+Origin: marian.ivanov@cern.ch
+
+Generate complex MC information - used for Comparison later on
+How to use it?
+
+gSystem->Load("libPWG1.so")
+AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",0,0)
+t->Exec();
+
+*/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <stdio.h>
+#include <string.h>
+//ROOT includes
+#include "TROOT.h"
+#include "Rtypes.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "TString.h"
+#include "TStopwatch.h"
+#include "TParticle.h"
+#include "TSystem.h"
+#include "TCanvas.h"
+#include "TGeometry.h"
+#include "TPolyLine3D.h"
+
+//ALIROOT includes
+#include "AliRun.h"
+#include "AliStack.h"
+#include "AliSimDigits.h"
+#include "AliTPCParam.h"
+#include "AliTPC.h"
+#include "AliTPCLoader.h"
+#include "AliDetector.h"
+#include "AliTrackReference.h"
+#include "AliTPCParamSR.h"
+#include "AliTracker.h"
+#include "AliMagF.h"
+#include "AliHelix.h"
+#include "AliTrackPointArray.h"
+
+#endif
+#include "AliGenInfo.h"
+#include "AliGenInfoMaker.h"
+//
+//
+
+ClassImp(AliGenInfoMaker)
+
+
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+AliGenInfoMaker::AliGenInfoMaker():
+ fDebug(0), //! debug flag
+ fEventNr(0), //! current event number
+ fLabel(0), //! track label
+ fNEvents(0), //! number of events to process
+ fFirstEventNr(0), //! first event to process
+ fNParticles(0), //! number of particles in TreeK
+ fTreeGenTracks(0), //! output tree with generated tracks
+ fTreeKinks(0), //! output tree with Kinks
+ fTreeV0(0), //! output tree with V0
+ fTreeHitLines(0), //! tree with hit lines
+ fFileGenTracks(0), //! output file with stored fTreeGenTracks
+ fLoader(0), //! pointer to the run loader
+ fTreeD(0), //! current tree with digits
+ fTreeTR(0), //! current tree with TR
+ fStack(0), //! current stack
+ fGenInfo(0), //! array with pointers to gen info
+ fNInfos(0), //! number of tracks with infos
+ fParamTPC(0), //! AliTPCParam
+ fTPCPtCut(0.03),
+ fITSPtCut(0.1),
+ fTRDPtCut(0.1),
+ fTOFPtCut(0.1)
+{
+}
+
+////////////////////////////////////////////////////////////////////////
+AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
+ Int_t nEvents, Int_t firstEvent):
+ fDebug(0), //! debug flag
+ fEventNr(0), //! current event number
+ fLabel(0), //! track label
+ fNEvents(0), //! number of events to process
+ fFirstEventNr(0), //! first event to process
+ fNParticles(0), //! number of particles in TreeK
+ fTreeGenTracks(0), //! output tree with generated tracks
+ fTreeKinks(0), //! output tree with Kinks
+ fTreeV0(0), //! output tree with V0
+ fTreeHitLines(0), //! tree with hit lines
+ fFileGenTracks(0), //! output file with stored fTreeGenTracks
+ fLoader(0), //! pointer to the run loader
+ fTreeD(0), //! current tree with digits
+ fTreeTR(0), //! current tree with TR
+ fStack(0), //! current stack
+ fGenInfo(0), //! array with pointers to gen info
+ fNInfos(0), //! number of tracks with infos
+ fParamTPC(0), //! AliTPCParam
+ fTPCPtCut(0.03),
+ fITSPtCut(0.1),
+ fTRDPtCut(0.1),
+ fTOFPtCut(0.1)
+
+{
+ //
+ //
+ //
+ fFirstEventNr = firstEvent;
+ fEventNr = firstEvent;
+ fNEvents = nEvents;
+ sprintf(fFnRes,"%s",fnRes);
+ //
+ fLoader = AliRunLoader::Open(fnGalice);
+ if (gAlice){
+ delete gAlice->GetRunLoader();
+ delete gAlice;
+ gAlice = 0x0;
+ }
+ if (fLoader->LoadgAlice()){
+ cerr<<"Error occured while l"<<endl;
+ }
+ Int_t nall = fLoader->GetNumberOfEvents();
+ if (nEvents==0) {
+ nEvents =nall;
+ fNEvents=nall;
+ fFirstEventNr=0;
+ }
+
+ if (nall<=0){
+ cerr<<"no events available"<<endl;
+ fEventNr = 0;
+ return;
+ }
+ if (firstEvent+nEvents>nall) {
+ fEventNr = nall-firstEvent;
+ cerr<<"restricted number of events availaible"<<endl;
+ }
+ AliMagF * magf = gAlice->Field();
+ AliTracker::SetFieldMap(magf,0);
+}
+
+
+AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
+{
+ //
+ if (i<fNParticles) {
+ if (fGenInfo[i]) return fGenInfo[i];
+ fGenInfo[i] = new AliMCInfo;
+ fNInfos++;
+ return fGenInfo[i];
+ }
+ else
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+AliGenInfoMaker::~AliGenInfoMaker()
+{
+
+ if (fLoader){
+ fLoader->UnloadgAlice();
+ gAlice = 0;
+ delete fLoader;
+ }
+}
+
+Int_t AliGenInfoMaker::SetIO()
+{
+ //
+ //
+ CreateTreeGenTracks();
+ if (!fTreeGenTracks) return 1;
+ // AliTracker::SetFieldFactor();
+
+ fParamTPC = GetTPCParam();
+ //
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
+{
+ //
+ //
+ // SET INPUT
+ fLoader->SetEventNumber(eventNr);
+ //
+ fLoader->LoadHeader();
+ fLoader->LoadKinematics();
+ fStack = fLoader->Stack();
+ //
+ fLoader->LoadTrackRefs();
+ fLoader->LoadHits();
+ fTreeTR = fLoader->TreeTR();
+ //
+ AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
+ tpcl->LoadDigits();
+ fTreeD = tpcl->TreeD();
+ return 0;
+}
+
+Int_t AliGenInfoMaker::CloseIOEvent()
+{
+ fLoader->UnloadHeader();
+ fLoader->UnloadKinematics();
+ fLoader->UnloadTrackRefs();
+ AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
+ tpcl->UnloadDigits();
+ return 0;
+}
+
+Int_t AliGenInfoMaker::CloseIO()
+{
+ fLoader->UnloadgAlice();
+ return 0;
+}
+
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
+{
+ fNEvents = nEvents;
+ fFirstEventNr = firstEventNr;
+ return Exec();
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::Exec()
+{
+ TStopwatch timer;
+ timer.Start();
+ Int_t status =SetIO();
+ if (status>0) return status;
+ //
+
+ for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
+ fEventNr++) {
+ SetIO(fEventNr);
+ fNParticles = fStack->GetNtrack();
+ //
+ fGenInfo = new AliMCInfo*[fNParticles];
+ for (UInt_t i = 0; i<fNParticles; i++) {
+ fGenInfo[i]=0;
+ }
+ //
+ cout<<"Start to process event "<<fEventNr<<endl;
+ cout<<"\tfNParticles = "<<fNParticles<<endl;
+ if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
+ if (TreeTRLoop()>0) return 1;
+ //
+ if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
+ if (TreeDLoop()>0) return 1;
+ //
+ if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
+ if (TreeKLoop()>0) return 1;
+ if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
+ //
+ if (BuildKinkInfo()>0) return 1;
+ if (BuildV0Info()>0) return 1;
+ //if (BuildHitLines()>0) return 1;
+ if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
+ //
+ for (UInt_t i = 0; i<fNParticles; i++) {
+ if (fGenInfo[i]) delete fGenInfo[i];
+ }
+ delete []fGenInfo;
+ CloseIOEvent();
+ }
+ //
+ CloseIO();
+ CloseOutputFile();
+
+ cerr<<"Exec finished"<<endl;
+
+ timer.Stop();
+ timer.Print();
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+void AliGenInfoMaker::CreateTreeGenTracks()
+{
+ fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
+ if (!fFileGenTracks) {
+ cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
+ return;
+ }
+ fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
+ AliMCInfo * info = new AliMCInfo;
+ fTreeGenTracks->Branch("MC","AliMCInfo",&info);
+ delete info;
+ //
+ AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
+ fTreeKinks = new TTree("genKinksTree","genKinksTree");
+ fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
+ delete kinkinfo;
+ //
+ AliGenV0Info *v0info = new AliGenV0Info;
+ fTreeV0 = new TTree("genV0Tree","genV0Tree");
+ fTreeV0->Branch("MC","AliGenV0Info",&v0info);
+ delete v0info;
+ //
+ //
+ AliTrackPointArray * points0 = new AliTrackPointArray;
+ AliTrackPointArray * points1 = new AliTrackPointArray;
+ AliTrackPointArray * points2 = new AliTrackPointArray;
+ fTreeHitLines = new TTree("HitLines","HitLines");
+ fTreeHitLines->Branch("TPC.","AliTrackPointArray",&points0,32000,10);
+ fTreeHitLines->Branch("TRD.","AliTrackPointArray",&points1,32000,10);
+ fTreeHitLines->Branch("ITS.","AliTrackPointArray",&points2,32000,10);
+ Int_t label=0;
+ fTreeHitLines->Branch("Label",&label,"label/I");
+ //
+ fTreeGenTracks->AutoSave();
+ fTreeKinks->AutoSave();
+ fTreeV0->AutoSave();
+ fTreeHitLines->AutoSave();
+}
+////////////////////////////////////////////////////////////////////////
+void AliGenInfoMaker::CloseOutputFile()
+{
+ if (!fFileGenTracks) {
+ cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
+ return;
+ }
+ fFileGenTracks->cd();
+ fTreeGenTracks->Write();
+ delete fTreeGenTracks;
+ fTreeKinks->Write();
+ delete fTreeKinks;
+ fTreeV0->Write();
+ delete fTreeV0;
+
+ fFileGenTracks->Close();
+ delete fFileGenTracks;
+ return;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::TreeKLoop()
+{
+//
+// open the file with treeK
+// loop over all entries there and save information about some tracks
+//
+
+ AliStack * stack = fStack;
+ if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
+
+ if (fDebug > 0) {
+ cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
+ <<fEventNr<<endl;
+ }
+ Int_t ipdg = 0;
+ TParticlePDG *ppdg = 0;
+ // not all generators give primary vertex position. Take the vertex
+ // of the particle 0 as primary vertex.
+ TDatabasePDG pdg; //get pdg table
+ //thank you very much root for this
+ TBranch * br = fTreeGenTracks->GetBranch("MC");
+ TParticle *particle = stack->ParticleFromTreeK(0);
+ fVPrim[0] = particle->Vx();
+ fVPrim[1] = particle->Vy();
+ fVPrim[2] = particle->Vz();
+ for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+ // load only particles with TR
+ AliMCInfo * info = GetInfo(iParticle);
+ if (!info) continue;
+ //////////////////////////////////////////////////////////////////////
+ info->fLabel = iParticle;
+ //
+ info->fParticle = *(stack->Particle(iParticle));
+ info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
+ info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
+ info->fVDist[2] = info->fParticle.Vz()-fVPrim[2];
+ info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
+ info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
+ //
+ //
+ ipdg = info->fParticle.GetPdgCode();
+ info->fPdg = ipdg;
+ ppdg = pdg.GetParticle(ipdg);
+ info->fEventNr = fEventNr;
+ info->Update();
+ //////////////////////////////////////////////////////////////////////
+ br->SetAddress(&info);
+ fTreeGenTracks->Fill();
+ }
+ fTreeGenTracks->AutoSave();
+ if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
+ return 0;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::BuildKinkInfo()
+{
+ //
+ // Build tree with Kink Information
+ //
+ AliStack * stack = fStack;
+ if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
+
+ if (fDebug > 0) {
+ cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
+ <<fEventNr<<endl;
+ }
+ // Int_t ipdg = 0;
+ //TParticlePDG *ppdg = 0;
+ // not all generators give primary vertex position. Take the vertex
+ // of the particle 0 as primary vertex.
+ TDatabasePDG pdg; //get pdg table
+ //thank you very much root for this
+ TBranch * br = fTreeKinks->GetBranch("MC");
+ //
+ AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
+ //
+ for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+ // load only particles with TR
+ AliMCInfo * info = GetInfo(iParticle);
+ if (!info) continue;
+ if (info->fCharge==0) continue;
+ if (info->fTRdecay.P()<0.13) continue; //momenta cut
+ if (info->fTRdecay.R()>500) continue; //R cut - decay outside barrel
+ TParticle & particle = info->fParticle;
+ Int_t first = particle.GetDaughter(0);
+ if (first ==0) continue;
+
+ Int_t last = particle.GetDaughter(1);
+ if (last ==0) last=first;
+ AliMCInfo * info2 = 0;
+ AliMCInfo * dinfo = 0;
+
+ for (Int_t id2=first;id2<=last;id2++){
+ info2 = GetInfo(id2);
+ if (!info2) continue;
+ TParticle &p2 = info2->fParticle;
+ Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
+ if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
+ if (!(p2.GetPDG())) continue;
+ dinfo =info2;
+
+ kinkinfo->SetInfoMother(*info);
+ kinkinfo->SetInfoDaughter(*dinfo);
+ if (kinkinfo->GetMinus().GetParticle().GetPDG()==0 || kinkinfo->GetPlus().GetParticle().GetPDG()==0) continue;
+ kinkinfo->Update();
+ br->SetAddress(&kinkinfo);
+ fTreeKinks->Fill();
+ }
+ /*
+ if (dinfo){
+ kinkinfo->fMCm = (*info);
+ kinkinfo->GetPlus() = (*dinfo);
+ kinkinfo->Update();
+ br->SetAddress(&kinkinfo);
+ fTreeKinks->Fill();
+ }
+ */
+ }
+ fTreeGenTracks->AutoSave();
+ if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
+ return 0;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::BuildV0Info()
+{
+ //
+ // Build tree with V0 Information
+ //
+ AliStack * stack = fStack;
+ if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
+
+ if (fDebug > 0) {
+ cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
+ <<fEventNr<<endl;
+ }
+ // Int_t ipdg = 0;
+ //TParticlePDG *ppdg = 0;
+ // not all generators give primary vertex position. Take the vertex
+ // of the particle 0 as primary vertex.
+ TDatabasePDG pdg; //get pdg table
+ //thank you very much root for this
+ TBranch * br = fTreeV0->GetBranch("MC");
+ //
+ AliGenV0Info * v0info = new AliGenV0Info;
+ //
+ //
+ for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+ // load only particles with TR
+ AliMCInfo * info = GetInfo(iParticle);
+ if (!info) continue;
+ if (info->fCharge==0) continue;
+ //
+ //
+ TParticle & particle = info->fParticle;
+ Int_t mother = particle.GetMother(0);
+ if (mother <=0) continue;
+ //
+ TParticle * motherparticle = stack->Particle(mother);
+ if (!motherparticle) continue;
+ //
+ Int_t last = motherparticle->GetDaughter(1);
+ if (last==(int)iParticle) continue;
+ AliMCInfo * info2 = info;
+ AliMCInfo * dinfo = GetInfo(last);
+ if (!dinfo) continue;
+ if (!dinfo->fParticle.GetPDG()) continue;
+ if (!info2->fParticle.GetPDG()) continue;
+ if (dinfo){
+ v0info->SetInfoP(*info);
+ v0info->SetInfoM(*dinfo);
+ v0info->SetMother(*motherparticle);
+ v0info->Update(fVPrim);
+ br->SetAddress(&v0info);
+ fTreeV0->Fill();
+ }
+ }
+ fTreeV0->AutoSave();
+ if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
+ return 0;
+}
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::BuildHitLines()
+{
+
+//
+// open the file with treeK
+// loop over all entries there and save information about some tracks
+//
+
+ AliStack * stack = fStack;
+ if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
+
+ if (fDebug > 0) {
+ cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
+ <<fEventNr<<endl;
+ }
+// Int_t ipdg = 0;
+// // TParticlePDG *ppdg = 0;
+// // not all generators give primary vertex position. Take the vertex
+// // of the particle 0 as primary vertex.
+// TDatabasePDG pdg; //get pdg table
+// //thank you very much root for this
+// AliTrackPointArray *tpcp = new AliTrackPointArray;
+// AliTrackPointArray *trdp = new AliTrackPointArray;
+// AliTrackPointArray *itsp = new AliTrackPointArray;
+// Int_t label =0;
+// //
+// TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
+// TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");
+// TBranch * brits = fTreeHitLines->GetBranch("ITS.");
+// TBranch * brlabel = fTreeHitLines->GetBranch("Label");
+// brlabel->SetAddress(&label);
+// brtpc->SetAddress(&tpcp);
+// brtrd->SetAddress(&trdp);
+// brits->SetAddress(&itsp);
+// //
+// AliDetector *dtpc = gAlice->GetDetector("TPC");
+// AliDetector *dtrd = gAlice->GetDetector("TRD");
+// AliDetector *dits = gAlice->GetDetector("ITS");
+
+// for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+// // load only particles with TR
+// AliMCInfo * info = GetInfo(iParticle);
+// if (!info) continue;
+// Int_t primpart = info->fPrimPart;
+// ipdg = info->fParticle.GetPdgCode();
+// label = iParticle;
+// //
+// gAlice->ResetHits();
+// tpcp->Reset();
+// itsp->Reset();
+// trdp->Reset();
+// tpcp->fLabel1 = ipdg;
+// trdp->fLabel1 = ipdg;
+// itsp->fLabel1 = ipdg;
+// if (dtpc->TreeH()->GetEvent(primpart)){
+// dtpc->LoadPoints(primpart);
+// tpcp->Reset(dtpc,iParticle);
+// }
+// if (dtrd->TreeH()->GetEvent(primpart)){
+// dtrd->LoadPoints(primpart);
+// trdp->Reset(dtrd,iParticle);
+// }
+// if (dits->TreeH()->GetEvent(primpart)){
+// dits->LoadPoints(primpart);
+// itsp->Reset(dits,iParticle);
+// }
+// //
+// fTreeHitLines->Fill();
+// dtpc->ResetPoints();
+// dtrd->ResetPoints();
+// dits->ResetPoints();
+// }
+// delete tpcp;
+// delete trdp;
+// delete itsp;
+// fTreeHitLines->AutoSave();
+// if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
+ return 0;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::TreeDLoop()
+{
+ //
+ // open the file with treeD
+ // loop over all entries there and save information about some tracks
+ //
+
+ Int_t nInnerSector = fParamTPC->GetNInnerSector();
+ Int_t rowShift = 0;
+ Int_t zero=fParamTPC->GetZeroSup()+6;
+ //
+ //
+ AliSimDigits digitsAddress, *digits=&digitsAddress;
+ fTreeD->GetBranch("Segment")->SetAddress(&digits);
+
+ Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
+ if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
+ for (Int_t i=0; i<sectorsByRows; i++) {
+ if (!fTreeD->GetEvent(i)) continue;
+ Int_t sec,row;
+ fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
+ if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
+ // here I expect that upper sectors follow lower sectors
+ if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
+ //
+ digits->ExpandTrackBuffer();
+ digits->First();
+ do {
+ Int_t iRow=digits->CurrentRow();
+ Int_t iColumn=digits->CurrentColumn();
+ Short_t digitValue = digits->CurrentDigit();
+ if (digitValue >= zero) {
+ Int_t label;
+ for (Int_t j = 0; j<3; j++) {
+ // label = digits->GetTrackID(iRow,iColumn,j);
+ label = digits->GetTrackIDFast(iRow,iColumn,j)-2;
+ if (label >= (Int_t)fNParticles) { //don't label from bakground event
+ continue;
+ }
+ if (label >= 0 && label <= (Int_t)fNParticles) {
+ if (fDebug > 6 ) {
+ cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
+ <<sec<<" "
+ <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
+ <<" "<<row<<endl;
+ }
+ AliMCInfo * info = GetInfo(label);
+ if (info){
+ info->fTPCRow.SetRow(row+rowShift);
+ }
+ }
+ }
+ }
+ } while (digits->Next());
+ }
+
+ if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
+ return 0;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::TreeTRLoop()
+{
+ //
+ // loop over TrackReferences and store the first one for each track
+ //
+ TTree * treeTR = fTreeTR;
+ Int_t nPrimaries = (Int_t) treeTR->GetEntries();
+ if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
+ //
+ //
+ //track references for TPC
+ TClonesArray* tpcArrayTR = new TClonesArray("AliTrackReference");
+ TClonesArray* itsArrayTR = new TClonesArray("AliTrackReference");
+ TClonesArray* trdArrayTR = new TClonesArray("AliTrackReference");
+ TClonesArray* tofArrayTR = new TClonesArray("AliTrackReference");
+ TClonesArray* runArrayTR = new TClonesArray("AliTrackReference");
+ //
+ if (treeTR->GetBranch("TPC")) treeTR->GetBranch("TPC")->SetAddress(&tpcArrayTR);
+ if (treeTR->GetBranch("ITS")) treeTR->GetBranch("ITS")->SetAddress(&itsArrayTR);
+ if (treeTR->GetBranch("TRD")) treeTR->GetBranch("TRD")->SetAddress(&trdArrayTR);
+ if (treeTR->GetBranch("TOF")) treeTR->GetBranch("TOF")->SetAddress(&tofArrayTR);
+ if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&runArrayTR);
+ //
+ //
+ //
+ for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
+ treeTR->GetEntry(iPrimPart);
+ //
+ // Loop over TPC references
+ //
+ for (Int_t iTrackRef = 0; iTrackRef < tpcArrayTR->GetEntriesFast(); iTrackRef++) {
+ AliTrackReference *trackRef = (AliTrackReference*)tpcArrayTR->At(iTrackRef);
+ //
+ if (trackRef->TestBit(BIT(2))){
+ //if decay
+ if (trackRef->P()<fTPCPtCut) continue;
+ Int_t label = trackRef->GetTrack();
+ AliMCInfo * info = GetInfo(label);
+ if (!info) info = MakeInfo(label);
+ info->fTRdecay = *trackRef;
+ }
+ //
+ if (trackRef->P()<fTPCPtCut) continue;
+ Int_t label = trackRef->GetTrack();
+ AliMCInfo * info = GetInfo(label);
+ if (!info) info = MakeInfo(label);
+ if (!info) continue;
+ info->fPrimPart = iPrimPart;
+ TClonesArray & arr = *(info->fTPCReferences);
+ new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
+ }
+ //
+ // Loop over ITS references
+ //
+ for (Int_t iTrackRef = 0; iTrackRef < itsArrayTR->GetEntriesFast(); iTrackRef++) {
+ AliTrackReference *trackRef = (AliTrackReference*)itsArrayTR->At(iTrackRef);
+ //
+ //
+ if (trackRef->P()<fTPCPtCut) continue;
+ Int_t label = trackRef->GetTrack();
+ AliMCInfo * info = GetInfo(label);
+ if ( (!info) && trackRef->Pt()<fITSPtCut) continue;
+ if (!info) info = MakeInfo(label);
+ if (!info) continue;
+ info->fPrimPart = iPrimPart;
+ TClonesArray & arr = *(info->fITSReferences);
+ new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
+ }
+ //
+ // Loop over TRD references
+ //
+ for (Int_t iTrackRef = 0; iTrackRef < trdArrayTR->GetEntriesFast(); iTrackRef++) {
+ AliTrackReference *trackRef = (AliTrackReference*)trdArrayTR->At(iTrackRef);
+ //
+ if (trackRef->P()<fTPCPtCut) continue;
+ Int_t label = trackRef->GetTrack();
+ AliMCInfo * info = GetInfo(label);
+ if ( (!info) && trackRef->Pt()<fTRDPtCut) continue;
+ if (!info) info = MakeInfo(label);
+ if (!info) continue;
+ info->fPrimPart = iPrimPart;
+ TClonesArray & arr = *(info->fTRDReferences);
+ new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
+ }
+ //
+ // Loop over TOF references
+ //
+ for (Int_t iTrackRef = 0; iTrackRef < tofArrayTR->GetEntriesFast(); iTrackRef++) {
+ AliTrackReference *trackRef = (AliTrackReference*)tofArrayTR->At(iTrackRef);
+ Int_t label = trackRef->GetTrack();
+ AliMCInfo * info = GetInfo(label);
+ if (!info){
+ if (trackRef->Pt()<fTPCPtCut) continue;
+ if ( (!info) && trackRef->Pt()<fTOFPtCut) continue;
+ }
+ if (!info) info = MakeInfo(label);
+ if (!info) continue;
+ info->fPrimPart = iPrimPart;
+ TClonesArray & arr = *(info->fTOFReferences);
+ new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
+ }
+ //
+ // get dacay position
+ for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
+ AliTrackReference *trackRef = (AliTrackReference*)runArrayTR->At(iTrackRef);
+ //
+ Int_t label = trackRef->GetTrack();
+ AliMCInfo * info = GetInfo(label);
+ if (!info) continue;
+ if (!trackRef->TestBit(BIT(2))) continue; //if not decay
+ // if (TMath::Abs(trackRef.X());
+ info->fTRdecay = *trackRef;
+ }
+ }
+ //
+ tpcArrayTR->Delete();
+ delete tpcArrayTR;
+ trdArrayTR->Delete();
+ delete trdArrayTR;
+ tofArrayTR->Delete();
+ delete tofArrayTR;
+ itsArrayTR->Delete();
+ delete itsArrayTR;
+ runArrayTR->Delete();
+ delete runArrayTR;
+ //
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
+ AliTPCParam *paramTPC) const {
+
+ Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
+ Int_t index[4];
+ paramTPC->Transform0to1(x,index);
+ paramTPC->Transform1to2(x,index);
+ return x[0];
+}
+////////////////////////////////////////////////////////////////////////
+
+
+
+AliTPCParam * AliGenInfoMaker::GetTPCParam(){
+ AliTPCParamSR * par = new AliTPCParamSR;
+ par->Update();
+ return par;
+}
+
+
--- /dev/null
+#ifndef ALIGENINFOMAKER_H
+#define ALIGENINFOMAKER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+
+
+//////////////////////////////////////////////////////////////////////////////
+// Class AliGenInfo //
+// collect together MC info for comparison purposes - effieciency studies and so on// //
+// marian.ivanov@cern.ch //
+//////////////////////////////////////////////////////////////////////////////
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class AliTPCdigitRow
+//
+////////////////////////////////////////////////////////////////////////
+
+#include <TParticle.h>
+#include "AliTrackReference.h"
+
+class TFile;
+class AliRunLoader;
+class AliStack;
+class AliTPCParam;
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class AliGenInfoMaker
+//
+////////////////////////////////////////////////////////////////////////
+
+class AliGenInfoMaker {
+
+public:
+ AliGenInfoMaker();
+ AliGenInfoMaker(const char * fnGalice, const char* fnRes ="genTracks.root",
+ Int_t nEvents=1, Int_t firstEvent=0);
+ virtual ~AliGenInfoMaker();
+ Int_t Exec();
+ Int_t Exec(Int_t nEvents, Int_t firstEventNr);
+ void CreateTreeGenTracks();
+ void CloseOutputFile();
+ Int_t TreeKLoop();
+ Int_t TreeTRLoop();
+ Int_t TreeDLoop();
+ Int_t BuildKinkInfo(); // build information about MC kinks
+ Int_t BuildV0Info(); // build information about MC kinks
+ Int_t BuildHitLines(); // build information about MC kinks
+ void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
+ void SetNEvents(Int_t i) {fNEvents = i;}
+ void SetDebug(Int_t level) {fDebug = level;}
+ Int_t SetIO(Int_t eventNr);
+ Int_t CloseIOEvent();
+ Int_t CloseIO();
+ Int_t SetIO();
+ Float_t TR2LocalX(AliTrackReference *trackRef,
+ AliTPCParam *paramTPC) const;
+ AliMCInfo * GetInfo(UInt_t i) const {return (i<fNParticles)? fGenInfo[i]:0;}
+ AliMCInfo * MakeInfo(UInt_t i);
+
+private:
+ AliTPCParam * GetTPCParam();
+ Float_t TPCBetheBloch(Float_t bg);
+ Int_t fDebug; //! debug flag
+ Int_t fEventNr; //! current event number
+ Int_t fLabel; //! track label
+ Int_t fNEvents; //! number of events to process
+ Int_t fFirstEventNr; //! first event to process
+ UInt_t fNParticles; //! number of particles in TreeK
+ TTree *fTreeGenTracks; //! output tree with generated tracks
+ TTree *fTreeKinks; //! output tree with Kinks
+ TTree *fTreeV0; //! output tree with V0
+ TTree *fTreeHitLines; //! tree with hit lines
+ char fFnRes[1000]; //! output file name with stored tracks
+ TFile *fFileGenTracks; //! output file with stored fTreeGenTracks
+ //
+ AliRunLoader * fLoader; //! pointer to the run loader
+ TTree * fTreeD; //! current tree with digits
+ TTree * fTreeTR; //! current tree with TR
+ AliStack *fStack; //! current stack
+ //
+ AliMCInfo ** fGenInfo; //! array with pointers to gen info
+ Int_t fNInfos; //! number of tracks with infos
+ //
+ AliTPCParam* fParamTPC; //! AliTPCParam
+ Float_t fVPrim[3]; //! primary vertex position
+ // the fVDist[3] contains size of the 3-vector
+ // cuts
+ //
+ Double_t fTPCPtCut; // do not store particles with generated pT less than this
+ Double_t fITSPtCut; // do not store particles with generated pT less than this
+ Double_t fTRDPtCut; // do not store particles with generated pT less than this
+ Double_t fTOFPtCut; // do not store particles with generated pT less than this
+
+ ClassDef(AliGenInfoMaker,0) // class which creates and fills tree with TPCGenTrack objects
+};
+
+
+
+
+
+
+#endif
static TH1F* CreateResHisto(TH2F* hRes2, TH1F **phMean,
Bool_t drawBinFits = kTRUE,Bool_t overflowBinFits = kFALSE);
- AliTreeDraw(const AliTreeDraw& t):fTree(0),fRes(0),fMean(0),fPoints(0){;}
+ AliTreeDraw(const AliTreeDraw& t):TObject(),fTree(0),fRes(0),fMean(0),fPoints(0){;}
AliTreeDraw & operator=(const AliTreeDraw & t){return *this;}
TTree * fTree; //the tree for visualization - NOT OWNER
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff