]>
Commit | Line | Data |
---|---|---|
94d260cf | 1 | /* $Id$ */ |
2 | ||
3 | // | |
4 | // Script to calculate PtvsEta correction map at different z of vtx. | |
5 | // | |
6 | ||
7 | makeCorrectionPtEtaVSz(Char_t* dataDir, Int_t nRuns=20) { | |
8 | ||
9 | Char_t str[256]; | |
10 | static const Int_t NZbin=5; | |
11 | Float_t Zbin[]={-10., -3.,0., 3., 10.}; | |
12 | ||
b0635849 | 13 | gSystem->Load("../libPWG0base"); |
94d260cf | 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 | ||
31 | CorrectionMatrix2D *PtEtaMap[NZbin]; | |
32 | for(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 | ||
154 | Int_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 | } |