2 // Script to make correction maps for dndeta measurements using the
3 // dNdEtaCorrection class.
6 makeCorrection(Char_t* dataDir, Int_t nRuns=20) {
10 gSystem->Load("../libPWG0base.so");
12 // ########################################################
13 // selection of esd tracks
14 AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();
15 esdTrackCuts->DefineHistograms(1);
17 esdTrackCuts->SetMinNClustersTPC(50);
18 esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
19 esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
20 esdTrackCuts->SetRequireTPCRefit(kTRUE);
22 esdTrackCuts->SetMinNsigmaToVertex(3);
23 esdTrackCuts->SetAcceptKingDaughters(kFALSE);
25 AliLog::SetClassDebugLevel("AliESDtrackCuts",1);
27 // ########################################################
28 // definition of dNdEta correction object
30 AlidNdEtaCorrection* dNdEtaMap = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
32 // ########################################################
35 sprintf(execDir,gSystem->pwd());
36 TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
37 TList* dirList = baseDir->GetListOfFiles();
38 Int_t nDirs = dirList->GetEntries();
39 // go back to the dir where this script is executed
42 // ########################################################
43 // definition of used pointers
50 // ########################################################
52 Int_t nRunCounter = 0;
53 for (Int_t r=0; r<nDirs; r++) {
55 TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
56 if (!presentDir || !presentDir->IsDirectory())
58 // check that the files are there
59 TString currentDataDir;
60 currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
61 cout << "Processing directory " << currentDataDir.Data() << endl;
62 if ((!gSystem->Which(currentDataDir,"galice.root")) ||
63 (!gSystem->Which(currentDataDir,"AliESDs.root")))
66 if (nRunCounter++ >= nRuns)
69 cout << "run #" << nRunCounter << endl;
71 // #########################################################
72 // setup galice and runloader
74 delete gAlice->GetRunLoader();
79 sprintf(str,"%s/galice.root",currentDataDir.Data());
80 AliRunLoader* runLoader = AliRunLoader::Open(str);
82 runLoader->LoadgAlice();
83 gAlice = runLoader->GetAliRun();
84 runLoader->LoadKinematics();
85 runLoader->LoadHeader();
87 // #########################################################
88 // open esd file and get the tree
90 sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
91 // close it first to avoid memory leak
93 if (esdFile->IsOpen())
96 esdFile = TFile::Open(str);
97 esdTree = (TTree*)esdFile->Get("esdTree");
100 esdBranch = esdTree->GetBranch("ESD");
101 esdBranch->SetAddress(&esd);
105 // ########################################################
107 AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
108 //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
110 // ########################################################
111 // getting number of events
112 Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
113 Int_t nESDEvents = esdBranch->GetEntries();
115 Int_t nEventsTriggered = 0;
116 Int_t nEventsAll = 0;
118 if (nEvents!=nESDEvents) {
119 cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
122 // ########################################################
123 // loop over number of events
124 cout << " looping over events..." << endl;
125 for(Int_t i=0; i<nEvents; i++) {
127 esdBranch->GetEntry(i);
128 runLoader->GetEvent(i);
130 // ########################################################
131 // get the EDS vertex
132 AliESDVertex* vtxESD = esd->GetVertex();
138 vtx_res[0] = vtxESD->GetXRes();
139 vtx_res[1] = vtxESD->GetYRes();
140 vtx_res[2] = vtxESD->GetZRes();
142 Bool_t vertexReconstructed = kTRUE;
144 // the vertex should be reconstructed
145 if (strcmp(vtxESD->GetName(),"default")==0)
146 vertexReconstructed = kFALSE;
148 // the resolution should be reasonable???
149 if (vtx_res[2]==0 || vtx_res[2]>0.01)
150 vertexReconstructed = kFALSE;
152 // ########################################################
153 // get the trigger info
155 Bool_t triggered = kFALSE;
158 // ITS_SPD_GFO_L0 : 32
160 // VZERO_OR_RIGHT : 2
162 ULong64_t triggerMask = esd->GetTriggerMask();
166 if (triggerMask&32 && ((triggerMask&1) || (triggerMask&2))) {
171 // ########################################################
173 AliGenPythiaEventHeader* genHeader =
174 (AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
177 genHeader->PrimaryVertex(vtxMC);
183 // ########################################################
184 // loop over mc particles
185 AliStack* particleStack = runLoader->Stack();
186 Int_t nPrim = particleStack->GetNprimary();
188 for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
190 TParticle* part = particleStack->Particle(i_mc);
191 if (!part || strcmp(part->GetName(),"XXX")==0)
194 TParticlePDG* pdgPart = part->GetPDG();
196 Bool_t prim = kFALSE;
197 // check if it's a primary particle - is there a better way ???
198 if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
199 if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
205 if (pdgPart->Charge()==0)
208 Float_t eta = part->Eta();
209 Float_t pt = part->Pt();
213 dNdEtaMap->FillParticleAllEvents(eta, pt);
216 dNdEtaMap->FillParticleWhenEventTriggered(eta, pt);
218 if (vertexReconstructed)
219 dNdEtaMap->FillParticle(vtx[2], eta, 1.);
222 }// end of mc particle
224 // ########################################################
225 // loop over esd tracks
226 Int_t nGoodTracks = 0;
228 Int_t nTracks = esd->GetNumberOfTracks();
229 for (Int_t t=0; t<nTracks; t++) {
230 AliESDtrack* esdTrack = esd->GetTrack(t);
232 // cut the esd track?
233 if (!esdTrackCuts->AcceptTrack(esdTrack))
239 esdTrack->GetPxPyPz(p);
240 Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
243 if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
244 eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
246 // using the eta of the mc particle
247 Int_t label = TMath::Abs(esdTrack->GetLabel());
249 cout << " cannot find corresponding mc part !!! " << label << endl;
252 TParticle* mcPart = particleStack->Particle(label);
254 Float_t pt = mcPart->Pt();
256 if (vertexReconstructed)
257 dNdEtaMap->FillParticleWhenMeasuredTrack(vtx[2], eta, pt);
259 } // end of track loop
261 dNdEtaMap->FillEvent(vtxMC[2], nGoodTracks);
263 if (vertexReconstructed)
264 dNdEtaMap->FillEventWithReconstructedVertex(vtxMC[2], nGoodTracks);
266 } // end of event loop
269 dNdEtaMap->Finish(nEventsAll, nEventsTriggered);
271 TFile* fout = new TFile("correction_map.root","RECREATE");
273 esdTrackCuts->SaveHistograms("esd_track_cuts");
274 dNdEtaMap->SaveHistograms();
279 dNdEtaMap->DrawHistograms();