PWGCF/EBYE/PIDFluctuation/task/AliEbyEParticleRatioFluctuationTask.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Satyajit Jena.                                                 *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 //=========================================================================//
  18 //             AliEbyE Analysis for Particle Ratio Fluctuation             //
  19 //                   Deepika Rathee  | Satyajit Jena                       //
  20 //                   drathee@cern.ch | sjena@cern.ch                       //
  21 //
  22 //             V0.1 2013/03/25 Using THnSparse
  23 //             V0.2 2013/04/03 Cleanup
  24 //             V1.0 2013/04/10 Cleanup Bins for Less Memory
  25 //             V1.1 2013/04/15 Bins Added
  26 //   Todo: pp and pA, Mix Events
  27 //=========================================================================//
  28
  29 #include "TChain.h"
  30 #include "TList.h"
  31 #include "TFile.h"
  32 #include "TTree.h"
  33 #include "TH1D.h"
  34 #include "TH2F.h"
  35 #include "TH3F.h"
  36 #include "TCanvas.h"
  37 #include "AliAnalysisTask.h"
  38 #include "AliAnalysisManager.h"
  39 #include "AliVEvent.h"
  40 #include "AliESD.h"
  41 #include "AliESDEvent.h"
  42 #include "AliAODEvent.h"
  43 #include "AliAODMCParticle.h"
  44 #include "AliAODMCHeader.h"
  45 #include "AliPIDResponse.h"
  46 #include "AliAODHeader.h"
  47 #include "AliAODpidUtil.h"
  48 #include "AliHelperPID.h"
  49 using std::endl;
  50 using std::cout;
  51 #include "AliEbyEParticleRatioFluctuationTask.h"
  52
  53 ClassImp(AliEbyEParticleRatioFluctuationTask)
  54
  55 //-----------------------------------------------------------------------
  56 AliEbyEParticleRatioFluctuationTask::AliEbyEParticleRatioFluctuationTask( const char *name ) : AliAnalysisTaskSE( name ),
  57   fThnList(NULL),
  58   fAnalysisType("AOD"),
  59   fAnalysisData("PbPb"),
  60   fCentralityEstimator("V0M"),
  61   fVxMax(3.),
  62   fVyMax(3.),
  63   fVzMax(10.),
  64   fDCAxy(2.4),
  65   fDCAz(3.2),
  66   fPtLowerLimit(0.2),
  67   fPtHigherLimit(5.),
  68   fEtaLowerLimit(-1.),
  69   fEtaHigherLimit(1.),
  70   fTPCNClus(80),
  71   fAODtrackCutBit(128),
  72   isQA(kFALSE),
  73   fDebug(kFALSE),
  74   fHelperPID(0x0),
  75   fEventCounter(NULL),
  76   fHistoCorrelation(NULL) {
  77   for(Int_t i = 0; i < 14; i++ ) fHistQA[i] = NULL;
  78   DefineOutput(1, TList::Class()); //! Connect Outpput....
  79 }
  80
  81 AliEbyEParticleRatioFluctuationTask::~AliEbyEParticleRatioFluctuationTask() {
  82   //!   Cleaning up
  83   if (fThnList)   delete fThnList;
  84   if (fHelperPID) delete fHelperPID;
  85 }
  86
  87 //---------------------------------------------------------------------------------
  88 void AliEbyEParticleRatioFluctuationTask::UserCreateOutputObjects() {
  89   fThnList = new TList();
  90   fThnList->SetOwner(kTRUE);
  91
  92   fEventCounter = new TH1D("fEventCounter","EventCounter", 300, 0.5,300.5);
  93   if (isQA) fThnList->Add(fEventCounter);
  94
  95   fHistQA[0] = new TH2F("fHistQAvx", "Histo Vx Selected;Centrality;Vx", 100,0,100, 5000, -5., 5.);
  96   fHistQA[1] = new TH2F("fHistQAvy", "Histo Vy Selected;Centrality;Vy", 100,0,100, 5000, -5., 5.);
  97   fHistQA[2] = new TH2F("fHistQAvz", "Histo Vz Selected;Centrality;Vz", 100,0,100, 5000, -25., 25.);
  98   fHistQA[3] = new TH2F("fHistQAvxA", "Histo Vx;Centrality;Vx", 100,0,100, 5000, -5., 5.);
  99   fHistQA[4] = new TH2F("fHistQAvyA", "Histo Vy;Centrality;Vy", 100,0,100, 5000, -5., 5.);
 100   fHistQA[5] = new TH2F("fHistQAvzA", "Histo Vz;Centrality;Vz", 100,0,100, 5000, -25., 25.);
 101
 102   fHistQA[6] = new TH2F("fHistQADcaXyA", "Histo DCAxy;Centrality;DCAxy",100,0,100, 600, -15., 15.);
 103   fHistQA[7] = new TH2F("fHistQADcaZA", "Histo DCAz;Centrality;DCAz ",100,0,100, 600, -15., 15.);
 104   fHistQA[8] = new TH2F("fHistQAPtA","p_{T} distribution;Centrality;p_{T}",100,0,100,1000,0,10);
 105   fHistQA[9] = new TH2F("fHistQAEtaA","#eta distribution;Centrality;#eta",100,0,100,240,-1.2,1.2);
 106
 107   fHistQA[10] = new TH2F("fHistQADcaXy", "Histo DCAxy after Selected;Centrality;DCAxy", 100,0,100,600, -15., 15.);
 108   fHistQA[11] = new TH2F("fHistQADcaZ", "Histo DCAz Selected;Centrality;DCAz", 100,0,100,600, -15., 15.);
 109   fHistQA[12] = new TH2F("fHistQAPt","p_{T} distribution Selected;Centrality;p_{T}",100,0,100,1000,0,10);
 110   fHistQA[13] = new TH2F("fHistQAEta","#eta distribution Selected;Centrality;#eta",100,0,100, 240,-1.2,1.2);
 111
 112   if (isQA) for(Int_t i = 0; i < 14; i++) fThnList->Add(fHistQA[i]);
 113
 114   Int_t fgSparseDataBins[kNSparseData]   = {100, 5000, 5000, 2500, 2500, 3000, 1500, 1500, 1000, 500, 500, 500, 250, 250};
 115   Double_t fgSparseDataMin[kNSparseData] = {0.,  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,  0.,  0.,  0.,  0.};
 116   Double_t fgSparseDataMax[kNSparseData] = {100.,5000.,5000.,2500.,2500.,3000.,1500.,1500.,1000.,500.,500.,500.,250.,250.};
 117
 118   const Char_t *fgkSparseDataTitle[] = {"centrality","RefMult","N_{ch}", "N_{+}","N_{-}","N_{#pi}", "N_{#pi^{+}}","N_{#pi^{-}}","N_{K}","N_{K^{+}}", "N_{K^{-}}","N_{p}","N_{p}","N_{#bar{p}}"};
 119
 120   fHistoCorrelation = new THnSparseI("fThnCorr", "", kNSparseData, fgSparseDataBins, fgSparseDataMin, fgSparseDataMax);
 121   for (Int_t iaxis = 0; iaxis < kNSparseData; iaxis++)
 122     fHistoCorrelation->GetAxis(iaxis)->SetTitle(fgkSparseDataTitle[iaxis]);
 123
 124   if(!isQA) fThnList->Add(fHistoCorrelation);
 125   if(isQA)
 126     if (fHelperPID)
 127       fThnList->Add(fHelperPID);
 128
 129   PostData(1, fThnList);
 130 }
 131
 132 //----------------------------------------------------------------------------------
 133 void AliEbyEParticleRatioFluctuationTask::UserExec( Option_t * ){
 134
 135   if (isQA) fEventCounter->Fill(1);
 136
 137   AliAODEvent* event = dynamic_cast<AliAODEvent*>(InputEvent());
 138   if (!event) {
 139     Printf("ERROR 01: AOD not found ");
 140     return;
 141   }
 142
 143   Int_t gCent   = -1;
 144   Float_t gRefMul = -1;
 145
 146   AliAODHeader *aodHeader = event->GetHeader();
 147   gCent = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 148   gRefMul = aodHeader->GetRefMultiplicity();
 149   if (gCent < 0 || gCent > 100) return;
 150   if (isQA) fEventCounter->Fill(2);
 151
 152   if (!AcceptEvent(event,gCent)) return;
 153
 154   Int_t icharge = -1;
 155   Int_t gCharge[2];
 156   Int_t gPid[3][2];
 157
 158   for (Int_t i = 0; i < 2; i++) {
 159     gCharge[i] = 0;
 160     for (Int_t ii = 0; ii < 3; ii++) {
 161       gPid[ii][i] = 0;
 162     }
 163   }
 164
 165
 166   if(fAnalysisType == "AOD") {
 167     if (isQA) {
 168         fEventCounter->Fill(5);
 169         fEventCounter->Fill(50+gCent);
 170     }
 171
 172     for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
 173       AliAODTrack* track = dynamic_cast<AliAODTrack *>(event->GetTrack(itrk));
 174       if (!track) continue;
 175
 176       if (!AcceptTrack(track, gCent)) continue;
 177
 178       Int_t a = fHelperPID->GetParticleSpecies((AliVTrack*) track,kTRUE);
 179
 180       if(a < 0 || a > 2) continue;
 181       icharge = track->Charge() > 0 ? 0 : 1;
 182       gCharge[icharge]++;
 183       gPid[a][icharge]++;
 184     }
 185   }
 186   else if(fAnalysisType == "MCAOD") {
 187     TClonesArray *arrayMC= 0;
 188     arrayMC = dynamic_cast<TClonesArray*> (event->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 189     if (!arrayMC) {
 190       Printf("Error: MC particles branch not found!\n");
 191       return;
 192     }
 193     AliAODMCHeader *mcHdr=0;
 194     mcHdr=(AliAODMCHeader*)event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 195     if(!mcHdr) {
 196       Printf("MC header branch not found!\n");
 197       return;
 198     }
 199
 200     if (isQA) {
 201       fEventCounter->Fill(5);
 202       fEventCounter->Fill(50+gCent);
 203     }
 204
 205     Int_t nMC = arrayMC->GetEntries();
 206     for (Int_t iMC = 0; iMC < nMC; iMC++) {
 207       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
 208       if(!AcceptMCTrack(partMC, gCent)) continue;
 209       Int_t a = fHelperPID->GetMCParticleSpecie((AliVEvent*) event,(AliVTrack*)partMC,1);
 210       if(a < 0 || a > 2) continue;
 211       icharge = partMC->Charge() > 0 ? 0 : 1;
 212       gCharge[icharge]++;
 213       gPid[a][icharge]++;
 214     }
 215   }
 216   else {
 217     printf(" No Event Type is Defined ");
 218     return;
 219   }
 220
 221   if( (gCharge[0] + gCharge[1]) != 0 ) {
 222     if (isQA) {
 223       fEventCounter->Fill(6);
 224       fEventCounter->Fill(160 + gCent);
 225     }
 226     else {
 227       Double_t vsparse[kNSparseData];
 228       vsparse[0]   = gCent;
 229       vsparse[1]   = gRefMul;
 230       vsparse[2]   = gCharge[0] + gCharge[1];
 231       vsparse[3]   = gCharge[0];
 232       vsparse[4]   = gCharge[1];
 233       vsparse[5]   = gPid[0][0] + gPid[0][1];
 234       vsparse[6]   = gPid[0][0];
 235       vsparse[7]   = gPid[0][1];
 236       vsparse[8]   = gPid[1][0] + gPid[1][0];
 237       vsparse[9]   = gPid[1][0];
 238       vsparse[10]  = gPid[1][1];
 239       vsparse[11]  = gPid[2][0] + gPid[2][1];
 240       vsparse[12]  = gPid[2][0];
 241       vsparse[13]  = gPid[2][1];
 242       fHistoCorrelation->Fill(vsparse);
 243     }
 244   }
 245
 246   if(fDebug && isQA)  Printf(" %6d  %6d %6d  %6d %6d  %6d %6d %6d %6d %6d %6d %6d %6d",
 247                              (Int_t)fEventCounter->GetBinContent(1),
 248                              (Int_t)fEventCounter->GetBinContent(2),
 249                              (Int_t)fEventCounter->GetBinContent(3),
 250                              (Int_t)gCent, (Int_t)gRefMul,
 251                              gCharge[0], gCharge[1],
 252                              gPid[0][0], gPid[0][1],  gPid[1][0],
 253                              gPid[1][1],  gPid[2][0], gPid[2][1]);
 254
 255   PostData(1, fThnList);
 256 }
 257
 258 void AliEbyEParticleRatioFluctuationTask::Terminate( Option_t * ){
 259   Info("AliEbyEParticleRatioFluctuationTask"," Task Successfully finished");
 260 }
 261
 262 //___________________________________________________________
 263 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptEvent(AliAODEvent *event, Int_t cent) const {
 264   Bool_t ver = kFALSE;
 265   const AliAODVertex *vertex = event->GetPrimaryVertex();
 266   if(vertex) {
 267     Double32_t fCov[6];
 268     vertex->GetCovarianceMatrix(fCov);
 269     if(vertex->GetNContributors() > 0) {
 270       if(fCov[5] != 0) {
 271
 272         if(isQA) {
 273           fEventCounter->Fill(3);
 274           fHistQA[3]->Fill(cent,vertex->GetX());
 275           fHistQA[4]->Fill(cent,vertex->GetY());
 276           fHistQA[5]->Fill(cent,vertex->GetZ());
 277         }
 278
 279         if(TMath::Abs(vertex->GetX()) < fVxMax) {
 280           if(TMath::Abs(vertex->GetY()) < fVyMax) {
 281             if(TMath::Abs(vertex->GetZ()) < fVzMax) {
 282               ver = kTRUE;
 283               if(isQA) {
 284                 fEventCounter->Fill(4);
 285                 fHistQA[0]->Fill(cent,vertex->GetX());
 286                 fHistQA[1]->Fill(cent,vertex->GetY());
 287                 fHistQA[2]->Fill(cent,vertex->GetZ());
 288               }
 289             }
 290           }
 291         }
 292       }
 293     }
 294   }
 295
 296   AliCentrality *centrality = event->GetCentrality();
 297   if (centrality->GetQuality() != 0) ver = kFALSE;
 298   return ver;
 299 }
 300
 301
 302 //___________________________________________________________
 303 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptTrack(AliAODTrack *track, Int_t cent) const {
 304   if(!track)                                  return kFALSE;
 305   if (track->Charge() == 0 )                  return kFALSE;
 306
 307   if(isQA) {
 308     fHistQA[6]->Fill(cent,track->DCA());
 309     fHistQA[7]->Fill(cent,track->ZAtDCA());
 310     fHistQA[8]->Fill(cent,track->Pt());
 311     fHistQA[9]->Fill(cent,track->Eta());
 312   }
 313
 314   if (!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;
 315
 316   if (track->Eta() < fEtaLowerLimit ||
 317       track->Eta() > fEtaHigherLimit)         return kFALSE;
 318   if (track->Pt() < fPtLowerLimit ||
 319       track->Pt() > fPtHigherLimit)           return kFALSE;
 320   if ( track->DCA() > fDCAxy )                return kFALSE;
 321   if ( track->ZAtDCA() > fDCAz )              return kFALSE;
 322
 323   if(isQA) {
 324     fHistQA[10]->Fill(cent,track->DCA());
 325     fHistQA[11]->Fill(cent,track->ZAtDCA());
 326     fHistQA[12]->Fill(cent,track->Pt());
 327     fHistQA[13]->Fill(cent,track->Eta());
 328   }
 329
 330   return kTRUE;
 331 }
 332
 333
 334 //___________________________________________________________
 335 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptMCTrack(AliAODMCParticle *track, Int_t cent) const {
 336   if(!track) return kFALSE;
 337   if(!track->IsPhysicalPrimary()) return kFALSE;
 338   if (track->Charge() == 0 ) return kFALSE;
 339   if(isQA) {
 340     fHistQA[8]->Fill(cent,track->Pt());
 341     fHistQA[9]->Fill(cent,track->Eta());
 342   }
 343
 344   if (track->Eta() < fEtaLowerLimit ||
 345       track->Eta() > fEtaHigherLimit) return kFALSE;
 346   if (track->Pt() < fPtLowerLimit ||
 347       track->Pt() > fPtHigherLimit) return kFALSE;
 348
 349   if(isQA) {
 350     fHistQA[12]->Fill(cent,track->Pt());
 351     fHistQA[13]->Fill(cent,track->Eta());
 352   }
 353
 354   return kTRUE;
 355 }