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
  11 #include "AliAnalysisTaskAODFilterBitQA.h"
  12
  13 // Analysis task for the QA of AOD track filter bits
  14 // Authors: m.weber@cern.ch
  15
  16 ClassImp(AliAnalysisTaskAODFilterBitQA)
  17
  18 //________________________________________________________________________
  19 AliAnalysisTaskAODFilterBitQA::AliAnalysisTaskAODFilterBitQA(const char *name)
  20   : AliAnalysisTaskSE(name),
  21   fHistTrackStats(0)
  22 {
  23
  24   for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
  25     fHistKinematics[iTrackBit] = NULL;
  26     fHistDCA[iTrackBit]        = NULL;
  27     fHistDCAprop[iTrackBit]    = NULL;
  28     fHistChiClus[iTrackBit]    = NULL;
  29   }
  30
  31   DefineInput(0, TChain::Class());
  32   // Output slot #0 writes into a TH1 container
  33   DefineOutput(1, TList::Class());
  34
  35 }
  36
  37 //________________________________________________________________________
  38 AliAnalysisTaskAODFilterBitQA::~AliAnalysisTaskAODFilterBitQA() {
  39   // Destructor
  40   // ... not implemented
  41 }
  42
  43 //________________________________________________________________________
  44 void AliAnalysisTaskAODFilterBitQA::UserCreateOutputObjects() {
  45   // Create histograms
  46   // Called once
  47
  48    // global switch disabling the reference
  49   // (to avoid "Replacing existing TH1" if several wagons are created in train)
  50   Bool_t oldStatus = TH1::AddDirectoryStatus();
  51   TH1::AddDirectory(kFALSE);
  52
  53   // QA list
  54   fListQA = new TList();
  55   fListQA->SetName("listQA");
  56   fListQA->SetOwner();
  57
  58   // QA histograms
  59   fHistTrackStats = new TH2D("fHistTrackStats","Track statistics;Centrality;TrackFilterBit;N_{events}",100,0,100,gBitMax,0,gBitMax);
  60   fListQA->Add(fHistTrackStats);
  61
  62   for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
  63     fHistKinematics[iTrackBit] = new TH3D(Form("Bit%d_Kinematics",iTrackBit),Form("Bit%d_Kinematics;#eta;#varphi (rad);p_{T} (GeV/c)",iTrackBit),100,-1.0,1.0,100,0,TMath::Pi()*2,100,0,10);
  64     fHistDCA[iTrackBit]        = new TH2D(Form("Bit%d_DCA",iTrackBit),Form("Bit%d_DCA;DCA XY [AODtrack] (cm);DCA Z [AODtrack] (cm)",iTrackBit),100,-5.0,5.0,220,-11.0,11.0);
  65     fHistDCAprop[iTrackBit]    = new TH2D(Form("Bit%d_DCAprop",iTrackBit),Form("Bit%d_DCAprop;DCA XY [Propagated] (cm);DCA Z [Propagated] (cm)",iTrackBit),100,-5.0,5.0,220,-11.0,11.0);
  66     fHistChiClus[iTrackBit]    = new TH2D(Form("Bit%d_ChiClus",iTrackBit),Form("Bit%d_ChiClus;#chi^{2} [Fit];N_{clus} [TPC]",iTrackBit),100,-1.0,5.0,160,0,160.0);
  67     fListQA->Add(fHistKinematics[iTrackBit]);
  68     fListQA->Add(fHistDCA[iTrackBit]);
  69     fListQA->Add(fHistDCAprop[iTrackBit]);
  70     fListQA->Add(fHistChiClus[iTrackBit]);
  71   }
  72
  73   // Post output data.
  74   PostData(1, fListQA);
  75
  76   AliInfo("Finished setting up the Output");
  77   TH1::AddDirectory(oldStatus);
  78 }
  79
  80 //________________________________________________________________________
  81 void AliAnalysisTaskAODFilterBitQA::UserExec(Option_t *) {
  82   // Main loop
  83   // Called for each event
  84
  85   AliVEvent* event = dynamic_cast<AliVEvent*>(InputEvent());
  86   if(!event) {
  87     AliError("event not available");
  88     return;
  89   }
  90
  91
  92
  93   // check event cuts
  94   Double_t lMultiplicityVar = -1;
  95   if((lMultiplicityVar = IsEventAccepted(event)) < 0){
  96     return;
  97   }
  98
  99
 100   // get the accepted tracks in main event
 101   GetAcceptedTracks(event,lMultiplicityVar);
 102
 103 }
 104
 105 //________________________________________________________________________
 106 void  AliAnalysisTaskAODFilterBitQA::FinishTaskOutput(){
 107   // Finish task output
 108   // not implemented ...
 109
 110 }
 111
 112 //________________________________________________________________________
 113 void AliAnalysisTaskAODFilterBitQA::Terminate(Option_t *) {
 114   // Draw result to the screen
 115   // Called once at the end of the query
 116   // not implemented ...
 117
 118 }
 119
 120
 121 //________________________________________________________________________
 122 Double_t AliAnalysisTaskAODFilterBitQA::IsEventAccepted(AliVEvent *event){
 123   // Checks the Event cuts
 124
 125   // still hard coded
 126   Double_t fCentralityPercentileMin = 0.;
 127   Double_t fCentralityPercentileMax = 80.;
 128   Double_t fVxMax = 0.5;
 129   Double_t fVyMax = 0.5;
 130   Double_t fVzMax = 10.0;
 131   TString fCentralityEstimator = "V0M";
 132
 133   Float_t gCentrality      = -1.;
 134   const AliVVertex *vertex = event->GetPrimaryVertex();
 135
 136   if(vertex) {
 137     Double32_t fCov[6];
 138     vertex->GetCovarianceMatrix(fCov);
 139     if(vertex->GetNContributors() > 0) {
 140       if(fCov[5] != 0) {
 141         if(TMath::Abs(vertex->GetX()) < fVxMax) {
 142           if(TMath::Abs(vertex->GetY()) < fVyMax) {
 143             if(TMath::Abs(vertex->GetZ()) < fVzMax) {
 144
 145               // get the reference multiplicty or centrality
 146               AliAODHeader *header = (AliAODHeader*) event->GetHeader();
 147               gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 148
 149               if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){
 150
 151                 return gCentrality;
 152
 153               }//centrality range
 154             }//Vz cut
 155           }//Vy cut
 156         }//Vx cut
 157       }//proper vertex resolution
 158     }//proper number of contributors
 159   }//vertex object valid
 160
 161   // in all other cases return -1 (event not accepted)
 162   return -1;
 163 }
 164
 165 //________________________________________________________________________
 166 void AliAnalysisTaskAODFilterBitQA::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality){
 167   // Checks track cuts (filter bits)
 168   // Fills QA histograms
 169
 170
 171   Short_t  vCharge;
 172   Double_t vEta;
 173   Double_t vY;
 174   Double_t vPhi;
 175   Double_t vPt;
 176   Double_t vDCAxy;
 177   Double_t vDCAz;
 178   Double_t vDCApropxy;
 179   Double_t vDCApropz;
 180   Double_t vChi2;
 181   Double_t vClus;
 182
 183   //for propagation to DCA
 184   const AliVVertex *vertex = event->GetPrimaryVertex();
 185   Double_t field = event->GetMagneticField();
 186   Double_t b[2];
 187   Double_t bCov[3];
 188
 189   // Loop over tracks in event
 190   for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
 191     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
 192     if (!aodTrack) {
 193       AliError(Form("Could not receive track %d", iTracks));
 194       continue;
 195     }
 196
 197     // track parameters
 198     vCharge = aodTrack->Charge();
 199     vEta    = aodTrack->Eta();
 200     vY      = aodTrack->Y();
 201     vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
 202     vPt     = aodTrack->Pt();
 203     vDCAxy  = aodTrack->DCA();
 204     vDCAz   = aodTrack->ZAtDCA();
 205     vChi2   = aodTrack->Chi2perNDF();
 206     vClus   = aodTrack->GetTPCNcls();
 207
 208     // propagate again to DCA, use copy to propagate (in order not to overwrite track parameters)
 209     AliAODTrack copy(*aodTrack);
 210     if (copy.PropagateToDCA(vertex, field, 100., b, bCov) )
 211       {
 212         vDCApropxy = b[0];
 213         vDCApropz  = b[1];
 214       }
 215     else
 216       {
 217         vDCApropxy = -999999;
 218         vDCApropz  = -999999;
 219       }
 220
 221     // AOD track cuts
 222     for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
 223       fHistTrackStats->Fill(gCentrality,iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
 224
 225       if(aodTrack->TestFilterBit(1<<iTrackBit)){
 226         fHistKinematics[iTrackBit]->Fill(vEta,vPhi,vPt);
 227         fHistDCA[iTrackBit]->Fill(vDCAxy,vDCAz);
 228         fHistDCAprop[iTrackBit]->Fill(vDCApropxy,vDCApropz);
 229         fHistChiClus[iTrackBit]->Fill(vChi2,vClus);
 230       }
 231
 232     }//bit loop
 233   }//track loop
 234
 235   return;
 236 }