Corrected cast
[u/mrichter/AliRoot.git] / PWG0 / ptSpectra / makeCorrectionPtEtaVSz.C
CommitLineData
94d260cf 1/* $Id$ */
2
3//
4// Script to calculate PtvsEta correction map at different z of vtx.
5//
6
7makeCorrectionPtEtaVSz(Char_t* dataDir, Int_t nRuns=20) {
8
9 Char_t str[256];
10static const Int_t NZbin=5;
11Float_t Zbin[]={-10., -3.,0., 3., 10.};
12
13 gSystem->Load("../libPWG0base.so");
14 // ########################################################
15 // selection of esd tracks
16 AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();
17 esdTrackCuts->DefineHistograms(1);
18 esdTrackCuts->SetMinNClustersTPC(50);
19 esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
20 esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
21 esdTrackCuts->SetRequireTPCRefit(kTRUE);
22 esdTrackCuts->SetMinNsigmaToVertex(3);
23 esdTrackCuts->SetAcceptKingDaughters(kFALSE);
24 AliLog::SetClassDebugLevel("AliESDtrackCuts",1);
25 AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
26
27 // ########################################################
28 // definition of PtEta correction objects
29
30
31CorrectionMatrix2D *PtEtaMap[NZbin];
32for(Int_t i=0; i<NZbin-1; i++)
33 PtEtaMap[i] = new CorrectionMatrix2D(Form("PtvsEta%d",i),Form("%f<z_{vtx}<%f",Zbin[i],Zbin[i+1]),80,0.,10.,120,-2.,2.);
34
35 PtEtaMap[NZbin-1] = new CorrectionMatrix2D("xxx","xxx",80,0.,10.,120,-2.,2.);
36
37
38 // ########################################################
39 // get the data dir
40 Char_t execDir[256];
41 sprintf(execDir,gSystem->pwd());
42 TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
43 TList* dirList = baseDir->GetListOfFiles();
44 Int_t nDirs = dirList->GetEntries();
45 // go back to the dir where this script is executed
46 gSystem->cd(execDir);
47
48 // ########################################################
49 // definition of used pointers
50 TFile* esdFile;
51 TTree* esdTree;
52 TBranch* esdBranch;
53
54 AliESD* esd =0;
55 // ########################################################
56 // loop over runs
57 Int_t nRunCounter = 0;
58 for (Int_t r=0; r<nDirs; r++) {
59
60 TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
61 if (!presentDir || !presentDir->IsDirectory())
62 continue;
63 // check that the files are there
64 TString currentDataDir;
65 currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
66 cout << "Processing directory " << currentDataDir.Data() << endl;
67 if ((!gSystem->Which(currentDataDir,"galice.root")) ||
68 (!gSystem->Which(currentDataDir,"AliESDs.root")))
69 continue;
70
71 if (nRunCounter++ >= nRuns)
72 break;
73
74 cout << "run #" << nRunCounter << endl;
75
76 // #########################################################
77 // setup galice and runloader
78 if (gAlice) {
79 delete gAlice->GetRunLoader();
80 delete gAlice;
81 gAlice=0;
82 }
83
84 sprintf(str,"%s/galice.root",currentDataDir.Data());
85 AliRunLoader* runLoader = AliRunLoader::Open(str);
86
87 runLoader->LoadgAlice();
88 gAlice = runLoader->GetAliRun();
89 runLoader->LoadKinematics();
90 runLoader->LoadHeader();
91
92
93 // #########################################################
94 // open esd file and get the tree
95
96 sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
97 // close it first to avoid memory leak
98 if (esdFile)
99 if (esdFile->IsOpen())
100 esdFile->Close();
101
102 esdFile = TFile::Open(str);
103 esdTree = (TTree*)esdFile->Get("esdTree");
104 if (!esdTree)
105 continue;
106 esdBranch = esdTree->GetBranch("ESD");
107 esdBranch->SetAddress(&esd);
108 if (!esdBranch)
109 continue;
110
111 // ########################################################
112 // Magnetic field
113 AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
114 //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
115
116 // ########################################################
117 // getting number of events
118 Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
119 Int_t nESDEvents = esdBranch->GetEntries();
120
121 if (nEvents!=nESDEvents) {
122 cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
123 return;
124 }
125 // ########################################################
126 // loop over number of events
127 cout << " looping over events..." << endl;
128 for(Int_t i=0; i<nEvents; i++) {
129
130 esdBranch->GetEntry(i);
131 runLoader->GetEvent(i);
132
133 // ########################################################
134 // get the EDS vertex
135 AliESDVertex* vtxESD = esd->GetVertex();
136
137 Double_t vtx[3];
138 Double_t vtx_res[3];
139 vtxESD->GetXYZ(vtx);
140
141
142 vtx_res[0] = vtxESD->GetXRes();
143 vtx_res[1] = vtxESD->GetYRes();
144 vtx_res[2] = vtxESD->GetZRes();
145
146 // the vertex should be reconstructed
147 if (strcmp(vtxESD->GetName(),"default")==0)
148 continue;
149
150 // the resolution should be reasonable???
151 if (vtx_res[2]==0 || vtx_res[2]>0.1)
152 continue;
153
154Int_t IZbin=NZbin-1;
155 for(Int_t l=0; l<NZbin-1; l++)
156 if((vtx[3]<Zbin[l+1])&&((vtx[3]>Zbin[l])))IZbin=l;
157
158
159
160 // ########################################################
161 // loop over mc particles
162 AliStack* particleStack = runLoader->Stack();
163 Int_t nPrim = particleStack->GetNprimary();
164
165 for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
166
167 TParticle* part = particleStack->Particle(i_mc);
168 if (!part || strcmp(part->GetName(),"XXX")==0)
169 continue;
170
171 TParticlePDG* pdgPart = part->GetPDG();
172
173 Bool_t prim = kFALSE;
174 // check if it's a primary particle - is there a better way ???
175 if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
176 if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
177 prim = kTRUE;
178 }
179 if (!prim)
180 continue;
181
182 if (pdgPart->Charge()==0)
183 continue;
184
185 if (prim)
186 PtEtaMap[IZbin]->FillGene(part->Pt(), part->Eta());
187
188 }// end of mc particle
189
190 // ########################################################
191 // loop over esd tracks
192 Int_t nTracks = esd->GetNumberOfTracks();
193 for (Int_t t=0; t<nTracks; t++) {
194 AliESDtrack* esdTrack = esd->GetTrack(t);
195
196 // cut the esd track?
197 if (!esdTrackCuts->AcceptTrack(esdTrack))
198 continue;
199
200 AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);
201 if (tpcTrack->GetAlpha()==0) {
202 cout << " Warning esd track alpha = 0" << endl;
203 continue;
204 }
205
206 Float_t theta = tpcTrack->Theta();
207 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
208 Float_t pt= = tpcTrack->Pt();
209
210
211 PtEtaMap[IZbin]->FillMeas(eta, pt);
212
213 } // end of track loop
214 } // end of event loop
215 } // end of run loop
216
217
218 for(Int_t l=0; l<NZbin-1; l++)
219 {
220 PtEtaMap[l]->Divide();
221 PtEtaMap[l]->RemoveEdges();
222 }
223 TFile* fout = new TFile("PtEtaCorrectionMap.root","RECREATE");
224
225 esdTrackCuts->SaveHistograms("EsdTrackCuts");
226 for(Int_t l=0; l<NZbin-1; l++) PtEtaMap[l]->SaveHistograms();
227
228 fout->Write();
229 fout->Close();
230
231}