PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidRatioPhy.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: ALICE Offline.                                                 *
   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 //                  Date: Wed Jul  9 18:38:30 CEST 2014                    //
  22 //          New approch to find particle ratio to reduce memory            //
  23 //                             (Test Only)                                 //
  24 //=========================================================================//
  25
  26 #include "TMath.h"
  27 #include "TAxis.h"
  28 #include "TProfile.h"
  29 #include "TProfile2D.h"
  30 #include "TH2D.h"
  31 #include "TH3D.h"
  32
  33 #include "AliStack.h"
  34 #include "AliMCEvent.h"
  35 #include "AliMCParticle.h"
  36 #include "AliESDEvent.h"
  37 #include "AliESDtrackCuts.h"
  38
  39 #include "AliCentrality.h"
  40 #include "AliTracker.h"
  41
  42 #include "AliAODEvent.h"
  43 #include "AliAODTrack.h"
  44 #include "AliAODMCParticle.h"
  45
  46 #include "AliEbyEPidRatioPhy.h"
  47
  48 using namespace std;
  49
  50 ClassImp(AliEbyEPidRatioPhy)
  51 //________________________________________________________________________
  52 AliEbyEPidRatioPhy::AliEbyEPidRatioPhy() :
  53   AliEbyEPidRatioBase("Dist", "Dist"),
  54   fOutList(NULL),
  55   fOrder(8),
  56   fNNp(6),
  57   fNp(NULL),
  58   fNMCNp(7),
  59   fMCNp(NULL),
  60   fRedFactp(NULL) {
  61   AliLog::SetClassDebugLevel("AliEbyEPidRatioPhy",10);
  62 }
  63
  64 //________________________________________________________________________
  65 AliEbyEPidRatioPhy::~AliEbyEPidRatioPhy() {
  66   // Destructor
  67
  68   for (Int_t ii = 0; ii < fNNp; ++ii) {
  69     for (Int_t kk = 0; kk < 2; ++kk)
  70       if (fNp[ii][kk]) delete[] fNp[ii][kk];
  71     if (fNp[ii]) delete[] fNp[ii];
  72   }
  73   if (fNp) delete[] fNp;
  74
  75  for (Int_t ii = 0; ii < fNMCNp; ++ii) {
  76     for (Int_t kk = 0; kk < 2; ++kk)
  77       if (fMCNp[ii][kk]) delete[] fMCNp[ii][kk];
  78     if (fMCNp[ii]) delete[] fMCNp[ii];
  79   }
  80   if (fMCNp) delete[] fMCNp;
  81
  82   for (Int_t ii = 0; ii <= fOrder; ++ii)
  83     if (fRedFactp[ii]) delete[] fRedFactp[ii];
  84   if (fRedFactp) delete[] fRedFactp;
  85
  86   return;
  87 }
  88
  89 //________________________________________________________________________
  90 void AliEbyEPidRatioPhy::Process() {
  91   ProcessTracks();
  92   if (fIsMC)
  93     ProcessParticles();
  94   return;
  95 }
  96
  97 //________________________________________________________________________
  98 void AliEbyEPidRatioPhy::Init() {
  99   fNp = new Int_t**[fNNp];
 100   for (Int_t ii = 0 ; ii < fNNp; ++ii) {
 101     fNp[ii] = new Int_t*[4];
 102     for (Int_t kk = 0 ; kk < 4; ++kk)
 103       fNp[ii][kk] = new Int_t[2];
 104   }
 105
 106   fMCNp = new Int_t**[fNMCNp];
 107   for (Int_t ii = 0 ; ii < fNMCNp; ++ii) {
 108     fMCNp[ii] = new Int_t*[4];
 109     for (Int_t kk = 0 ; kk < 4; ++kk)
 110       fMCNp[ii][kk] = new Int_t[2];
 111   }
 112
 113   fRedFactp = new Double_t*[fOrder+1];
 114   for (Int_t ii = 0 ; ii <= fOrder; ++ii)
 115     fRedFactp[ii] = new Double_t[2];
 116 }
 117
 118 //________________________________________________________________________
 119 void AliEbyEPidRatioPhy::Reset() {
 120   for (Int_t ii = 0; ii < fNNp; ++ii)
 121     for (Int_t kk = 0; kk < 4; ++kk)
 122       for (Int_t jj = 0; jj < 2; ++jj)
 123         fNp[ii][kk][jj] = 0;
 124   for (Int_t ii = 0; ii < fNMCNp; ++ii)
 125     for (Int_t kk = 0; kk < 4; ++kk)
 126       for (Int_t jj = 0; jj < 2; ++jj)
 127         fMCNp[ii][kk][jj] = 0;
 128   for (Int_t ii = 0; ii <= fOrder; ++ii)
 129     for (Int_t jj = 0; jj < 2; ++jj)
 130       fRedFactp[ii][jj] = 1.;
 131 }
 132
 133 //________________________________________________________________________
 134 void AliEbyEPidRatioPhy::CreateHistograms() {
 135   Float_t etaRange[2];
 136   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 137
 138   Float_t ptRange[2];
 139   fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 140   TString sTitle("");
 141
 142  if (fIsRatio) AddHistSetRatio("Ratio",       Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], ptRange[1]));
 143
 144   AddHistSetCent("Phy",       Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], ptRange[1]));
 145   //  AddHistSetCent("PhyTPC",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired()));
 146   // AddHistSetCent("PhyTOF",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1]));
 147
 148 #if USE_PHI
 149   AddHistSetCent("Phyphi",    Form("%s,#it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", sTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 150   //  AddHistSetCent("PhyTPCphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired(), fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 151   // AddHistSetCent("PhyTOFphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 152 #endif
 153
 154   if (fIsMC) {
 155     TString sMCTitle("");
 156
 157     if (fIsRatio) AddHistSetRatio("MCRatio",       Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], ptRange[1]));
 158
 159     AddHistSetCent("MC",      Form("%s", sTitle.Data()));
 160     AddHistSetCent("MCpt",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], ptRange[1]));
 161     // AddHistSetCent("MCTPC",   Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired()));
 162     //AddHistSetCent("MCTOF",   Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1]));
 163
 164 #if USE_PHI
 165     AddHistSetCent("MCphi",   Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sMCTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 166     //AddHistSetCent("MCTPCphi",Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sMCTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired(), fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 167     //AddHistSetCent("MCTOFphi",Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sMCTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 168 #endif
 169
 170   }
 171
 172   return;
 173 }
 174
 175 //________________________________________________________________________
 176 Int_t AliEbyEPidRatioPhy::ProcessTracks() {
 177   Float_t etaRange[2];
 178   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 179
 180   Float_t ptRange[2];
 181   fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 182   // -- Track Loop
 183   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 184     AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
 185     // -- Check if track is accepted for basic parameters
 186     if (!fHelper->IsTrackAcceptedBasicCharged(track))
 187       continue;
 188     // -- Check if accepted - ESD
 189     if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
 190       continue;
 191     // -- Check if accepted - AOD
 192     if (fAOD){
 193       AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
 194
 195       if (!trackAOD) {
 196         AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
 197         continue;
 198       }
 199       if (!trackAOD->TestFilterBit(fAODtrackCutBit))
 200         continue;
 201       if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
 202         continue;
 203     }
 204
 205     if (!fHelper->IsTrackAcceptedDCA(track))
 206       continue;
 207
 208     Int_t iPid = 0;
 209     Double_t pid[3];
 210     if      (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion)))   iPid = 1;
 211     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon)))   iPid = 2;
 212     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))) iPid = 3;
 213     else iPid = 0;
 214
 215     Double_t yP;
 216     if (iPid != 0 && !fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
 217       continue;
 218
 219     Int_t idxPart = (track->Charge() < 0) ? 0 : 1;
 220     // -- in pt Range
 221     fNp[0][0][idxPart] += 1;
 222     if(iPid != 0) fNp[0][iPid][idxPart] += 1;
 223     // -- in TPC pt Range
 224     /*    if (track->Pt() <= fHelper->GetMinPtForTOFRequired()) {
 225       fNp[1][0][idxPart] += 1;
 226       if(iPid != 0) fNp[1][iPid][idxPart] += 1;
 227     }
 228     // -- in TPC+TOF pt Range
 229     if (track->Pt() > fHelper->GetMinPtForTOFRequired()){
 230       fNp[2][0][idxPart] += 1;
 231       if(iPid != 0) fNp[2][iPid][idxPart] += 1;
 232       }*/
 233
 234 #if USE_PHI
 235     if(!fHelper->IsTrackAcceptedPhi(track))
 236       continue;
 237
 238     // -- in pt Range
 239     fNp[3][iPid][idxPart] += 1;
 240     if(iPid != 0) fNp[3][0][idxPart] += 1;
 241     // -- in TPC pt Range
 242     /*   if (track->Pt() <= fHelper->GetMinPtForTOFRequired()) {
 243       fNp[4][0][idxPart] += 1;
 244       if(iPid != 0)fNp[4][iPid][idxPart] += 1;
 245     }
 246     // -- in TPC+TOF pt Range
 247     if (track->Pt() > fHelper->GetMinPtForTOFRequired()) {
 248       fNp[5][0][idxPart] += 1;
 249       if(iPid != 0) fNp[5][iPid][idxPart] += 1;*/
 250     }
 251 #endif
 252   } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
 253
 254 FillHistSetCent("Phy",        0, kFALSE);
 255 if (fIsRatio) FillHistSetRatio("Ratio",        0, kFALSE);
 256 //  FillHistSetCent("PhyTPC",     1, kFALSE);
 257 //  FillHistSetCent("PhyTOF",     2, kFALSE);
 258
 259 #if USE_PHI
 260   FillHistSetCent("Phyphi",     1, kFALSE);
 261 //  FillHistSetCent("PhyTPCphi",  4, kFALSE);
 262 // FillHistSetCent("PhyTOFphi",  5, kFALSE);
 263
 264  #endif
 265   /*  Printf("<<<<<<<<<< Inside Loop >>>>>>>>>>");
 266   for (Int_t id = 0; id < 6; ++id) {
 267     printf(" == %2d == ", id);
 268     for (Int_t iPid = 0; iPid < 4; ++iPid) {
 269       printf("%7d %7d ", fNp[id][iPid][0] , fNp[id][iPid][1]);
 270     }
 271     Printf("");
 272     }  */
 273
 274
 275   return 0;
 276 }
 277
 278 //________________________________________________________________________
 279 Int_t AliEbyEPidRatioPhy::ProcessParticles() {
 280   // -- Process primary particles from the stack and fill histograms
 281
 282   Float_t etaRange[2];
 283   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 284
 285   Float_t ptRange[2];
 286   fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 287
 288   for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
 289     AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(idxMC) : NULL;
 290
 291     if (!particle)
 292       continue;
 293     if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
 294       continue;
 295
 296     Int_t iPid = 0;
 297     if      (TMath::Abs(particle->PdgCode()) ==  211) iPid = 1; // pion
 298     else if (TMath::Abs(particle->PdgCode()) ==  321) iPid = 2; // kaon
 299     else if (TMath::Abs(particle->PdgCode()) == 2212) iPid = 3; // proton
 300     else    iPid = 0;
 301
 302     Double_t yMC;
 303     if ((iPid != 0) && !fHelper->IsParticleAcceptedRapidity(particle, yMC, iPid))
 304       continue;
 305
 306     // -- Check eta window -- for charged particles
 307     if ((iPid == 0) && TMath::Abs(particle->Eta()) > etaRange[1])
 308       continue;
 309
 310     Int_t idxPart = (particle->PdgCode() < 0) ? 0 : 1;
 311
 312     fMCNp[0][0][idxPart]    += 1.;
 313     if(iPid != 0) fMCNp[0][iPid][idxPart] += 1.;
 314
 315     // -- Check main pt window
 316     if (!(particle->Pt() > ptRange[0] && particle->Pt() <= ptRange[1]))
 317       continue;
 318
 319     // -- in pt Range
 320     fMCNp[1][0][idxPart]    += 1.;
 321     if(iPid != 0)fMCNp[1][iPid][idxPart] += 1.;
 322
 323     // -- in TPC pt Range
 324     /*   if (particle->Pt() <= fHelper->GetMinPtForTOFRequired()) {
 325       fMCNp[2][0][idxPart]    += 1;
 326       if(iPid != 0)fMCNp[2][iPid][idxPart] += 1;
 327     }
 328     // -- in TPC+TOF pt Range
 329     if (particle->Pt() > fHelper->GetMinPtForTOFRequired()) {
 330       fMCNp[3][0][idxPart]    += 1;
 331       if(iPid != 0) fMCNp[3][iPid][idxPart] += 1;
 332       }*/
 333
 334 #if USE_PHI
 335     if(!fHelper->IsParticleAcceptedPhi(particle))
 336       continue;
 337     // idxPhi = 1;
 338
 339     // -- in pt Range
 340     fMCNp[4][0][idxPart]    += 1;
 341     if(iPid != 0)fMCNp[4][iPid][idxPart] += 1;
 342
 343     /*  // -- in TPC pt Range
 344     if (particle->Pt() <= fHelper->GetMinPtForTOFRequired()) {
 345       fMCNp[5][0][idxPart]    += 1;
 346       if(iPid != 0)fMCNp[5][iPid][idxPart] += 1;
 347     }
 348
 349     // -- in TPC+TOF pt Range
 350     if (particle->Pt() > fHelper->GetMinPtForTOFRequired()){
 351       fMCNp[6][0][idxPart]    += 1;
 352       if(iPid != 0)fMCNp[6][iPid][idxPart] += 1;
 353       }*/
 354 #endif
 355   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
 356
 357   FillHistSetCent("MC",        0, kTRUE);
 358   if (fIsRatio) FillHistSetRatio("MCRatio",     0, kTRUE);
 359   FillHistSetCent("MCpt",      1, kTRUE);
 360   //  FillHistSetCent("MCTPC",     2, kTRUE);
 361   // FillHistSetCent("MCTOF",     3, kTRUE);
 362
 363 #if USE_PHI
 364   FillHistSetCent("MCphi",     2, kTRUE);
 365
 366   // FillHistSetCent("MCTPCphi",  5, kTRUE);
 367   // FillHistSetCent("MCTOFphi",  6, kTRUE);
 368
 369
 370 #endif
 371
 372   return 0;
 373 }
 374
 375 //________________________________________________________________________
 376 void  AliEbyEPidRatioPhy::AddHistSetCent(const Char_t *name, const Char_t *title)  {
 377   TString sName(name);
 378   TString sTitle(title);
 379
 380   Float_t etaRange[2];
 381   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 382
 383   //TList *list[4];
 384   fOutList->Add(new TList);
 385   TList *list =  static_cast<TList*>(fOutList->Last());
 386   list->SetName(Form("f%s", name));
 387   list->SetOwner(kTRUE);
 388
 389
 390   Int_t nBinsCent         =  AliEbyEPidRatioHelper::fgkfHistNBinsCent;
 391   Double_t centBinRange[] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], AliEbyEPidRatioHelper::fgkfHistRangeCent[1]};
 392
 393   TString xyz = Form("|y| < %.1f",fHelper->GetRapidityMax());
 394
 395   list->Add(new TH2F(Form("fHistRatioKPi%s",name),
 396                      Form("(%s %s) : K/#pi;Centrality(11);K/#pi", xyz.Data(), sTitle.Data()),
 397                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 398
 399   list->Add(new TH2F(Form("fHistRatioKpPip%s",name),
 400                      Form("(%s %s) : K^{+}/#pi^{+};Centrality(11);K^{+}/#pi^{+}", xyz.Data(), sTitle.Data()),
 401                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 402
 403   list->Add(new TH2F(Form("fHistRatioKmPip%s",name),
 404                      Form("(%s %s) : K^{-}/#pi^{+};Centrality(11);K^{-}/#pi^{+}", xyz.Data(), sTitle.Data()),
 405                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 406
 407   list->Add(new TH2F(Form("fHistRatioKmPim%s",name),
 408                      Form("(%s %s) : K^{-}/#pi^{-};Centrality(11);K^{-}/#pi^{-}", xyz.Data(), sTitle.Data()),
 409                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 410
 411
 412
 413   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 414     // fOutList->Add(new TList);
 415     // list[iPid] = static_cast<TList*>(fOutList->Last());
 416     // list[iPid]->SetName(Form("f%s_%s", name,AliEbyEPidRatioHelper::fgkPidName[iPid]));
 417     // list[iPid]->SetOwner(kTRUE);
 418
 419     TString sNetTitle(Form("%s - %s", AliEbyEPidRatioHelper::fgkPidLatex[iPid][1], AliEbyEPidRatioHelper::fgkPidLatex[iPid][0]));
 420
 421     sTitle = (iPid != 0 ) ? Form("|y| < %.1f", fHelper->GetRapidityMax()) : Form(" |#eta|<%.1f", etaRange[1]);
 422
 423     for (Int_t idx = 1; idx <= fOrder; ++idx) {
 424       list->Add(new TProfile(Form("fProf%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid],name, idx),
 425                              Form("(%s)^{%d} : %s;Centrality(100);(%s)^{%d}",sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
 426                              100,-0.5,99.5));
 427     }
 428
 429     for (Int_t ii = 0; ii <= fOrder; ++ii) {
 430       for (Int_t kk = 0; kk <= fOrder; ++kk) {
 431         list->Add(new TProfile(Form("fProf%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk),
 432                                Form("f_{%02d%02d} : %s;Centrality(100);f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
 433                                100,-0.5,99.5));
 434       }
 435     }
 436
 437   }
 438
 439   /* fOutList->Add(new TList);
 440      TList *list_nu = static_cast<TList*>(fOutList->Last());
 441      list_nu->SetName(Form("f%s_nu", name));
 442      list_nu->SetOwner(kTRUE);
 443   */
 444
 445   for (Int_t iPhy = 0; iPhy < 46; ++iPhy) {
 446     list->Add(new TProfile(Form("fProf%sNu%02d",name,iPhy),Form("Physics Variable for index %d | %s ; Centrality;",iPhy,name),100,-0.5,99.5));
 447   }
 448
 449
 450   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 451     TString sNetTitle(Form("%s - %s", AliEbyEPidRatioHelper::fgkPidLatex[iPid][1], AliEbyEPidRatioHelper::fgkPidLatex[iPid][0]));
 452     sTitle = (iPid != 0 ) ? Form(" |y|<%.1f", fHelper->GetRapidityMax()) : Form(" |#eta| < %.1f", etaRange[1]);
 453
 454     for (Int_t idx = 1; idx <= fOrder; ++idx) {
 455       list->Add(new TProfile(Form("fProfBin%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid],name, idx),
 456                              Form("(%s)^{%d} : %s;Centrality(11);(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
 457                              nBinsCent, centBinRange[0], centBinRange[1]));
 458     }
 459
 460     for (Int_t ii = 0; ii <= fOrder; ++ii) {
 461       for (Int_t kk = 0; kk <= fOrder; ++kk) {
 462         list->Add(new TProfile(Form("fProfBin%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk),
 463                                Form("f_{%02d%02d} : %s;Centrality(11);f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
 464                                nBinsCent, centBinRange[0], centBinRange[1]));
 465       }
 466     }
 467
 468   }
 469
 470
 471   for (Int_t iPhy = 0; iPhy < 46; ++iPhy) {
 472     list->Add(new TProfile(Form("fProfBin%sNu%02d",name,iPhy),Form("Physics Variable for index %d | %s ; Centrality;",iPhy,name),nBinsCent, centBinRange[0], centBinRange[1]));
 473   }
 474
 475   return;
 476 }
 477
 478
 479
 480 //________________________________________________________________________
 481 void  AliEbyEPidRatioPhy::AddHistSetRatio(const Char_t *name, const Char_t *title)  {
 482   TString sName(name);
 483   TString sTitle(title);
 484
 485   Float_t etaRange[2];
 486   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 487
 488   //TList *list[4];
 489   fOutList->Add(new TList);
 490   TList *list =  static_cast<TList*>(fOutList->Last());
 491   list->SetName(Form("f%s", name));
 492   list->SetOwner(kTRUE);
 493
 494
 495   Int_t nBinsCent         =  AliEbyEPidRatioHelper::fgkfHistNBinsCent;
 496   Double_t centBinRange[] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], AliEbyEPidRatioHelper::fgkfHistRangeCent[1]};
 497
 498   TString xyz = Form("|y| < %.1f",fHelper->GetRapidityMax());
 499
 500   list->Add(new TH2F(Form("fHistRatioKPi%s",name),
 501                      Form("(%s %s) : K/#pi;Centrality(11);K/#pi", xyz.Data(), sTitle.Data()),
 502                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 503
 504   list->Add(new TH2F(Form("fHistRatioKpPip%s",name),
 505                      Form("(%s %s) : K^{+}/#pi^{+};Centrality(11);K^{+}/#pi^{+}", xyz.Data(), sTitle.Data()),
 506                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 507
 508   list->Add(new TH2F(Form("fHistRatioKmPip%s",name),
 509                      Form("(%s %s) : K^{-}/#pi^{+};Centrality(11);K^{-}/#pi^{+}", xyz.Data(), sTitle.Data()),
 510                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 511
 512   list->Add(new TH2F(Form("fHistRatioKmPim%s",name),
 513                      Form("(%s %s) : K^{-}/#pi^{-};Centrality(11);K^{-}/#pi^{-}", xyz.Data(), sTitle.Data()),
 514                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 515
 516
 517
 518  list->Add(new TH2F(Form("fHistRatioPK%s",name),
 519                      Form("(%s %s) : P/K;Centrality(11);P/K", xyz.Data(), sTitle.Data()),
 520                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 521
 522   list->Add(new TH2F(Form("fHistRatioPpKp%s",name),
 523                      Form("(%s %s) : P/K^{+};Centrality(11);P/K^{+}", xyz.Data(), sTitle.Data()),
 524                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 525
 526   list->Add(new TH2F(Form("fHistRatioPmKp%s",name),
 527                      Form("(%s %s) : #bar{P}/K^{+};Centrality(11);#bar{P}/K^{+}", xyz.Data(), sTitle.Data()),
 528                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 529
 530   list->Add(new TH2F(Form("fHistRatioPmKm%s",name),
 531                      Form("(%s %s) : #bar{P}/K^{-};Centrality(11);#bar{P}/K^{-}", xyz.Data(), sTitle.Data()),
 532                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 533
 534
 535
 536  list->Add(new TH2F(Form("fHistRatioPPi%s",name),
 537                      Form("(%s %s) : P/#pi;Centrality(11);K/#pi", xyz.Data(), sTitle.Data()),
 538                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 539
 540   list->Add(new TH2F(Form("fHistRatioPpPip%s",name),
 541                      Form("(%s %s) : P/#pi^{+};Centrality(11);P/#pi^{+}", xyz.Data(), sTitle.Data()),
 542                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 543
 544   list->Add(new TH2F(Form("fHistRatioPmPip%s",name),
 545                      Form("(%s %s) : #bar{P}/#pi^{+};Centrality(11);#bar{P}/#pi^{+}", xyz.Data(), sTitle.Data()),
 546                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 547
 548   list->Add(new TH2F(Form("fHistRatioPmPim%s",name),
 549                      Form("(%s %s) : #bar{P}/#pi^{-};Centrality(11);#bar{P}/#pi^{-}", xyz.Data(), sTitle.Data()),
 550                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 551
 552
 553
 554   return;
 555 }
 556
 557
 558
 559
 560 //________________________________________________________________________
 561 void AliEbyEPidRatioPhy::FillHistSetCent(const Char_t *name, Int_t idx, Bool_t isMC)  {
 562   /*
 563     printf(" !! %2d !! ", idx);
 564     for (Int_t iPid = 0; iPid < 4; ++iPid) {
 565     printf("%7d %7d ", fNp[idx][iPid][0] , fNp[idx][iPid][1]);
 566     }
 567     Printf("");
 568   */
 569
 570   Int_t ***np = (isMC) ? fMCNp : fNp;
 571
 572     /*
 573           printf(" ** %2d ** ", idx);
 574           for (Int_t iPid = 0; iPid < 4; ++iPid) {
 575           printf("%7d %7d ", np[idx][iPid][0] , np[idx][iPid][1]);
 576           }
 577           Printf("");
 578     */
 579
 580     // TList *list[4];
 581
 582   Float_t centralityBin = fHelper->GetCentralityBin();
 583   Float_t centralityPer = fHelper->GetCentralityPercentile();//fHelper->GetCentralityBin();
 584
 585   TList *list = static_cast<TList*>(fOutList->FindObject(Form("f%s",name)));
 586
 587
 588
 589   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 590     Int_t deltaNp = np[idx][iPid][1]-np[idx][iPid][0];
 591     Double_t delta = 1.;
 592     for (Int_t idxOrder = 1; idxOrder <= fOrder; ++idxOrder) {
 593       delta *= deltaNp;
 594       (static_cast<TProfile*>(list->FindObject(Form("fProfBin%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid], name, idxOrder))))->Fill(centralityBin, delta);
 595       (static_cast<TProfile*>(list->FindObject(Form("fProf%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid], name, idxOrder))))->Fill(centralityPer, delta);
 596     }
 597
 598     for (Int_t idxOrder = 0; idxOrder <= fOrder; ++ idxOrder) {
 599       fRedFactp[idxOrder][0]  = 1.;
 600       fRedFactp[idxOrder][1]  = 1.;
 601     }
 602
 603     for (Int_t idxOrder = 1; idxOrder <= fOrder; ++ idxOrder) {
 604       fRedFactp[idxOrder][0]  = fRedFactp[idxOrder-1][0]  * Double_t(np[idx][iPid][0]-(idxOrder-1));
 605       fRedFactp[idxOrder][1]  = fRedFactp[idxOrder-1][1]  * Double_t(np[idx][iPid][1]-(idxOrder-1));
 606     }
 607
 608     for (Int_t ii = 0; ii <= fOrder; ++ii) {   // ii -> p    -> n1
 609       for (Int_t kk = 0; kk <= fOrder; ++kk) { // kk -> pbar -> n2
 610         Double_t fik = fRedFactp[ii][1] * fRedFactp[kk][0];   // n1 *n2 -> p * pbar
 611         (static_cast<TProfile*>(list->FindObject(Form("fProfBin%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk))))->Fill(centralityBin, fik);
 612         (static_cast<TProfile*>(list->FindObject(Form("fProf%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk))))->Fill(centralityPer, fik);
 613       }
 614     }
 615   }
 616
 617   Double_t KPi   = -1; if(np[idx][1][1]+np[idx][1][0] != 0 ) KPi = (np[idx][2][1]+np[idx][2][0])/(np[idx][1][1]+np[idx][1][0]);
 618   Double_t KpPip = -1;
 619   Double_t KmPip = -1;if (np[idx][1][1] != 0 ) { KpPip  = (np[idx][2][1])/(np[idx][1][1]); KmPip =  (np[idx][2][0])/(np[idx][1][1]); }
 620   Double_t KmPim = -1;if (np[idx][1][0] != 0) KmPim = (np[idx][2][0])/(np[idx][1][0]);
 621
 622   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKPi%s",name))))->Fill(centralityBin, KPi);
 623   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKpPip%s",name))))->Fill(centralityBin, KpPip);
 624   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKmPip%s",name))))->Fill(centralityBin, KmPip);
 625   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKmPim%s",name))))->Fill(centralityBin, KmPim);
 626
 627
 628   Double_t a[6][4]; Double_t b[22];
 629   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 630     a[0][iPid] = np[idx][iPid][1]+np[idx][iPid][0];       // 0  n+ + n-
 631     a[1][iPid] = np[idx][iPid][1];                        // 1  n+
 632     a[2][iPid] = np[idx][iPid][0];                        // 2  n-
 633     a[3][iPid] = np[idx][iPid][1]*np[idx][iPid][0];       // 3  n+ . n-
 634     a[4][iPid] = np[idx][iPid][1]*(np[idx][iPid][1]-1);   // 4  n+ (n+ - 1)
 635     a[5][iPid] = np[idx][iPid][0]*(np[idx][iPid][0]-1);   // 5  n- (n- - 1)
 636   }
 637
 638   b[0]  = a[0][0]*a[0][2];       // 24 N   K
 639   b[1]  = a[0][1]*a[0][2];       // 25 Pi  K
 640   b[2]  = a[1][1]*a[1][2];       // 26 pi+ k+
 641   b[3]  = a[1][1]*a[2][2];       // 27 pi+ k-
 642   b[4]  = a[2][1]*a[1][2];       // 28 pi- k+
 643   b[5]  = a[2][1]*a[2][2];       // 29 pi- k-
 644
 645   b[6]  = a[0][0]*a[0][3];       // 30 N   P
 646   b[7]  = a[0][2]*a[0][3];       // 31 K   P
 647   b[8]  = a[1][2]*a[1][3];       // 32 k+  p+
 648   b[9]  = a[1][2]*a[2][3];       // 33 k+  p-
 649   b[10] = a[2][2]*a[1][3];       // 34 k-  p+
 650   b[11] = a[2][2]*a[2][3];       // 35 k-  p-
 651
 652   b[12] = a[0][0]*a[0][1];       // 36 N  Pi
 653   b[13] = a[0][3]*a[0][1];       // 37 P  Pi
 654   b[14] = a[1][3]*a[1][1];       // 38 p+ pi+
 655   b[15] = a[1][3]*a[2][1];       // 39 p+ pi-
 656   b[16] = a[2][3]*a[1][1];       // 40 p- pi+
 657   b[17] = a[2][3]*a[2][1];       // 41 p- pi-
 658
 659   b[18] = a[0][0]*(a[0][0] - 1); // 42 N ( N - 1 )
 660   b[19] = a[0][1]*(a[0][1] - 1); // 43 Pi( Pi- 1 )
 661   b[20] = a[0][2]*(a[0][1] - 1); // 44 K ( K - 1 )
 662   b[21] = a[0][3]*(a[0][3] - 1); // 45 P ( P - 1 )
 663   // TList *list_nu = static_cast<TList*>(fOutList->FindObject(Form("f%s_nu",name)));
 664   Int_t k = 0;
 665   for (Int_t i = 0; i < 6; i++) {
 666     for (Int_t j = 0; j < 4; j++) {
 667       (static_cast<TProfile*>(list->FindObject(Form("fProfBin%sNu%02d", name,k))))->Fill(centralityBin,a[i][j]);
 668       (static_cast<TProfile*>(list->FindObject(Form("fProf%sNu%02d", name,k))))->Fill(centralityPer,a[i][j]);
 669       k++;
 670     }
 671   }
 672
 673   for (Int_t j = 0; j < 22; j++) {
 674     (static_cast<TProfile*>(list->FindObject(Form("fProfBin%sNu%02d", name,j+24))))->Fill(centralityBin,b[j]);
 675     (static_cast<TProfile*>(list->FindObject(Form("fProf%sNu%02d", name,j+24))))->Fill(centralityPer,b[j]);
 676   }
 677
 678   return;
 679 }
 680
 681
 682 //________________________________________________________________________
 683 void AliEbyEPidRatioPhy::FillHistSetRatio(const Char_t *name, Int_t idx, Bool_t isMC)  {
 684
 685   Int_t ***np = (isMC) ? fMCNp : fNp;
 686   Float_t centralityBin = fHelper->GetCentralityBin();
 687
 688   TList *list = static_cast<TList*>(fOutList->FindObject(Form("f%s",name)));
 689
 690   Double_t KPi   = -1; if(np[idx][1][1]+np[idx][1][0] != 0 ) KPi = (np[idx][2][1]+np[idx][2][0])/(np[idx][1][1]+np[idx][1][0]);
 691   Double_t KpPip = -1;
 692   Double_t KmPip = -1;if (np[idx][1][1] != 0 ) { KpPip  = (np[idx][2][1])/(np[idx][1][1]); KmPip =  (np[idx][2][0])/(np[idx][1][1]); }
 693   Double_t KmPim = -1;if (np[idx][1][0] != 0) KmPim = (np[idx][2][0])/(np[idx][1][0]);
 694
 695   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKPi%s",name))))->Fill(centralityBin, KPi);
 696   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKpPip%s",name))))->Fill(centralityBin, KpPip);
 697   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKmPip%s",name))))->Fill(centralityBin, KmPip);
 698   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKmPim%s",name))))->Fill(centralityBin, KmPim);
 699
 700   Double_t PK   = -1; if(np[idx][2][1]+np[idx][2][0] != 0 ) PK = (np[idx][3][1]+np[idx][3][0])/(np[idx][2][1]+np[idx][2][0]);
 701   Double_t PpKp = -1;
 702   Double_t PmKp = -1;if (np[idx][2][1] != 0 ) { PpKp  = (np[idx][3][1])/(np[idx][2][1]); PmKp =  (np[idx][3][0])/(np[idx][2][1]); }
 703   Double_t PmKm = -1;if (np[idx][2][0] != 0) PmKm = (np[idx][3][0])/(np[idx][2][0]);
 704
 705   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioPK%s",name))))->Fill(centralityBin, PK);
 706   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioPpKp%s",name))))->Fill(centralityBin, PpKp);
 707   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioPmKp%s",name))))->Fill(centralityBin, PmKp);
 708   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioPmKm%s",name))))->Fill(centralityBin, PmKm);
 709
 710   Double_t PPi   = -1; if(np[idx][1][1]+np[idx][1][0] != 0 ) PPi = (np[idx][3][1]+np[idx][3][0])/(np[idx][1][1]+np[idx][1][0]);
 711   Double_t PpPip = -1;
 712   Double_t PmPip = -1;if (np[idx][1][1] != 0 ) { PpPip  = (np[idx][3][1])/(np[idx][1][1]); PmPip =  (np[idx][3][0])/(np[idx][1][1]); }
 713   Double_t PmPim = -1;if (np[idx][1][0] != 0) PmPim = (np[idx][3][0])/(np[idx][1][0]);
 714
 715   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioPPi%s",name))))->Fill(centralityBin, PPi);
 716   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioPpPip%s",name))))->Fill(centralityBin, PpPip);
 717   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioPmPip%s",name))))->Fill(centralityBin, PmPip);
 718   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioPmPim%s",name))))->Fill(centralityBin, PmPim);
 719
 720   return;
 721 }
 722
 723
 724
 725