next50 trigger mask in AliHLTGlobalEsdConverterComponent
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDhistTRDpid.C
CommitLineData
8470c96e 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
36Bool_t AliTRDhistTRDpid(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