64ba207a4c7c29ac60680cfbee76fdd95edcdf90
[u/mrichter/AliRoot.git] / HLT / QA / tasks / AliAnalysisTaskHLTCentralBarrel.cxx
1 // $Id$
2 //**************************************************************************
3 //* This file is property of and copyright by the ALICE HLT Project        *
4 //* ALICE Experiment at CERN, All rights reserved.                         *
5 //*                                                                        *
6 //* Primary Authors: Kalliopi Kanaki <Kalliopi.Kanaki@ift.uib.no>          *
7 //*                  for The ALICE HLT Project.                            *
8 //*                                                                        *
9 //* Permission to use, copy, modify and distribute this software and its   *
10 //* documentation strictly for non-commercial purposes is hereby granted   *
11 //* without fee, provided that the above copyright notice appears in all   *
12 //* copies and that both the copyright notice and this permission notice   *
13 //* appear in the supporting documentation. The authors make no claims     *
14 //* about the suitability of this software for any purpose. It is          *
15 //* provided "as is" without express or implied warranty.                  *
16 //**************************************************************************
17
18
19 /** @file  AliAnalysisTaskHLTCentralBarrel.cxx  
20     @author Per Ivar L√łnne, Hege Erdal, Kalliopi Kanaki
21     @date 
22     @brief An analysis task containing
23     loops over HLT and offline ESD trees for comparing
24     event and track properties
25     
26     After the task has been run, HLT/QA/tasks/macros/drawTHnSparse.C 
27     can be used to visualize the contents of the THnSparse objects stored
28     in the output file. Documentation about how to run it is included 
29     in the macro.
30 */
31
32 #include <iostream>
33
34 #include "AliAnalysisTaskHLTCentralBarrel.h"
35 #include "AliESDEvent.h"
36 #include "AliESDtrack.h"
37 #include "AliESDInputHandler.h"
38 #include "AliTracker.h"
39 #include "AliESDVZERO.h"
40 //#include "AliHLTGlobalTriggerDecision.h"
41 #include "AliCentrality.h"
42
43 #include "TAxis.h"
44 #include "TString.h"
45 #include "TText.h"
46 #include "TTimeStamp.h"
47 #include "TH1.h"
48
49 ClassImp(AliAnalysisTaskHLTCentralBarrel)
50 //_______________________________________________________________________________________________//
51
52 AliAnalysisTaskHLTCentralBarrel::AliAnalysisTaskHLTCentralBarrel()
53 :AliAnalysisTaskSE()
54   ,fUseHLTTrigger(kFALSE)
55   ,fCentrality()
56   ,fBeamType()
57   ,fOutputList(0)
58   ,fEventOFF(0)
59   ,fEventHLT(0)
60   ,fTrackOFF(0) 
61   ,fTrackHLT(0)
62   ,fOptions()
63   ,fTextBox(0)
64   ,fSwitch(kTRUE)
65 {
66   // Constructor
67   // Define input and output slots here
68   // Input slot #0 works with a TChain
69   // DefineInput(0, TChain::Class());
70   // Output slot #0 writes into a TH1 container
71   
72   //DefineOutput(1, TList::Class());
73 }
74
75 AliAnalysisTaskHLTCentralBarrel::AliAnalysisTaskHLTCentralBarrel(const char *name)
76 :AliAnalysisTaskSE(name)    
77   ,fUseHLTTrigger(kFALSE)   
78   ,fCentrality()   
79   ,fBeamType()
80   ,fOutputList(0)
81   ,fEventOFF(0)
82   ,fEventHLT(0)
83   ,fTrackOFF(0) 
84   ,fTrackHLT(0) 
85   ,fOptions()
86   ,fTextBox(0)
87   ,fSwitch(kTRUE)
88 {
89   // Constructor
90   // Define input and output slots here
91   // Input slot #0 works with a TChain
92   // DefineInput(0, TChain::Class());
93   // Output slot #0 writes into a TH1 container
94
95   DefineOutput(1, TList::Class());
96 }
97
98 AliAnalysisTaskHLTCentralBarrel::~AliAnalysisTaskHLTCentralBarrel(){
99   // destructor
100 }
101
102 void AliAnalysisTaskHLTCentralBarrel::UserCreateOutputObjects(){
103   // Create THnSparse objects
104   
105   OpenFile(1);  
106   fOutputList = new TList();
107   fOutputList->SetOwner();
108   fOutputList->SetName(GetName());
109   
110   if(fBeamType.Contains("Pb")){ // in case of a Pb+Pb run the V0 centrality is added to the THnSparse  
111      static const int sizeEvent = 10;
112      //                               0    1     2    3     4     5      6       7       8               9
113      //                               x    y     z   spdx spdy  spdz  #contr   mult  vertexStatus  V0centrality
114      int    binsEvent[sizeEvent] = { 100, 100,  60, 100,  100,   60,    200,    200,    2,            100 }; // binning
115      double minEvent [sizeEvent] = {  -1,  -1, -20,  -1,   -1,  -20,      0,      0,    0,              0 }; // min x
116      double maxEvent [sizeEvent] = {   1,   1,  20,   1,    1,   20,   2000,   2000,    2,            100 }; // max x    
117      fEventHLT = CreateEventTHnSparse("fEventHLT",sizeEvent,binsEvent,minEvent,maxEvent);
118      fEventOFF = CreateEventTHnSparse("fEventOFF",sizeEvent,binsEvent,minEvent,maxEvent);
119   }
120   else {
121      static const int sizeEvent = 9;
122      //                               0    1     2    3     4     5      6       7       8       
123      //                               x    y     z   spdx spdy  spdz  #contr   mult  vertexStatus
124      int    binsEvent[sizeEvent] = { 100, 100,  60, 100,  100,   60,    200,    200,    2 }; // binning
125      double minEvent [sizeEvent] = {  -1,  -1, -20,  -1,   -1,  -20,      0,      0,    0 }; // min x
126      double maxEvent [sizeEvent] = {   1,   1,  20,   1,    1,   20,   2000,   2000,    2 }; // max x   
127      fEventHLT = CreateEventTHnSparse("fEventHLT",sizeEvent,binsEvent,minEvent,maxEvent);
128      fEventOFF = CreateEventTHnSparse("fEventOFF",sizeEvent,binsEvent,minEvent,maxEvent);  
129   }
130   
131   if(fBeamType.Contains("Pb")){ // in case of a Pb+Pb run the V0 centrality is added to the THnSparse
132      static const int sizeTrack = 11;
133      //                                 0    1     2   3      4    5     6       7     8             9         10          
134      //                                pt  TPCcl  eta phi   DCAr  DCAz charge  ITScl vertex status  vertexZ  V0centrality
135      Int_t    binsTrack[sizeTrack] = { 400, 200, 200, 200,  100,  100,    3,   10,    2,            60,      100 }; // binning
136      Double_t minTrack [sizeTrack] = {   0,   0,  -2,  -1,  -10,  -10, -1.5,    0,    0,           -20,        0 }; // min x
137      Double_t maxTrack [sizeTrack] = {  10, 200,   2,   7,   10,  -10,  1.5,   10,    2,            20,      100 }; // max x      
138      fTrackHLT = CreateTrackTHnSparse("fTrackHLT",sizeTrack,binsTrack,minTrack,maxTrack);
139      fTrackOFF = CreateTrackTHnSparse("fTrackOFF",sizeTrack,binsTrack,minTrack,maxTrack);
140   }
141   else {    
142      static const int sizeTrack = 10;
143      //                                 0    1     2    3      4    5     6         7           8            9     
144      //                                pt  TPCcl  eta  phi   DCAr  DCAz  charge   ITScl   vertex status   vertexZ  
145      Int_t    binsTrack[sizeTrack] = {400, 200,  200,  200,  100,  100,      3,    10,       2,           60 }; // binning
146      Double_t minTrack [sizeTrack] = {  0,   0,   -2,   -1,  -10,  -10,   -1.5,     0,       0,          -20 }; // min x
147      Double_t maxTrack [sizeTrack] = { 10 ,200,    2,    7,   10,   10,    1.5,    10,       2,           20 }; // max x        
148      fTrackHLT = CreateTrackTHnSparse("fTrackHLT",sizeTrack,binsTrack,minTrack,maxTrack);
149      fTrackOFF = CreateTrackTHnSparse("fTrackOFF",sizeTrack,binsTrack,minTrack,maxTrack);     
150   }
151   
152   fTextBox = new TText();  
153   fOutputList->Add(fEventOFF);
154   fOutputList->Add(fEventHLT);
155   fOutputList->Add(fTrackOFF);
156   fOutputList->Add(fTrackHLT);  
157   fOutputList->Add(fTextBox);
158   
159   PostData(1, fOutputList);
160 }
161
162 void AliAnalysisTaskHLTCentralBarrel::UserExec(Option_t *){
163 // see header for documentation
164     
165   AliESDEvent *esdOFF = dynamic_cast<AliESDEvent*>(InputEvent());  
166   if(!esdOFF){
167      printf("ERROR: offline ESD is not available.\n");
168      return;
169   }
170    
171   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(fInputHandler);
172   AliESDEvent *esdHLT = NULL;   
173
174   if(esdH) esdHLT = esdH->GetHLTEvent();  
175   if(!esdHLT){
176      printf("ERROR: HLT ESD is not available.\n");
177      return;
178   }
179   
180   if(fSwitch==kTRUE){  
181      TTimeStamp *timestamp = new TTimeStamp(esdHLT->GetTimeStamp());
182      fTextBox->SetName("text");
183      TString s = Form("Run %d, Date %d", esdHLT->GetRunNumber(), timestamp->GetDate());
184      printf("You are analyzing run %d from date %d\n\n", esdHLT->GetRunNumber(), timestamp->GetDate());
185      fTextBox->SetTitle(s);
186      fSwitch=kFALSE;
187   }
188     
189   // if(fUseHLTTrigger && !((AliHLTGlobalTriggerDecision*)esdHLT->GetHLTTriggerDecision())->Result()) return;  
190   
191   //============================ OFFLINE =============================//
192
193   const AliESDVertex *vertOFF = esdOFF->GetPrimaryVertexTracks();  
194   Double_t bfield = esdOFF->GetMagneticField();
195
196   Int_t nr_tracksOFF = 0;  
197   if(esdOFF->GetEventSpecie()==16) return; // skip calibration events
198
199   if(fBeamType.Contains("Pb")){
200      fCentrality = esdOFF->GetCentrality(); 
201      // this information is only available from the offline ESD for 2010 PbPb data, the V0 info was not stored in the HLTesd (23.03.11, Kelly)
202      if(!fCentrality){
203         printf("Centrality pointer is empty\n");
204         return;
205      }
206   }
207                
208   for(Int_t i=0; i<esdOFF->GetNumberOfTracks(); i++){     
209       AliESDtrack *esdTrackOFF = esdOFF->GetTrack(i);
210       if (!esdTrackOFF) continue;
211       if(!(esdTrackOFF->GetStatus()&AliESDtrack::kTPCin)) continue;
212       nr_tracksOFF++;
213
214       //DCA calculations(from offline)
215       Double_t x[3];
216       Double_t b[3]; 
217       esdTrackOFF->GetXYZ(x);
218      
219       AliTracker::GetBxByBz(x,b);
220       Bool_t isOK = esdTrackOFF->RelateToVertexTPCBxByBz(vertOFF, b, kVeryBig);
221       if(!isOK) continue;
222         
223       const AliExternalTrackParam *track = esdTrackOFF->GetTPCInnerParam();
224       if(!track) continue;
225         
226       Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
227       if(vertOFF->GetX()==0 && vertOFF->GetY()==0 && vertOFF->GetZ()==0 ){
228          dca[0]=-99;
229          dca[1]=-99;
230       }
231       esdTrackOFF->GetImpactParametersTPC(dca,cov);
232
233       if(fBeamType.Contains("Pb")){       
234          Double_t trackOFF[] = {
235                                    TMath::Abs(esdTrackOFF->Pt()) 
236                                   ,esdTrackOFF->GetTPCNcls()      
237                                   ,esdTrackOFF->Eta()             
238                                   ,esdTrackOFF->Phi()             
239                                   ,dca[0]                         
240                                   ,dca[1]                         
241                                   ,esdTrackOFF->Charge() 
242                                   ,esdTrackOFF->GetNcls(0)
243                                   ,nr_tracksOFF
244                                   ,vertOFF->GetStatus()
245                                   ,vertOFF->GetZ()
246                                   ,fCentrality->GetCentralityPercentile("V0M")
247                                }; 
248         if(fOptions.Contains("track-off")) fTrackOFF->Fill(trackOFF);
249       } else {
250         Double_t trackOFF[] = {
251                                    TMath::Abs(esdTrackOFF->Pt()) 
252                                   ,esdTrackOFF->GetTPCNcls()      
253                                   ,esdTrackOFF->Eta()             
254                                   ,esdTrackOFF->Phi()             
255                                   ,dca[0]                         
256                                   ,dca[1]                         
257                                   ,esdTrackOFF->Charge() 
258                                   ,esdTrackOFF->GetNcls(0)
259                                   ,nr_tracksOFF
260                                   ,vertOFF->GetStatus()
261                                   ,vertOFF->GetZ()                                
262                               };
263         if(fOptions.Contains("track-off")) fTrackOFF->Fill(trackOFF);
264       }
265   } // end of track loop
266     
267   if(fBeamType.Contains("Pb")){
268      Double_t eventOFF[] = { vertOFF->GetX(), vertOFF->GetY(), vertOFF->GetZ(), 
269                              esdOFF->GetPrimaryVertexSPD()->GetX(), esdOFF->GetPrimaryVertexSPD()->GetY(), esdOFF->GetPrimaryVertexSPD()->GetZ(),  
270                              vertOFF->GetNContributors(), nr_tracksOFF, vertOFF->GetStatus(),fCentrality->GetCentralityPercentile("V0M")};
271      if(fOptions.Contains("event-off")) fEventOFF->Fill(eventOFF);  
272   }
273   else {
274      Double_t eventOFF[] = { vertOFF->GetX(), vertOFF->GetY(), vertOFF->GetZ(), 
275                              esdOFF->GetPrimaryVertexSPD()->GetX(), esdOFF->GetPrimaryVertexSPD()->GetY(), esdOFF->GetPrimaryVertexSPD()->GetZ(),      
276                              vertOFF->GetNContributors(), nr_tracksOFF, vertOFF->GetStatus()};                      
277      if(fOptions.Contains("event-off")) fEventOFF->Fill(eventOFF);
278   }
279   // Inspite of the different options to fill event or track properties, all the loops over tracks are being executed. 
280   // The options influence only whether the respective THnSparse is filled or not.
281   // Can definitely be improved to save processing time in unnecessary loops that won't fill anything at the end. 
282   // One thing to consider is the track multiplicity which is calculated as the number of tracks that have the kTPCin flag.
283   // This is calculated inside the track loop.  
284   
285   //======================================== HLT ==========================================//
286
287   Int_t nr_tracksHLT = 0; 
288   // the following line does not have an effect yet, since HLT does not fill the relevant information (14.04.11, Kelly)
289   // It does not harm to keep it in, in case this changes in the future 
290   if(esdHLT->GetEventSpecie()==16) return; // skip calibration events
291
292   const AliESDVertex *vertHLT = esdHLT->GetPrimaryVertexTracks();
293   for(Int_t i=0; i<esdHLT->GetNumberOfTracks(); i++){
294       
295       AliESDtrack *esdTrackHLT = esdHLT->GetTrack(i);
296       if(!esdTrackHLT) continue;
297       if(!(esdTrackHLT->GetStatus()&AliESDtrack::kTPCin)) continue;
298       nr_tracksHLT++; 
299       
300       Float_t dca[2];
301       if(vertHLT->GetX()==0 && vertHLT->GetY()==0 && vertHLT->GetZ() ==0 ){
302         dca[0]=-99;
303         dca[1]=-99;
304       }
305       else{
306         //Calculating DCA "old" fashion
307         esdTrackHLT->GetDZ(esdHLT->GetPrimaryVertex()->GetXv(), esdHLT->GetPrimaryVertex()->GetYv(), esdHLT->GetPrimaryVertex()->GetZv(), bfield, dca);
308         // plotting the DCA calculated by Sergey 
309         //esdTrackHLT->GetImpactParametersTPC(DCAr,DCAz);
310       }
311       
312       if(fBeamType.Contains("Pb")){
313          Double_t trackHLT[] = {
314                                  TMath::Abs(esdTrackHLT->Pt())
315                                 ,esdTrackHLT->GetTPCNcls()    
316                                 ,esdTrackHLT->Eta()           
317                                 ,esdTrackHLT->Phi()
318                                 ,dca[0]                       
319                                 ,dca[1]                       
320                                 ,esdTrackHLT->Charge()        
321                                 ,esdTrackHLT->GetNcls(0)
322                                 ,vertHLT->GetStatus()
323                                 ,vertHLT->GetZ()
324                                 ,fCentrality->GetCentralityPercentile("V0M")
325                                };
326         if(fOptions.Contains("track-hlt")) fTrackHLT->Fill(trackHLT);   
327       } else {
328          Double_t trackHLT[] = {
329                                  TMath::Abs(esdTrackHLT->Pt())
330                                 ,esdTrackHLT->GetTPCNcls()    
331                                 ,esdTrackHLT->Eta()           
332                                 ,esdTrackHLT->Phi()
333                                 ,dca[0]                       
334                                 ,dca[1]                       
335                                 ,esdTrackHLT->Charge()  
336                                 ,esdTrackHLT->GetNcls(0)
337                                 ,vertHLT->GetStatus()
338                                 ,vertHLT->GetZ()                       
339                                };
340         if(fOptions.Contains("track-hlt")) fTrackHLT->Fill(trackHLT);
341       }   
342   }  // end of track loop  
343   
344   if(fBeamType.Contains("Pb")){
345      Double_t eventHLT[] = { vertHLT->GetX(), vertHLT->GetY(), vertHLT->GetZ(), 
346                              esdHLT->GetPrimaryVertexSPD()->GetX(), esdHLT->GetPrimaryVertexSPD()->GetY(), esdHLT->GetPrimaryVertexSPD()->GetZ(),  
347                              vertHLT->GetNContributors(), nr_tracksHLT, vertHLT->GetStatus(),fCentrality->GetCentralityPercentile("V0M")};
348      if(fOptions.Contains("event-hlt")) fEventHLT->Fill(eventHLT);  
349   }
350   else{
351      Double_t eventHLT[] = { vertHLT->GetX(), vertHLT->GetY(), vertHLT->GetZ(),
352                              esdHLT->GetPrimaryVertexSPD()->GetX(), esdHLT->GetPrimaryVertexSPD()->GetY(), esdHLT->GetPrimaryVertexSPD()->GetZ(),  
353                              vertHLT->GetNContributors(), nr_tracksHLT, vertHLT->GetStatus()};
354      if(fOptions.Contains("event-hlt")) fEventHLT->Fill(eventHLT);
355   }
356              
357   // Post output data.
358   PostData(1, fOutputList);
359 }
360
361 void AliAnalysisTaskHLTCentralBarrel::Terminate(Option_t *){
362   // see header file of AliAnalysisTask for documentation 
363 }
364
365 THnSparseF* AliAnalysisTaskHLTCentralBarrel::CreateEventTHnSparse(const char* name, Int_t size, const Int_t* bins, Double_t* min, Double_t* max){
366 //see header for documentation                     
367   
368   THnSparseF *thn = new THnSparseF(name,"",size,bins,min,max);
369   thn->GetAxis(0)->SetTitle("primary vertex x");
370   thn->GetAxis(1)->SetTitle("primary vertex y");
371   thn->GetAxis(2)->SetTitle("primary vertex z"); 
372   thn->GetAxis(3)->SetTitle("SPD primary vertex x");
373   thn->GetAxis(4)->SetTitle("SPD primary vertex y");
374   thn->GetAxis(5)->SetTitle("SPD primary vertex z");
375   thn->GetAxis(6)->SetTitle("number of contributors");
376   thn->GetAxis(7)->SetTitle("track multiplicity");
377   thn->GetAxis(8)->SetTitle("vertex status"); 
378   if(fBeamType.Contains("Pb")) thn->GetAxis(9)->SetTitle("V0 centrality");
379   return thn;
380 }
381
382 THnSparseF* AliAnalysisTaskHLTCentralBarrel::CreateTrackTHnSparse(const char* name, Int_t size, const Int_t* bins, Double_t* min, Double_t* max){
383 //see header for documentation                     
384   
385   THnSparseF *thn = new THnSparseF(name,"",size,bins,min,max);
386   thn->GetAxis(0)->SetTitle("p_{T}");
387   thn->GetAxis(1)->SetTitle("TPC clusters/track");
388   thn->GetAxis(2)->SetTitle("#eta");
389   thn->GetAxis(3)->SetTitle("#phi");
390   thn->GetAxis(4)->SetTitle("DCAr");
391   thn->GetAxis(5)->SetTitle("DCAz");
392   thn->GetAxis(6)->SetTitle("polarity");
393   thn->GetAxis(7)->SetTitle("ITS clusters/track"); 
394   return thn;
395 }
396
397 // void AliAnalysisTaskHLTCentralBarrel::Fill(AliESDevent *esd, THnSparseF* thn){
398 // //see header for documentation                     
399 //  
400 //   int nTracks = 0;  
401 //   
402 //   for(int i=0; i<esdHLT->GetNumberOfTracks(); i++){
403 //       
404 //       AliESDtrack *esdTrack = esd->GetTrack(i);
405 //       if(!esdTrack) continue;
406 //       if(!(esdTrack->GetStatus()&AliESDtrack::kTPCin)) continue;
407 //       
408 //       //DCA calculations
409 //       Float_t DCAr=-99;
410 //       Float_t DCAz=-99;
411 //       Float_t dca[2];
412 //       if(vertHLT->GetX()==0 && vertHLT->GetY()==0 && vertHLT->GetZ() ==0 ){
413 //      dca[0]=-99;
414 //      dca[1]=-99;
415 //       }
416 //       else{
417 //      //Calculating DCA "old" fashion
418 //      esdTrackHLT->GetDZ(esdHLT->GetPrimaryVertex()->GetXv(), esdHLT->GetPrimaryVertex()->GetYv(), esdHLT->GetPrimaryVertex()->GetZv(), bfield, dca);
419 //      // plotting the DCA calculated by Sergey 
420 //      esdTrackHLT->GetImpactParametersTPC(DCAr,DCAz);
421 //       }
422 //      
423 //       Double_t trackHLT[] = {
424 //                              TMath::Abs(esdTrackHLT->Pt())
425 //                             ,esdTrackHLT->GetTPCNcls()    
426 //                             ,esdTrackHLT->Eta()           
427 //                             ,esdTrackHLT->Phi()           
428 //                             ,DCAr                         
429 //                             ,DCAz                         
430 //                             ,esdTrackHLT->Charge()        
431 //                             ,dca[0]                       
432 //                             ,dca[1]                       
433 //                             ,esdTrackHLT->GetStatus()
434 //                             ,esdTrackHLT->GetNcls(0)   
435 //                           };
436 //       fTrackHLT->Fill(trackHLT);      
437 //       nr_tracksHLT++;    
438 //   }  
439 //   Double_t eventHLT[] = { vertHLT->GetX(), vertHLT->GetY(), vertHLT->GetZ(), vertHLT->GetNContributors(), nr_tracksHLT, vertHLT->GetStatus()};
440 //   fEventHLT->Fill(eventHLT);  
441 // }