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