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 "AliAODTrack.h"
   8 #include "AliLog.h"
   9
  10 #include "AliAnalysisTaskSE.h"
  11 #include "AliMCEvent.h"
  12 #include "AliAODMCParticle.h"
  13
  14
  15
  16 #include "AliAnalysisTaskAODFilterBitQA.h"
  17
  18 // Analysis task for the QA of AOD track filter bits
  19 // Authors: m.weber@cern.ch
  20
  21 ClassImp(AliAnalysisTaskAODFilterBitQA)
  22
  23 //________________________________________________________________________
  24 AliAnalysisTaskAODFilterBitQA::AliAnalysisTaskAODFilterBitQA(const char *name)
  25   : AliAnalysisTaskSE(name),
  26   fArrayMC(0x0),
  27   fListQA(0x0),
  28   fillOnlySecondaries(kFALSE),
  29   fCentralityPercentileMin(0.),
  30   fCentralityPercentileMax(80.),
  31   fPtMin(0),
  32   fPtMax(1000),
  33   fEtaMin(-10),
  34   fEtaMax(10),
  35   fHistTrackStats(0)
  36 {
  37
  38   for(Int_t iCharge = 0; iCharge < gNCharge; iCharge++){
  39     for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
  40       fHistKinematics[iCharge][iTrackBit] = NULL;
  41       fHistDCAconstrained[iCharge][iTrackBit] = NULL;
  42       fHistDCAglobal[iCharge][iTrackBit]  = NULL;
  43       fHistChiClus[iCharge][iTrackBit]    = NULL;
  44     }
  45   }
  46
  47   DefineInput(0, TChain::Class());
  48   // Output slot #0 writes into a TH1 container
  49   DefineOutput(1, TList::Class());
  50
  51 }
  52
  53 //________________________________________________________________________
  54 AliAnalysisTaskAODFilterBitQA::~AliAnalysisTaskAODFilterBitQA() {
  55   // Destructor
  56   // ... not implemented
  57 }
  58
  59 //________________________________________________________________________
  60 void AliAnalysisTaskAODFilterBitQA::UserCreateOutputObjects() {
  61   // Create histograms
  62   // Called once
  63
  64    // global switch disabling the reference
  65   // (to avoid "Replacing existing TH1" if several wagons are created in train)
  66   Bool_t oldStatus = TH1::AddDirectoryStatus();
  67   TH1::AddDirectory(kFALSE);
  68
  69   // QA list
  70   fListQA = new TList();
  71   fListQA->SetName("listQA");
  72   fListQA->SetOwner();
  73
  74   // QA histograms
  75   fHistTrackStats = new TH2D("fHistTrackStats","Track statistics;Centrality;TrackFilterBit;N_{events}",100,0,100,gBitMax,0,gBitMax);
  76   fListQA->Add(fHistTrackStats);
  77
  78   TString sCharge[gNCharge] = {"Plus","Minus"};
  79
  80   for(Int_t iCharge = 0; iCharge < gNCharge; iCharge++){
  81     for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
  82       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);
  83       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);
  84       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);
  85 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);
  86       fListQA->Add(fHistKinematics[iCharge][iTrackBit]);
  87       fListQA->Add(fHistDCAconstrained[iCharge][iTrackBit]);
  88       fListQA->Add(fHistDCAglobal[iCharge][iTrackBit]);
  89       fListQA->Add(fHistChiClus[iCharge][iTrackBit]);
  90     }
  91   }
  92
  93   // Post output data.
  94   PostData(1, fListQA);
  95
  96   AliInfo("Finished setting up the Output");
  97   TH1::AddDirectory(oldStatus);
  98 }
  99
 100 //________________________________________________________________________
 101 void AliAnalysisTaskAODFilterBitQA::UserExec(Option_t *) {
 102   // Main loop
 103   // Called for each event
 104
 105   AliVEvent* event = dynamic_cast<AliVEvent*>(InputEvent());
 106   if(!event) {
 107     AliError("event not available");
 108     return;
 109   }
 110
 111   // MC information (set if available)
 112   fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
 113
 114   // check event cuts
 115   Double_t lMultiplicityVar = -1;
 116   if((lMultiplicityVar = IsEventAccepted(event)) < 0){
 117     return;
 118   }
 119
 120
 121   // get the accepted tracks in main event
 122   GetAcceptedTracks(event,lMultiplicityVar);
 123
 124 }
 125
 126 //________________________________________________________________________
 127 void  AliAnalysisTaskAODFilterBitQA::FinishTaskOutput(){
 128   // Finish task output
 129   // not implemented ...
 130
 131 }
 132
 133 //________________________________________________________________________
 134 void AliAnalysisTaskAODFilterBitQA::Terminate(Option_t *) {
 135   // Draw result to the screen
 136   // Called once at the end of the query
 137   // not implemented ...
 138
 139 }
 140
 141
 142 //________________________________________________________________________
 143 Double_t AliAnalysisTaskAODFilterBitQA::IsEventAccepted(AliVEvent *event){
 144   // Checks the Event cuts
 145
 146   // still hard coded
 147   Double_t fVxMax = 0.5;
 148   Double_t fVyMax = 0.5;
 149   Double_t fVzMax = 10.0;
 150   TString fCentralityEstimator = "V0M";
 151
 152   Float_t gCentrality      = -1.;
 153   const AliVVertex *vertex = event->GetPrimaryVertex();
 154
 155   if(vertex) {
 156     Double32_t fCov[6];
 157     vertex->GetCovarianceMatrix(fCov);
 158     if(vertex->GetNContributors() > 0) {
 159       if(fCov[5] != 0) {
 160         if(TMath::Abs(vertex->GetX()) < fVxMax) {
 161           if(TMath::Abs(vertex->GetY()) < fVyMax) {
 162             if(TMath::Abs(vertex->GetZ()) < fVzMax) {
 163
 164               // get the reference multiplicty or centrality
 165               AliAODHeader *header = (AliAODHeader*) event->GetHeader();
 166               gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 167
 168               if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){
 169
 170                 return gCentrality;
 171
 172               }//centrality range
 173             }//Vz cut
 174           }//Vy cut
 175         }//Vx cut
 176       }//proper vertex resolution
 177     }//proper number of contributors
 178   }//vertex object valid
 179
 180   // in all other cases return -1 (event not accepted)
 181   return -1;
 182 }
 183
 184 //________________________________________________________________________
 185 void AliAnalysisTaskAODFilterBitQA::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality){
 186   // Checks track cuts (filter bits)
 187   // Fills QA histograms
 188
 189
 190   Short_t  vCharge;
 191   Double_t vEta;
 192   Double_t vY;
 193   Double_t vPhi;
 194   Double_t vPt;
 195   Double_t vDCAconstrainedxy;
 196   Double_t vDCAconstrainedz;
 197   Double_t vDCAglobalx;
 198   Double_t vDCAglobaly;
 199   Double_t vDCAglobalz;
 200   Double_t vChi2;
 201   Double_t vClus;
 202
 203   Double_t pos[3];
 204   Double_t v[3];
 205
 206   const AliVVertex *vertex = event->GetPrimaryVertex();
 207   vertex->GetXYZ(v);
 208
 209   // Loop over tracks in event
 210   for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
 211     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
 212     if (!aodTrack) {
 213       AliError(Form("Could not receive track %d", iTracks));
 214       continue;
 215     }
 216
 217     // get MC information (if available)
 218     if(fArrayMC && fillOnlySecondaries){
 219
 220       Int_t label = aodTrack->GetLabel();
 221       AliAODMCParticle *mcTrack = (AliAODMCParticle *)fArrayMC->At(TMath::Abs(label));
 222
 223       if(mcTrack->IsPhysicalPrimary())
 224         continue;
 225     }
 226
 227     // track parameters
 228     vCharge = aodTrack->Charge();
 229     vEta    = aodTrack->Eta();
 230     vY      = aodTrack->Y();
 231     vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
 232     vPt     = aodTrack->Pt();
 233     vDCAconstrainedxy  = aodTrack->DCA();
 234     vDCAconstrainedz   = aodTrack->ZAtDCA();
 235     vChi2   = aodTrack->Chi2perNDF();
 236     vClus   = aodTrack->GetTPCNcls();
 237
 238     // kinematic cuts
 239     if( vPt > fPtMax || vPt < fPtMin )
 240       continue;
 241     if( vEta > fEtaMax || vEta < fEtaMin )
 242       continue;
 243
 244     // if not constrained track the position is stored (primary vertex to be subtracted)
 245     aodTrack->GetXYZ(pos);
 246     vDCAglobalx  = pos[0] - v[0];
 247     vDCAglobaly  = pos[1] - v[1];
 248     vDCAglobalz  = pos[2] - v[2];
 249
 250     // fill for separately for positive and negative charges
 251     Int_t iCharge = -1;
 252     // positive
 253     if(aodTrack->Charge() > 0)
 254       iCharge = 0;
 255     else if(aodTrack->Charge() < 0)
 256       iCharge = 1;
 257     else{
 258       AliError("Charge==0?");
 259       iCharge = -1;
 260     }
 261
 262
 263     // AOD track cuts
 264     if(iCharge > -1){
 265       for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
 266         fHistTrackStats->Fill(gCentrality,iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
 267
 268         if(aodTrack->TestFilterBit(1<<iTrackBit)){
 269           fHistKinematics[iCharge][iTrackBit]->Fill(vEta,vPhi,vPt);
 270           fHistDCAconstrained[iCharge][iTrackBit]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
 271           fHistDCAglobal[iCharge][iTrackBit]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
 272           fHistChiClus[iCharge][iTrackBit]->Fill(vChi2,vClus);
 273         }
 274       }//bit loop
 275     }//charge positive or negative
 276   }//track loop
 277
 278   return;
 279 }