Update of calibration task and replace AliTRDrunCalib.C by AddTaskTRDCalib.C (Raphaelle)
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDhistTRDpid.C
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 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