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