PWGCF/EBYE/PIDFluctuation/task/AliEbyEParticleRatioFluctuationTask.cxx

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