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   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 382     // fOutList->Add(new TList);
 383     // list[iPid] = static_cast<TList*>(fOutList->Last());
 384     // list[iPid]->SetName(Form("f%s_%s", name,AliEbyEPidRatioHelper::fgkPidName[iPid]));
 385     // list[iPid]->SetOwner(kTRUE);
 386
 387     TString sNetTitle(Form("%s - %s", AliEbyEPidRatioHelper::fgkPidLatex[iPid][1], AliEbyEPidRatioHelper::fgkPidLatex[iPid][0]));
 388
 389     sTitle = (iPid != 0 ) ? Form("|y| < %.1f", fHelper->GetRapidityMax()) : Form(" |#eta|<%.1f", etaRange[1]);
 390
 391     for (Int_t idx = 1; idx <= fOrder; ++idx) {
 392       list->Add(new TProfile(Form("fProf%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid],name, idx),
 393                              Form("(%s)^{%d} : %s;Centrality(100);(%s)^{%d}",sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
 394                              100,-0.5,99.5));
 395     }
 396
 397     for (Int_t ii = 0; ii <= fOrder; ++ii) {
 398       for (Int_t kk = 0; kk <= fOrder; ++kk) {
 399         list->Add(new TProfile(Form("fProf%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk),
 400                                Form("f_{%02d%02d} : %s;Centrality(100);f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
 401                                100,-0.5,99.5));
 402       }
 403     }
 404
 405   }
 406
 407   /* fOutList->Add(new TList);
 408      TList *list_nu = static_cast<TList*>(fOutList->Last());
 409      list_nu->SetName(Form("f%s_nu", name));
 410      list_nu->SetOwner(kTRUE);
 411   */
 412
 413   for (Int_t iPhy = 0; iPhy < 46; ++iPhy) {
 414     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));
 415   }
 416
 417   Int_t nBinsCent         =  AliEbyEPidRatioHelper::fgkfHistNBinsCent;
 418   Double_t centBinRange[] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], AliEbyEPidRatioHelper::fgkfHistRangeCent[1]};
 419
 420   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 421     TString sNetTitle(Form("%s - %s", AliEbyEPidRatioHelper::fgkPidLatex[iPid][1], AliEbyEPidRatioHelper::fgkPidLatex[iPid][0]));
 422     sTitle = (iPid != 0 ) ? Form(" |y|<%.1f", fHelper->GetRapidityMax()) : Form(" |#eta| < %.1f", etaRange[1]);
 423
 424     for (Int_t idx = 1; idx <= fOrder; ++idx) {
 425       list->Add(new TProfile(Form("fProfBin%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid],name, idx),
 426                              Form("(%s)^{%d} : %s;Centrality(11);(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
 427                              nBinsCent, centBinRange[0], centBinRange[1]));
 428     }
 429
 430     for (Int_t ii = 0; ii <= fOrder; ++ii) {
 431       for (Int_t kk = 0; kk <= fOrder; ++kk) {
 432         list->Add(new TProfile(Form("fProfBin%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk),
 433                                Form("f_{%02d%02d} : %s;Centrality(11);f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
 434                                nBinsCent, centBinRange[0], centBinRange[1]));
 435       }
 436     }
 437
 438   }
 439
 440
 441   for (Int_t iPhy = 0; iPhy < 46; ++iPhy) {
 442     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]));
 443   }
 444
 445   return;
 446 }
 447
 448 //________________________________________________________________________
 449 void AliEbyEPidRatioPhy::FillHistSetCent(const Char_t *name, Int_t idx, Bool_t isMC)  {
 450   /*
 451     printf(" !! %2d !! ", idx);
 452     for (Int_t iPid = 0; iPid < 4; ++iPid) {
 453     printf("%7d %7d ", fNp[idx][iPid][0] , fNp[idx][iPid][1]);
 454     }
 455     Printf("");
 456   */
 457
 458   Int_t ***np = (isMC) ? fMCNp : fNp;
 459
 460     /*
 461           printf(" ** %2d ** ", idx);
 462           for (Int_t iPid = 0; iPid < 4; ++iPid) {
 463           printf("%7d %7d ", np[idx][iPid][0] , np[idx][iPid][1]);
 464           }
 465           Printf("");
 466     */
 467
 468     // TList *list[4];
 469
 470   Float_t centralityBin = fHelper->GetCentralityBin();
 471   Float_t centralityPer = fHelper->GetCentralityPercentile();//fHelper->GetCentralityBin();
 472
 473   TList *list = static_cast<TList*>(fOutList->FindObject(Form("f%s",name)));
 474
 475   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 476     Int_t deltaNp = np[idx][iPid][1]-np[idx][iPid][0];
 477     Double_t delta = 1.;
 478     for (Int_t idxOrder = 1; idxOrder <= fOrder; ++idxOrder) {
 479       delta *= deltaNp;
 480       (static_cast<TProfile*>(list->FindObject(Form("fProfBin%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid], name, idxOrder))))->Fill(centralityBin, delta);
 481       (static_cast<TProfile*>(list->FindObject(Form("fProf%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid], name, idxOrder))))->Fill(centralityPer, delta);
 482     }
 483
 484     for (Int_t idxOrder = 0; idxOrder <= fOrder; ++ idxOrder) {
 485       fRedFactp[idxOrder][0]  = 1.;
 486       fRedFactp[idxOrder][1]  = 1.;
 487     }
 488
 489     for (Int_t idxOrder = 1; idxOrder <= fOrder; ++ idxOrder) {
 490       fRedFactp[idxOrder][0]  = fRedFactp[idxOrder-1][0]  * Double_t(np[idx][iPid][0]-(idxOrder-1));
 491       fRedFactp[idxOrder][1]  = fRedFactp[idxOrder-1][1]  * Double_t(np[idx][iPid][1]-(idxOrder-1));
 492     }
 493
 494     for (Int_t ii = 0; ii <= fOrder; ++ii) {   // ii -> p    -> n1
 495       for (Int_t kk = 0; kk <= fOrder; ++kk) { // kk -> pbar -> n2
 496         Double_t fik = fRedFactp[ii][1] * fRedFactp[kk][0];   // n1 *n2 -> p * pbar
 497         (static_cast<TProfile*>(list->FindObject(Form("fProfBin%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk))))->Fill(centralityBin, fik);
 498         (static_cast<TProfile*>(list->FindObject(Form("fProf%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk))))->Fill(centralityPer, fik);
 499       }
 500     }
 501   }
 502   Double_t a[6][4]; Double_t b[22];
 503   for (Int_t iPid = 0; iPid < 4; ++iPid) {
 504     a[0][iPid] = np[idx][iPid][1]+np[idx][iPid][0];       // 0  n+ + n-
 505     a[1][iPid] = np[idx][iPid][1];                        // 1  n+
 506     a[2][iPid] = np[idx][iPid][0];                        // 2  n-
 507     a[3][iPid] = np[idx][iPid][1]*np[idx][iPid][0];       // 3  n+ . n-
 508     a[4][iPid] = np[idx][iPid][1]*(np[idx][iPid][1]-1);   // 4  n+ (n+ - 1)
 509     a[5][iPid] = np[idx][iPid][0]*(np[idx][iPid][0]-1);   // 5  n- (n- - 1)
 510   }
 511
 512   b[0]  = a[0][0]*a[0][2];       // 24 N   K
 513   b[1]  = a[0][1]*a[0][2];       // 25 Pi  K
 514   b[2]  = a[1][1]*a[1][2];       // 26 pi+ k+
 515   b[3]  = a[1][1]*a[2][2];       // 27 pi+ k-
 516   b[4]  = a[2][1]*a[1][2];       // 28 pi- k+
 517   b[5]  = a[2][1]*a[2][2];       // 29 pi- k-
 518
 519   b[6]  = a[0][0]*a[0][3];       // 30 N   P
 520   b[7]  = a[0][2]*a[0][3];       // 31 K   P
 521   b[8]  = a[1][2]*a[1][3];       // 32 k+  p+
 522   b[9]  = a[1][2]*a[2][3];       // 33 k+  p-
 523   b[10] = a[2][2]*a[1][3];       // 34 k-  p+
 524   b[11] = a[2][2]*a[2][3];       // 35 k-  p-
 525
 526   b[12] = a[0][0]*a[0][1];       // 36 N  Pi
 527   b[13] = a[0][3]*a[0][1];       // 37 P  Pi
 528   b[14] = a[1][3]*a[1][1];       // 38 p+ pi+
 529   b[15] = a[1][3]*a[2][1];       // 39 p+ pi-
 530   b[16] = a[2][3]*a[1][1];       // 40 p- pi+
 531   b[17] = a[2][3]*a[2][1];       // 41 p- pi-
 532
 533   b[18] = a[0][0]*(a[0][0] - 1); // 42 N ( N - 1 )
 534   b[19] = a[0][1]*(a[0][1] - 1); // 43 Pi( Pi- 1 )
 535   b[20] = a[0][2]*(a[0][1] - 1); // 44 K ( K - 1 )
 536   b[21] = a[0][3]*(a[0][3] - 1); // 45 P ( P - 1 )
 537   // TList *list_nu = static_cast<TList*>(fOutList->FindObject(Form("f%s_nu",name)));
 538   Int_t k = 0;
 539   for (Int_t i = 0; i < 6; i++) {
 540     for (Int_t j = 0; j < 4; j++) {
 541       (static_cast<TProfile*>(list->FindObject(Form("fProfBin%sNu%02d", name,k))))->Fill(centralityBin,a[i][j]);
 542       (static_cast<TProfile*>(list->FindObject(Form("fProf%sNu%02d", name,k))))->Fill(centralityPer,a[i][j]);
 543       k++;
 544     }
 545   }
 546
 547   for (Int_t j = 0; j < 22; j++) {
 548     (static_cast<TProfile*>(list->FindObject(Form("fProfBin%sNu%02d", name,j+24))))->Fill(centralityBin,b[j]);
 549     (static_cast<TProfile*>(list->FindObject(Form("fProf%sNu%02d", name,j+24))))->Fill(centralityPer,b[j]);
 550   }
 551
 552   return;
 553 }