PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidRatioTaskOnFly.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 Fluctuaion Analysis for PID                     //
  19 //                   Deepika Jena | drathee@cern.ch                        //
  20 //                   Satyajit Jena | sjena@cern.ch                         //
  21 //=========================================================================//
  22
  23 #include "TChain.h"
  24 #include "TList.h"
  25 #include "TFile.h"
  26 #include "TTree.h"
  27 #include "TH1D.h"
  28 #include "TH2F.h"
  29 #include "TH3F.h"
  30 #include "TProfile.h"
  31
  32 #include "AliAnalysisTask.h"
  33 #include "AliAnalysisManager.h"
  34
  35 #include "AliVEvent.h"
  36 #include "AliESDEvent.h"
  37 #include "AliMCEvent.h"
  38 #include "AliAODEvent.h"
  39
  40 #include "AliStack.h"
  41 #include "AliGenEventHeader.h"
  42 #include "AliGenPythiaEventHeader.h"
  43 #include "AliGenHijingEventHeader.h"
  44 #include "AliGenDPMjetEventHeader.h"
  45 #include "TDatabasePDG.h"
  46
  47 #include "AliEbyEPidRatioTaskOnFly.h"
  48
  49 ClassImp(AliEbyEPidRatioTaskOnFly)
  50
  51 //-----------------------------------------------------------------------
  52 AliEbyEPidRatioTaskOnFly::AliEbyEPidRatioTaskOnFly( const char *name )
  53 : AliAnalysisTaskSE( name ),
  54   fThnList(0),
  55   fPtLowerLimit(0.2),
  56   fPtHigherLimit(5.),
  57   fEtaLowerLimit(-1.),
  58   fEtaHigherLimit(1.),
  59   fCentrality(-1),
  60   fOrder(8),
  61   fRedFactp(NULL)
  62
  63
  64 {
  65
  66   DefineOutput(1, TList::Class());
  67 }
  68
  69 AliEbyEPidRatioTaskOnFly::~AliEbyEPidRatioTaskOnFly()
  70 {
  71   if (fThnList) delete fThnList;
  72
  73   for (Int_t ii = 0; ii <= fOrder; ++ii)
  74     if (fRedFactp[ii]) delete[] fRedFactp[ii];
  75   if (fRedFactp) delete[] fRedFactp;
  76
  77 }
  78
  79
  80
  81
  82 //---------------------------------------------------------------------------------
  83 void AliEbyEPidRatioTaskOnFly::UserCreateOutputObjects() {
  84   fThnList = new TList();
  85   fThnList->SetOwner(kTRUE);
  86
  87   const Char_t *name = "Mc";
  88   const Char_t *title = Form(" #eta [%2.1f-%2.1f]",fEtaLowerLimit,fEtaHigherLimit);
  89
  90
  91   TString sName(name);
  92   TString sTitle(title);
  93
  94
  95   fRedFactp = new Double_t*[fOrder+1];
  96   for (Int_t ii = 0 ; ii <= fOrder; ++ii)
  97     fRedFactp[ii] = new Double_t[2];
  98
  99   //TList *list[4];
 100   fThnList->Add(new TList);
 101   TList *list =  static_cast<TList*>(fThnList->Last());
 102   list->SetName(Form("f%s",name));
 103   list->SetOwner(kTRUE);
 104
 105   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 106     TString sNetTitle(Form("%s - %s",fgkPidLatex[iPid][1],fgkPidLatex[iPid][0]));
 107
 108     list->Add(new TProfile(Form("fProfTot%sPlus%s",fgkPidName[iPid],name),
 109                            Form("(%s) : %s;Centrality(100);(%s)",fgkPidName[iPid], sTitle.Data(), sNetTitle.Data()),
 110                            20,-0.5,19.5));
 111
 112     list->Add(new TProfile(Form("fProfTot%sMinus%s",fgkPidName[iPid],name),
 113                            Form("(%s) : %s;Centrality(100);(%s)",fgkPidName[iPid], sTitle.Data(), sNetTitle.Data()),
 114                            20,-0.5,19.5));
 115
 116     for (Int_t idx = 1; idx <= fOrder; ++idx) {
 117       list->Add(new TProfile(Form("fProf%s%sNet%dM",fgkPidName[iPid],name, idx),
 118                              Form("(%s)^{%d} : %s;Centrality(100);(%s)^{%d}",sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
 119                              20,-0.5,19.5));
 120     }
 121
 122     for (Int_t ii = 0; ii <= fOrder; ++ii) {
 123       for (Int_t kk = 0; kk <= fOrder; ++kk) {
 124         list->Add(new TProfile(Form("fProf%s%sNetF%02d%02d",fgkPidName[iPid], name, ii, kk),
 125                                Form("f_{%02d%02d} : %s;Centrality(100);f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
 126                                20,-0.5,19.5));
 127       }
 128     }
 129
 130   }
 131
 132
 133   for (Int_t iPhy = 0; iPhy < 46; ++iPhy) {
 134     list->Add(new TProfile(Form("fProf%sNu%02d",name,iPhy),Form("Physics Variable for index %d | %s ; Centrality;",iPhy,name),20,-0.5,19.5));
 135   }
 136
 137
 138   PostData(1, fThnList);
 139 }
 140
 141
 142 //----------------------------------------------------------------------------------
 143 void AliEbyEPidRatioTaskOnFly::UserExec( Option_t * ){
 144   AliMCEvent* mcEvent = MCEvent();
 145   if (!mcEvent) {
 146     Printf("ERROR: Could not retrieve MC event");
 147     return;
 148   }
 149
 150   for (Int_t ii = 0 ; ii < 4; ++ii)
 151     for (Int_t kk = 0 ; kk < 2; ++kk)
 152     fNp[ii][kk] = 0;
 153
 154
 155   const AliVVertex *vtxMC = mcEvent->GetPrimaryVertex();
 156   if (vtxMC->GetZ() > 10.) return;
 157
 158   AliStack *stack = mcEvent->Stack();
 159
 160   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 161   if(!genHeader){
 162     printf("  Event generator header not available!!!\n");
 163     return;
 164   }
 165
 166   Double_t imp = 21;
 167
 168   if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
 169     imp = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
 170   }
 171   if (imp < 20) fCentrality = Int_t(imp);
 172   else return;
 173
 174   Int_t nParticles = stack->GetNtrack();
 175   for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
 176     TParticle* part = stack->Particle(iParticle);
 177     if (!part) {
 178       Printf(" No Particle Available ");
 179       continue;
 180     }
 181
 182     Float_t pt = part->Pt();
 183     if(pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
 184     Float_t gEta = part->Eta();
 185     if(gEta < fEtaLowerLimit || gEta > fEtaHigherLimit ) continue;
 186
 187     //  Float_t gCharge = (pdgPart ? pdgPart->Charge() : 0);
 188
 189     Int_t pid = part->GetPdgCode();
 190     if(pid == 211)         { fNp[1][1]++;  fNp[0][1]++; }
 191     else if(pid ==  -211)  { fNp[1][0]++;  fNp[0][0]++; }
 192     else if(pid ==   321)  { fNp[2][1]++;  fNp[0][1]++; }
 193     else if(pid ==  -321)  { fNp[2][0]++;  fNp[0][0]++; }
 194     else if(pid ==  2212)  { fNp[3][1]++;  fNp[0][1]++; }
 195     else if(pid == -2212)  { fNp[3][0]++;  fNp[0][0]++; }
 196
 197   }
 198
 199   /*Printf("%6d %6d %6d %6d %6d %6d %6d %6d %6d", fCentrality,
 200          fNp[0][0], fNp[0][1],
 201          fNp[1][0], fNp[1][1],
 202          fNp[2][0], fNp[2][1],
 203          fNp[3][0], fNp[3][1]);*/
 204
 205
 206   FillHistSetCent();
 207   PostData(1, fThnList);
 208 }
 209
 210 void AliEbyEPidRatioTaskOnFly::Terminate( Option_t * ){
 211   Info("AliEbyEPidRatioTaskOnFly"," Task Successfully finished");
 212 }
 213 //________________________________________________________________________
 214 const Char_t* AliEbyEPidRatioTaskOnFly::fgkPidName[4]      = {"Nch","Npi","Nka","Npr"};
 215 //________________________________________________________________________
 216 const Char_t* AliEbyEPidRatioTaskOnFly::fgkPidLatex[4][2]  = {{"N_{-}","N_{+}"}, {"N_{#pi^{-}}","N_{#pi^{+}}"},{"N_{K^{-}}","N_{K^{+}}"}, {"N_{#bar{p}}","N_{p}"}};
 217 //________________________________________________________________________
 218 const Char_t* AliEbyEPidRatioTaskOnFly::fgkPidTitles[4][2] = {{"Negative","Positive"},{"Anti-Pions","Pions"},{"Anti-Kaons","Kaons"}, {"Anti-Protons","Protons"}};
 219 //________________________________________________________________________
 220 void AliEbyEPidRatioTaskOnFly::FillHistSetCent()  {
 221
 222   const Char_t *name = "Mc";
 223
 224   TList *list = static_cast<TList*>(fThnList->FindObject(Form("f%s",name)));
 225
 226   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 227     Int_t deltaNp = fNp[iPid][1]-fNp[iPid][0];
 228     Double_t delta = 1.;
 229     (static_cast<TProfile*>(list->FindObject(Form("fProfTot%sPlus%s", fgkPidName[iPid], name))))->Fill(fCentrality, fNp[iPid][1]);
 230     (static_cast<TProfile*>(list->FindObject(Form("fProfTot%sMinus%s", fgkPidName[iPid], name))))->Fill(fCentrality, fNp[iPid][0]);
 231     for (Int_t idxOrder = 1; idxOrder <= fOrder; ++idxOrder) {
 232       delta *= deltaNp;
 233       (static_cast<TProfile*>(list->FindObject(Form("fProf%s%sNet%dM", fgkPidName[iPid], name, idxOrder))))->Fill(fCentrality, delta);
 234     }
 235
 236     for (Int_t idxOrder = 0; idxOrder <= fOrder; ++ idxOrder) {
 237       fRedFactp[idxOrder][0]  = 1.;
 238       fRedFactp[idxOrder][1]  = 1.;
 239     }
 240
 241     for (Int_t idxOrder = 1; idxOrder <= fOrder; ++ idxOrder) {
 242       fRedFactp[idxOrder][0]  = fRedFactp[idxOrder-1][0]  * Double_t(fNp[iPid][0]-(idxOrder-1));
 243       fRedFactp[idxOrder][1]  = fRedFactp[idxOrder-1][1]  * Double_t(fNp[iPid][1]-(idxOrder-1));
 244     }
 245
 246     for (Int_t ii = 0; ii <= fOrder; ++ii) {
 247       for (Int_t kk = 0; kk <= fOrder; ++kk) {
 248         Double_t fik = fRedFactp[ii][1] * fRedFactp[kk][0];
 249         (static_cast<TProfile*>(list->FindObject(Form("fProf%s%sNetF%02d%02d", fgkPidName[iPid], name, ii, kk))))->Fill(fCentrality, fik);
 250       }
 251     }
 252   }
 253
 254   //Printf("%6d %20s %6.2f %6d %6d %6d %6d  %6d %6d %6d %6d", idx, name, centralityBin,
 255   //     fNp[0][1],  fNp[0][0],
 256   //     fNp[1][1],  fNp[1][0],
 257   ///    fNp[2][1],  fNp[2][0],
 258   //     fNp[3][1],  fNp[3][0]);
 259   //
 260
 261    Int_t a[6][4]; Int_t b[22];
 262    for (Int_t iPid = 0; iPid < 4; ++iPid) {
 263      a[0][iPid] = fNp[iPid][1]+fNp[iPid][0];       // 0  n+ + n-
 264      a[1][iPid] = fNp[iPid][1];                        // 1  n+
 265      a[2][iPid] = fNp[iPid][0];                        // 2  n-
 266      a[3][iPid] = fNp[iPid][1]*fNp[iPid][0];       // 3  n+ . n-
 267      a[4][iPid] = fNp[iPid][1]*(fNp[iPid][1]-1);   // 4  n+ (n+ - 1)
 268      a[5][iPid] = fNp[iPid][0]*(fNp[iPid][0]-1);   // 5  n- (n- - 1)
 269
 270      // Printf("%6d %20s %6.2f %6d %6d %6d ", idx, name, centralityBin,
 271      //    a[0][iPid], a[1][iPid], a[2][iPid]);
 272
 273   }
 274
 275   b[0]  = a[0][0]*a[0][2];       // 24 N   K
 276   b[1]  = a[0][1]*a[0][2];       // 25 Pi  K
 277   b[2]  = a[1][1]*a[1][2];       // 26 pi+ k+
 278   b[3]  = a[1][1]*a[2][2];       // 27 pi+ k-
 279   b[4]  = a[2][1]*a[1][2];       // 28 pi- k+
 280   b[5]  = a[2][1]*a[2][2];       // 29 pi- k-
 281
 282   b[6]  = a[0][0]*a[0][3];       // 30 N   P
 283   b[7]  = a[0][2]*a[0][3];       // 31 K   P
 284   b[8]  = a[1][2]*a[1][3];       // 32 k+  p+
 285   b[9]  = a[1][2]*a[2][3];       // 33 k+  p-
 286   b[10] = a[2][2]*a[1][3];       // 34 k-  p+
 287   b[11] = a[2][2]*a[2][3];       // 35 k-  p-
 288
 289   b[12] = a[0][0]*a[0][1];       // 36 N  Pi
 290   b[13] = a[0][3]*a[0][1];       // 37 P  Pi
 291   b[14] = a[1][3]*a[1][1];       // 38 p+ pi+
 292   b[15] = a[1][3]*a[2][1];       // 39 p+ pi-
 293   b[16] = a[2][3]*a[1][1];       // 40 p- pi+
 294   b[17] = a[2][3]*a[2][1];       // 41 p- pi-
 295
 296   b[18] = a[0][0]*(a[0][0] - 1); // 42 N ( N - 1 )
 297   b[19] = a[0][1]*(a[0][1] - 1); // 43 Pi( Pi- 1 )
 298   b[20] = a[0][2]*(a[0][1] - 1); // 44 K ( K - 1 )
 299   b[21] = a[0][3]*(a[0][3] - 1); // 45 P ( P - 1 )
 300   // TList *list_nu = static_cast<TList*>(fOutList->FindObject(Form("f%s_nu",name)));
 301   Int_t k = 0;
 302   for (Int_t j = 0; j < 4; j++) {
 303     for (Int_t i = 0; i < 6; i++) {
 304       (static_cast<TProfile*>(list->FindObject(Form("fProf%sNu%02d", name,k))))->Fill(fCentrality,a[i][j]);
 305       k++;
 306     }
 307   }
 308
 309   for (Int_t j = 0; j < 22; j++) {
 310     (static_cast<TProfile*>(list->FindObject(Form("fProf%sNu%02d", name,j+24))))->Fill(fCentrality,b[j]);
 311   }
 312
 313   /*
 314   Printf("%6d %6d %6d %6d %6d %6d %6d %6d %6d", fCentrality,
 315          fNp[0][0], fNp[0][1],
 316          fNp[1][0], fNp[1][1],
 317          fNp[2][0], fNp[2][1],
 318          fNp[3][0], fNp[3][1]);
 319   */
 320
 321   return;
 322 }
 323