TRD module
[u/mrichter/AliRoot.git] / TRD / TRDqaAnalysis / AliTRDqaJPsi.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as iLs" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //
19 // This class is a part of a package of high level QA monitoring for TRD.
20 //
21 // S. Radomski
22 // radomski@physi.uni-heidelberg.de
23 // March 2008
24 //
25
26 #include "AliTRDqaJPsi.h"
27 #include "AliTRDqaAT.h"
28
29 #include "TFile.h"
30 #include "TH1D.h"
31 #include "TH2D.h"
32 #include "TChain.h"
33
34 #include "AliESDEvent.h" 
35 #include "AliESDtrack.h"
36 #include "AliKFParticle.h"
37
38 #include "TLorentzVector.h"
39
40 //______________________________________________________________________________
41
42 AliTRDqaJPsi::AliTRDqaJPsi() 
43   : AliAnalysisTask("",""),  
44     fChain(0),
45     fESD(0),
46     fOutputContainer(0),
47     fnKFtracks(0)
48 {
49   //
50   // default dummy constructor
51   //
52
53   for (Int_t k = 0; k < fgknSteps; k++) {
54     fStatus[k]      = 0x0;
55     fAngleSM[k]     = 0x0;
56     fInvMass[k]     = 0x0;
57     fInvMassVec[k]  = 0x0;
58     fInvMassDiff[k] = 0x0;
59     fPtAngle[k]     = 0x0;
60   }
61
62   for (Int_t l = 0; l < 2*fgknSteps; l++) {
63     fnTracks[l]     = 0x0;
64     fPt[l]          = 0x0;
65     fPID[l]         = 0x0;
66   }
67
68   for (Int_t i = 0; i < 1000; i++) {
69     fSM[i]          = 0;
70     fTracks[i]      = 0x0;
71     fVec[i]         = 0x0;
72     for (Int_t j = 0; j < fgknSteps; j++) {
73       fInSample[i][j] = 0;
74     }
75   }
76
77 }
78
79 //______________________________________________________________________________
80 AliTRDqaJPsi:: AliTRDqaJPsi(const AliTRDqaJPsi & /*trd*/)
81   : AliAnalysisTask("",""),  
82     fChain(0),
83     fESD(0),
84     fOutputContainer(0),
85     fnKFtracks(0)
86 {
87   //
88   // Copy constructor
89   //
90
91   for (Int_t k = 0; k < fgknSteps; k++) {
92     fStatus[k]      = 0x0;
93     fAngleSM[k]     = 0x0;
94     fInvMass[k]     = 0x0;
95     fInvMassVec[k]  = 0x0;
96     fInvMassDiff[k] = 0x0;
97     fPtAngle[k]     = 0x0;
98   }
99
100   for (Int_t l = 0; l < 2*fgknSteps; l++) {
101     fnTracks[l]     = 0x0;
102     fPt[l]          = 0x0;
103     fPID[l]         = 0x0;
104   }
105
106   for (Int_t i = 0; i < 1000; i++) {
107     fSM[i]          = 0;
108     fTracks[i]      = 0x0;
109     fVec[i]         = 0x0;
110     for (Int_t j = 0; j < fgknSteps; j++) {
111       fInSample[i][j] = 0;
112     }
113   }
114
115 }
116
117 //______________________________________________________________________________
118 AliTRDqaJPsi::AliTRDqaJPsi(const char *name) 
119   : AliAnalysisTask(name,""),  
120     fChain(0),
121     fESD(0),
122     fOutputContainer(0),
123     fnKFtracks(0)
124     
125 {
126   // Constructor.
127   // Input slot #0 works with an Ntuple
128   DefineInput(0, TChain::Class());
129   // Output slot #0 writes into a TH1 container
130   DefineOutput(0,  TObjArray::Class()) ; 
131
132   for (Int_t k = 0; k < fgknSteps; k++) {
133     fStatus[k]      = 0x0;
134     fAngleSM[k]     = 0x0;
135     fInvMass[k]     = 0x0;
136     fInvMassVec[k]  = 0x0;
137     fInvMassDiff[k] = 0x0;
138     fPtAngle[k]     = 0x0;
139   }
140
141   for (Int_t l = 0; l < 2*fgknSteps; l++) {
142     fnTracks[l]     = 0x0;
143     fPt[l]          = 0x0;
144     fPID[l]         = 0x0;
145   }
146
147   for (Int_t i = 0; i < 1000; i++) {
148     fSM[i]          = 0;
149     fTracks[i]      = 0x0;
150     fVec[i]         = 0x0;
151     for (Int_t j = 0; j < fgknSteps; j++) {
152       fInSample[i][j] = 0;
153     }
154   }
155
156 }
157
158 //______________________________________________________________________________
159 void AliTRDqaJPsi::ConnectInputData(const Option_t *)
160 {
161   // Initialisation of branch container and histograms 
162
163   //AliInfo(Form("*** Initialization of %s", GetName())) ; 
164
165   fChain = (TChain*)GetInputData(0);
166   fESD = new AliESDEvent();
167   fESD->ReadFromTree(fChain);
168 }
169
170 //________________________________________________________________________
171 void AliTRDqaJPsi::CreateOutputObjects()
172 {
173   //
174   // Create the objects that contain the analysis output
175   //
176   
177   Int_t c = 0;
178   fOutputContainer = new TObjArray(100);
179   
180   const char *charge[2] = {"Neg", "Pos"};
181   
182   // build histograms
183   
184   for(Int_t i=0; i<fgknSteps; i++) {
185
186     fStatus[i] = new TH1D(Form("status_%d", i), "status", 32, -0.5, 31.5);
187     fOutputContainer->AddAt(fStatus[i], c++);
188     
189     fInvMass[i] = new TH1D(Form("mass_%d", i), ";m_{inv} (GeV);", 100, 0, 5);
190     fOutputContainer->AddAt(fInvMass[i], c++);
191     
192     fInvMassVec[i] = new TH1D(Form("massVec_%d", i), ";m_{inv} (GeV);", 100, 0, 5);
193     fOutputContainer->AddAt(fInvMassVec[i], c++);
194     
195     fInvMassDiff[i] = new TH1D(Form("massDiff_%d", i), ";m_{inv} (GeV);", 100, -1, 1);
196     fOutputContainer->AddAt(fInvMassDiff[i], c++);
197     
198     fAngleSM[i] = new TH1D(Form("angleSM_%d", i), ";#delta SM", 19, -0.5, 18.5);        
199     fOutputContainer->AddAt(fAngleSM[i], c++);
200     
201     fPtAngle[i] = new TH2D(Form("ptAngle_%d", i), ";p_{T} (GeV/c);#delta SM",
202                            20, 0, 5, 10, -0.5, 9.5);
203     fOutputContainer->AddAt(fPtAngle[i], c++);
204     
205
206     for(Int_t j=0; j<2; j++) {
207       fnTracks[j*fgknSteps+i] = 
208         new TH1D(Form("nTracks%s_%d", charge[j],i), Form("%s;number of tracks",charge[j]), 100, -0.5, 99.5);
209       fPt[j*fgknSteps+i] = new TH1D(Form("pt%s_%d", charge[j], i), Form("%s;p_{T} (GeV/c)", charge[j]), 100, 0, 5);
210       fPID[j*fgknSteps+i] = new TH1D(Form("pid%s_%d", charge[j], i), ";electron LQ", 100, 0, 1);
211
212       fOutputContainer->AddAt(fnTracks[j*fgknSteps+i], c++);
213       fOutputContainer->AddAt(fPt[j*fgknSteps+i], c++); 
214       fOutputContainer->AddAt(fPID[j*fgknSteps+i], c++); 
215     }
216   }
217
218   //TH2D *fnGoodTracks;  
219
220   printf("n hist = %d\n", c);
221 }
222 //______________________________________________________________________________
223 void AliTRDqaJPsi::Exec(Option_t *) 
224 {
225   /*
226     Selection steps:
227     - Parameters In and Out
228     - TRDrefit bit
229     - TRDpid and quality
230     - pt > 0.8 GeV/c
231     - PID > 0.9
232    */
233
234
235   // Process one event
236   Long64_t entry = fChain->GetReadEntry() ;
237   if (!(entry%100)) Info("Exec", "Entry = %lld", entry);
238
239   // Processing of one event 
240    
241   if (!fESD) {
242     //AliError("fESD is not connected to the input!") ; 
243     return ; 
244   }
245   
246   
247   Int_t nTracks = fESD->GetNumberOfTracks();
248   Int_t cTracks[2*fgknSteps] = {0,0,0,0,0,0,0,0,0,0};
249   fnKFtracks = 0;
250
251   // track loop
252   for(Int_t i=0; i<nTracks; i++) {
253     
254     //
255     // track selection 
256     //
257     // param in and Out
258     // TRDrefit and TRDPid bit
259     //
260  
261     AliESDtrack *track = fESD->GetTrack(i);
262     const AliExternalTrackParam *paramOut = track->GetOuterParam();
263     const AliExternalTrackParam *paramIn = track->GetInnerParam();
264
265     // long track ..
266     if (!paramIn) continue;
267     if (!paramOut) continue;
268
269     Int_t step    = 0;
270     Int_t charge  = (track->Charge() > 0) ? 1 : 0;
271     UInt_t status = track->GetStatus();
272     Double_t pt   = track->Pt();
273     Double_t pid  = track->GetTRDpid(AliPID::kElectron);
274     
275     Double_t esdPid[5];
276     track->GetESDpid(esdPid); 
277
278     // create a kalman particle
279     Int_t pdg = (charge == 0)? -11 : 11;
280     for(Int_t k=0; k<fgknSteps; k++) fInSample[fnKFtracks][k] = 0;  
281  
282     fVec[fnKFtracks] = CreateVector(track);
283     fTracks[fnKFtracks] = new AliKFParticle(*track, pdg);
284     fSM[fnKFtracks] = AliTRDqaAT::GetSector(paramOut->GetAlpha());
285     fnKFtracks++;
286
287     //AliTRDqaAT::PrintPID(track);
288
289     // apply the cuts
290
291     cTracks[fgknSteps *charge + step]++;
292     FillHist(track, step++);
293
294     if (!(status & AliESDtrack::kTRDrefit)) continue;
295
296     cTracks[fgknSteps *charge + step]++;    
297     FillHist(track, step++);    
298
299     if (!(status & AliESDtrack::kTRDpid)) continue;
300     if (track->GetTRDntracklets() < 6) continue;
301
302     cTracks[fgknSteps *charge + step]++;
303     FillHist(track, step++);  
304
305     if (pt < 0.8) continue;
306
307     cTracks[fgknSteps *charge + step]++;
308     FillHist(track, step++);      
309
310     if (pid < 0.3) continue; // 
311     //if (esdPid[AliPID::kElectron] < 0.5) continue;
312     
313     cTracks[fgknSteps *charge + step]++;
314     FillHist(track, step);  
315
316     for(Int_t k=0; k<2*fgknSteps; k++) fnTracks[k]->Fill(cTracks[k]);
317   }
318
319   // calculate invariant mass
320
321   for(Int_t k=0; k<fgknSteps; k++) {
322     for(Int_t i=0; i<fnKFtracks; i++) {
323       if (!fInSample[i][k]) continue;
324       for(Int_t j=i+1; j<fnKFtracks; j++) {
325         if (!fInSample[j][k]) continue;
326         AliKFParticle jpsi(*(fTracks[i]), *(fTracks[j]));
327         TLorentzVector jpsiVec = (*(fVec[i])) + (*fVec[j]);
328         fInvMass[k]->Fill(jpsi.GetMass());
329         fInvMassVec[k]->Fill(jpsiVec.M());
330         fInvMassDiff[k]->Fill(jpsiVec.M() - jpsi.GetMass());
331         
332         if (jpsi.GetMass() > 2.5 && jpsi.GetMass() < 3.5) {
333           Int_t dSM = TMath::Abs(fSM[i] - fSM[j]);
334           if (dSM > 9) dSM = (18-dSM);
335           fAngleSM[k]->Fill(dSM);
336           fPtAngle[k]->Fill(TMath::Hypot(jpsi.GetPx(), jpsi.GetPy()), dSM);
337         }
338         
339       }
340     }
341   }
342
343   PostData(0, fOutputContainer);
344 }
345
346 //______________________________________________________________________________
347 void AliTRDqaJPsi::Terminate(Option_t *)
348 {
349   // save histograms
350   fOutputContainer = (TObjArray*)GetOutputData(0);
351   
352   TFile *file = new TFile("outJPsi.root", "RECREATE");
353   fOutputContainer->Write();
354   
355   file->Flush();
356   file->Close();
357   delete file;  
358
359   //for(Int_t i=0; i<fOutputContainer->GetEntries(); i++) {
360   //  TObject *obj = fOu
361   // }
362 }
363
364 //______________________________________________________________________________
365
366 void AliTRDqaJPsi::FillHist(AliESDtrack *track, Int_t step) {
367   //
368   // Fill the histograms
369   //
370
371   Int_t charge  = (track->Charge() > 0) ? 1 : 0;
372   UInt_t status = track->GetStatus();
373   Double_t pt   = track->Pt();
374   Double_t pid  = track->GetTRDpid(AliPID::kElectron);
375   
376   Double_t esdPid[5];
377   track->GetESDpid(esdPid);
378
379   Int_t id = charge * fgknSteps + step;
380   AliTRDqaAT::FillStatus(fStatus[step], status);
381   fPt[id]->Fill(pt);
382   fPID[id]->Fill(pid);
383   //fPID[id]->Fill(esdPid[AliPID::kElectron]);
384
385   fInSample[fnKFtracks-1][step] = 1;
386 }
387
388 //______________________________________________________________________________
389
390 TLorentzVector *AliTRDqaJPsi::CreateVector(AliESDtrack *track) {
391
392   TLorentzVector *vec = new TLorentzVector();
393   vec->SetXYZM(track->Px(), track->Py(), track->Pz(), 0);
394   return vec;
395 }
396
397 //______________________________________________________________________________