- fix names and units of axes
[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
27 #include <iostream>
28
29 #include "AliAnalysisTaskHLTCentralBarrel.h"
30 #include "AliESDEvent.h"
31 #include "AliESDtrack.h"
32 #include "AliESDInputHandler.h"
33 #include "AliTracker.h"
34 #include "AliESDVZERO.h"
35 #include "AliHLTGlobalTriggerDecision.h"
36
37 #include "TAxis.h"
38 #include "TString.h"
39 #include "TText.h"
40 #include "TTimeStamp.h"
41 #include "TH1.h"
42
43 ClassImp(AliAnalysisTaskHLTCentralBarrel)
44 //_______________________________________________________________________________________________//
45
46 AliAnalysisTaskHLTCentralBarrel::AliAnalysisTaskHLTCentralBarrel()
47 :AliAnalysisTaskSE()
48   ,fUseHLTTrigger(kFALSE)   
49   ,fOutputList(0)
50   ,fEventOFF(0)
51   ,fEventHLT(0)
52   ,fTrackOFF(0) 
53   ,fTrackHLT(0)
54   ,fTextBox(0)
55 {
56   // Constructor
57   // Define input and output slots here
58   // Input slot #0 works with a TChain
59   // DefineInput(0, TChain::Class());
60   // Output slot #0 writes into a TH1 container
61   
62   //DefineOutput(1, TList::Class());
63 }
64
65 AliAnalysisTaskHLTCentralBarrel::AliAnalysisTaskHLTCentralBarrel(const char *name)
66 :AliAnalysisTaskSE(name)    
67   ,fUseHLTTrigger(kFALSE)   
68   ,fOutputList(0)
69   ,fEventOFF(0)
70   ,fEventHLT(0)
71   ,fTrackOFF(0) 
72   ,fTrackHLT(0) 
73   ,fTextBox(0)
74 {
75   // Constructor
76   // Define input and output slots here
77   // Input slot #0 works with a TChain
78   //DefineInput(0, TChain::Class());
79   // Output slot #0 writes into a TH1 container
80   
81   DefineOutput(1, TList::Class());
82 }
83
84 AliAnalysisTaskHLTCentralBarrel::~AliAnalysisTaskHLTCentralBarrel(){
85   // destructor
86 }
87
88 void AliAnalysisTaskHLTCentralBarrel::UserCreateOutputObjects(){
89   // Create THnSpare objects
90   
91   OpenFile(1);  
92   fOutputList = new TList();
93   fOutputList->SetName(GetName());
94   
95   static const int sizeEvent = 6;
96  
97   int    binsEvent[sizeEvent] = { 200, 200, 250,  2000,  2000, 2 };
98   double minEvent [sizeEvent] = {  -2,  -2, -30,     0,     0, 0 };
99   double maxEvent [sizeEvent] = {   2,   2,  30, 20000, 20000, 2 };
100   
101   fEventHLT = CreateEventTHnSparse("fEventHLT",sizeEvent,binsEvent,minEvent,maxEvent);
102   fEventOFF = CreateEventTHnSparse("fEventOFF",sizeEvent,binsEvent,minEvent,maxEvent);
103   
104   static const int sizeTrack = 13;
105   
106   Int_t    binsTrack[sizeTrack] = {1500, 200, 200, 200, 200,  400,  400,    3,  400,  400, 6,  2000, 2 };
107   Double_t minTrack [sizeTrack] = {   0,   0,  -1,  -3,  -1, -100, -100, -1.5, -100, -100, 0,     0, 0 };
108   Double_t maxTrack [sizeTrack] = { 150, 200,   4,   3,   7,  100,  100,  1.5,  100,  100, 6, 20000, 2 };
109    
110   fTrackHLT = CreateTrackTHnSparse("fTrackHLT",sizeTrack,binsTrack,minTrack,maxTrack);
111   fTrackOFF = CreateTrackTHnSparse("fTrackOFF",sizeTrack,binsTrack,minTrack,maxTrack);
112      
113   fTextBox = new TText();  
114   fOutputList->Add(fEventOFF);
115   fOutputList->Add(fEventHLT);
116   fOutputList->Add(fTrackOFF);
117   fOutputList->Add(fTrackHLT);
118   fOutputList->Add(fTextBox);
119 }
120
121 void AliAnalysisTaskHLTCentralBarrel::NotifyRun(){
122  
123   AliESDEvent *esdOFF = dynamic_cast<AliESDEvent*>(InputEvent());
124   TTimeStamp* timestamp = new TTimeStamp(esdOFF->GetTimeStamp());
125   fTextBox->SetName("text");
126   
127   TString s = Form("Run %d, Date %d", esdOFF->GetRunNumber(), timestamp->GetDate());
128   fTextBox->SetTitle(s);
129 }
130
131 void AliAnalysisTaskHLTCentralBarrel::UserExec(Option_t *){
132 // see header for documentation
133     
134   AliESDEvent *esdOFF = dynamic_cast<AliESDEvent*>(InputEvent());  
135   if(!esdOFF){
136      printf("ERROR: offline ESD is not available.\n");
137      return;
138   }
139       
140   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(fInputHandler);
141   AliESDEvent *esdHLT = NULL;   
142
143   if(esdH) esdHLT = esdH->GetHLTEvent();  
144   if(!esdHLT){
145      printf("ERROR: HLT ESD is not available.\n");
146      return;
147   }
148     
149   if(fUseHLTTrigger && !((AliHLTGlobalTriggerDecision*)esdHLT->GetHLTTriggerDecision())->Result()) return;
150   
151   if(fUseCentrality){
152     Int_t centbin = CalculateCentrality(esdHLT);
153     Printf("Centrality bin = %d", centbin);
154   }
155
156  
157   
158   //============================ OFFLINE =============================//
159
160   const AliESDVertex *vertOFF = esdOFF->GetPrimaryVertexTracks();
161   
162   Double_t bfield = esdOFF->GetMagneticField();
163   Int_t nr_tracksOFF = 0;
164   
165   if(esdOFF->GetEventSpecie()==16) return;
166    
167   for(Int_t i=0; i<esdOFF->GetNumberOfTracks(); i++){
168       AliESDtrack *esdTrackOFF = esdOFF->GetTrack(i);
169       if (!esdTrackOFF) continue;
170       if(!(esdTrackOFF->GetStatus()&AliESDtrack::kTPCin)) continue;
171       nr_tracksOFF++; 
172   }
173          
174   for(Int_t i=0; i<esdOFF->GetNumberOfTracks(); i++){
175       
176       AliESDtrack *esdTrackOFF = esdOFF->GetTrack(i);
177       if (!esdTrackOFF) continue;
178       if(!(esdTrackOFF->GetStatus()&AliESDtrack::kTPCin)) continue;
179
180       //DCA calculations(from offline)
181       Double_t x[3];
182       Double_t b[3]; 
183       esdTrackOFF->GetXYZ(x);
184      
185       AliTracker::GetBxByBz(x,b);
186       Bool_t isOK = esdTrackOFF->RelateToVertexTPCBxByBz(vertOFF, b, kVeryBig);
187       if(!isOK) continue;
188         
189       const AliExternalTrackParam *track = esdTrackOFF->GetTPCInnerParam();
190       if(!track) continue;
191         
192       Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
193       if(vertOFF->GetX()==0 && vertOFF->GetY()==0 && vertOFF->GetZ()==0 ){
194          dca[0]=-99;
195          dca[1]=-99;
196       }
197       else esdTrackOFF->GetImpactParametersTPC(dca,cov);
198
199       Float_t DCAr =-99, DCAz = -99.;   
200       Double_t trackOFF[] = {
201                                 TMath::Abs(esdTrackOFF->Pt()) 
202                                ,esdTrackOFF->GetTPCNcls()      
203                                ,esdTrackOFF->Theta()           
204                                ,esdTrackOFF->Eta()             
205                                ,esdTrackOFF->Phi()             
206                                ,dca[0]                         
207                                ,dca[1]                         
208                                ,esdTrackOFF->Charge() 
209                                ,DCAr                           
210                                ,DCAz                           
211                                ,esdTrackOFF->GetNcls(0)
212                                ,nr_tracksOFF
213                                ,vertOFF->GetStatus()
214                             };
215       fTrackOFF->Fill(trackOFF);
216     }
217     
218     Double_t eventOFF[] = { vertOFF->GetX(), vertOFF->GetY(), vertOFF->GetZ(), vertOFF->GetNContributors(), nr_tracksOFF, vertOFF->GetStatus()};
219     fEventOFF->Fill(eventOFF);  
220     
221   
222   
223   //======================================== HLT ==========================================//
224
225   Int_t nr_tracksHLT = 0;  
226   if(esdHLT->GetEventSpecie()==16) return;
227   const AliESDVertex *vertHLT = esdHLT->GetPrimaryVertexTracks();
228   
229   for(Int_t i=0; i<esdHLT->GetNumberOfTracks(); i++){
230       AliESDtrack *esdTrackHLT = esdHLT->GetTrack(i);
231       if (!esdTrackHLT) continue;
232       if(!(esdTrackHLT->GetStatus()&AliESDtrack::kTPCin)) continue;
233       nr_tracksHLT++; 
234   }
235   
236   Double_t eventHLT[] = { vertHLT->GetX(), vertHLT->GetY(), vertHLT->GetZ(), vertHLT->GetNContributors(), nr_tracksHLT, vertHLT->GetStatus()};
237   fEventHLT->Fill(eventHLT);  
238
239   for(Int_t i=0; i<esdHLT->GetNumberOfTracks(); i++){
240       
241       AliESDtrack *esdTrackHLT = esdHLT->GetTrack(i);
242       if(!esdTrackHLT) continue;
243       if(!(esdTrackHLT->GetStatus()&AliESDtrack::kTPCin)) continue;
244       
245       //DCA calculations
246       Float_t DCAr=-99;
247       Float_t DCAz=-99;
248       Float_t dca[2];
249       if(vertHLT->GetX()==0 && vertHLT->GetY()==0 && vertHLT->GetZ() ==0 ){
250         dca[0]=-99;
251         dca[1]=-99;
252       }
253       else{
254         //Calculating DCA "old" fashion
255         esdTrackHLT->GetDZ(esdHLT->GetPrimaryVertex()->GetXv(), esdHLT->GetPrimaryVertex()->GetYv(), esdHLT->GetPrimaryVertex()->GetZv(), bfield, dca);
256         // plotting the DCA calculated by Sergey 
257         esdTrackHLT->GetImpactParametersTPC(DCAr,DCAz);
258       }
259         
260       Double_t trackHLT[] = {
261                                TMath::Abs(esdTrackHLT->Pt())
262                               ,esdTrackHLT->GetTPCNcls()    
263                               ,esdTrackHLT->Theta()
264                               ,esdTrackHLT->Eta()           
265                               ,esdTrackHLT->Phi()
266                               ,dca[0]                       
267                               ,dca[1]                       
268                               ,esdTrackHLT->Charge()        
269                               ,DCAr                         
270                               ,DCAz                         
271                               ,esdTrackHLT->GetNcls(0)
272                               ,nr_tracksHLT
273                               ,vertHLT->GetStatus() 
274                             };
275       fTrackHLT->Fill(trackHLT);      
276   }               
277   // Post output data.
278   PostData(1, fOutputList);
279 }
280
281 void AliAnalysisTaskHLTCentralBarrel::Terminate(Option_t *){
282   // see header file of AliAnalysisTask for documentation 
283 }
284
285 Int_t AliAnalysisTaskHLTCentralBarrel::CalculateCentrality(AliESDEvent* esd){
286 //see header for documentation
287
288   Int_t centrality =-1;
289
290   AliESDVZERO* esdV0 = esd->GetVZEROData();
291   //AliESDZDC* esdZDC = esd->GetZDCData();
292   //Int_t partZDC = esdZDC->GetZDCParticipants();
293
294   Float_t multV0 = esdV0->GetMTotV0A() + esdV0->GetMTotV0C();
295     
296   if (      multV0 >=    0.  && multV0 <=   124.5 ) centrality = 90;
297   else if ( multV0 >   124.5 && multV0 <=   274.5 ) centrality = 80;
298   else if ( multV0 >   274.5 && multV0 <=   574.5 ) centrality = 70;
299   else if ( multV0 >   574.5 && multV0 <=  1224.5 ) centrality = 60;
300   else if ( multV0 >  1224.5 && multV0 <=  2174.5 ) centrality = 50;
301   else if ( multV0 >  2174.5 && multV0 <=  3624.5 ) centrality = 40;
302   else if ( multV0 >  3624.5 && multV0 <=  5574.5 ) centrality = 30;
303   else if ( multV0 >  5574.5 && multV0 <=  8274.5 ) centrality = 20;
304   else if ( multV0 >  8274.5 && multV0 <= 12024.5 ) centrality = 10;
305   else if ( multV0 > 12024.5 && multV0 <= 14674.5 ) centrality = 5;
306   else if ( multV0 > 14674.5 && multV0 <= 19449.5 ) centrality = 0;
307
308   return centrality;
309 }
310
311 THnSparseF* AliAnalysisTaskHLTCentralBarrel::CreateEventTHnSparse(const char* name, Int_t size, const Int_t* bins, Double_t* min, Double_t* max){
312 //see header for documentation                     
313   
314   THnSparseF *thn = new THnSparseF(name,"",size,bins,min,max);
315   thn->GetAxis(0)->SetTitle("primary vertex x");
316   thn->GetAxis(1)->SetTitle("primary vertex y");
317   thn->GetAxis(2)->SetTitle("primary vertex z");
318   thn->GetAxis(3)->SetTitle("number of contributors");
319   thn->GetAxis(4)->SetTitle("track multiplicity");
320   thn->GetAxis(5)->SetTitle("vertex status"); 
321   return thn;
322 }
323
324 THnSparseF* AliAnalysisTaskHLTCentralBarrel::CreateTrackTHnSparse(const char* name, Int_t size, const Int_t* bins, Double_t* min, Double_t* max){
325 //see header for documentation                     
326   
327   THnSparseF *thn = new THnSparseF(name,"",size,bins,min,max);
328   thn->GetAxis(0)->SetTitle("p_{T}");
329   thn->GetAxis(1)->SetTitle("TPC clusters/track");
330   thn->GetAxis(2)->SetTitle("#theta");
331   thn->GetAxis(3)->SetTitle("#eta");
332   thn->GetAxis(4)->SetTitle("#phi");
333   thn->GetAxis(5)->SetTitle("DCAr");
334   thn->GetAxis(6)->SetTitle("DCAz");
335   thn->GetAxis(7)->SetTitle("polarity");
336   thn->GetAxis(8)->SetTitle("DCArSG");
337   thn->GetAxis(9)->SetTitle("DCAzSG");
338   thn->GetAxis(10)->SetTitle("ITS clusters/track");  
339   return thn;
340 }
341
342 // void AliAnalysisTaskHLTCentralBarrel::Fill(AliESDevent *esd, THnSparseF* thn){
343 // //see header for documentation                     
344 //  
345 //   int nTracks = 0;  
346 //   
347 //   for(int i=0; i<esdHLT->GetNumberOfTracks(); i++){
348 //       
349 //       AliESDtrack *esdTrack = esd->GetTrack(i);
350 //       if(!esdTrack) continue;
351 //       if(!(esdTrack->GetStatus()&AliESDtrack::kTPCin)) continue;
352 //       
353 //       //DCA calculations
354 //       Float_t DCAr=-99;
355 //       Float_t DCAz=-99;
356 //       Float_t dca[2];
357 //       if(vertHLT->GetX()==0 && vertHLT->GetY()==0 && vertHLT->GetZ() ==0 ){
358 //      dca[0]=-99;
359 //      dca[1]=-99;
360 //       }
361 //       else{
362 //      //Calculating DCA "old" fashion
363 //      esdTrackHLT->GetDZ(esdHLT->GetPrimaryVertex()->GetXv(), esdHLT->GetPrimaryVertex()->GetYv(), esdHLT->GetPrimaryVertex()->GetZv(), bfield, dca);
364 //      // plotting the DCA calculated by Sergey 
365 //      esdTrackHLT->GetImpactParametersTPC(DCAr,DCAz);
366 //       }
367 //      
368 //       Double_t trackHLT[] = {
369 //                              TMath::Abs(esdTrackHLT->Pt())
370 //                             ,esdTrackHLT->GetTPCNcls()    
371 //                             ,esdTrackHLT->Theta()         
372 //                             ,esdTrackHLT->Eta()           
373 //                             ,esdTrackHLT->Phi()           
374 //                             ,DCAr                         
375 //                             ,DCAz                         
376 //                             ,esdTrackHLT->Charge()        
377 //                             ,dca[0]                       
378 //                             ,dca[1]                       
379 //                             ,esdTrackHLT->GetStatus()
380 //                             ,esdTrackHLT->GetNcls(0)   
381 //                           };
382 //       fTrackHLT->Fill(trackHLT);      
383 //       nr_tracksHLT++;    
384 //   }  
385 //   Double_t eventHLT[] = { vertHLT->GetX(), vertHLT->GetY(), vertHLT->GetZ(), vertHLT->GetNContributors(), nr_tracksHLT, vertHLT->GetStatus()};
386 //   fEventHLT->Fill(eventHLT);  
387 // }