PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskAODFilterBitQA.cxx

   1 #include "TChain.h"
   2 #include "TH2D.h"
   3 #include "TH3D.h"
   4
   5
   6 #include "AliAODEvent.h"
   7 #include "AliAODRecoDecay.h"
   8 #include "AliAODTrack.h"
   9 #include "AliLog.h"
  10
  11 #include "AliAnalysisTaskSE.h"
  12 #include "AliMCEvent.h"
  13 #include "AliAODMCParticle.h"
  14
  15
  16
  17 #include "AliAnalysisTaskAODFilterBitQA.h"
  18
  19 // Analysis task for the QA of AOD track filter bits
  20 // Authors: m.weber@cern.ch
  21
  22 ClassImp(AliAnalysisTaskAODFilterBitQA)
  23
  24 //________________________________________________________________________
  25 AliAnalysisTaskAODFilterBitQA::AliAnalysisTaskAODFilterBitQA(const char *name)
  26   : AliAnalysisTaskSE(name),
  27   fArrayMC(0x0),
  28   fListQA(0x0),
  29   useCentrality(kFALSE),
  30   fillOnlySecondaries(kFALSE),
  31   fillHFVertexingTracks(kFALSE),
  32   fHFBranchName("D0toKpi"),
  33   fBitIgnore1(-1),
  34   fBitIgnore2(-1),
  35   fCentralityPercentileMin(0.),
  36   fCentralityPercentileMax(80.),
  37   fPtMin(0),
  38   fPtMax(1000),
  39   fEtaMin(-10),
  40   fEtaMax(10),
  41   fHistTrackStats(0)
  42 {
  43
  44   for(Int_t iCharge = 0; iCharge < gNCharge; iCharge++){
  45     for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
  46       fHistKinematics[iCharge][iTrackBit] = NULL;
  47       fHistDCAconstrained[iCharge][iTrackBit] = NULL;
  48       fHistDCAglobal[iCharge][iTrackBit]  = NULL;
  49       fHistChiClus[iCharge][iTrackBit]    = NULL;
  50     }
  51   }
  52
  53   DefineInput(0, TChain::Class());
  54   // Output slot #0 writes into a TH1 container
  55   DefineOutput(1, TList::Class());
  56
  57 }
  58
  59 //________________________________________________________________________
  60 AliAnalysisTaskAODFilterBitQA::~AliAnalysisTaskAODFilterBitQA() {
  61   // Destructor
  62   // ... not implemented
  63 }
  64
  65 //________________________________________________________________________
  66 void AliAnalysisTaskAODFilterBitQA::UserCreateOutputObjects() {
  67   // Create histograms
  68   // Called once
  69
  70    // global switch disabling the reference
  71   // (to avoid "Replacing existing TH1" if several wagons are created in train)
  72   Bool_t oldStatus = TH1::AddDirectoryStatus();
  73   TH1::AddDirectory(kFALSE);
  74
  75   // QA list
  76   fListQA = new TList();
  77   fListQA->SetName("listQA");
  78   fListQA->SetOwner();
  79
  80   // QA histograms
  81   fHistTrackStats = new TH2D("fHistTrackStats","Track statistics;Centrality;TrackFilterBit;N_{events}",100,0,100,gBitMax,0,gBitMax);
  82   fListQA->Add(fHistTrackStats);
  83
  84   TString sCharge[gNCharge] = {"Plus","Minus"};
  85
  86   for(Int_t iCharge = 0; iCharge < gNCharge; iCharge++){
  87     for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
  88       fHistKinematics[iCharge][iTrackBit] = new TH3D(Form("Bit%d_%s_Kinematics",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_Kinematics;#eta;#varphi (rad);p_{T} (GeV/c)",iTrackBit,sCharge[iCharge].Data()),100,-1.0,1.0,100,0,TMath::Pi()*2,100,0,10);
  89       fHistDCAconstrained[iCharge][iTrackBit] = new TH2D(Form("Bit%d_%s_DCAconstrained",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_DCAconstrained;DCA XY [Constrained] (cm);DCA Z [Constrained] (cm)",iTrackBit,sCharge[iCharge].Data()),100,-5.0,5.0,100,-5.0,5.0);
  90       fHistDCAglobal[iCharge][iTrackBit]  = new TH3D(Form("Bit%d_%s_DCAglobal",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_DCAglobal;DCA X [Global] (cm);DCA Y [Global] (cm);DCA Z [Global] (cm)",iTrackBit,sCharge[iCharge].Data()),100,-5.0,5.0,100,-5.0,5.0,100,-5.0,5.0);
  91 fHistChiClus[iCharge][iTrackBit]    = new TH2D(Form("Bit%d_%s_ChiClus",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_ChiClus;#chi^{2} [Fit];N_{clus} [TPC]",iTrackBit,sCharge[iCharge].Data()),100,-1.0,5.0,160,0,160.0);
  92       fListQA->Add(fHistKinematics[iCharge][iTrackBit]);
  93       fListQA->Add(fHistDCAconstrained[iCharge][iTrackBit]);
  94       fListQA->Add(fHistDCAglobal[iCharge][iTrackBit]);
  95       fListQA->Add(fHistChiClus[iCharge][iTrackBit]);
  96     }
  97   }
  98
  99   // Post output data.
 100   PostData(1, fListQA);
 101
 102   AliInfo("Finished setting up the Output");
 103   TH1::AddDirectory(oldStatus);
 104 }
 105
 106 //________________________________________________________________________
 107 void AliAnalysisTaskAODFilterBitQA::UserExec(Option_t *) {
 108   // Main loop
 109   // Called for each event
 110
 111   AliVEvent* event = dynamic_cast<AliVEvent*>(InputEvent());
 112   if(!event) {
 113     AliError("event not available");
 114     return;
 115   }
 116
 117   // MC information (set if available)
 118   fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
 119
 120   // check event cuts
 121   Double_t lMultiplicityVar = -1;
 122   if((lMultiplicityVar = IsEventAccepted(event)) < 0){
 123     return;
 124   }
 125
 126   // fill HF vertexing tracks
 127   if(fillHFVertexingTracks){
 128     GetAcceptedHFVertexingTracks(event,lMultiplicityVar);
 129   }
 130   else{
 131     // get the accepted tracks in main event
 132     GetAcceptedTracks(event,lMultiplicityVar);
 133   }
 134 }
 135
 136 //________________________________________________________________________
 137 void  AliAnalysisTaskAODFilterBitQA::FinishTaskOutput(){
 138   // Finish task output
 139   // not implemented ...
 140
 141 }
 142
 143 //________________________________________________________________________
 144 void AliAnalysisTaskAODFilterBitQA::Terminate(Option_t *) {
 145   // Draw result to the screen
 146   // Called once at the end of the query
 147   // not implemented ...
 148
 149 }
 150
 151
 152 //________________________________________________________________________
 153 Double_t AliAnalysisTaskAODFilterBitQA::IsEventAccepted(AliVEvent *event){
 154   // Checks the Event cuts
 155
 156   // still hard coded
 157   Double_t fVxMax = 0.5;
 158   Double_t fVyMax = 0.5;
 159   Double_t fVzMax = 10.0;
 160   TString fCentralityEstimator = "V0M";
 161
 162   Float_t gCentrality      = -1.;
 163   const AliVVertex *vertex = event->GetPrimaryVertex();
 164
 165   if(vertex) {
 166     Double32_t fCov[6];
 167     vertex->GetCovarianceMatrix(fCov);
 168     if(vertex->GetNContributors() > 0) {
 169       if(fCov[5] != 0) {
 170         if(TMath::Abs(vertex->GetX()) < fVxMax) {
 171           if(TMath::Abs(vertex->GetY()) < fVyMax) {
 172             if(TMath::Abs(vertex->GetZ()) < fVzMax) {
 173
 174               // get the reference multiplicty or centrality (if required)
 175               if(useCentrality){
 176
 177                 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
 178                 gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 179
 180                 if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){
 181                   return gCentrality;
 182                 }
 183                 else{
 184                   return -1;
 185                 }//centrality range
 186               }//use centrality
 187
 188               // if not using centrality/multiplicty, return 1
 189               return 1;
 190
 191             }//Vz cut
 192           }//Vy cut
 193         }//Vx cut
 194       }//proper vertex resolution
 195     }//proper number of contributors
 196   }//vertex object valid
 197
 198   // in all other cases return -1 (event not accepted)
 199   return -1;
 200 }
 201
 202 //________________________________________________________________________
 203 void AliAnalysisTaskAODFilterBitQA::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality){
 204   // Checks track cuts (filter bits)
 205   // Fills QA histograms
 206
 207
 208   Double_t vEta;
 209   Double_t vPhi;
 210   Double_t vPt;
 211   Double_t vDCAconstrainedxy;
 212   Double_t vDCAconstrainedz;
 213   Double_t vDCAglobalx;
 214   Double_t vDCAglobaly;
 215   Double_t vDCAglobalz;
 216   Double_t vChi2;
 217   Double_t vClus;
 218
 219   Double_t pos[3];
 220   Double_t v[3];
 221
 222   const AliVVertex *vertex = event->GetPrimaryVertex();
 223   vertex->GetXYZ(v);
 224
 225   // Loop over tracks in event
 226   for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
 227     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
 228     if (!aodTrack) {
 229       AliError(Form("Could not receive track %d", iTracks));
 230       continue;
 231     }
 232
 233     // get MC information (if available)
 234     if(fArrayMC && fillOnlySecondaries){
 235
 236       Int_t label = aodTrack->GetLabel();
 237       AliAODMCParticle *mcTrack = (AliAODMCParticle *)fArrayMC->At(TMath::Abs(label));
 238
 239       if(mcTrack->IsPhysicalPrimary())
 240         continue;
 241     }
 242
 243     // track parameters
 244     vEta    = aodTrack->Eta();
 245     vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
 246     vPt     = aodTrack->Pt();
 247     vDCAconstrainedxy  = aodTrack->DCA();
 248     vDCAconstrainedz   = aodTrack->ZAtDCA();
 249     vChi2   = aodTrack->Chi2perNDF();
 250     vClus   = aodTrack->GetTPCNcls();
 251
 252     // kinematic cuts
 253     if( vPt > fPtMax || vPt < fPtMin )
 254       continue;
 255     if( vEta > fEtaMax || vEta < fEtaMin )
 256       continue;
 257
 258     // if not constrained track the position is stored (primary vertex to be subtracted)
 259     aodTrack->GetXYZ(pos);
 260     vDCAglobalx  = pos[0] - v[0];
 261     vDCAglobaly  = pos[1] - v[1];
 262     vDCAglobalz  = pos[2] - v[2];
 263
 264     // fill for separately for positive and negative charges
 265     Int_t iCharge = -1;
 266     // positive
 267     if(aodTrack->Charge() > 0)
 268       iCharge = 0;
 269     else if(aodTrack->Charge() < 0)
 270       iCharge = 1;
 271     else{
 272       AliError("Charge==0?");
 273       iCharge = -1;
 274     }
 275
 276
 277     // AOD track cuts
 278     if(iCharge > -1){
 279       for(Int_t iTrackBit = 0; iTrackBit < gBitMax-1; iTrackBit++){
 280         fHistTrackStats->Fill(gCentrality,iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
 281
 282         if(aodTrack->TestFilterBit(1<<iTrackBit)){
 283           fHistKinematics[iCharge][iTrackBit]->Fill(vEta,vPhi,vPt);
 284           fHistDCAconstrained[iCharge][iTrackBit]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
 285           fHistDCAglobal[iCharge][iTrackBit]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
 286           fHistChiClus[iCharge][iTrackBit]->Fill(vChi2,vClus);
 287         }
 288       }//bit loop
 289
 290       // fill all tracks in last bit
 291       fHistTrackStats->Fill(gCentrality,gBitMax-1,1);
 292       fHistKinematics[iCharge][gBitMax-1]->Fill(vEta,vPhi,vPt);
 293       fHistDCAconstrained[iCharge][gBitMax-1]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
 294       fHistDCAglobal[iCharge][gBitMax-1]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
 295       fHistChiClus[iCharge][gBitMax-1]->Fill(vChi2,vClus);
 296
 297     }//charge positive or negative
 298   }//track loop
 299
 300   return;
 301 }
 302
 303
 304 //________________________________________________________________________
 305 void AliAnalysisTaskAODFilterBitQA::GetAcceptedHFVertexingTracks(AliVEvent *event, Double_t gCentrality){
 306   // Checks track cuts (filter bits)
 307   // from daughters of HF candidates
 308   // Fills QA histograms
 309
 310   Double_t vEta;
 311   Double_t vPhi;
 312   Double_t vPt;
 313   Double_t vDCAconstrainedxy;
 314   Double_t vDCAconstrainedz;
 315   Double_t vDCAglobalx;
 316   Double_t vDCAglobaly;
 317   Double_t vDCAglobalz;
 318   Double_t vChi2;
 319   Double_t vClus;
 320
 321   Double_t pos[3];
 322   Double_t v[3];
 323
 324   Int_t IDList[1000];
 325   Int_t IDListLength = 0;
 326   for(Int_t i = 0; i < 1000; i++){
 327     IDList[i] = -5;
 328   }
 329
 330   const AliVVertex *vertex = event->GetPrimaryVertex();
 331   vertex->GetXYZ(v);
 332
 333   // =================================================================================
 334   // HF part (taken from AliAnalysisTaskSEDmesonsFilterCJ)
 335
 336   TClonesArray *arrayDStartoD0pi = (TClonesArray*)event->GetList()->FindObject(fHFBranchName.Data());
 337
 338   if (!arrayDStartoD0pi) {
 339     AliInfo(Form("Could not find array %s, skipping the event",fHFBranchName.Data()));
 340     return;
 341   } else {
 342     AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
 343   }
 344
 345   // loop over tracks
 346    const Int_t NVertices = arrayDStartoD0pi->GetEntriesFast();
 347
 348    for (Int_t iVertex = 0; iVertex < NVertices; ++iVertex) {
 349
 350      AliAODRecoDecay *HFvertex = static_cast<AliAODRecoDecay*>(arrayDStartoD0pi->At(iVertex));
 351
 352      // Loop over tracks (daughters of D candidates)
 353      for (Int_t iProng = 0; iProng<HFvertex->GetNProngs(); iProng++) {
 354
 355        AliAODTrack *daughter = (AliAODTrack*)HFvertex->GetDaughter(iProng);
 356        if (!daughter){
 357          AliError(Form("Could not receive track %d %d", iVertex, iProng));
 358          continue;
 359        }
 360
 361
 362        // get MC information (if available)
 363        if(fArrayMC && fillOnlySecondaries){
 364
 365          Int_t label = daughter->GetLabel();
 366          AliAODMCParticle *mcTrack = (AliAODMCParticle *)fArrayMC->At(TMath::Abs(label));
 367
 368          if(mcTrack->IsPhysicalPrimary())
 369            continue;
 370        }
 371
 372        // track parameters
 373        vEta    = daughter->Eta();
 374        vPhi    = daughter->Phi();// * TMath::RadToDeg();
 375        vPt     = daughter->Pt();
 376        vDCAconstrainedxy  = daughter->DCA();
 377        vDCAconstrainedz   = daughter->ZAtDCA();
 378        vChi2   = daughter->Chi2perNDF();
 379        vClus   = daughter->GetTPCNcls();
 380
 381        // kinematic cuts
 382        if( vPt > fPtMax || vPt < fPtMin )
 383          continue;
 384        if( vEta > fEtaMax || vEta < fEtaMin )
 385          continue;
 386
 387        // avoid double counting (can be optimized)
 388        Bool_t doubleCount = kFALSE;
 389        Int_t  daughterID  = daughter->GetID();
 390        for(Int_t idx = 0; idx < IDListLength; idx++){
 391          if(IDList[idx]==daughterID){
 392            doubleCount = kTRUE;
 393            break;
 394          }
 395        }
 396        if(!doubleCount){
 397          IDList[IDListLength] = daughterID;
 398          IDListLength++;
 399        }
 400        else{
 401          continue;
 402        }
 403
 404
 405        // if not constrained track the position is stored (primary vertex to be subtracted)
 406        daughter->GetXYZ(pos);
 407        vDCAglobalx  = pos[0] - v[0];
 408        vDCAglobaly  = pos[1] - v[1];
 409        vDCAglobalz  = pos[2] - v[2];
 410
 411        // fill for separately for positive and negative charges
 412        Int_t iCharge = -1;
 413        // positive
 414        if(daughter->Charge() > 0)
 415          iCharge = 0;
 416        else if(daughter->Charge() < 0)
 417          iCharge = 1;
 418        else{
 419          AliError("Charge==0?");
 420          iCharge = -1;
 421        }
 422
 423
 424        // AOD track cuts
 425        if(iCharge > -1){
 426
 427          // if some filter bits should be ignored, skip them here
 428          if(fBitIgnore1 > -1 && daughter->TestFilterBit(1<<fBitIgnore1))
 429            continue;
 430          if(fBitIgnore2 > -1 && daughter->TestFilterBit(1<<fBitIgnore2))
 431            continue;
 432
 433          for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
 434            fHistTrackStats->Fill(gCentrality,iTrackBit,daughter->TestFilterBit(1<<iTrackBit));
 435
 436            if(daughter->TestFilterBit(1<<iTrackBit)){
 437              fHistKinematics[iCharge][iTrackBit]->Fill(vEta,vPhi,vPt);
 438              fHistDCAconstrained[iCharge][iTrackBit]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
 439              fHistDCAglobal[iCharge][iTrackBit]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
 440              fHistChiClus[iCharge][iTrackBit]->Fill(vChi2,vClus);
 441            }
 442          }//bit loop
 443
 444          // fill all tracks in last bit
 445          fHistTrackStats->Fill(gCentrality,gBitMax-1,1);
 446          fHistKinematics[iCharge][gBitMax-1]->Fill(vEta,vPhi,vPt);
 447          fHistDCAconstrained[iCharge][gBitMax-1]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
 448          fHistDCAglobal[iCharge][gBitMax-1]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
 449          fHistChiClus[iCharge][gBitMax-1]->Fill(vChi2,vClus);
 450
 451        }//charge positive or negative
 452
 453      }//prong loop
 454    }//HF vertex loop
 455
 456   return;
 457 }