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   AddHistSetCent("Phy",       Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], ptRange[1]));
 142   //  AddHistSetCent("PhyTPC",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired()));
 143   // AddHistSetCent("PhyTOF",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1]));
 144
 145 #if USE_PHI
 146   AddHistSetCent("Phyphi",    Form("%s,#it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", sTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 147   //  AddHistSetCent("PhyTPCphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired(), fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 148   // AddHistSetCent("PhyTOFphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 149 #endif
 150
 151   if (fIsMC) {
 152     TString sMCTitle("");
 153
 154     AddHistSetCent("MC",      Form("%s", sTitle.Data()));
 155     AddHistSetCent("MCpt",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], ptRange[1]));
 156     // AddHistSetCent("MCTPC",   Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired()));
 157     //AddHistSetCent("MCTOF",   Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1]));
 158
 159 #if USE_PHI
 160     AddHistSetCent("MCphi",   Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sMCTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 161     //AddHistSetCent("MCTPCphi",Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sMCTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired(), fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 162     //AddHistSetCent("MCTOFphi",Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sMCTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
 163 #endif
 164
 165   }
 166
 167   return;
 168 }
 169
 170 //________________________________________________________________________
 171 Int_t AliEbyEPidRatioPhy::ProcessTracks() {
 172   Float_t etaRange[2];
 173   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 174
 175   Float_t ptRange[2];
 176   fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 177   // -- Track Loop
 178   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 179     AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
 180     // -- Check if track is accepted for basic parameters
 181     if (!fHelper->IsTrackAcceptedBasicCharged(track))
 182       continue;
 183     // -- Check if accepted - ESD
 184     if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
 185       continue;
 186     // -- Check if accepted - AOD
 187     if (fAOD){
 188       AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
 189
 190       if (!trackAOD) {
 191         AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
 192         continue;
 193       }
 194       if (!trackAOD->TestFilterBit(fAODtrackCutBit))
 195         continue;
 196       if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
 197         continue;
 198     }
 199
 200     if (!fHelper->IsTrackAcceptedDCA(track))
 201       continue;
 202
 203     Int_t iPid = 0;
 204     Double_t pid[3];
 205     if      (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion)))   iPid = 1;
 206     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon)))   iPid = 2;
 207     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))) iPid = 3;
 208     else iPid = 0;
 209
 210     Double_t yP;
 211     if (iPid != 0 && !fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
 212       continue;
 213
 214     Int_t idxPart = (track->Charge() < 0) ? 0 : 1;
 215     // -- in pt Range
 216     fNp[0][0][idxPart] += 1;
 217     if(iPid != 0) fNp[0][iPid][idxPart] += 1;
 218     // -- in TPC pt Range
 219     /*    if (track->Pt() <= fHelper->GetMinPtForTOFRequired()) {
 220       fNp[1][0][idxPart] += 1;
 221       if(iPid != 0) fNp[1][iPid][idxPart] += 1;
 222     }
 223     // -- in TPC+TOF pt Range
 224     if (track->Pt() > fHelper->GetMinPtForTOFRequired()){
 225       fNp[2][0][idxPart] += 1;
 226       if(iPid != 0) fNp[2][iPid][idxPart] += 1;
 227       }*/
 228
 229 #if USE_PHI
 230     if(!fHelper->IsTrackAcceptedPhi(track))
 231       continue;
 232
 233     // -- in pt Range
 234     fNp[3][iPid][idxPart] += 1;
 235     if(iPid != 0) fNp[3][0][idxPart] += 1;
 236     // -- in TPC pt Range
 237     /*   if (track->Pt() <= fHelper->GetMinPtForTOFRequired()) {
 238       fNp[4][0][idxPart] += 1;
 239       if(iPid != 0)fNp[4][iPid][idxPart] += 1;
 240     }
 241     // -- in TPC+TOF pt Range
 242     if (track->Pt() > fHelper->GetMinPtForTOFRequired()) {
 243       fNp[5][0][idxPart] += 1;
 244       if(iPid != 0) fNp[5][iPid][idxPart] += 1;*/
 245     }
 246 #endif
 247   } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
 248
 249   FillHistSetCent("Phy",        0, kFALSE);
 250 //  FillHistSetCent("PhyTPC",     1, kFALSE);
 251 //  FillHistSetCent("PhyTOF",     2, kFALSE);
 252
 253 #if USE_PHI
 254   FillHistSetCent("Phyphi",     1, kFALSE);
 255 //  FillHistSetCent("PhyTPCphi",  4, kFALSE);
 256 // FillHistSetCent("PhyTOFphi",  5, kFALSE);
 257
 258  #endif
 259   /*  Printf("<<<<<<<<<< Inside Loop >>>>>>>>>>");
 260   for (Int_t id = 0; id < 6; ++id) {
 261     printf(" == %2d == ", id);
 262     for (Int_t iPid = 0; iPid < 4; ++iPid) {
 263       printf("%7d %7d ", fNp[id][iPid][0] , fNp[id][iPid][1]);
 264     }
 265     Printf("");
 266     }  */
 267
 268
 269   return 0;
 270 }
 271
 272 //________________________________________________________________________
 273 Int_t AliEbyEPidRatioPhy::ProcessParticles() {
 274   // -- Process primary particles from the stack and fill histograms
 275
 276   Float_t etaRange[2];
 277   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 278
 279   Float_t ptRange[2];
 280   fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 281
 282   for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
 283     AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(idxMC) : NULL;
 284
 285     if (!particle)
 286       continue;
 287     if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
 288       continue;
 289
 290     Int_t iPid = 0;
 291     if      (TMath::Abs(particle->PdgCode()) ==  211) iPid = 1; // pion
 292     else if (TMath::Abs(particle->PdgCode()) ==  321) iPid = 2; // kaon
 293     else if (TMath::Abs(particle->PdgCode()) == 2212) iPid = 3; // proton
 294     else    iPid = 0;
 295
 296     Double_t yMC;
 297     if ((iPid != 0) && !fHelper->IsParticleAcceptedRapidity(particle, yMC, iPid))
 298       continue;
 299
 300     // -- Check eta window -- for charged particles
 301     if ((iPid == 0) && TMath::Abs(particle->Eta()) > etaRange[1])
 302       continue;
 303
 304     Int_t idxPart = (particle->PdgCode() < 0) ? 0 : 1;
 305
 306     fMCNp[0][0][idxPart]    += 1.;
 307     if(iPid != 0) fMCNp[0][iPid][idxPart] += 1.;
 308
 309     // -- Check main pt window
 310     if (!(particle->Pt() > ptRange[0] && particle->Pt() <= ptRange[1]))
 311       continue;
 312
 313     // -- in pt Range
 314     fMCNp[1][0][idxPart]    += 1.;
 315     if(iPid != 0)fMCNp[1][iPid][idxPart] += 1.;
 316
 317     // -- in TPC pt Range
 318     /*   if (particle->Pt() <= fHelper->GetMinPtForTOFRequired()) {
 319       fMCNp[2][0][idxPart]    += 1;
 320       if(iPid != 0)fMCNp[2][iPid][idxPart] += 1;
 321     }
 322     // -- in TPC+TOF pt Range
 323     if (particle->Pt() > fHelper->GetMinPtForTOFRequired()) {
 324       fMCNp[3][0][idxPart]    += 1;
 325       if(iPid != 0) fMCNp[3][iPid][idxPart] += 1;
 326       }*/
 327
 328 #if USE_PHI
 329     if(!fHelper->IsParticleAcceptedPhi(particle))
 330       continue;
 331     // idxPhi = 1;
 332
 333     // -- in pt Range
 334     fMCNp[4][0][idxPart]    += 1;
 335     if(iPid != 0)fMCNp[4][iPid][idxPart] += 1;
 336
 337     /*  // -- in TPC pt Range
 338     if (particle->Pt() <= fHelper->GetMinPtForTOFRequired()) {
 339       fMCNp[5][0][idxPart]    += 1;
 340       if(iPid != 0)fMCNp[5][iPid][idxPart] += 1;
 341     }
 342
 343     // -- in TPC+TOF pt Range
 344     if (particle->Pt() > fHelper->GetMinPtForTOFRequired()){
 345       fMCNp[6][0][idxPart]    += 1;
 346       if(iPid != 0)fMCNp[6][iPid][idxPart] += 1;
 347       }*/
 348 #endif
 349   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
 350
 351   FillHistSetCent("MC",        0, kTRUE);
 352   FillHistSetCent("MCpt",      1, kTRUE);
 353   //  FillHistSetCent("MCTPC",     2, kTRUE);
 354   // FillHistSetCent("MCTOF",     3, kTRUE);
 355
 356 #if USE_PHI
 357   FillHistSetCent("MCphi",     2, kTRUE);
 358   // FillHistSetCent("MCTPCphi",  5, kTRUE);
 359   // FillHistSetCent("MCTOFphi",  6, kTRUE);
 360
 361
 362 #endif
 363
 364   return 0;
 365 }
 366
 367 //________________________________________________________________________
 368 void  AliEbyEPidRatioPhy::AddHistSetCent(const Char_t *name, const Char_t *title)  {
 369   TString sName(name);
 370   TString sTitle(title);
 371
 372   Float_t etaRange[2];
 373   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 374
 375   //TList *list[4];
 376   fOutList->Add(new TList);
 377   TList *list =  static_cast<TList*>(fOutList->Last());
 378   list->SetName(Form("f%s", name));
 379   list->SetOwner(kTRUE);
 380
 381
 382   Int_t nBinsCent         =  AliEbyEPidRatioHelper::fgkfHistNBinsCent;
 383   Double_t centBinRange[] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], AliEbyEPidRatioHelper::fgkfHistRangeCent[1]};
 384
 385   TString xyz = Form("|y| < %.1f",fHelper->GetRapidityMax());
 386
 387   list->Add(new TH2F(Form("fHistRatioKPi%s",name),
 388                      Form("(%s %s) : K/#pi;Centrality(11);K/#pi", xyz.Data(), sTitle.Data()),
 389                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 390
 391   list->Add(new TH2F(Form("fHistRatioKpPip%s",name),
 392                      Form("(%s %s) : K^{+}/#pi^{+};Centrality(11);K^{+}/#pi^{+}", xyz.Data(), sTitle.Data()),
 393                      nBinsCent, centBinRange[0], centBinRange[1], 2500,0,0.25));
 394
 395   list->Add(new TH2F(Form("fHistRatioKmPip%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("fHistRatioKmPim%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
 404
 405   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 406     // fOutList->Add(new TList);
 407     // list[iPid] = static_cast<TList*>(fOutList->Last());
 408     // list[iPid]->SetName(Form("f%s_%s", name,AliEbyEPidRatioHelper::fgkPidName[iPid]));
 409     // list[iPid]->SetOwner(kTRUE);
 410
 411     TString sNetTitle(Form("%s - %s", AliEbyEPidRatioHelper::fgkPidLatex[iPid][1], AliEbyEPidRatioHelper::fgkPidLatex[iPid][0]));
 412
 413     sTitle = (iPid != 0 ) ? Form("|y| < %.1f", fHelper->GetRapidityMax()) : Form(" |#eta|<%.1f", etaRange[1]);
 414
 415     for (Int_t idx = 1; idx <= fOrder; ++idx) {
 416       list->Add(new TProfile(Form("fProf%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid],name, idx),
 417                              Form("(%s)^{%d} : %s;Centrality(100);(%s)^{%d}",sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
 418                              100,-0.5,99.5));
 419     }
 420
 421     for (Int_t ii = 0; ii <= fOrder; ++ii) {
 422       for (Int_t kk = 0; kk <= fOrder; ++kk) {
 423         list->Add(new TProfile(Form("fProf%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk),
 424                                Form("f_{%02d%02d} : %s;Centrality(100);f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
 425                                100,-0.5,99.5));
 426       }
 427     }
 428
 429   }
 430
 431   /* fOutList->Add(new TList);
 432      TList *list_nu = static_cast<TList*>(fOutList->Last());
 433      list_nu->SetName(Form("f%s_nu", name));
 434      list_nu->SetOwner(kTRUE);
 435   */
 436
 437   for (Int_t iPhy = 0; iPhy < 46; ++iPhy) {
 438     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));
 439   }
 440
 441
 442   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 443     TString sNetTitle(Form("%s - %s", AliEbyEPidRatioHelper::fgkPidLatex[iPid][1], AliEbyEPidRatioHelper::fgkPidLatex[iPid][0]));
 444     sTitle = (iPid != 0 ) ? Form(" |y|<%.1f", fHelper->GetRapidityMax()) : Form(" |#eta| < %.1f", etaRange[1]);
 445
 446     for (Int_t idx = 1; idx <= fOrder; ++idx) {
 447       list->Add(new TProfile(Form("fProfBin%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid],name, idx),
 448                              Form("(%s)^{%d} : %s;Centrality(11);(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
 449                              nBinsCent, centBinRange[0], centBinRange[1]));
 450     }
 451
 452     for (Int_t ii = 0; ii <= fOrder; ++ii) {
 453       for (Int_t kk = 0; kk <= fOrder; ++kk) {
 454         list->Add(new TProfile(Form("fProfBin%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk),
 455                                Form("f_{%02d%02d} : %s;Centrality(11);f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
 456                                nBinsCent, centBinRange[0], centBinRange[1]));
 457       }
 458     }
 459
 460   }
 461
 462
 463   for (Int_t iPhy = 0; iPhy < 46; ++iPhy) {
 464     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]));
 465   }
 466
 467   return;
 468 }
 469
 470 //________________________________________________________________________
 471 void AliEbyEPidRatioPhy::FillHistSetCent(const Char_t *name, Int_t idx, Bool_t isMC)  {
 472   /*
 473     printf(" !! %2d !! ", idx);
 474     for (Int_t iPid = 0; iPid < 4; ++iPid) {
 475     printf("%7d %7d ", fNp[idx][iPid][0] , fNp[idx][iPid][1]);
 476     }
 477     Printf("");
 478   */
 479
 480   Int_t ***np = (isMC) ? fMCNp : fNp;
 481
 482     /*
 483           printf(" ** %2d ** ", idx);
 484           for (Int_t iPid = 0; iPid < 4; ++iPid) {
 485           printf("%7d %7d ", np[idx][iPid][0] , np[idx][iPid][1]);
 486           }
 487           Printf("");
 488     */
 489
 490     // TList *list[4];
 491
 492   Float_t centralityBin = fHelper->GetCentralityBin();
 493   Float_t centralityPer = fHelper->GetCentralityPercentile();//fHelper->GetCentralityBin();
 494
 495   TList *list = static_cast<TList*>(fOutList->FindObject(Form("f%s",name)));
 496
 497
 498
 499   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 500     Int_t deltaNp = np[idx][iPid][1]-np[idx][iPid][0];
 501     Double_t delta = 1.;
 502     for (Int_t idxOrder = 1; idxOrder <= fOrder; ++idxOrder) {
 503       delta *= deltaNp;
 504       (static_cast<TProfile*>(list->FindObject(Form("fProfBin%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid], name, idxOrder))))->Fill(centralityBin, delta);
 505       (static_cast<TProfile*>(list->FindObject(Form("fProf%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid], name, idxOrder))))->Fill(centralityPer, delta);
 506     }
 507
 508     for (Int_t idxOrder = 0; idxOrder <= fOrder; ++ idxOrder) {
 509       fRedFactp[idxOrder][0]  = 1.;
 510       fRedFactp[idxOrder][1]  = 1.;
 511     }
 512
 513     for (Int_t idxOrder = 1; idxOrder <= fOrder; ++ idxOrder) {
 514       fRedFactp[idxOrder][0]  = fRedFactp[idxOrder-1][0]  * Double_t(np[idx][iPid][0]-(idxOrder-1));
 515       fRedFactp[idxOrder][1]  = fRedFactp[idxOrder-1][1]  * Double_t(np[idx][iPid][1]-(idxOrder-1));
 516     }
 517
 518     for (Int_t ii = 0; ii <= fOrder; ++ii) {   // ii -> p    -> n1
 519       for (Int_t kk = 0; kk <= fOrder; ++kk) { // kk -> pbar -> n2
 520         Double_t fik = fRedFactp[ii][1] * fRedFactp[kk][0];   // n1 *n2 -> p * pbar
 521         (static_cast<TProfile*>(list->FindObject(Form("fProfBin%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk))))->Fill(centralityBin, fik);
 522         (static_cast<TProfile*>(list->FindObject(Form("fProf%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk))))->Fill(centralityPer, fik);
 523       }
 524     }
 525   }
 526
 527   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]);
 528   Double_t KpPip = -1;
 529   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]); }
 530   Double_t KmPim = -1;if (np[idx][1][0] != 0) KmPim = (np[idx][2][0])/(np[idx][1][0]);
 531
 532   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKPi%s",name))))->Fill(centralityBin, KPi);
 533   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKpPip%s",name))))->Fill(centralityBin, KpPip);
 534   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKmPip%s",name))))->Fill(centralityBin, KmPip);
 535   (static_cast<TProfile*>(list->FindObject(Form("fHistRatioKmPim%s",name))))->Fill(centralityBin, KmPim);
 536
 537
 538   Double_t a[6][4]; Double_t b[22];
 539   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 540     a[0][iPid] = np[idx][iPid][1]+np[idx][iPid][0];       // 0  n+ + n-
 541     a[1][iPid] = np[idx][iPid][1];                        // 1  n+
 542     a[2][iPid] = np[idx][iPid][0];                        // 2  n-
 543     a[3][iPid] = np[idx][iPid][1]*np[idx][iPid][0];       // 3  n+ . n-
 544     a[4][iPid] = np[idx][iPid][1]*(np[idx][iPid][1]-1);   // 4  n+ (n+ - 1)
 545     a[5][iPid] = np[idx][iPid][0]*(np[idx][iPid][0]-1);   // 5  n- (n- - 1)
 546   }
 547
 548   b[0]  = a[0][0]*a[0][2];       // 24 N   K
 549   b[1]  = a[0][1]*a[0][2];       // 25 Pi  K
 550   b[2]  = a[1][1]*a[1][2];       // 26 pi+ k+
 551   b[3]  = a[1][1]*a[2][2];       // 27 pi+ k-
 552   b[4]  = a[2][1]*a[1][2];       // 28 pi- k+
 553   b[5]  = a[2][1]*a[2][2];       // 29 pi- k-
 554
 555   b[6]  = a[0][0]*a[0][3];       // 30 N   P
 556   b[7]  = a[0][2]*a[0][3];       // 31 K   P
 557   b[8]  = a[1][2]*a[1][3];       // 32 k+  p+
 558   b[9]  = a[1][2]*a[2][3];       // 33 k+  p-
 559   b[10] = a[2][2]*a[1][3];       // 34 k-  p+
 560   b[11] = a[2][2]*a[2][3];       // 35 k-  p-
 561
 562   b[12] = a[0][0]*a[0][1];       // 36 N  Pi
 563   b[13] = a[0][3]*a[0][1];       // 37 P  Pi
 564   b[14] = a[1][3]*a[1][1];       // 38 p+ pi+
 565   b[15] = a[1][3]*a[2][1];       // 39 p+ pi-
 566   b[16] = a[2][3]*a[1][1];       // 40 p- pi+
 567   b[17] = a[2][3]*a[2][1];       // 41 p- pi-
 568
 569   b[18] = a[0][0]*(a[0][0] - 1); // 42 N ( N - 1 )
 570   b[19] = a[0][1]*(a[0][1] - 1); // 43 Pi( Pi- 1 )
 571   b[20] = a[0][2]*(a[0][1] - 1); // 44 K ( K - 1 )
 572   b[21] = a[0][3]*(a[0][3] - 1); // 45 P ( P - 1 )
 573   // TList *list_nu = static_cast<TList*>(fOutList->FindObject(Form("f%s_nu",name)));
 574   Int_t k = 0;
 575   for (Int_t i = 0; i < 6; i++) {
 576     for (Int_t j = 0; j < 4; j++) {
 577       (static_cast<TProfile*>(list->FindObject(Form("fProfBin%sNu%02d", name,k))))->Fill(centralityBin,a[i][j]);
 578       (static_cast<TProfile*>(list->FindObject(Form("fProf%sNu%02d", name,k))))->Fill(centralityPer,a[i][j]);
 579       k++;
 580     }
 581   }
 582
 583   for (Int_t j = 0; j < 22; j++) {
 584     (static_cast<TProfile*>(list->FindObject(Form("fProfBin%sNu%02d", name,j+24))))->Fill(centralityBin,b[j]);
 585     (static_cast<TProfile*>(list->FindObject(Form("fProf%sNu%02d", name,j+24))))->Fill(centralityPer,b[j]);
 586   }
 587
 588   return;
 589 }
 590
 591
 592
 593