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()->GetX(), esdHLT->GetPrimaryVertex()->GetY(), esdHLT->GetPrimaryVertex()->GetZ(), 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()->GetX(), esdHLT->GetPrimaryVertex()->GetY(), esdHLT->GetPrimaryVertex()->GetZ(), 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 // }