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 ), \r
  56   fThnList(NULL), \r
  57   fAnalysisType("AOD"), \r
  58   fAnalysisData("PbPb"), \r
  59   fCentralityEstimator("V0M"), \r
  60   fVxMax(3.), \r
  61   fVyMax(3.), \r
  62   fVzMax(10.), \r
  63   fDCAxy(2.4),  \r
  64   fDCAz(3.2), \r
  65   fPtLowerLimit(0.2), \r
  66   fPtHigherLimit(5.), \r
  67   fEtaLowerLimit(-1.), \r
  68   fEtaHigherLimit(1.), \r
  69   fTPCNClus(80),\r
  70   fAODtrackCutBit(128), \r
  71   isQA(kFALSE), \r
  72   fDebug(kFALSE), \r
  73   fHelperPID(0x0),\r
  74   fEventCounter(NULL), \r
  75   fHistoCorrelation(NULL) { \r
  76   for(Int_t i = 0; i < 14; i++ ) fHistQA[i] = NULL;\r
  77   DefineOutput(1, TList::Class()); //! Connect Outpput....\r
  78 }\r
  79 \r
  80 AliEbyEParticleRatioFluctuationTask::~AliEbyEParticleRatioFluctuationTask() {\r
  81   //!   Cleaning up\r
  82   if (fThnList)   delete fThnList;\r
  83   if (fHelperPID) delete fHelperPID;\r
  84 }\r
  85 \r
  86 //---------------------------------------------------------------------------------\r
  87 void AliEbyEParticleRatioFluctuationTask::UserCreateOutputObjects() {\r
  88   fThnList = new TList();\r
  89   fThnList->SetOwner(kTRUE);\r
  90 \r
  91   fEventCounter = new TH1D("fEventCounter","EventCounter", 300, 0.5,300.5);\r
  92   if (isQA) fThnList->Add(fEventCounter);\r
  93   \r
  94   fHistQA[0] = new TH2F("fHistQAvx", "Histo Vx Selected;Centrality;Vx", 100,0,100, 5000, -5., 5.);\r
  95   fHistQA[1] = new TH2F("fHistQAvy", "Histo Vy Selected;Centrality;Vy", 100,0,100, 5000, -5., 5.);\r
  96   fHistQA[2] = new TH2F("fHistQAvz", "Histo Vz Selected;Centrality;Vz", 100,0,100, 5000, -25., 25.);  \r
  97   fHistQA[3] = new TH2F("fHistQAvxA", "Histo Vx;Centrality;Vx", 100,0,100, 5000, -5., 5.);\r
  98   fHistQA[4] = new TH2F("fHistQAvyA", "Histo Vy;Centrality;Vy", 100,0,100, 5000, -5., 5.);\r
  99   fHistQA[5] = new TH2F("fHistQAvzA", "Histo Vz;Centrality;Vz", 100,0,100, 5000, -25., 25.);\r
 100 \r
 101   fHistQA[6] = new TH2F("fHistQADcaXyA", "Histo DCAxy;Centrality;DCAxy",100,0,100, 600, -15., 15.);\r
 102   fHistQA[7] = new TH2F("fHistQADcaZA", "Histo DCAz;Centrality;DCAz ",100,0,100, 600, -15., 15.);   \r
 103   fHistQA[8] = new TH2F("fHistQAPtA","p_{T} distribution;Centrality;p_{T}",100,0,100,1000,0,10);\r
 104   fHistQA[9] = new TH2F("fHistQAEtaA","#eta distribution;Centrality;#eta",100,0,100,240,-1.2,1.2);\r
 105 \r
 106   fHistQA[10] = new TH2F("fHistQADcaXy", "Histo DCAxy after Selected;Centrality;DCAxy", 100,0,100,600, -15., 15.);\r
 107   fHistQA[11] = new TH2F("fHistQADcaZ", "Histo DCAz Selected;Centrality;DCAz", 100,0,100,600, -15., 15.);   \r
 108   fHistQA[12] = new TH2F("fHistQAPt","p_{T} distribution Selected;Centrality;p_{T}",100,0,100,1000,0,10);\r
 109   fHistQA[13] = new TH2F("fHistQAEta","#eta distribution Selected;Centrality;#eta",100,0,100, 240,-1.2,1.2);\r
 110   \r
 111   if (isQA) for(Int_t i = 0; i < 14; i++) fThnList->Add(fHistQA[i]);\r
 112   \r
 113   Int_t fgSparseDataBins[kNSparseData]   = {100, 5000, 5000, 2500, 2500, 3000, 1500, 1500, 1000, 500, 500, 500, 250, 250};\r
 114   Double_t fgSparseDataMin[kNSparseData] = {0.,  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,  0.,  0.,  0.,  0.};\r
 115   Double_t fgSparseDataMax[kNSparseData] = {100.,5000.,5000.,2500.,2500.,3000.,1500.,1500.,1000.,500.,500.,500.,250.,250.};\r
 116   \r
 117   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
 118     \r
 119   fHistoCorrelation = new THnSparseI("fThnCorr", "", kNSparseData, fgSparseDataBins, fgSparseDataMin, fgSparseDataMax);\r
 120   for (Int_t iaxis = 0; iaxis < kNSparseData; iaxis++)\r
 121     fHistoCorrelation->GetAxis(iaxis)->SetTitle(fgkSparseDataTitle[iaxis]);\r
 122   \r
 123   if(!isQA) fThnList->Add(fHistoCorrelation);\r
 124   if(isQA)  \r
 125     if (fHelperPID)\r
 126       fThnList->Add(fHelperPID);\r
 127   \r
 128   PostData(1, fThnList);\r
 129 }\r
 130 \r
 131 //----------------------------------------------------------------------------------\r
 132 void AliEbyEParticleRatioFluctuationTask::UserExec( Option_t * ){\r
 133 \r
 134   if (isQA) fEventCounter->Fill(1);\r
 135 \r
 136   AliAODEvent* event = dynamic_cast<AliAODEvent*>(InputEvent());\r
 137   if (!event) {\r
 138     cout<< "ERROR 01: AOD not found " <<endl;\r
 139     return;\r
 140   }\r
 141 \r
 142   Int_t gCent   = -1;\r
 143   Float_t gRefMul = -1;\r
 144   \r
 145   AliAODHeader *aodHeader = event->GetHeader();\r
 146   gCent = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 147   gRefMul = aodHeader->GetRefMultiplicity();\r
 148   if (gCent < 0 || gCent > 100) return;\r
 149   if (isQA) fEventCounter->Fill(2);  \r
 150 \r
 151   if (!AcceptEvent(event,gCent)) return;\r
 152   \r
 153   Int_t icharge = -1;\r
 154   Int_t gCharge[2];\r
 155   Int_t gPid[3][2];\r
 156   \r
 157   for (Int_t i = 0; i < 2; i++) {\r
 158     gCharge[i] = 0;\r
 159     for (Int_t ii = 0; ii < 3; ii++) {\r
 160       gPid[ii][i] = 0;\r
 161     }\r
 162   }\r
 163 \r
 164 \r
 165   if(fAnalysisType == "AOD") {\r
 166     if (isQA) {\r
 167         fEventCounter->Fill(5);\r
 168         fEventCounter->Fill(50+gCent);\r
 169     }\r
 170     \r
 171     for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {\r
 172       AliAODTrack* track = dynamic_cast<AliAODTrack *>(event->GetTrack(itrk));\r
 173       if (!track) continue;\r
 174       \r
 175       if (!AcceptTrack(track, gCent)) continue;\r
 176             \r
 177       Int_t a = fHelperPID->GetParticleSpecies((AliVTrack*) track,kTRUE);\r
 178 \r
 179       if(a < 0 || a > 2) continue;\r
 180       icharge = track->Charge() > 0 ? 0 : 1;\r
 181       gCharge[icharge]++;\r
 182       gPid[a][icharge]++;\r
 183     }\r
 184   }\r
 185   else if(fAnalysisType == "MCAOD") {\r
 186     TClonesArray *arrayMC= 0; \r
 187     arrayMC = dynamic_cast<TClonesArray*> (event->GetList()->FindObject(AliAODMCParticle::StdBranchName()));\r
 188     if (!arrayMC) {\r
 189       Printf("Error: MC particles branch not found!\n");\r
 190       return;\r
 191     }\r
 192     AliAODMCHeader *mcHdr=0;\r
 193     mcHdr=(AliAODMCHeader*)event->GetList()->FindObject(AliAODMCHeader::StdBranchName());  \r
 194     if(!mcHdr) {\r
 195       Printf("MC header branch not found!\n");\r
 196       return;\r
 197     }\r
 198     \r
 199     if (isQA) {\r
 200       fEventCounter->Fill(5);\r
 201       fEventCounter->Fill(50+gCent);\r
 202     }\r
 203     \r
 204     Int_t nMC = arrayMC->GetEntries();\r
 205     for (Int_t iMC = 0; iMC < nMC; iMC++) {\r
 206       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
 207       if(!AcceptMCTrack(partMC, gCent)) continue;\r
 208       Int_t a = fHelperPID->GetMCParticleSpecie((AliVEvent*) event,(AliVTrack*)partMC,1);\r
 209       if(a < 0 || a > 2) continue;\r
 210       icharge = partMC->Charge() > 0 ? 0 : 1;\r
 211       gCharge[icharge]++;\r
 212       gPid[a][icharge]++;\r
 213     }\r
 214   }\r
 215   else {\r
 216     printf(" No Event Type is Defined ");\r
 217     return;\r
 218   }\r
 219   \r
 220   if( (gCharge[0] + gCharge[1]) != 0 ) {\r
 221     if (isQA) {\r
 222       fEventCounter->Fill(6); \r
 223       fEventCounter->Fill(160 + gCent);\r
 224     }\r
 225     else { \r
 226       Double_t vsparse[kNSparseData];\r
 227       vsparse[0]   = gCent;\r
 228       vsparse[1]   = gRefMul;\r
 229       vsparse[2]   = gCharge[0] + gCharge[1];\r
 230       vsparse[3]   = gCharge[0];\r
 231       vsparse[4]   = gCharge[1];\r
 232       vsparse[5]   = gPid[0][0] + gPid[0][1];\r
 233       vsparse[6]   = gPid[0][0];\r
 234       vsparse[7]   = gPid[0][1];\r
 235       vsparse[8]   = gPid[1][0] + gPid[1][0];\r
 236       vsparse[9]   = gPid[1][0];\r
 237       vsparse[10]  = gPid[1][1];\r
 238       vsparse[11]  = gPid[2][0] + gPid[2][1];\r
 239       vsparse[12]  = gPid[2][0];\r
 240       vsparse[13]  = gPid[2][1];\r
 241       fHistoCorrelation->Fill(vsparse);\r
 242     }\r
 243   }\r
 244   \r
 245   if(fDebug && isQA)  Printf(" %6d  %6d %6d  %6d %6d  %6d %6d %6d %6d %6d %6d %6d %6d", \r
 246                              (Int_t)fEventCounter->GetBinContent(1), \r
 247                              (Int_t)fEventCounter->GetBinContent(2), \r
 248                              (Int_t)fEventCounter->GetBinContent(3), \r
 249                              (Int_t)gCent, (Int_t)gRefMul, \r
 250                              gCharge[0], gCharge[1], \r
 251                              gPid[0][0], gPid[0][1],  gPid[1][0], \r
 252                              gPid[1][1],  gPid[2][0], gPid[2][1]);\r
 253   \r
 254   PostData(1, fThnList);\r
 255 }\r
 256 \r
 257 void AliEbyEParticleRatioFluctuationTask::Terminate( Option_t * ){\r
 258   Info("AliEbyEParticleRatioFluctuationTask"," Task Successfully finished");\r
 259 }\r
 260 \r
 261 //___________________________________________________________\r
 262 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptEvent(AliAODEvent *event, Int_t cent) const {\r
 263   Bool_t ver = kFALSE;\r
 264   const AliAODVertex *vertex = event->GetPrimaryVertex();\r
 265   if(vertex) {\r
 266     Double32_t fCov[6];\r
 267     vertex->GetCovarianceMatrix(fCov);\r
 268     if(vertex->GetNContributors() > 0) {\r
 269       if(fCov[5] != 0) {\r
 270         \r
 271         if(isQA) {      \r
 272           fEventCounter->Fill(3);\r
 273           fHistQA[3]->Fill(cent,vertex->GetX());\r
 274           fHistQA[4]->Fill(cent,vertex->GetY());\r
 275           fHistQA[5]->Fill(cent,vertex->GetZ());\r
 276         }\r
 277         \r
 278         if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
 279           if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
 280             if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
 281               ver = kTRUE;\r
 282               if(isQA) {        \r
 283                 fEventCounter->Fill(4);\r
 284                 fHistQA[0]->Fill(cent,vertex->GetX());\r
 285                 fHistQA[1]->Fill(cent,vertex->GetY());\r
 286                 fHistQA[2]->Fill(cent,vertex->GetZ());\r
 287               }\r
 288             }\r
 289           }\r
 290         }\r
 291       }\r
 292     }\r
 293   }\r
 294   \r
 295   AliCentrality *centrality = event->GetCentrality();\r
 296   if (centrality->GetQuality() != 0) ver = kFALSE;\r
 297   return ver;\r
 298 }\r
 299 \r
 300 \r
 301 //___________________________________________________________\r
 302 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptTrack(AliAODTrack *track, Int_t cent) const {\r
 303   if(!track)                                  return kFALSE;\r
 304   if (track->Charge() == 0 )                  return kFALSE;\r
 305   \r
 306   if(isQA) {\r
 307     fHistQA[6]->Fill(cent,track->DCA());\r
 308     fHistQA[7]->Fill(cent,track->ZAtDCA());\r
 309     fHistQA[8]->Fill(cent,track->Pt());\r
 310     fHistQA[9]->Fill(cent,track->Eta());\r
 311   }\r
 312   \r
 313   if (!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;\r
 314 \r
 315   if (track->Eta() < fEtaLowerLimit ||\r
 316       track->Eta() > fEtaHigherLimit)         return kFALSE;\r
 317   if (track->Pt() < fPtLowerLimit ||\r
 318       track->Pt() > fPtHigherLimit)           return kFALSE;  \r
 319   if ( track->DCA() > fDCAxy )                return kFALSE; \r
 320   if ( track->ZAtDCA() > fDCAz )              return kFALSE;\r
 321    \r
 322   if(isQA) {\r
 323     fHistQA[10]->Fill(cent,track->DCA());\r
 324     fHistQA[11]->Fill(cent,track->ZAtDCA());\r
 325     fHistQA[12]->Fill(cent,track->Pt());\r
 326     fHistQA[13]->Fill(cent,track->Eta());\r
 327   }\r
 328  \r
 329   return kTRUE;\r
 330 }\r
 331 \r
 332 \r
 333 //___________________________________________________________\r
 334 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptMCTrack(AliAODMCParticle *track, Int_t cent) const {\r
 335   if(!track) return kFALSE;\r
 336   if(!track->IsPhysicalPrimary()) return kFALSE;\r
 337   if (track->Charge() == 0 ) return kFALSE;\r
 338   if(isQA) {\r
 339     fHistQA[8]->Fill(cent,track->Pt());\r
 340     fHistQA[9]->Fill(cent,track->Eta());\r
 341   }\r
 342 \r
 343   if (track->Eta() < fEtaLowerLimit ||\r
 344       track->Eta() > fEtaHigherLimit) return kFALSE;\r
 345   if (track->Pt() < fPtLowerLimit ||\r
 346       track->Pt() > fPtHigherLimit) return kFALSE;  \r
 347     \r
 348   if(isQA) {\r
 349     fHistQA[12]->Fill(cent,track->Pt());\r
 350     fHistQA[13]->Fill(cent,track->Eta());\r
 351   }\r
 352  \r
 353   return kTRUE;\r
 354 }\r