]>
Commit | Line | Data |
---|---|---|
defc9040 | 1 | // //////////////////////////////////////////////////////////////////////////// |
2 | // Macro to get TRD signal distribution from all six layers of TRD | |
3 | // This function takes the data from following directories | |
4 | // "rmom.6", "rmom.8", "rmom1", "rmom1.5", "rmom2", "rmom3", "rmom4", "rmom5",// "rmom6", "rmom8", "rmom10" | |
5 | // corresponding to momenta | |
6 | // 0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0 | |
7 | // Running this macro will generate the root file TRDdEdxHistogramsV1.root | |
8 | // containing histograms for all momenta. | |
9 | // Prashant Shukla (shukla@physi.uni-heidelberg.de) | |
10 | // //////////////////////////////////////////////////////////////////////////// | |
11 | ||
12 | #if !defined( __CINT__) || defined(__MAKECINT__) | |
13 | #include <TFile.h> | |
14 | #include <TH1.h> | |
15 | #include <TH2.h> | |
16 | #include <TCanvas.h> | |
17 | #include <TVector3.h> | |
18 | #include <TPDGCode.h> | |
19 | #include <TParticle.h> | |
20 | #include <TTree.h> | |
21 | ||
22 | #include "AliTRDtrack.h" | |
23 | #include "AliRunLoader.h" | |
24 | #include "AliLoader.h" | |
25 | #include "AliESD.h" | |
26 | #include "AliRun.h" | |
27 | #include "AliStack.h" | |
28 | #include "AliHeader.h" | |
29 | #include "AliGenEventHeader.h" | |
30 | #else | |
31 | #endif | |
32 | ||
33 | ||
34 | //_____________________________________________________________________________ | |
35 | ||
36 | Bool_t histTRDpid(Int_t oFile=1) | |
37 | { | |
38 | // | |
39 | // | |
40 | ||
41 | const Int_t kNMom=11; | |
42 | ||
43 | TH1F *h1dEdxEL[kNMom]; | |
44 | TH1F *h1dEdxPI[kNMom]; | |
45 | TH1F *h1dEdxMU[kNMom]; | |
46 | TH1F *h1dEdxKA[kNMom]; | |
47 | TH1F *h1dEdxPR[kNMom]; | |
48 | TH1F *h1MaxTimBinEL[kNMom]; | |
49 | TH1F *h1MaxTimBinPI[kNMom]; | |
50 | Char_t text[200]; | |
51 | ||
52 | // Histograms for dedx of e and pi track in TRD | |
53 | ||
54 | for (Int_t imom = 0; imom < kNMom; imom++) { | |
55 | sprintf(text,"h1dEdxEL%01d",imom+1); | |
56 | h1dEdxEL[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000); | |
57 | h1dEdxEL[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)"); | |
58 | h1dEdxEL[imom]->GetYaxis()->SetTitle("Entries"); | |
59 | h1dEdxEL[imom]->SetLineColor(kRed); | |
60 | ||
61 | sprintf(text,"h1dEdxPI%01d",imom+1); | |
62 | h1dEdxPI[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000); | |
63 | // TH1F *h1dEdxPI[imom] = new TH1F("","",250,0,2500); | |
64 | h1dEdxPI[imom]->GetXaxis()->SetTitle("dE/dx_{TRD} (arbitrary units)"); | |
65 | h1dEdxPI[imom]->GetYaxis()->SetTitle("Entries"); | |
66 | h1dEdxPI[imom]->SetLineColor(kBlue); | |
67 | ||
68 | sprintf(text,"h1dEdxMU%01d",imom+1); | |
69 | h1dEdxMU[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000); | |
70 | h1dEdxMU[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)"); | |
71 | h1dEdxMU[imom]->GetYaxis()->SetTitle("Entries"); | |
72 | h1dEdxMU[imom]->SetLineColor(kGreen); | |
73 | ||
74 | sprintf(text,"h1dEdxKA%01d",imom+1); | |
75 | h1dEdxKA[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000); | |
76 | h1dEdxKA[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)"); | |
77 | h1dEdxKA[imom]->GetYaxis()->SetTitle("Entries"); | |
78 | h1dEdxKA[imom]->SetLineColor(7); | |
79 | ||
80 | sprintf(text,"h1dEdxPR%01d",imom+1); | |
81 | h1dEdxPR[imom] = new TH1F(text,"dE/dx Dist.",200,0,4000); | |
82 | h1dEdxPR[imom]->GetXaxis()->SetTitle("dE/dx_{TRD}(arbitrary units)"); | |
83 | h1dEdxPR[imom]->GetYaxis()->SetTitle("Entries"); | |
84 | h1dEdxPR[imom]->SetLineColor(6); | |
85 | ||
86 | sprintf(text,"h1MaxTimBinEL%01d",imom+1); | |
87 | h1MaxTimBinEL[imom] = new TH1F(text,"Dist. Time Bin of max. Cluster for e(Red) and pi(Blue)",30,0,30); | |
88 | h1MaxTimBinEL[imom]->GetXaxis()->SetTitle("Time Bin of Maximum Cluster"); | |
89 | h1MaxTimBinEL[imom]->GetYaxis()->SetTitle("Entries"); | |
90 | h1MaxTimBinEL[imom]->SetLineColor(2); | |
91 | ||
92 | sprintf(text,"h1MaxTimBinPI%01d",imom+1); | |
93 | h1MaxTimBinPI[imom] = new TH1F(text,"Time Bin of max. Cluster Pion",30,0,30); | |
94 | h1MaxTimBinPI[imom]->SetLineColor(4); | |
95 | } | |
96 | ||
97 | // PID | |
98 | Int_t partCode[5] = | |
99 | {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; | |
100 | ||
101 | // Files | |
102 | ||
103 | Char_t *FileDir[11] = {"rmom.6", "rmom.8", "rmom1", "rmom1.5", "rmom2", "rmom3", "rmom4", "rmom5", "rmom6", "rmom8", "rmom10"}; | |
104 | ||
105 | Float_t mome[11]={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0}; | |
106 | ||
107 | Char_t gAliceFileName[200]; | |
108 | Char_t esdFileName[200]; | |
109 | for (Int_t imom = 0; imom < kNMom; imom++) { | |
110 | ||
111 | sprintf(gAliceFileName,"~/work/data/%s/galice.root", FileDir[imom]); | |
112 | sprintf(esdFileName,"~/work/data/%s/AliESDs.root", FileDir[imom]); | |
113 | ||
114 | // open run loader and load gAlice, kinematics and header | |
115 | AliRunLoader* runLoader = AliRunLoader::Open(gAliceFileName); | |
116 | if (!runLoader) { | |
117 | Error("CheckESD", "getting run loader from file %s failed",gAliceFileName); | |
118 | return kFALSE; | |
119 | } | |
120 | runLoader->LoadgAlice(); | |
121 | gAlice = runLoader->GetAliRun(); | |
122 | if (!gAlice) { | |
123 | Error("CheckESD", "no galice object found"); | |
124 | return kFALSE; | |
125 | } | |
126 | runLoader->LoadKinematics(); | |
127 | runLoader->LoadHeader(); | |
128 | ||
129 | // AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField()); | |
130 | // open the ESD file | |
131 | TFile* esdFile = TFile::Open(esdFileName); | |
132 | if (!esdFile || !esdFile->IsOpen()) { | |
133 | Error("CheckESD", "opening ESD file %s failed",esdFileName); | |
134 | return kFALSE; | |
135 | } | |
136 | AliESD* esd = new AliESD; | |
137 | TTree* tree = (TTree*) esdFile->Get("esdTree"); | |
138 | if (!tree) { | |
139 | Error("CheckESD", "no ESD tree found"); | |
140 | return kFALSE; | |
141 | } | |
142 | tree->SetBranchAddress("ESD", &esd); | |
143 | ||
144 | // loop over events | |
145 | for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { | |
146 | runLoader->GetEvent(iEvent); | |
147 | AliStack* stack = gAlice->Stack(); | |
148 | Int_t nParticles = stack->GetNtrack(); | |
149 | TArrayF vertex(3); | |
150 | runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex); | |
151 | Int_t nGene=0, nGenpi=0; | |
152 | for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) { | |
153 | TParticle* particle = stack->Particle(iParticle); | |
154 | if (!particle) continue; | |
155 | if (particle->Pt() < 0.001) continue; | |
156 | // if (TMath::Abs(particle->Eta()) > 0.3) continue; | |
157 | TVector3 dVertex(particle->Vx() - vertex[0], | |
158 | particle->Vy() - vertex[1], | |
159 | particle->Vz() - vertex[2]); | |
160 | if (dVertex.Mag() > 0.0001) continue; | |
161 | if (TMath::Abs(particle->GetPdgCode()) == partCode[0]) nGene++; | |
162 | if (TMath::Abs(particle->GetPdgCode()) == partCode[2]) nGenpi++; | |
163 | } // end loop over particles | |
164 | printf("Gen Prim. el= %d, Gen Prim. pi= %d in Event No.= %d in File = %s \n", nGene, nGenpi, iEvent, FileDir[imom]); | |
165 | ||
166 | // get the event summary data | |
167 | tree->GetEvent(iEvent); | |
168 | if (!esd) { | |
169 | Error("CheckESD", "no ESD object found for event %d", iEvent); | |
170 | return kFALSE; | |
171 | } | |
172 | ||
173 | /////////////////////////////////////////// | |
174 | // AliTRDpidESD::MakePID(esd); | |
175 | // AliESDpid::MakePID(esd); | |
176 | /////////////////////////////////////////// | |
177 | ||
178 | // Initializing number of electrons and pions | |
179 | Int_t nne=0, nnpi=0; | |
180 | // loop over tracks | |
181 | Int_t nTracks = esd->GetNumberOfTracks(); | |
182 | printf("Found %d tracks. \n",nTracks); | |
183 | ||
184 | for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { | |
185 | AliESDtrack* track = esd->GetTrack(iTrack); | |
186 | // select tracks of selected particles | |
187 | Int_t label = TMath::Abs(track->GetLabel()); | |
188 | if (label > stack->GetNtrack()) continue; // background | |
189 | TParticle* particle = stack->Particle(label); | |
190 | // if (TMath::Abs(particle->Eta()) > 0.3) continue; | |
191 | TVector3 dVertex(particle->Vx() - vertex[0], | |
192 | particle->Vy() - vertex[1], | |
193 | particle->Vz() - vertex[2]); | |
194 | if (dVertex.Mag() > 0.0001) continue; | |
195 | if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue; | |
196 | if (track->GetConstrainedChi2() > 1e9) continue; | |
197 | Double_t p[3]; | |
198 | track->GetConstrainedPxPyPz(p); | |
199 | TVector3 pTrack(p); | |
200 | // Double_t mom = track->GetP(); | |
201 | ||
202 | // PID///////////////////////////////////////////////////// | |
203 | if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue; | |
204 | Int_t iGen = 5; | |
205 | for (Int_t i = 0; i < 5; i++) { | |
206 | if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i; | |
207 | } | |
208 | ||
209 | // Filling muon kan and proton////////////////////////////// | |
210 | Int_t Itot=0; | |
211 | if(track->GetTRDsignal()) Itot=1; | |
212 | if(!Itot) continue; | |
213 | ||
214 | // dE/dx and TEST TRD | |
215 | Double_t TRDsignal = track->GetTRDsignal(); | |
216 | // if (iGen == 0) h1dEdxEL[imom]->Fill(TRDsignal); | |
217 | // if (iGen == 2) h1dEdxPI[imom]->Fill(TRDsignal); | |
218 | // if (iGen == 1) h1dEdxMU[imom]->Fill(TRDsignal); | |
219 | // if (iGen == 3) h1dEdxKA[imom]->Fill(TRDsignal); | |
220 | // if (iGen == 4) h1dEdxPR[imom]->Fill(TRDsignal); | |
221 | ||
222 | //Filling elcron and pions//////////////////////////////// | |
223 | const Int_t kNPlane = 6; | |
224 | // Track Selection | |
225 | Double_t Ise=1.; | |
226 | for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) { | |
227 | Ise*=track->GetTRDsignals(iPlane); | |
228 | } | |
229 | if(Ise==0) continue; | |
230 | ||
231 | for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) { | |
232 | Double_t TRDsignal = track->GetTRDsignals(iPlane); | |
233 | if (iGen == 0) h1dEdxEL[imom]->Fill(TRDsignal); | |
234 | if (iGen == 2) h1dEdxPI[imom]->Fill(TRDsignal); | |
235 | if (iGen == 1) h1dEdxMU[imom]->Fill(TRDsignal); | |
236 | if (iGen == 3) h1dEdxKA[imom]->Fill(TRDsignal); | |
237 | if (iGen == 4) h1dEdxPR[imom]->Fill(TRDsignal); | |
238 | ||
239 | if (iGen == 0) h1MaxTimBinEL[imom]->Fill(track->GetTRDTimBin(iPlane)); | |
240 | if (iGen == 2) h1MaxTimBinPI[imom]->Fill(track->GetTRDTimBin(iPlane)); | |
241 | } | |
242 | ||
243 | // Int_t TRDncls = track->GetTRDncls(); | |
244 | if (iGen == 0) nne++; | |
245 | if (iGen == 2) nnpi++; | |
246 | } // end of loop over tracks | |
247 | printf("Electron Tracks = %d, Pion Tracks = %d \n", nne, nnpi); | |
248 | }// end of loop over events | |
249 | delete esd; | |
250 | esdFile->Close(); | |
251 | delete esdFile; | |
252 | ||
253 | runLoader->UnloadHeader(); | |
254 | runLoader->UnloadKinematics(); | |
255 | delete runLoader; | |
256 | } // Loop over files | |
257 | ||
258 | ||
259 | Char_t fileprint[200]; | |
260 | ||
261 | // draw the histograms if not in batch mode | |
262 | ||
263 | new TCanvas; | |
264 | gPad->SetBorderMode(0); | |
265 | gPad->SetFillColor(0); | |
266 | // gStyle->SetStatColor(0); | |
267 | h1dEdxPI[0]->DrawCopy(); | |
268 | h1dEdxEL[0]->DrawCopy("SAME"); | |
269 | // h1dEdxMU[0]->DrawCopy("SAME"); | |
270 | // h1dEdxKA[0]->DrawCopy("SAME"); | |
271 | // h1dEdxPR[0]->DrawCopy("SAME"); | |
272 | ||
273 | new TCanvas; | |
274 | gPad->SetBorderMode(0); | |
275 | gPad->SetFillColor(0); | |
276 | h1MaxTimBinEL[0]->DrawCopy(); | |
277 | h1MaxTimBinPI[0]->DrawCopy("SAME"); | |
278 | ||
279 | ||
280 | // write the output histograms to a file | |
281 | ||
282 | Char_t fileresponse[200]; | |
283 | sprintf(fileresponse,"TRDdEdxHistogramsV%d.root",oFile); | |
284 | ||
285 | TFile* outputFile = TFile::Open(fileresponse, "recreate"); | |
286 | if (!outputFile || !outputFile->IsOpen()) { | |
287 | Error("CheckESD", "opening output file check.root failed"); | |
288 | return kFALSE; | |
289 | } | |
290 | ||
291 | for (Int_t imom = 0; imom < kNMom; imom++) { | |
292 | h1dEdxPI[imom]->Write(); | |
293 | h1dEdxEL[imom]->Write(); | |
294 | h1dEdxMU[imom]->Write(); | |
295 | h1dEdxKA[imom]->Write(); | |
296 | h1dEdxPR[imom]->Write(); | |
297 | h1MaxTimBinEL[imom]->Write(); | |
298 | h1MaxTimBinPI[imom]->Write(); | |
299 | } | |
300 | ||
301 | delete [] h1dEdxPI; | |
302 | delete [] h1dEdxEL; | |
303 | delete [] h1dEdxMU; | |
304 | delete [] h1dEdxKA; | |
305 | delete [] h1dEdxPR; | |
306 | delete [] h1MaxTimBinEL; | |
307 | delete [] h1MaxTimBinPI; | |
308 | delete outputFile; | |
309 | ||
310 | // result of check | |
311 | Info("GetPIDsignals", "GetPIDsignals was successfull"); | |
312 | return kTRUE; | |
313 | ||
314 | } // end of main program | |
315 |