]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/SPECTRA/Nuclei/masses/AliAnalysisNucleiMass.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / masses / AliAnalysisNucleiMass.cxx
index ce71cd2daeb64763dcf732e59f00bdea499b4d10..c91b4770899898c911b125d25583639912412c9c 100644 (file)
 #include "TH2D.h"
 #include "TH1F.h"
 #include "TF1.h"
+#include "TF2.h"
 #include "TGraph.h"
 #include "TProfile.h"
 #include "AliESDtrackCuts.h"
 #include "AliAnalysisManager.h"
+#include "TFile.h"
 
 ClassImp(AliAnalysisNucleiMass)
 
 //_____________________________________________________________________________
 AliAnalysisNucleiMass::AliAnalysisNucleiMass():
   AliAnalysisTaskSE(),
-  fMC(kFALSE),
-  FilterBit(16),
-  NminTPCcluster(0),
-  DCAzCUT(100),
-  DCAxyCUT(0.1),
-  kTPCcut(kTRUE),
-  kTPC(0),
+  Centrality(),                         
+  FilterBit(16),                                
+  EtaLimit(),                           
+  DCAxyCut(0.1),                               
+  DCAzCut(1000.),                                
+  NsigmaTpcCut(2.0),                           
+  NminTpcCluster(0),                           
+  iTrdCut(0),
+  kSignalCheck(1),
+  iMtof(1),
+  kPvtxCorr(1),
+  iBconf(0),                                  
   kTOF(0),
-  iBconf(0),
-  fAOD(NULL),
+//iTriggerSel(-99),
+  fAOD(NULL), 
   fESD(NULL),
   fEvent(NULL),
   fPIDResponse(NULL)
 {
-   // Default constructor (should not be used)
-  fList1[0]=new TList();
-  fList1[0]->SetName("results");
+  Centrality[0]=0.0;
+  Centrality[1]=100.0;
+  
+  EtaLimit[0]=-99.0;
+  EtaLimit[1]=99.0;
+  
+  fList[0]=new TList();
+  fList[0]->SetName("results");
   
-  fList1[1]=new TList();
-  fList1[1]->SetName("results2");
+  fList[1]=new TList();
+  fList[1]->SetName("results2");
 }
 //______________________________________________________________________________
 AliAnalysisNucleiMass::AliAnalysisNucleiMass(const char *name):
   AliAnalysisTaskSE(name),
-  fMC(kFALSE),
-  FilterBit(16),
-  NminTPCcluster(0),
-  DCAzCUT(100),
-  DCAxyCUT(0.1),
-  kTPCcut(kTRUE),
-  kTPC(0),
+  Centrality(),                         
+  FilterBit(16),                                
+  EtaLimit(),                           
+  DCAxyCut(0.1),                               
+  DCAzCut(1000.),                                
+  NsigmaTpcCut(2.0),                           
+  NminTpcCluster(0),
+  iTrdCut(0),
+  kSignalCheck(1),
+  iMtof(1),
+  kPvtxCorr(1),
+  iBconf(0),                                  
   kTOF(0),
-  iBconf(0),
+  //iTriggerSel(-99),
   fAOD(NULL), 
   fESD(NULL),
   fEvent(NULL),
   fPIDResponse(NULL)
 {
-  fList1[0]=new TList();
+
+  Centrality[0]=0.0;
+  Centrality[1]=100.0;
+
+  EtaLimit[0]=-99.0;
+  EtaLimit[1]=99.0;
+
+  fList[0]=new TList();
   DefineOutput(1, TList::Class());
-  fList1[0]->SetName("results");
+  fList[0]->SetName("results");
   
-  fList1[1]=new TList();
+  fList[1]=new TList();
   DefineOutput(2, TList::Class());
-  fList1[1]->SetName("results2");
+  fList[1]->SetName("results2");
 }
 //_____________________________________________________________________________
 AliAnalysisNucleiMass::~AliAnalysisNucleiMass()
 {
-  if(fList1[0]) delete fList1[0];
-  if(fList1[1]) delete fList1[1];
+  if(fList[0]) delete fList[0];
+  if(fList[1]) delete fList[1];
 }
 //______________________________________________________________________________
 void AliAnalysisNucleiMass::UserCreateOutputObjects()
 {
+  Char_t namePart[nPart][30];
+  snprintf(namePart[0],30,"e");
+  snprintf(namePart[1],30,"#mu");
+  snprintf(namePart[2],30,"#pi");
+  snprintf(namePart[3],30,"K");
+  snprintf(namePart[4],30,"p");
+  snprintf(namePart[5],30,"d");
+  snprintf(namePart[6],30,"t");
+  snprintf(namePart[7],30,"He3");
+  snprintf(namePart[8],30,"He4");
   
-  for(Int_t iB=0;iB<2;iB++) {
-
-    hNevent[iB] = new TH1F("hNevent_Analyzed","Centrality(analyzed)",20,0,100);
-    
-    hNeventSelected[iB] = new TH1F("hNevent_Selected","Centrality(selected)",20,0,100);
-
-    hZvertex[iB] = new TH1F("hZvertex","Vertex distribution of selected events; z vertex (cm)",240,-30,30);
-    
-    hTOFSignalPion[iB] = new TH1F("hTOFSignalPion","TOF signal 0.9<p_{T}<1.0; t-t_{exp}^{#pi} (ps)",1500,-1500,1500);
-    
-    fdEdxVSp[iB][0] = new TH2F("fdEdxVSp","dE/dx vs p; p/|Z| (GeV/c); dE/dx_{TPC} (a.u.)",500,0,5,2000,0,1000);
-    fdEdxVSp[iB][1] = new TH2F("fdEdxVSp_pos","dE/dx vs p positive charge; p/|Z| (GeV/c); dE/dx_{TPC} (a.u.)",500,0,5,2000,0,1000);
-    fdEdxVSp[iB][2] = new TH2F("fdEdxVSp_neg","dE/dx vs p negative charge; p/|Z| (GeV/c); dE/dx_{TPC} (a.u.)",500,0,5,2000,0,1000);
-    
-    fBetaTofVSp[iB] = new TH2F("fBetaTofVSp","#beta_{TOF} vs p; p(GeV/c); #beta_{TOF}",1000,0,5,1300,0.4,1.05);
-    
-    fM2vsP_NoTpcCut[iB][0] = new TH2F("fM2vsP_NoTpcCut","M_{TOF}^{2} vs p; M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4}); p/|Z| (GeV/c)",8000,0,10,200,0,10);
-    fM2vsP_NoTpcCut[iB][1] = new TH2F("fM2vsP_NoTpcCut_Positive","M_{TOF}^{2} vs p Pos Part; M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4}); p/|Z| (GeV/c)",8000,0,10,200,0,10);
-    fM2vsP_NoTpcCut[iB][2] = new TH2F("fM2vsP_NoTpcCut_Negative","M_{TOF}^{2} vs p Neg Part; M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4}); p/|Z| (GeV/c)",8000,0,10,200,0,10);
-    
-    fM2vsP_NoTpcCut_DCAxyCut[iB][0] = new TH2F("fM2vsP_NoTpcCut_DCAxycut","M_{TOF}^{2} vs p with DCAxy cut; M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4}); p/|Z| (GeV/c)",8000,0,10,200,0,10);
-    fM2vsP_NoTpcCut_DCAxyCut[iB][1] = new TH2F("fM2vsP_NoTpcCut_Positive_DCAxycut","M_{TOF}^{2} vs p Pos Part with DCAxy cut; M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4}); p/|Z| (GeV/c)",8000,0,10,200,0,10);
-    fM2vsP_NoTpcCut_DCAxyCut[iB][2] = new TH2F("fM2vsP_NoTpcCut_Negative_DCAxycut","M_{TOF}^{2} vs p Neg Part with DCAxy cut; M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4}); p/|Z| (GeV/c)",8000,0,10,200,0,10);
-    
-    fM2vsZ[iB][0] = new TH2F("fM2vsZ","M_{TOF}^{2} vs Z^{2} Integrated p_{T};Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][1] = new TH2F("fM2vsZ_0.3pT0.5","M_{TOF}^{2} vs Z^{2} 0.3<pT<0.5;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][2] = new TH2F("fM2vsZ_0.5pT1.0","M_{TOF}^{2} vs Z^{2} 0.5<pT<1.0;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][3] = new TH2F("fM2vsZ_1.0pT1.5","M_{TOF}^{2} vs Z^{2} 1.0<pT<1.5;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][4] = new TH2F("fM2vsZ_1.5pT2.0","M_{TOF}^{2} vs Z^{2} 1.5<pT<2.0;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][5] = new TH2F("fM2vsZ_2.0pT2.5","M_{TOF}^{2} vs Z^{2} 2.0<pT<2.5;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][6] = new TH2F("fM2vsZ_2.5pT3.0","M_{TOF}^{2} vs Z^{2} 2.5<pT<3.0;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][7] = new TH2F("fM2vsZ_3.0pT3.5","M_{TOF}^{2} vs Z^{2} 3.0<pT<3.5;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][8] = new TH2F("fM2vsZ_3.5pT4.0","M_{TOF}^{2} vs Z^{2} 3.5<pT<4.0;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][9] = new TH2F("fM2vsZ_4.0pT4.5","M_{TOF}^{2} vs Z^{2} 4.0<pT<4.5;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][10] = new TH2F("fM2vsZ_0.0pT1.0","M_{TOF}^{2} vs Z^{2} 2.0<pT<2.5;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][11] = new TH2F("fM2vsZ_1.0pT2.0","M_{TOF}^{2} vs Z^{2} 2.5<pT<3.0;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][12] = new TH2F("fM2vsZ_2.0pT3.0","M_{TOF}^{2} vs Z^{2} 3.0<pT<3.5;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][13] = new TH2F("fM2vsZ_3.0pT4.0","M_{TOF}^{2} vs Z^{2} 3.5<pT<4.0;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZ[iB][14] = new TH2F("fM2vsZ_4.0pT4.5","M_{TOF}^{2} vs Z^{2} 4.0<pT<4.5;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    
-    fM2vsZwithTPC[iB][0] = new TH2F("fM2vsZwithTPC","M_{TOF}^{2} vs Z^{2} Integrated p_{T} withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][1] = new TH2F("fM2vsZwithTPC_0.3pT0.5","M_{TOF}^{2} vs Z^{2} 0.3<pT<0.5 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][2] = new TH2F("fM2vsZwithTPC_0.5pT1.0","M_{TOF}^{2} vs Z^{2} 0.5<pT<1.0 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][3] = new TH2F("fM2vsZwithTPC_1.0pT1.5","M_{TOF}^{2} vs Z^{2} 1.0<pT<1.5 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][4] = new TH2F("fM2vsZwithTPC_1.5pT2.0","M_{TOF}^{2} vs Z^{2} 1.5<pT<2.0 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][5] = new TH2F("fM2vsZwithTPC_2.0pT2.5","M_{TOF}^{2} vs Z^{2} 2.0<pT<2.5 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][6] = new TH2F("fM2vsZwithTPC_2.5pT3.0","M_{TOF}^{2} vs Z^{2} 2.5<pT<3.0 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][7] = new TH2F("fM2vsZwithTPC_3.0pT3.5","M_{TOF}^{2} vs Z^{2} 3.0<pT<3.5 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][8] = new TH2F("fM2vsZwithTPC_3.5pT4.0","M_{TOF}^{2} vs Z^{2} 3.5<pT<4.0 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][9] = new TH2F("fM2vsZwithTPC_4.0pT4.5","M_{TOF}^{2} vs Z^{2} 4.0<pT<4.5 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][10] = new TH2F("fM2vsZwithTPC_0.0pT1.0","M_{TOF}^{2} vs Z^{2} 2.0<pT<2.5 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][11] = new TH2F("fM2vsZwithTPC_1.0pT2.0","M_{TOF}^{2} vs Z^{2} 2.5<pT<3.0 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][12] = new TH2F("fM2vsZwithTPC_2.0pT3.0","M_{TOF}^{2} vs Z^{2} 3.0<pT<3.5 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][13] = new TH2F("fM2vsZwithTPC_3.0pT4.0","M_{TOF}^{2} vs Z^{2} 3.5<pT<4.0 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    fM2vsZwithTPC[iB][14] = new TH2F("fM2vsZwithTPC_4.0pT4.5","M_{TOF}^{2} vs Z^{2} 4.0<pT<4.5 withTPCcut;Z;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4})",4000,-4,4,1000,0,10);
-    
-    char namePart[9][30];
-    char namePart_par_TPC[9][40];
-    char namePart_title_TPC[9][120];
-    
-    char namePart_par_TOF[9][40];
-    char namePart_title_TOF[9][120];
-    
-    char namePart_par_ProfileTPC[9][40];
-    char namePart_title_ProfileTPC[9][80];
-    
-    char namePart_par_ProfileTOF[9][40];
-    char namePart_title_ProfileTOF[9][80];
-    
-    
-    snprintf(namePart[0],20,"e");
-    snprintf(namePart[1],20,"#mu");
-    snprintf(namePart[2],20,"#pi");
-    snprintf(namePart[3],20,"K");
-    snprintf(namePart[4],20,"p");
-    snprintf(namePart[5],20,"d");
-    snprintf(namePart[6],20,"t");
-    snprintf(namePart[7],20,"He3");
-    snprintf(namePart[8],20,"He4");
-    
-    for(Int_t i=0;i<9;i++) {
-      snprintf(namePart_par_TPC[i],40,"NsigmaTPC_%s",namePart[i]);
-      snprintf(namePart_title_TPC[i],120,"NsigmaTPC_%s;p_{T} (GeV/c);n_{#sigma_{TPC}}^{%s}",namePart[i],namePart[i]);
-      
-      snprintf(namePart_par_TOF[i],40,"NsigmaTOF_%s",namePart[i]);
-      snprintf(namePart_title_TOF[i],120,"NsigmaTOF_%s;p_{T} (GeV/c);n_{#sigma_{TOF}}^{%s}",namePart[i],namePart[i]);
-      
-      snprintf(namePart_par_ProfileTPC[i],40,"hDeDxExp_%s",namePart[i]);
-      snprintf(namePart_title_ProfileTPC[i],80,"hDeDxExp_%s;p (GeV/c);dE/dx_{TPC} (a.u.)",namePart[i]);
-      
-      snprintf(namePart_par_ProfileTOF[i],40,"hBetaVsP_Exp_%s",namePart[i]);
-      snprintf(namePart_title_ProfileTOF[i],80,"hBetaVsP_Exp%s;p (GeV/c); #beta_{TOF}",namePart[i]);
-    }
-    
-    for(Int_t i=0;i<9;i++) {
-      fNsigmaTPC[iB][i] = new TH2F(namePart_par_TPC[i],namePart_title_TPC[i],250,0,5,200,-5,5);
-      fNsigmaTPC[iB][i]->GetYaxis()->CenterTitle();
-      fNsigmaTOF[iB][i] = new TH2F(namePart_par_TOF[i],namePart_title_TOF[i],250,0,5,200,-5,5);
-      fNsigmaTOF[iB][i]->GetYaxis()->CenterTitle();
-      hDeDxExp[iB][i] = new TProfile(namePart_par_ProfileTPC[i],namePart_title_ProfileTPC[i],500,0,5,0,1000,"");
-      hBetaExp[iB][i] = new TProfile(namePart_par_ProfileTOF[i],namePart_title_ProfileTOF[i],400,0,5,0.4,1.05,"");
-    }
-    
-    char namePart_par_TPCvsP_kTOFtrue[18][80];
-    char namePart_title_TPCvsP_kTOFtrue[18][120];
-    
-    char name[18][30];
-    
-    snprintf(name[0],20,"e^{+}");
-    snprintf(name[1],20,"#mu^{+}");
-    snprintf(name[2],20,"#pi^{+}");
-    snprintf(name[3],20,"K^{+}");
-    snprintf(name[4],20,"p");
-    snprintf(name[5],20,"d");
-    snprintf(name[6],20,"t");
-    snprintf(name[7],20,"He3");
-    snprintf(name[8],20,"He4");
-    
-    snprintf(name[9],20,"e^{-}");
-    snprintf(name[10],20,"#mu^{-}");
-    snprintf(name[11],20,"#pi^{-}");
-    snprintf(name[12],20,"K^{-}");
-    snprintf(name[13],20,"#bar{p}");
-    snprintf(name[14],20,"#bar{d}");
-    snprintf(name[15],20,"#bar{t}");
-    snprintf(name[16],20,"#bar{He3}");
-    snprintf(name[17],20,"#bar{He4}");
-    
-    for(Int_t iS=0;iS<18;iS++) {
-      snprintf(namePart_par_TPCvsP_kTOFtrue[iS],40,"NsigmaTPCvsP_kTOFout&&kTIME_%s",name[iS]);
-      snprintf(namePart_title_TPCvsP_kTOFtrue[iS],120,"NsigmaTPCvsP_kTOFout&&kTIME_%s;p (GeV/c);n_{#sigma_{TPC}}^{%s}",name[iS],name[iS]);
-    }
-    
-    for(Int_t iS=0;iS<18;iS++) {
-      fNsigmaTPCvsP_kTOFtrue[iB][iS] = new TH2F(namePart_par_TPCvsP_kTOFtrue[iS],namePart_title_TPCvsP_kTOFtrue[iS],250,0,5,200,-5,5);
-      fNsigmaTPCvsP_kTOFtrue[iB][iS]->GetYaxis()->CenterTitle();
-    }
-    
-    char name_par_MvsP[18][60];
-    char name_title_MvsP[18][150];
-    
-    for(Int_t i=0;i<18;i++) {
-      snprintf(name_par_MvsP[i],60,"M2vsP_%s",name[i]);
-      snprintf(name_title_MvsP[i],150,"M_{TOF}^{2}_%s_2#sigma_TPCcut;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4});p/|Z| (GeV/c)",name[i]);
-    }
+  Char_t name[nSpec][30];
+  snprintf(name[0],20,"e^{+}");
+  snprintf(name[1],20,"#mu^{+}");
+  snprintf(name[2],20,"#pi^{+}");
+  snprintf(name[3],20,"K^{+}");
+  snprintf(name[4],20,"p");
+  snprintf(name[5],20,"d");
+  snprintf(name[6],20,"t");
+  snprintf(name[7],20,"He3");
+  snprintf(name[8],20,"He4");
+  snprintf(name[9],20,"e^{-}");
+  snprintf(name[10],20,"#mu^{-}");
+  snprintf(name[11],20,"#pi^{-}");
+  snprintf(name[12],20,"K^{-}");
+  snprintf(name[13],20,"#bar{p}");
+  snprintf(name[14],20,"#bar{d}");
+  snprintf(name[15],20,"#bar{t}");
+  snprintf(name[16],20,"#bar{He3}");
+  snprintf(name[17],20,"#bar{He4}");
+  
+  Double_t binP[nbin+1];
+  for(Int_t i=0;i<nbin+1;i++) {
+    binP[i]=0.4+0.1*i;
+  }
+  
+  Char_t name_nbin[nbin][200];
+  for(Int_t j=0;j<nbin;j++) {
+    snprintf(name_nbin[j],200,"%.1f<P<%.1f",binP[j],binP[j+1]);
+  }
+  
+  for(Int_t iB=0;iB<nBconf;iB++) {
     
-    for (Int_t i=0;i<18;i++) fM2vsP[iB][i] = new TH2F(name_par_MvsP[i],name_title_MvsP[i],8000,0,10,200,0,10);
+    htemp[iB] = new TH1F("htemp","htemp (avoid the problem with the empty list...);B field",20,-10,10);
+
+    //htriggerbits[iB] = new TH1I("htriggerbits","htriggerbits; bits",10,-5,5);
+    htriggerbits[iB][0] = new TH1I("htriggerbits_0","trigger mask; bits",45,-5,40);
+    htriggerbits[iB][1] = new TH1I("htriggerbits_1","trigger bits (exclusive); bits",45,-5,40);
     
-    char name_par_MvsP_DCAxyCut[18][60];
-    char name_title_MvsP_DCAxyCut[18][150];
+    hCentrality[iB][0] = new TH1F("hCentrality_Selected","Centrality (selected events);centrality(%)",20,0,100);//20,0,100
+    hCentrality[iB][1] = new TH1F("hCentrality_Analyzed","Centrality (analyzed events);centrality (%)",20,0,100);//20,0,100
     
-    for(Int_t i=0;i<18;i++) {
-      snprintf(name_par_MvsP_DCAxyCut[i],60,"M2vsP_DCAxyCut_%s",name[i]);
-      snprintf(name_title_MvsP_DCAxyCut[i],150,"M_{TOF}^{2}_%s_2#sigma_TPCcut_DCAxyCut;M_{TOF}^{2}/Z^{2} (GeV^{2}/c^{4});p/|Z| (GeV/c)",name[i]);
-    }
+    hZvertex[iB][0] = new TH1F("hZvertex_Selected","Vertex distribution of selected events;z vertex (cm)",240,-30,30);
+    hZvertex[iB][1] = new TH1F("hZvertex_Analyzed","Vertex distribution of analyzed events;z vertex (cm)",240,-30,30);
     
-    for (Int_t i=0;i<18;i++) fM2vsP_DCAxyCut[iB][i] = new TH2F(name_par_MvsP_DCAxyCut[i],name_title_MvsP_DCAxyCut[i],8000,0,10,200,0,10);
-    
-    
-    Char_t par_name_nbin[nbin][30];
-    
-    snprintf(par_name_nbin[0],30,"0.4<Pt<0.5");
-    snprintf(par_name_nbin[1],30,"0.5<Pt<0.6");
-    snprintf(par_name_nbin[2],30,"0.6<Pt<0.7");
-    snprintf(par_name_nbin[3],30,"0.7<Pt<0.8");
-    snprintf(par_name_nbin[4],30,"0.8<Pt<0.9");
-    snprintf(par_name_nbin[5],30,"0.9<Pt<1.0");
-    snprintf(par_name_nbin[6],30,"1.0<Pt<1.1");
-    snprintf(par_name_nbin[7],30,"1.1<Pt<1.2");
-    snprintf(par_name_nbin[8],30,"1.2<Pt<1.3");
-    snprintf(par_name_nbin[9],30,"1.3<Pt<1.4");
-    snprintf(par_name_nbin[10],30,"1.4<Pt<1.5");
-    snprintf(par_name_nbin[11],30,"1.5<Pt<1.6");
-    snprintf(par_name_nbin[12],30,"1.6<Pt<1.7");
-    snprintf(par_name_nbin[13],30,"1.7<Pt<1.8");
-    snprintf(par_name_nbin[14],30,"1.8<Pt<1.9");
-    snprintf(par_name_nbin[15],30,"1.9<Pt<2.0");
-    snprintf(par_name_nbin[16],30,"2.0<Pt<2.1");
-    snprintf(par_name_nbin[17],30,"2.1<Pt<2.2");
-    snprintf(par_name_nbin[18],30,"2.2<Pt<2.3");
-    snprintf(par_name_nbin[19],30,"2.3<Pt<2.4");
-    snprintf(par_name_nbin[20],30,"2.4<Pt<2.5");
-    snprintf(par_name_nbin[21],30,"2.5<Pt<2.6");
-    snprintf(par_name_nbin[22],30,"2.6<Pt<2.7");
-    snprintf(par_name_nbin[23],30,"2.7<Pt<2.8");
-    snprintf(par_name_nbin[24],30,"2.8<Pt<2.9");
-    snprintf(par_name_nbin[25],30,"2.9<Pt<3.0");
-    snprintf(par_name_nbin[26],30,"3.0<Pt<3.1");
-    snprintf(par_name_nbin[27],30,"3.1<Pt<3.2");
-    snprintf(par_name_nbin[28],30,"3.2<Pt<3.3");
-    snprintf(par_name_nbin[29],30,"3.3<Pt<3.4");
-    snprintf(par_name_nbin[30],30,"3.4<Pt<3.5");
-    snprintf(par_name_nbin[31],30,"3.5<Pt<3.6");
-    snprintf(par_name_nbin[32],30,"3.6<Pt<3.7");
-    snprintf(par_name_nbin[33],30,"3.7<Pt<3.8");
-    snprintf(par_name_nbin[34],30,"3.8<Pt<3.9");
-    snprintf(par_name_nbin[35],30,"3.9<Pt<4.0");
-    snprintf(par_name_nbin[36],30,"4.0<Pt<4.1");
-    snprintf(par_name_nbin[37],30,"4.1<Pt<4.2");
-    snprintf(par_name_nbin[38],30,"4.2<Pt<4.3");
-    snprintf(par_name_nbin[39],30,"4.3<Pt<4.4");
-    snprintf(par_name_nbin[40],30,"4.4<Pt<4.5");
-    snprintf(par_name_nbin[41],30,"4.5<Pt<4.6");
-    snprintf(par_name_nbin[42],30,"4.6<Pt<4.7");
-    snprintf(par_name_nbin[43],30,"4.7<Pt<4.8");
-    snprintf(par_name_nbin[44],30,"4.8<Pt<4.9");
-    snprintf(par_name_nbin[45],30,"4.9<Pt<5.0");
-    
-    
-    Char_t nameDCAxy[18][nbin][120];
-    Char_t titleDCAxy[18][nbin][120];
-    
-    Char_t nameDCAz[18][nbin][120];
-    Char_t titleDCAz[18][nbin][120];
-    
-    Char_t nameM2CutDCAxy[18][nbin][120];
-    Char_t titleM2CutDCAxy[18][nbin][120];
-    
-    Char_t nameM2CutGroundDCAxy[18][nbin][120];
-    Char_t titleM2CutGroundDCAxy[18][nbin][120];
-    
-    
-    for(Int_t iS=0;iS<18;iS++) {
-      for(Int_t j=0;j<nbin;j++) {
-       snprintf(nameDCAxy[iS][j],120,"hDCAxy_%s_%s",name[iS],par_name_nbin[j]);
-       snprintf(titleDCAxy[iS][j],120,"hDCAxy_%s_%s;DCA_{xy} (cm)",name[iS],par_name_nbin[j]);
-       
-       snprintf(nameDCAz[iS][j],120,"hDCAz_%s_%s",name[iS],par_name_nbin[j]);
-       snprintf(titleDCAz[iS][j],120,"hDCAz_%s_%s;DCA_{z} (cm)",name[iS],par_name_nbin[j]);
-       
-       snprintf(nameM2CutDCAxy[iS][j],120,"hM2_CutDCAxy_%s_%s",name[iS],par_name_nbin[j]);
-       snprintf(titleM2CutDCAxy[iS][j],120,"hM2_CutDCAxy_%s_%s;M^{2}_{TOF} (GeV^{2}/c^{4})",name[iS],par_name_nbin[j]);
-       
-       snprintf(nameM2CutGroundDCAxy[iS][j],120,"hM2_GroundCatDCAxy_%s_%s",name[iS],par_name_nbin[j]);
-       snprintf(titleM2CutGroundDCAxy[iS][j],120,"hM2_GroundCatDCAxy_%s_%s;M^{2}_{TOF} (GeV^{2}/c^{4})",name[iS],par_name_nbin[j]);
-      }
-    }
+    hEta[iB] = new TH1F("hEta_Analyzed","|#eta| distribution after the track cuts;#eta",200,-1.0,1.0);
     
-    for(Int_t iS=0;iS<18;iS++) {
-      for(Int_t j=0;j<nbin;j++) {
-       hDCAxy[iB][iS][j] = new TH1D(nameDCAxy[iS][j],titleDCAxy[iS][j],875,-3.5,3.5);//125 bins
-       hDCAxy[iB][iS][j]->GetXaxis()->CenterTitle();
-       hDCAz[iB][iS][j] = new TH1D(nameDCAz[iS][j],titleDCAz[iS][j],875,-3.5,3.5);//125 bins
-       hDCAz[iB][iS][j]->GetXaxis()->CenterTitle();
-      }
-    }
+    hPhi[iB] = new TH1F("hPhi_Analyzed","#phi distribution after the track cuts;#phi (rad.)",90,0,6.3);//Each TRD supermodule is divided for 5 (DeltaPhi(TRD)=0.35 theoretical)
     
-    //for e,#mu,#pi and antiparticle (e and #mu will not be drawn)
-    //the binning is chosen for #pi distribution:
-    for(Int_t iSp=0;iSp<3;iSp++) {
-      for(Int_t j=0;j<nbin;j++) {
-       hM2CutDCAxy[iB][iSp][j] = new TH1D(nameM2CutDCAxy[iSp][j],titleM2CutDCAxy[iSp][j],600,-0.1,0.5);
-       hM2CutGroundDCAxy[iB][iSp][j] = new TH1D(nameM2CutGroundDCAxy[iSp][j],titleM2CutGroundDCAxy[iSp][j],600,-0.1,0.5);
-       hM2CutDCAxy[iB][iSp+9][j] = new TH1D(nameM2CutDCAxy[iSp+9][j],titleM2CutDCAxy[iSp+9][j],600,-0.1,0.5);
-       hM2CutGroundDCAxy[iB][iSp+9][j] = new TH1D(nameM2CutGroundDCAxy[iSp+9][j],titleM2CutGroundDCAxy[iSp+9][j],600,-0.1,0.5);
-      }
+    Int_t hbins[2];
+    if(kSignalCheck!=0) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=100; hbins[1]=90;} to reduce RAM consuming (toram)
+    else {hbins[0]=1; hbins[1]=1;}
+    fEtaPhi[iB] = new TH2F("fEtaPhi_Analyzed","#eta vs. #phi after the track cuts;#eta;#phi (rad.)",hbins[0],-1.0,1.0,hbins[1],0,6.3);
+
+    hNTpcCluster[iB] = new TH1F("hNTpcCluster","Number of the TPC clusters after the track cuts;n_{cl}^{TPC}",300,0,300);
+
+    hNTrdSlices[iB] = new TH1F("hNTrdSlices","Number of the TRD slices after the track cuts;n_{slices}^{TRD}",40,0,40);
+
+    if(kSignalCheck==1) {hbins[0]=500; hbins[1]=2000;}
+    else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
+    else if(kSignalCheck==2) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=100; hbins[1]=500;} toram
+    fdEdxVSp[iB][0] = new TH2F("fdEdxVSp_pos","dE/dx vs p (positive charge); p/|z| (GeV/c); dE/dx_{TPC} (a.u.)",hbins[0],0,5,hbins[1],0,1000);
+    fdEdxVSp[iB][1] = new TH2F("fdEdxVSp_neg","dE/dx vs p (negative charge); p/|z| (GeV/c); dE/dx_{TPC} (a.u.)",hbins[0],0,5,hbins[1],0,1000);
+
+    Char_t name_hDeDxExp[nPart][200];
+    Char_t title_hDeDxExp[nPart][200];
+    for(Int_t i=0;i<nPart;i++) {
+      snprintf(name_hDeDxExp[i],200,"hDeDxExp_%s",namePart[i]);
+      snprintf(title_hDeDxExp[i],200,"Expected dE/dx of %s in the TPC;p/|z| (GeV/c);dE/dx_{TPC} (a.u.)",namePart[i]);
+      hDeDxExp[iB][i] = new TProfile(name_hDeDxExp[i],title_hDeDxExp[i],1,0,5,0,1,"");//,500,0,5,0,1000,""); toram
     }
-    
-    for(Int_t j=0;j<nbin;j++) {
-      hM2CutDCAxy[iB][3][j] = new TH1D(nameM2CutDCAxy[3][j],titleM2CutDCAxy[3][j],400,0,1);
-      hM2CutGroundDCAxy[iB][3][j] = new TH1D(nameM2CutGroundDCAxy[3][j],titleM2CutGroundDCAxy[3][j],400,0,1);
-      hM2CutDCAxy[iB][3+9][j] = new TH1D(nameM2CutDCAxy[3+9][j],titleM2CutDCAxy[3+9][j],400,0,1);
-      hM2CutGroundDCAxy[iB][3+9][j] = new TH1D(nameM2CutGroundDCAxy[3+9][j],titleM2CutGroundDCAxy[3+9][j],400,0,1);
+
+    Char_t name_fNsigmaTpc[nSpec][200];
+    Char_t title_fNsigmaTpc[nSpec][200];
+    if(kSignalCheck==1) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=100; hbins[1]=100;} toram
+    else {hbins[0]=100; hbins[1]=100;}//temp!
+    for(Int_t i=0;i<nSpec;i++) {
+      snprintf(name_fNsigmaTpc[i],200,"NsigmaTpc_%s",name[i]);
+      snprintf(title_fNsigmaTpc[i],200,"NsigmaTpc_%s;p_{TPC}/|z| (GeV/c);n_{#sigma_{TPC}}^{%s}",name[i],name[i]);
+      fNsigmaTpc[iB][i] = new TH2F(name_fNsigmaTpc[i],title_fNsigmaTpc[i],hbins[0],0,5,hbins[1],-5,5);
     }
     
-    for(Int_t j=0;j<nbin;j++) {
-      hM2CutDCAxy[iB][4][j] = new TH1D(nameM2CutDCAxy[4][j],titleM2CutDCAxy[4][j],500,0,4);
-      hM2CutGroundDCAxy[iB][4][j] = new TH1D(nameM2CutGroundDCAxy[4][j],titleM2CutGroundDCAxy[4][j],500,0,4);
-      hM2CutDCAxy[iB][4+9][j] = new TH1D(nameM2CutDCAxy[4+9][j],titleM2CutDCAxy[4+9][j],500,0,4);
-      hM2CutGroundDCAxy[iB][4+9][j] = new TH1D(nameM2CutGroundDCAxy[4+9][j],titleM2CutGroundDCAxy[4+9][j],500,0,4);
+    if(kSignalCheck>0) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=100; hbins[1]=100;} toram
+    else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
+    Char_t name_fNsigmaTpc_kTOF[nSpec][200];
+    Char_t title_fNsigmaTpc_kTOF[nSpec][200];
+    for(Int_t i=0;i<nSpec;i++) {
+      snprintf(name_fNsigmaTpc_kTOF[i],200,"NsigmaTpc_%s_kTOF",name[i]);
+      snprintf(title_fNsigmaTpc_kTOF[i],200,"NsigmaTpc_kTOF_%s;p/|z| (GeV/c);n_{#sigma_{TPC}}^{%s}",name[i],name[i]);
+      fNsigmaTpc_kTOF[iB][i] = new TH2F(name_fNsigmaTpc_kTOF[i],title_fNsigmaTpc_kTOF[i],hbins[0],0,5,hbins[1],-5,5);
     }
-    
-    for(Int_t j=0;j<nbin;j++) {
-      hM2CutDCAxy[iB][5][j] = new TH1D(nameM2CutDCAxy[5][j],titleM2CutDCAxy[5][j],500,0,6);
-      hM2CutGroundDCAxy[iB][5][j] = new TH1D(nameM2CutGroundDCAxy[5][j],titleM2CutGroundDCAxy[5][j],500,0,6);
-      hM2CutDCAxy[iB][5+9][j] = new TH1D(nameM2CutDCAxy[5+9][j],titleM2CutDCAxy[5+9][j],500,0,6);
-      hM2CutGroundDCAxy[iB][5+9][j] = new TH1D(nameM2CutGroundDCAxy[5+9][j],titleM2CutGroundDCAxy[5+9][j],500,0,6);
+
+    if(kSignalCheck==1) {hbins[0]=1000; hbins[1]=1300;}
+    else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
+    else if(kSignalCheck==2) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=100; hbins[1]=260;}
+    fBetaTofVSp[iB][0] = new TH2F("fBetaTofVSp_pos","#beta_{TOF} vs p/|z| (positive charge);p(GeV/c);#beta_{TOF}",hbins[0],0,5,hbins[1],0.4,1.05);
+    fBetaTofVSp[iB][1] = new TH2F("fBetaTofVSp_neg","#beta_{TOF} vs p/|z| (negative charge);p(GeV/c);#beta_{TOF}",hbins[0],0,5,hbins[1],0.4,1.05);
+    
+    Char_t name_hBetaExp[nPart][200];
+    Char_t title_hBetaExp[nPart][200];
+    for(Int_t i=0;i<nPart;i++) {
+      snprintf(name_hBetaExp[i],200,"hBetaTofVsP_Exp_%s",namePart[i]);
+      snprintf(title_hBetaExp[i],200,"Expected #beta_{TOF} vs p/|z| of %s;p/|z| (GeV/c); #beta_{TOF}",namePart[i]);
+      hBetaExp[iB][i] = new TProfile(name_hBetaExp[i],title_hBetaExp[i],1,0,5,0.4,1.05,"");//,400,0,5,0.4,1.05,""); toram
     }
     
-    for(Int_t j=0;j<nbin;j++) {
-      hM2CutDCAxy[iB][6][j] = new TH1D(nameM2CutDCAxy[6][j],titleM2CutDCAxy[6][j],1000,0,12);
-      hM2CutGroundDCAxy[iB][6][j] = new TH1D(nameM2CutGroundDCAxy[6][j],titleM2CutGroundDCAxy[6][j],1000,0,12);
-      hM2CutDCAxy[iB][6+9][j] = new TH1D(nameM2CutDCAxy[6+9][j],titleM2CutDCAxy[6+9][j],1000,0,12);
-      hM2CutGroundDCAxy[iB][6+9][j] = new TH1D(nameM2CutGroundDCAxy[6+9][j],titleM2CutGroundDCAxy[6+9][j],1000,0,12);
+    Char_t name_fNsigmaTof[nPart][200];
+    Char_t title_fNsigmaTof[nPart][200];    
+    if(kSignalCheck==1) {hbins[0]=100; hbins[1]=100;}
+    else {hbins[0]=1; hbins[1]=1;}
+    for(Int_t i=0;i<nPart;i++) {
+      snprintf(name_fNsigmaTof[i],200,"NsigmaTof_%s",namePart[i]);
+      snprintf(title_fNsigmaTof[i],200,"NsigmaTof_%s;p_{T}/|z| (GeV/c);n_{#sigma_{TOF}}^{%s}",namePart[i],namePart[i]);
+      fNsigmaTof[iB][i] = new TH2F(name_fNsigmaTof[i],title_fNsigmaTof[i],hbins[0],0,5,hbins[1],-5,5);
     }
-    
-    for(Int_t j=0;j<nbin;j++) {
-      hM2CutDCAxy[iB][7][j] = new TH1D(nameM2CutDCAxy[7][j],titleM2CutDCAxy[7][j],200,0,4);
-      hM2CutGroundDCAxy[iB][7][j] = new TH1D(nameM2CutGroundDCAxy[7][j],titleM2CutGroundDCAxy[7][j],200,0,4);
-      hM2CutDCAxy[iB][7+9][j] = new TH1D(nameM2CutDCAxy[7+9][j],titleM2CutDCAxy[7+9][j],200,0,4);
-      hM2CutGroundDCAxy[iB][7+9][j] = new TH1D(nameM2CutGroundDCAxy[7+9][j],titleM2CutGroundDCAxy[7+9][j],200,0,4);
+
+    if(kSignalCheck==1) {hbins[0]=8000; hbins[1]=100;}
+    else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
+    else if(kSignalCheck==2) {hbins[0]=1; hbins[1]=1;}// {hbins[0]=1000; hbins[1]=100;} toram
+    fM2vsP_NoTpcCut[iB][0][0] = new TH2F("fM2vsP_NoTpcCut_pos","m^{2}/z^{2}_{TOF} vs p/|z| (positive charge);m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4});p/|z| (GeV/c)",hbins[0],0,10,hbins[1],0,5);
+    fM2vsP_NoTpcCut[iB][0][1] = new TH2F("fM2vsP_NoTpcCut_neg","m^{2}/z^{2}_{TOF} vs p/|z| (negative charge);m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4});p/|z| (GeV/c)",hbins[0],0,10,hbins[1],0,5);
+
+    Char_t name_fM2vsP[1][18][300]; 
+    Char_t title_fM2vsP[1][18][300]; 
+
+    for(Int_t i=0;i<nSpec;i++) {
+      snprintf(name_fM2vsP[0][i],300,"fM2vsPc_%s",name[i]);
+      snprintf(title_fM2vsP[0][i],300,"m^{2}/z^{2}_{TOF} vs p/|z| of %s with a NsigmaTpcCut (pReco->pTrue for nuclei);m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4});p/|z| (GeV/c)",name[i]);
+      
+      if(kSignalCheck==1) {hbins[0]=8000; hbins[1]=100;}
+      else {hbins[0]=1; hbins[1]=1;}
+      fM2vsP[iB][0][i] = new TH2F(name_fM2vsP[0][i],title_fM2vsP[0][i],hbins[0],0,10,hbins[1],0,5);
     }
     
-    for(Int_t j=0;j<nbin;j++) {
-      hM2CutDCAxy[iB][8][j] = new TH1D(nameM2CutDCAxy[8][j],titleM2CutDCAxy[8][j],600,0,6);
-      hM2CutGroundDCAxy[iB][8][j] = new TH1D(nameM2CutGroundDCAxy[8][j],titleM2CutGroundDCAxy[8][j],600,0,6);
-      hM2CutDCAxy[iB][8+9][j] = new TH1D(nameM2CutDCAxy[8+9][j],titleM2CutDCAxy[8+9][j],600,0,6);
-      hM2CutGroundDCAxy[iB][8+9][j] = new TH1D(nameM2CutGroundDCAxy[8+9][j],titleM2CutGroundDCAxy[8+9][j],600,0,6);
+    if(kSignalCheck==1) {hbins[0]=4000; hbins[1]=1000;}
+    else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
+    else if(kSignalCheck==2) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=1000 oppure 500; hbins[1]=100;} toram
+    fM2vsZ[iB][0] = new TH2F("fM2vsZ","m^{2}/z^{2}_{TOF} vs z_{TPC} Integrated p_{T};z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
+    fM2vsZ[iB][1] = new TH2F("fM2vsZ_0.5pT1.0","m^{2}/z^{2}_{TOF} vs z_{TPC} 0.5<pT<1.0;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
+    fM2vsZ[iB][2] = new TH2F("fM2vsZ_1.0pT1.5","m^{2}/z^{2}_{TOF} vs z_{TPC} 1.0<pT<1.5;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
+    fM2vsZ[iB][3] = new TH2F("fM2vsZ_1.5pT2.0","m^{2}/z^{2}_{TOF} vs z_{TPC} 1.5<pT<2.0;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
+    fM2vsZ[iB][4] = new TH2F("fM2vsZ_2.0pT2.5","m^{2}/z^{2}_{TOF} vs z_{TPC} 2.0<pT<2.5;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
+    fM2vsZ[iB][5] = new TH2F("fM2vsZ_2.5pT3.0","m^{2}/z^{2}_{TOF} vs z_{TPC} 2.5<pT<3.0;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
+    fM2vsZ[iB][6] = new TH2F("fM2vsZ_3.0pT3.5","m^{2}/z^{2}_{TOF} vs z_{TPC} 3.0<pT<3.5;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
+    fM2vsZ[iB][7] = new TH2F("fM2vsZ_3.5pT4.0","m^{2}/z^{2}_{TOF} vs z_{TPC} 3.5<pT<4.0;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
+    fM2vsZ[iB][8] = new TH2F("fM2vsZ_4.0pT4.5","m^{2}/z^{2}_{TOF} vs z_{TPC} 4.0<pT<4.5;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
+    fM2vsZ[iB][9] = new TH2F("fM2vsZ_4.5pT5.0","m^{2}/z^{2}_{TOF} vs z_{TPC} 2.0<pT<2.5;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
+  
+    Char_t name_h2DCAap[18][200];
+    Char_t title_h2DCAap[18][200];
+    
+    for(Int_t iS=0;iS<nSpec;iS++) {
+      snprintf(name_h2DCAap[iS],200,"h2DCAap_%s",name[iS]);
+      snprintf(title_h2DCAap[iS],200,"h2DCA_%s in for p/z<1.5GeV;DCA_{xy} (cm);DCA_{z} (cm)",name[iS]);
+      if(iS==5 || iS==7 || iS==5+9 || iS==7+9) h2DCAap[iB][iS] = new TH2F(name_h2DCAap[iS],title_h2DCAap[iS],1,-3.5,3.5,1,-3.5,3.5);//1750,-3.5,3.5,1750,-3.5,3.5
+      else h2DCAap[iB][iS] = new TH2F(name_h2DCAap[iS],title_h2DCAap[iS],1,-3.5,3.5,1,-3.5,3.5);//1750,-3.5,3.5,1750,-3.5,3.5
     }
+        
+    Char_t name_hDCAxy[18][nbin][200];
+    Char_t title_hDCAxy[18][nbin][200];
+    Char_t name_hDCAz[18][nbin][200];
+    Char_t title_hDCAz[18][nbin][200];
     
-    fList1[iB]->Add(hNeventSelected[iB]);
-    fList1[iB]->Add(hNevent[iB]);
-    fList1[iB]->Add(hZvertex[iB]);
-    fList1[iB]->Add(hTOFSignalPion[iB]);
-    for(Int_t iS=0;iS<18;iS++) fList1[iB]->Add(fNsigmaTPCvsP_kTOFtrue[iB][iS]);
-    
-    for(Int_t i=0;i<3;i++) fList1[iB]->Add(fdEdxVSp[iB][i]);
-    for(Int_t i=0;i<9;i++) fList1[iB]->Add(hDeDxExp[iB][i]);
-    fList1[iB]->Add(fBetaTofVSp[iB]);
-    for(Int_t i=0;i<9;i++) fList1[iB]->Add(hBetaExp[iB][i]);
-    for(Int_t i=0;i<9;i++) fList1[iB]->Add(fNsigmaTPC[iB][i]);
-    for(Int_t i=0;i<9;i++) fList1[iB]->Add(fNsigmaTOF[iB][i]);
-    
-    for(Int_t i=0;i<15;i++) fList1[iB]->Add(fM2vsZ[iB][i]);
-    for(Int_t i=0;i<15;i++) fList1[iB]->Add(fM2vsZwithTPC[iB][i]);
-    
-    for(Int_t i=0;i<3;i++) fList1[iB]->Add(fM2vsP_NoTpcCut[iB][i]);
-    for(Int_t i=0;i<18;i++) fList1[iB]->Add(fM2vsP[iB][i]);
+    //Char_t name_h2DCA[18][nbin][200];
+    //Char_t title_h2DCA[18][nbin][200];
     
-    for(Int_t i=0;i<3;i++) fList1[iB]->Add(fM2vsP_NoTpcCut_DCAxyCut[iB][i]);
-    for(Int_t i=0;i<18;i++) fList1[iB]->Add(fM2vsP_DCAxyCut[iB][i]);
-    
-    /*
-      for(Int_t j=0;j<nbin;j++) {//electron
-      fList1[iB]->Add(hDCAxy[iB][0][j]);
-      fList1[iB]->Add(hDCAz[iB][0][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][0][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][0][j]);
-      fList1[iB]->Add(hDCAxy[iB][9][j]);
-      fList1[iB]->Add(hDCAz[iB][9][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][9][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][9][j]);
+    for(Int_t iS=0;iS<nSpec;iS++) {
+      for(Int_t j=0;j<nbin;j++) {
+       snprintf(name_hDCAxy[iS][j],200,"hDCAxy_%s_%s",name[iS],name_nbin[j]);
+       snprintf(title_hDCAxy[iS][j],200,"hDCAxy_%s_%s in DCAzCut;DCA_{xy} (cm)",name[iS],name_nbin[j]);
+       if(iS==5 || iS==7 || iS==5+9 || iS==7+9) hDCAxy[iB][iS][j] = new TH1D(name_hDCAxy[iS][j],title_hDCAxy[iS][j],875,-3.5,3.5);
+       else hDCAxy[iB][iS][j] = new TH1D(name_hDCAxy[iS][j],title_hDCAxy[iS][j],1,-3.5,3.5);
+
+       snprintf(name_hDCAz[iS][j],200,"hDCAz_%s_%s",name[iS],name_nbin[j]);
+       snprintf(title_hDCAz[iS][j],200,"hDCAz_%s_%s in DCAxyCut;DCA_{z} (cm)",name[iS],name_nbin[j]);
+       if(iS==5 || iS==7 || iS==5+9 || iS==7+9) hDCAz[iB][iS][j] = new TH1D(name_hDCAz[iS][j],title_hDCAz[iS][j],875,-3.5,3.5);
+       else hDCAz[iB][iS][j] = new TH1D(name_hDCAz[iS][j],title_hDCAz[iS][j],1,-3.5,3.5);
+      
+       //snprintf(name_h2DCA[iS][j],200,"h2DCA_%s_%s",name[iS],name_nbin[j]);
+       //snprintf(title_h2DCA[iS][j],200,"h2DCA_%s_%s;DCA_{xy} (cm);DCA_{z} (cm)",name[iS],name_nbin[j]);
+       //if(iS==2 || iS==5 || iS==7 || iS==2+9 || iS==5+9 || iS==7+9) h2DCA[iB][iS][j] = new TH2F(name_h2DCA[iS][j],title_h2DCA[iS][j],1,-4,4,1,-4,4);//,160,-4,4,160,-4,4);
+       //else h2DCA[iB][iS][j] = new TH2F(name_h2DCA[iS][j],title_h2DCA[iS][j],1,-4,4,1,-4,4);
       }
+    }
       
-      for(Int_t j=0;j<nbin;j++) {//muon
-      fList1[iB]->Add(hDCAxy[iB][1][j]);
-      fList1[iB]->Add(hDCAz[iB][1][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][1][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][1][j]);
-      fList1[iB]->Add(hDCAxy[iB][10][j]);
-      fList1[iB]->Add(hDCAz[iB][10][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][10][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][10][j]);
+    Char_t name_hM2CutDCAxy[18][nbin][200];
+    Char_t title_hM2CutDCAxy[18][nbin][200];
+    for(Int_t iS=0;iS<nSpec;iS++) {
+      for(Int_t j=0;j<nbin;j++) {
+       snprintf(name_hM2CutDCAxy[iS][j],200,"hM2_CutDCAxy_%s_%s",name[iS],name_nbin[j]);
+       snprintf(title_hM2CutDCAxy[iS][j],200,"m^{2}/z^{2} Tof distribution of %s and in %s;m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",name[iS],name_nbin[j]);
       }
-    */
-    
-    for(Int_t j=0;j<nbin;j++) {//pion
-      fList1[iB]->Add(hDCAxy[iB][2][j]);
-      fList1[iB]->Add(hDCAz[iB][2][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][2][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][2][j]);
-      fList1[iB]->Add(hDCAxy[iB][11][j]);
-      fList1[iB]->Add(hDCAz[iB][11][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][11][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][11][j]);
     }
 
-    for(Int_t j=0;j<nbin;j++) {//kaon
-      fList1[iB]->Add(hDCAxy[iB][3][j]);
-      fList1[iB]->Add(hDCAz[iB][3][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][3][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][3][j]);
-      fList1[iB]->Add(hDCAxy[iB][12][j]);
-      fList1[iB]->Add(hDCAz[iB][12][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][12][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][12][j]);
+    const Int_t BinM2pT[nPart]={1,1,1,1000,500,500,1,400,1};//1,1,600,250,500,500,1000,400,600//1,1,1,250,500,500,1,400,1//1,1,1,1000,500,500,1,400,1
+    const Double_t RangeM2min[nPart]={0.0,0.0,-0.1,0.0,0.0,0.0,0.0,0.0,0.0};
+    const Double_t RangeM2max[nPart]={1.0,1.0,0.5,2.0,4.0,6.0,12.0,4.0,6.0};
+
+    for(Int_t iS=0;iS<nPart;iS++) {
+      for(Int_t j=0;j<nbin;j++) {
+       hM2CutDCAxy[iB][iS][j] = new TH1D(name_hM2CutDCAxy[iS][j],title_hM2CutDCAxy[iS][j],BinM2pT[iS],RangeM2min[iS],RangeM2max[iS]);
+       hM2CutDCAxy[iB][iS+nPart][j] = new TH1D(name_hM2CutDCAxy[iS+nPart][j],title_hM2CutDCAxy[iS+nPart][j],BinM2pT[iS],RangeM2min[iS],RangeM2max[iS]);
+      }
     }
-    
-    for(Int_t j=0;j<nbin;j++) {//proton
-      fList1[iB]->Add(hDCAxy[iB][4][j]);
-      fList1[iB]->Add(hDCAz[iB][4][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][4][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][4][j]);
-      fList1[iB]->Add(hDCAxy[iB][13][j]);
-      fList1[iB]->Add(hDCAz[iB][13][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][13][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][13][j]);
+
+    Char_t name_fPmeanVsBetaGamma[18][200];
+    Char_t title_fPmeanVsBetaGamma[18][200];
+    
+    if(iMtof==2) {hbins[0]=1; hbins[1]=1;}//if(iMtof==2) {hbins[0]=200; hbins[1]=200;}
+    else {hbins[0]=1; hbins[1]=1;}
+    for(Int_t iS=0;iS<nSpec;iS++) {
+      snprintf(name_fPmeanVsBetaGamma[iS],200,"fPmeanVsPvtx_%s",name[iS]);
+      snprintf(title_fPmeanVsBetaGamma[iS],200,"<p>/p_{vtx} vs #beta#gamma of %s;p_{vtx}/m_{%s};<p>_{%s}/p_{vtx}",name[iS],name[iS],name[iS]);
+      fPmeanVsBetaGamma[iB][iS]=new TH2F(name_fPmeanVsBetaGamma[iS],title_fPmeanVsBetaGamma[iS],hbins[0],0,10,hbins[1],0.8,1.2);
+    }  
+    
+    Char_t name_prPmeanVsBetaGamma[18][200];
+    Char_t title_prPmeanVsBetaGamma[18][200];
+    
+    if(iMtof==2) {hbins[0]=1; hbins[1]=1;}//if(iMtof==2) {hbins[0]=200; hbins[1]=200;}
+    else {hbins[0]=1; hbins[1]=1;}
+    for(Int_t iS=0;iS<nSpec;iS++) {
+      snprintf(name_prPmeanVsBetaGamma[iS],200,"prPmeanVsPvtx_%s",name[iS]);
+      snprintf(title_prPmeanVsBetaGamma[iS],200,"<p>/p_{vtx} vs #beta#gamma of %s;p_{vtx}/m_{%s};<p>_{%s}/p_{vtx}",name[iS],name[iS],name[iS]);
+      prPmeanVsBetaGamma[iB][iS]=new TProfile(name_prPmeanVsBetaGamma[iS],title_prPmeanVsBetaGamma[iS],hbins[0],0,10,0.8,1.2,"");
+    }  
+    
+    SetPvtxCorrections();
+
+    prPvtxTrueVsReco[iB][0]=new TProfile("prPvtxTrueVsReco_d","p_{true} vs p_{reco} of d and dbar;p_{reco} (GeV/c); p_{true}/p_{reco} (d)",1,0,10);//,100,0,10
+    prPvtxTrueVsReco[iB][1]=new TProfile("prPvtxTrueVsReco_t","p_{true} vs p_{reco} of t and tbar;p_{reco} (GeV/c);p_{true}/p_{reco} (t)",1,0,10);//,100,0,10
+    prPvtxTrueVsReco[iB][2]=new TProfile("prPvtxTrueVsReco_He3","p_{true} vs p_{reco} of He3 and He3bar;p_{reco} (GeV/c);p_{true}/p_{reco} (He3)",1,0,10);//,100,0,10
+    prPvtxTrueVsReco[iB][3]=new TProfile("prPvtxTrueVsReco_He4","p_{true} vs p_{reco} of He4 and He4bar;p_{reco} (GeV/c);p_{true}/p_{reco} (He4)",1,0,10);//,100,0,10
+
+    SetPmeanCorrections();
+       
+    Char_t nameTemp[14][200];
+    snprintf(nameTemp[0],200,"#pi^{+}");
+    snprintf(nameTemp[1],200,"K^{+}");
+    snprintf(nameTemp[2],200,"p");
+    snprintf(nameTemp[3],200,"d");
+    snprintf(nameTemp[4],200,"t");
+    snprintf(nameTemp[5],200,"He3");
+    snprintf(nameTemp[6],200,"He4");
+    snprintf(nameTemp[7],200,"#pi^{-}");
+    snprintf(nameTemp[8],200,"K^{-}");
+    snprintf(nameTemp[9],200,"#bar{p}");
+    snprintf(nameTemp[10],200,"#bar{d}");
+    snprintf(nameTemp[11],200,"#bar{t}");
+    snprintf(nameTemp[12],200,"#bar{He3}");
+    snprintf(nameTemp[13],200,"#bar{He4}");
+    Char_t name_prPmeanVsBGcorr[14][200];
+    Char_t title_prPmeanVsBGcorr[14][200];
+   
+    hbins[0]=200;
+    for(Int_t iS=0;iS<14;iS++) {
+      snprintf(name_prPmeanVsBGcorr[iS],200,"prPmeanVsBGcorr_%s",nameTemp[iS]);
+      snprintf(title_prPmeanVsBGcorr[iS],200,"<p>/p_{vtx} vs #beta#gamma of %s as parameterized in input TF1;p_{vtx}/m_{%s};<p>_{%s}/p_{vtx}",nameTemp[iS],nameTemp[iS],nameTemp[iS]);
+      prPmeanVsBGcorr[iB][iS]=new TProfile(name_prPmeanVsBGcorr[iS],title_prPmeanVsBGcorr[iS],hbins[0],0,20,0.8,1.2,"");
+    }  
+
+    fList[iB]->Add(htemp[iB]);
+    for(Int_t i=0;i<2;i++) fList[iB]->Add(htriggerbits[iB][i]);
+    for(Int_t i=0;i<2;i++) fList[iB]->Add(hCentrality[iB][i]);
+    for(Int_t i=0;i<2;i++) fList[iB]->Add(hZvertex[iB][i]);
+    fList[iB]->Add(hEta[iB]);
+    fList[iB]->Add(hPhi[iB]);
+    //fList[iB]->Add(fEtaPhi[iB]);
+    fList[iB]->Add(hNTpcCluster[iB]);
+    fList[iB]->Add(hNTrdSlices[iB]);
+    //for(Int_t i=0;i<2;i++) fList[iB]->Add(fdEdxVSp[iB][i]);
+    //for(Int_t i=0;i<nPart;i++) fList[iB]->Add(hDeDxExp[iB][i]);
+    for(Int_t i=0;i<nSpec;i++) fList[iB]->Add(fNsigmaTpc[iB][i]);
+    for(Int_t i=0;i<nPart;i++) {
+      if(kSignalCheck!=1) 
+       if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
+      //fList[iB]->Add(fNsigmaTpc_kTOF[iB][i]);
+      //fList[iB]->Add(fNsigmaTpc_kTOF[iB][i+nPart]);
     }
-    
-    for(Int_t j=0;j<nbin;j++) {//deuteron
-      fList1[iB]->Add(hDCAxy[iB][5][j]);
-      fList1[iB]->Add(hDCAz[iB][5][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][5][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][5][j]);
-      fList1[iB]->Add(hDCAxy[iB][14][j]);
-      fList1[iB]->Add(hDCAz[iB][14][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][14][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][14][j]);
+    //for(Int_t i=0;i<2;i++) fList[iB]->Add(fBetaTofVSp[iB][i]);
+    //for(Int_t i=0;i<nPart;i++) fList[iB]->Add(hBetaExp[iB][i]);
+    //for(Int_t i=0;i<nPart;i++) fList[iB]->Add(fNsigmaTof[iB][i]);
+    //for(Int_t i=0;i<2;i++) fList[iB]->Add(fM2vsP_NoTpcCut[iB][0][i]);
+    for(Int_t i=0;i<nPart;i++) {
+      if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
+      //fList[iB]->Add(fM2vsP[iB][0][i]);
+      //fList[iB]->Add(fM2vsP[iB][0][i+nPart]);
     }
-    
-    /*
-      for(Int_t j=0;j<nbin;j++) {//triton
-      fList1[iB]->Add(hDCAxy[iB][6][j]);
-      fList1[iB]->Add(hDCAz[iB][6][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][6][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][6][j]);
-      fList1[iB]->Add(hDCAxy[iB][15][j]);
-      fList1[iB]->Add(hDCAz[iB][15][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][15][j]);
-    fList1[iB]->Add(hM2CutGroundDCAxy[iB][15][j]);
+  
+    for(Int_t i=0;i<2;i++){
+      //fList[iB]->Add(fPvtxTrueVsReco[i]);
+      //fList[iB]->Add(prPvtxTrueVsReco[iB][i]);
     }
-    */
-    
-    for(Int_t j=0;j<nbin;j++) {//He3
-      fList1[iB]->Add(hDCAxy[iB][7][j]);
-      fList1[iB]->Add(hDCAz[iB][7][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][7][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][7][j]);
-      fList1[iB]->Add(hDCAxy[iB][16][j]);
-      fList1[iB]->Add(hDCAz[iB][16][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][16][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][16][j]);
+    if(iMtof==2) {
+      for(Int_t i=0;i<nPart;i++){
+       if(i<2) continue;//e,mu excluded
+       //fList[iB]->Add(fPmeanVsBetaGamma[iB][i]);
+       //fList[iB]->Add(prPmeanVsBetaGamma[iB][i]);
+       //fList[iB]->Add(fPmeanVsBetaGamma[iB][i+nPart]);
+       //fList[iB]->Add(prPmeanVsBetaGamma[iB][i+nPart]);
+      }
+    }
+    if(iMtof>2) {
+      //for(Int_t i=0;i<14;i++)fList[iB]->Add(fPmeanVsBGcorr[i]);
+      //for(Int_t i=0;i<14;i++)fList[iB]->Add(prPmeanVsBGcorr[iB][i]);
+    }
+    for(Int_t i=0;i<nPart;i++) {
+      if(i<5 || i==6 || i==8) continue;//e,mu,pi,K,p,t,he4 excluded//i<5 || i==6 || i==8
+      fList[iB]->Add(h2DCAap[iB][i]);
+      fList[iB]->Add(h2DCAap[iB][i+nPart]);
     }
-    
     /*
-      for(Int_t j=0;j<nbin;j++) {//4He
-      fList1[iB]->Add(hDCAxy[iB][8][j]);
-      fList1[iB]->Add(hDCAz[iB][8][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][8][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][8][j]);
-      fList1[iB]->Add(hDCAxy[iB][17][j]);
-      fList1[iB]->Add(hDCAz[iB][17][j]);
-      fList1[iB]->Add(hM2CutDCAxy[iB][17][j]);
-      fList1[iB]->Add(hM2CutGroundDCAxy[iB][17][j]);
+      for(Int_t i=0;i<nPart;i++) {
+      if(i<5 || i==6 || i==8) continue;//e,mu,pi,K,p,t,he4 excluded//i<5 || i==6 || i==8
+      for(Int_t j=0;j<nbin;j++){
+      fList[iB]->Add(h2DCA[iB][i][j]);
+      fList[iB]->Add(h2DCA[iB][i+nPart][j]);
+      }
       }
     */
+    //for(Int_t i=0;i<10;i++) fList[iB]->Add(fM2vsZ[iB][i]);
+    for(Int_t i=0;i<nPart;i++){
+      if(kSignalCheck!=1) 
+       if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
+      for(Int_t j=0;j<nbin;j++){
+       fList[iB]->Add(hDCAxy[iB][i][j]);
+       fList[iB]->Add(hDCAz[iB][i][j]);
+       fList[iB]->Add(hM2CutDCAxy[iB][i][j]);
+       fList[iB]->Add(hDCAxy[iB][i+nPart][j]);
+       fList[iB]->Add(hDCAz[iB][i+nPart][j]);
+       fList[iB]->Add(hM2CutDCAxy[iB][i+nPart][j]);
+      }
+    }
     
     // Post output data.
-    PostData(1, fList1[0]);
-    PostData(2, fList1[1]);
-  }//close the iB loop
+    PostData(1, fList[0]);
+    PostData(2, fList[1]);
+        
+  }//end iB loop
 }
 //______________________________________________________________________________
 void AliAnalysisNucleiMass::UserExec(Option_t *) 
@@ -532,262 +473,319 @@ void AliAnalysisNucleiMass::UserExec(Option_t *)
 
   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
-  fPIDResponse=inputHandler->GetPIDResponse(); // data member di tipo "const AliPIDResponse *fPIDResponse;"
-
-  //Centrality
-  Float_t v0Centr  = -10.;
-  //Float_t trkCentr  = -10.;
-  AliCentrality *centrality = fEvent->GetCentrality();
-  if (centrality){
-    v0Centr  = centrality->GetCentralityPercentile("V0M"); // VZERO
-    //trkCentr = centrality->GetCentralityPercentile("TRK"); // TPC
-  }
+  fPIDResponse=inputHandler->GetPIDResponse();
   
+  //--------------------------Magnetic field polarity--------------------
   Double_t fBfield=fEvent->GetMagneticField();
   if(fBfield<0.0) iBconf=0;//B--
   else iBconf=1;//B++
+  for(Int_t i=0;i<nBconf;i++) htemp[i]->Fill(fBfield);
+    
+  //--------------------------Centrality--------------------------------
+  Double_t v0Centr  = -10.;
+  AliCentrality *centrality = fEvent->GetCentrality();
+  if (centrality){
+    v0Centr=centrality->GetCentralityPercentile("V0M"); // VZERO
+  }
+  hCentrality[iBconf][0]->Fill(v0Centr);
 
-  hNeventSelected[iBconf]->Fill(v0Centr);//selected events
+  //-------------------------zVertex determination of event----------------
+  Double_t zvtx = 9999.9;
+  const AliVVertex* vtxEVENT = fEvent->GetPrimaryVertex();
+  if(vtxEVENT->GetNContributors()>0) zvtx = vtxEVENT->GetZ();
   
-  const AliAODVertex* vtxAOD = fAOD->GetPrimaryVertex();
+  hZvertex[iBconf][0]->Fill(zvtx);
   
-  Float_t zvtx = 10000.0;
+  //---------------------------EVENT CUTS-----------------------------
+  if(TMath::Abs(zvtx) < 10.0 && v0Centr>Centrality[0] && v0Centr<Centrality[1]){
 
-  if(vtxAOD->GetNContributors()>0)
-    zvtx = vtxAOD->GetZ();
-  
-   hZvertex[iBconf]->Fill(zvtx);
+    //TRIGGER SELECTION
+    Int_t iTrigger=-2;
+
+    if(inputHandler->IsEventSelected() & AliVEvent::kMB) iTrigger = 0;
+    if(inputHandler->IsEventSelected() & AliVEvent::kCentral) iTrigger = 16;
+    if(inputHandler->IsEventSelected() & AliVEvent::kSemiCentral) iTrigger = 17;
+    //if((((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()) & AliVEvent::kAny) iTrigger = 35;
     
-  if(TMath::Abs(zvtx) < 10.0){ // consistency cut on centrality selection AND d(zPrimaryVertez;NominalPointInteraction)<10cm
+    if(iTriggerSel!=-99) {//if a dedicated trigger is required
+      if(iTrigger!=iTriggerSel) return;
+    }
     
-    //Bool_t isTrack=1;
+    for(Int_t i=0;i<32;i++) {
+      Int_t bit=(1<<i);
+      if(inputHandler->IsEventSelected() & bit) htriggerbits[iBconf][0]->Fill(i);
+    }
+    if(inputHandler->IsEventSelected() & AliVEvent::kAny) htriggerbits[iBconf][0]->Fill(35);
+    if(inputHandler->IsEventSelected() & AliVEvent::kAnyINT) htriggerbits[iBconf][0]->Fill(36);
+    
+    htriggerbits[iBconf][1]->Fill(iTrigger);
+    
+    hCentrality[iBconf][1]->Fill(v0Centr);
+    hZvertex[iBconf][1]->Fill(zvtx);
     
     Int_t nTracks = fEvent->GetNumberOfTracks();
     
-    if(v0Centr>=fCentrality[0] && v0Centr<=fCentrality[1]) {//window cut centrality open
+    //----------------------loop on the TRACKS-----------------------------
+    for(Int_t iT = 0; iT < nTracks; iT++) { 
+      AliVTrack* track = (AliVTrack *) fEvent->GetTrack(iT);
       
-      hNevent[iBconf]->Fill(v0Centr);//analyzed events
-                  
-      for(Int_t iT = 0; iT < nTracks; iT++) { // loop on the tracks
-       AliVTrack* track = (AliVTrack *) fEvent->GetTrack(iT);
-       
-       if (!track){
-         continue;
-       }
-       
-       Bool_t trkFlag = 0;
-       trkFlag = ((AliAODTrack *) track)->TestFilterBit(FilterBit);
-       //TestFilterBit(16) -- Standard Cuts with very loose DCA: GetStandardITSTPCTrackCuts2011(kFALSE) && SetMaxDCAToVertexXY(2.4) && SetMaxDCAToVertexZ(3.2) && SetDCaToVertex2D(kTRUE)
-       //TestFilterBit(32) (STARDARD) -- Standard Cuts with very tight DCA cut ( 7sigma^primaries: 7*(0.0015+0.0050/pt^1.1) ) : GetStandardITSTPCTrackCuts2011(). 
-               
-               Int_t NTpcCls=track->GetTPCNcls();
-       if(NTpcCls>NminTPCcluster) kTPC=kTRUE;
-       else kTPC=kFALSE;
+      if (!track){
+       continue;
+      }
+      
+     //For the geometrical cuts
+      Double_t eta = track->Eta();
+      
+      Bool_t trkFlag = 0;
+      trkFlag = ((AliAODTrack *) track)->TestFilterBit(FilterBit);
+      //TestFilterBit(16) -- Standard Cuts with very loose DCA: GetStandardITSTPCTrackCuts2011(kFALSE) && SetMaxDCAToVertexXY(2.4) && SetMaxDCAToVertexZ(3.2) && SetDCaToVertex2D(kTRUE)
+      //TestFilterBit(32) (STARDARD) -- Standard Cuts with very tight DCA cut ( 7sigma^primaries: 7*(0.0015+0.0050/pt^1.1) ) : GetStandardITSTPCTrackCuts2011(). 
+      
+      //Cut on the Minumum Number of the TPC clusters
+      Bool_t isMinTpcCluster=kFALSE;
+      Int_t nTpcCluster=0;
+      nTpcCluster=track->GetTPCNcls();
+      if(nTpcCluster>NminTpcCluster) isMinTpcCluster=kTRUE;
 
-       if ((TMath::Abs(track->Eta()) > 0.8) || (track->Pt() < 0.2) || !trkFlag || !kTPC){
-         continue;
-       }
-       
-       Double_t b[2] = {-99., -99.};
-       Double_t bCov[3] = {-99., -99., -99.};
-       if (!track->PropagateToDCA(fEvent->GetPrimaryVertex(), fEvent->GetMagneticField(), 100., b, bCov))
-         continue;
-       
-//     Float_t Eta = TMath::Abs(track->Eta());
-       Float_t charge = (Float_t)track->Charge();
-       Float_t p = track->P();
-       Float_t pt = track->Pt();
-       Float_t dedx = track->GetTPCsignal();
-       Float_t tof  = track->GetTOFsignal()-fPIDResponse->GetTOFResponse().GetStartTime(p);
-       Float_t pTPC = track->GetTPCmomentum();
-       Float_t beta = 0.0;
-       Float_t M2 = 1000.0;
-       Float_t M = 1000.0;
-       Float_t Z2 = 1000.0;
-       Float_t DCAxy = b[0];
-       Float_t DCAz = b[1];
-
-       if(TMath::Abs(DCAz)>DCAzCUT)//CUT ON DCAz
-         continue;
-       
-       Bool_t kTpcPure;
-       kTpcPure = track->GetTPCsignal()>10;
-       if(kTpcPure==kFALSE) continue;
+      //-------------------------------------start TRACK CUTS (I): for (II) see below--------
+      if ((track->Pt() < 0.2) || (eta<EtaLimit[0]) || (eta>EtaLimit[1]) || !trkFlag || !isMinTpcCluster)
+       continue;
+      
+      //For the Tpc purity cut
+      Double_t dedx = track->GetTPCsignal();
+      if(dedx<10) continue;
 
-       
-       kTOF = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
+      Int_t nTrdSlices = track->GetNumberOfTRDslices();
+      if(nTrdSlices<2 && iTrdCut==1) continue; 
+      if(nTrdSlices>0 && iTrdCut==2) continue;
+      
+      //-------------------------------------end TRACK CUTS (I)----------------------------------
 
-               Float_t nsigmaTPC[9];
-       Float_t nsigmaTOF[9];
-       
-       Double_t expdedx[9];
-       
-       for(Int_t iS=0;iS < 9;iS++){ //TPC expected signal
-         expdedx[iS] = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, (AliPID::EParticleType) iS, AliTPCPIDResponse::kdEdxDefault, kTRUE);
+      //-------------------------------------Track info--------------------------------------
+      Double_t phi= track->Phi();
+      Double_t charge = (Double_t)track->Charge();
+      Double_t p = track->P();
+      Double_t pt = track->Pt();
+      Double_t tof  = track->GetTOFsignal()-fPIDResponse->GetTOFResponse().GetStartTime(p);
+      Double_t pTPC = track->GetTPCmomentum();
+      Double_t beta = 0.0;
+      Double_t M2 = 999.9;
+      Double_t Z2 = 999.9;
+      
+       //Vertex determination
+      Double_t b[2] = {-99., -99.};
+      Double_t bCov[3] = {-99., -99., -99.};
+      if (!track->PropagateToDCA(fEvent->GetPrimaryVertex(), fEvent->GetMagneticField(), 100., b, bCov))
+       continue;
+      
+      Double_t DCAxy = b[0];
+      Double_t DCAz = b[1];
+     
+      kTOF = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
+      
+      //-----------------------------TPC info------------------------------
+      Double_t nsigmaTPC[nPart];
+      Double_t expdedx[nPart];
+      
+      Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
+      Int_t FlagPid = 0;
+      
+      for(Int_t iS=0;iS<9;iS++){
+       nsigmaTPC[iS] = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType) iS);
+       //TPC identification:
+       if(TMath::Abs(nsigmaTPC[iS])<NsigmaTpcCut) {
+         FlagPid += ((Int_t)TMath::Power(2,iS));
        }
-               
-       for(Int_t iS=0;iS < 9;iS++){
-         nsigmaTPC[iS] = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType) iS);
-         fNsigmaTPC[iBconf][iS]->Fill(pt,nsigmaTPC[iS]);
-         hDeDxExp[iBconf][iS]->Fill(pTPC,expdedx[iS]);
+      }
+      //Correction of the momentum to the vertex for (anti)nuclei
+      Double_t pC[9];
+      for(Int_t iS=0;iS<9;iS++)        pC[iS]=p;
+      this->MomVertexCorrection(p,pC,eta,FlagPid);
+      
+      this->FillDCAdist(DCAxy,DCAz,charge,FlagPid,stdFlagPid,pC);
+      
+      //-------------------------------------start TRACK CUTS (II)-------------------------------------
+      //Cut on the DCAxy
+      Bool_t isDCAxyCut=kFALSE;
+      if(TMath::Abs(DCAxy)<DCAxyCut) isDCAxyCut=kTRUE;
+      
+      //Cut on the DCAz
+      Bool_t isDCAzCut=kFALSE;
+      if(TMath::Abs(DCAz)<DCAzCut) isDCAzCut=kTRUE;
+      
+      if (!isDCAxyCut || !isDCAzCut)
+       continue;
+          
+      //-------------------------------------end TRACK CUTS (II)----------------------------------
+      
+      hEta[iBconf]->Fill(eta);
+      hPhi[iBconf]->Fill(phi);
+      fEtaPhi[iBconf]->Fill(eta,phi);
+      hNTpcCluster[iBconf]->Fill(nTpcCluster);
+      hNTrdSlices[iBconf]->Fill(nTrdSlices);
+      
+      //More TPC info:
+      for(Int_t iS=0;iS<9;iS++){
+       expdedx[iS] = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, (AliPID::EParticleType) iS, AliTPCPIDResponse::kdEdxDefault, kTRUE);
+       hDeDxExp[iBconf][iS]->Fill(pTPC,expdedx[iS]);
+       nsigmaTPC[iS] = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType) iS);
+       //fNsigmaTpc[iBconf][iS]->Fill(pTPC,nsigmaTPC[iS]);
+       if(charge>0) {//positive particle
+         fNsigmaTpc[iBconf][iS]->Fill(pTPC,nsigmaTPC[iS]);
+         if(kTOF) fNsigmaTpc_kTOF[iBconf][iS]->Fill(p,nsigmaTPC[iS]);
        }
-       fdEdxVSp[iBconf][0]->Fill(pTPC,dedx);
-       if(charge>0) fdEdxVSp[iBconf][1]->Fill(pTPC,dedx);
-       else fdEdxVSp[iBconf][2]->Fill(pTPC,dedx);
-               
-       Float_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.877837,2.817402,1.408701,1.877837};
-       
-       Double_t exptimes[9]; // TOF expected times
+       else {//negative particle
+         fNsigmaTpc[iBconf][iS+nPart]->Fill(pTPC,nsigmaTPC[iS]);
+         if(kTOF) fNsigmaTpc_kTOF[iBconf][iS+nPart]->Fill(p,nsigmaTPC[iS]);
+       }
+       /*
+         if(TMath::Abs(nsigmaTPC[iS])<NsigmaTpcCut) {
+         FlagPid += ((Int_t)TMath::Power(2,iS));
+         }*/
+      }
+          
+      if(charge>0) fdEdxVSp[iBconf][0]->Fill(pTPC,dedx);
+      else fdEdxVSp[iBconf][1]->Fill(pTPC,dedx);
+
+      //-----------------------------TOF info------------------------------
+      
+      Double_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.875612859,2.808921005,1.404195741,1.863689620};
+
+      //----------------------------------------kTOF available-----------------------------
+           
+      if(kTOF) {
+       Double_t exptimes[9];
        track->GetIntegratedTimes(exptimes);
-       exptimes[5] = exptimes[0] / p * massOverZ[5] * TMath::Sqrt(1+p*p/massOverZ[5]/massOverZ[5]);
-       exptimes[6] = exptimes[0] / p * massOverZ[6] * TMath::Sqrt(1+p*p/massOverZ[6]/massOverZ[6]);
-       exptimes[7] = exptimes[0] / p * massOverZ[7] * TMath::Sqrt(1+p*p/massOverZ[7]/massOverZ[7]);
-       exptimes[8] = exptimes[0] / p * massOverZ[8] * TMath::Sqrt(1+p*p/massOverZ[8]/massOverZ[8]);
+       //Integrated times of the Nuclei:
+       for(Int_t iN=5;iN<9;iN++) {
+         exptimes[iN] = exptimes[4]*exptimes[4]*(massOverZ[iN]*massOverZ[iN]/p/p+1)/(massOverZ[4]*massOverZ[4]/p/p+1);
+         exptimes[iN] = TMath::Sqrt(exptimes[iN]);
+       }  
        
-       beta=exptimes[0];//expected times of the electron (it will be diveded for the T.o.f.)
+       beta=exptimes[0];
+       beta=beta/tof;//beta = L/tof/c = t_e/tof
        
-       if(kTOF) {
-         for(Int_t iS=0;iS < 9;iS++){
-           nsigmaTOF[iS] = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType) iS);
-           fNsigmaTOF[iBconf][iS]->Fill(pt,nsigmaTOF[iS]);
-           hBetaExp[iBconf][iS]->Fill(p,beta/exptimes[iS]);
+       Int_t FlagPidTof = 0;
+       Double_t NsigmaTofCut = 2.0;
+       
+       Double_t nsigmaTOF[9];
+       for(Int_t iS=0;iS<9;iS++){
+         nsigmaTOF[iS] = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType) iS);
+         fNsigmaTof[iBconf][iS]->Fill(pt,nsigmaTOF[iS]);
+         if(charge>0) {
+           hBetaExp[iBconf][iS]->Fill(p,exptimes[0]/exptimes[iS]);
+         }
+         else {
+           hBetaExp[iBconf][iS+nPart]->Fill(p,exptimes[0]/exptimes[iS]);
+         }
+
+         //TOF identification:
+         if(TMath::Abs(nsigmaTOF[iS])<NsigmaTofCut) {
+           FlagPidTof += ((Int_t)TMath::Power(2,iS));
          }
-         if(pt>0.9 && pt<1.0) hTOFSignalPion[iBconf]->Fill(tof-exptimes[2]);
-         beta=beta/tof;
-         fBetaTofVSp[iBconf]->Fill(p,beta);
        }
        
-       Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
-       Int_t FlagPid = 0;
-       Float_t binPt[nbin+1];
+       if(charge>0) fBetaTofVSp[iBconf][0]->Fill(p,beta);
+       else fBetaTofVSp[iBconf][1]->Fill(p,beta);
+               
+       this->GetMassFromPvertex(beta,p,M2);
+       this->GetZTpc(dedx,pTPC,M2,Z2);
+       
+       Double_t Mass2[9];
+       //-----------------------------M2 as a function of momentum to the primary vertex if iMtof==1---------------------------------
+       if(iMtof==1) this->GetMassFromPvertexCorrected(beta,pC,Mass2);
 
-       for(Int_t i=0;i<nbin+1;i++) {
-         binPt[i]=0.4+i*0.1;
-       }
+       if(iMtof==2) this->GetPmeanVsBetaGamma(exptimes,pC,FlagPid,FlagPidTof,charge);
        
-       Float_t binCutPt[10] = {0.3,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5};
-
-       Float_t binCutLargePt[6] = {0.0,1.0,2.0,3.0,4.0,5.0};
-
-       if(kTOF) {
-         
-         M2 = (p*p*(1-beta*beta))/(beta*beta);
-         
-         fM2vsP_NoTpcCut[iBconf][0]->Fill(M2,p);
-         if(TMath::Abs(DCAxy)<DCAxyCUT) fM2vsP_NoTpcCut_DCAxyCut[iBconf][0]->Fill(M2,p);
-         
-         if(M2>0.0) {
-           M=TMath::Sqrt(M2);
-           Z2 = TMath::Power(dedx/fPIDResponse->GetTPCResponse().GetExpectedSignal(pTPC*massOverZ[4]/M, AliPID::kProton),0.862);
-           fM2vsZ[iBconf][0]->Fill(charge*TMath::Sqrt(Z2),M2);
-           
-           for(Int_t i=0;i<9;i++) {
-             if(pt>binCutPt[i] && pt<binCutPt[i+1]){
-               fM2vsZ[iBconf][i+1]->Fill(charge*TMath::Sqrt(Z2),M2);
-               break;
-             }
-           }
-           for(Int_t i=0;i<5;i++) {
-             if(pt>binCutLargePt[i] && pt<binCutLargePt[i+1]){
-               fM2vsZ[iBconf][10+i]->Fill(charge*TMath::Sqrt(Z2),M2);
-               break;          
-             }       
-           }
-         }
-         
-         if(charge>0) {
-           fM2vsP_NoTpcCut[iBconf][1]->Fill(M2,p);
-           if(TMath::Abs(DCAxy)<DCAxyCUT) fM2vsP_NoTpcCut_DCAxyCut[iBconf][1]->Fill(M2,p);
-           for(Int_t iS=0;iS<9;iS++){
-             fNsigmaTPCvsP_kTOFtrue[iBconf][iS]->Fill(p,nsigmaTPC[iS]);
-           }
-         }
-         else {//else charge<0
-           fM2vsP_NoTpcCut[iBconf][2]->Fill(M2,p);
-           if(TMath::Abs(DCAxy)<DCAxyCUT) fM2vsP_NoTpcCut_DCAxyCut[iBconf][2]->Fill(M2,p);
-           
-           for(Int_t iS=0;iS < 9;iS++){
-             fNsigmaTPCvsP_kTOFtrue[iBconf][iS+9]->Fill(p,nsigmaTPC[iS]);
+       //-----------------------------M2 as a function of expected times---------------------------------
+       if(iMtof==2) this->GetMassFromExpTimes(beta,exptimes,Mass2);
+        
+       //-----------------------------M2 as a function of mean momentum calculated from expected time and extrapolated to the (anti)nuclei---------------------------------
+       if(iMtof>2) this->GetMassFromMeanMom(beta,exptimes,pC,eta,charge,Mass2,FlagPid,FlagPidTof);
+
+       //-------------------------------Squared Mass TH2 distributions-----------------------
+       if(charge>0) {
+         //without TPC
+         fM2vsP_NoTpcCut[iBconf][0][0]->Fill(M2,p);
+         //with TPC
+         for(Int_t iS=0;iS<9;iS++) {
+           M2=999.9;
+           M2=Mass2[iS];
+           //-----------------
+           if(FlagPid & stdFlagPid[iS]) {
+             fM2vsP[iBconf][0][iS]->Fill(M2,pC[iS]);
            }
          }
-         
+       }
+       else {//charge<0
+         //without TPC
+         fM2vsP_NoTpcCut[iBconf][0][1]->Fill(M2,p);
+         //with TPC
          for(Int_t iS=0;iS<9;iS++) {
-           if(TMath::Abs(nsigmaTPC[iS])<2.0) {
-             FlagPid += ((Int_t)TMath::Power(2,iS));
+           M2=999.9;
+           M2=Mass2[iS];
+           //-----------------
+           if(FlagPid & stdFlagPid[iS]) {
+             fM2vsP[iBconf][0][iS+nPart]->Fill(M2,pC[iS]);
            }
          }
+       }
+       
+       //------------------------------start Squared Mass TH1 distributions-------------------------
+       Double_t binP[nbin+1];
+       for(Int_t i=0;i<nbin+1;i++) {
+         binP[i]=0.4+i*0.1;
+       }
 
-         if(M2>0.0) {
-           for(Int_t iS=0;iS<9;iS++) {
-             if(FlagPid & stdFlagPid[iS]) {
-               fM2vsZwithTPC[iBconf][0]->Fill(charge*TMath::Sqrt(Z2),M2);
-               for(Int_t i=0;i<9;i++) {
-                 if(pt>binCutPt[i] && pt<binCutPt[i+1]) {
-                   fM2vsZwithTPC[iBconf][i+1]->Fill(charge*TMath::Sqrt(Z2),M2);
-                   break;
-                 }
-               }
-               for(Int_t i=0;i<5;i++) {
-                 if(pt>binCutLargePt[i] && pt<binCutLargePt[i+1]){
-                   fM2vsZwithTPC[iBconf][10+i]->Fill(charge*TMath::Sqrt(Z2),M2);
-                   break;              
-                 }           
+       if(charge>0) {
+         for(Int_t iS=0;iS<9;iS++) {
+           M2=999.9;
+           M2=Mass2[iS];
+           
+           if(FlagPid & stdFlagPid[iS]) {
+             for(Int_t j=0;j<nbin;j++) {
+               if(pC[iS]>binP[j] && pC[iS]<binP[j+1]) {
+                 hM2CutDCAxy[iBconf][iS][j]->Fill(M2);
+                 break;
                }
-             }
+             }//end loop on the p bins (j)
            }
-         }       
-
+         }//end loop on the particle species (iS)
+       }
+       else {//charge<0
          for(Int_t iS=0;iS<9;iS++) {
-           if(FlagPid & stdFlagPid[iS] || !kTPCcut) {
-             if(charge>0) {
-               fM2vsP[iBconf][iS]->Fill(M2,p);
-               if(TMath::Abs(DCAxy)<DCAxyCUT) {
-                 fM2vsP_DCAxyCut[iBconf][iS]->Fill(M2,p);
-               }
-               for(Int_t j=0;j<nbin;j++) {
-                 if(pt>binPt[j] && pt<binPt[j+1]) {
-                   hDCAxy[iBconf][iS][j]->Fill(DCAxy);
-                   hDCAxy[iBconf][iS][j]->Fill(-DCAxy);
-                   hDCAz[iBconf][iS][j]->Fill(DCAz);
-                   hDCAz[iBconf][iS][j]->Fill(-DCAz);
-                   if(TMath::Abs(DCAxy)<DCAxyCUT) {
-                     hM2CutDCAxy[iBconf][iS][j]->Fill(M2);
-                   }
-                   if(TMath::Abs(DCAxy+0.5)<DCAxyCUT) hM2CutGroundDCAxy[iBconf][iS][j]->Fill(M2);
-                   break;
-                 }
-               }
-             }
-             else {//if(charge<0)
-               fM2vsP[iBconf][iS+9]->Fill(M2,p);
-               if(TMath::Abs(DCAxy)<DCAxyCUT) {
-                 fM2vsP_DCAxyCut[iBconf][iS+9]->Fill(M2,p);
+           M2=999.9;
+           M2=Mass2[iS];
+                   
+           if(FlagPid & stdFlagPid[iS]) {
+             for(Int_t j=0;j<nbin;j++) {
+               if(pC[iS]>binP[j] && pC[iS]<binP[j+1]) {
+                 hM2CutDCAxy[iBconf][iS+nPart][j]->Fill(M2);
+                 break;
                }
-               for(Int_t j=0;j<nbin;j++) {
-                 if(pt>binPt[j] && pt<binPt[j+1]) {
-                   hDCAxy[iBconf][iS+9][j]->Fill(DCAxy);
-                   hDCAxy[iBconf][iS+9][j]->Fill(-DCAxy);
-                   hDCAz[iBconf][iS+9][j]->Fill(DCAz);
-                   hDCAz[iBconf][iS+9][j]->Fill(-DCAz);
-                   if(TMath::Abs(DCAxy)<DCAxyCUT) {
-                     hM2CutDCAxy[iBconf][iS+9][j]->Fill(M2);
-                   }
-                   if(TMath::Abs(DCAxy+0.5)<DCAxyCUT) hM2CutGroundDCAxy[iBconf][iS+9][j]->Fill(M2);
-                   break;
-                 }
-               }
-             }
+             }//end loop on the p bins (j)
            }
-         }
-       }// close (KTOF request)
+         }//end loop on the particle species (iS)
+       }
+               
+       //-------------------------------------------------M2/Z2 vs Z-------------------------
        
+
+       Double_t binCutPt[10] = {0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0};
+       Double_t Z=999.9;
+       if(Z2>0) Z=TMath::Sqrt(Z2);
        
-      } // end track loop
-      
-    } //window cut centrality close
-    
-  }  
-  
+       fM2vsZ[iBconf][0]->Fill(charge*TMath::Sqrt(Z2),M2);
+       for(Int_t i=1;i<10;i++) {
+         if(pt>binCutPt[i-1] && pt<binCutPt[i]){
+           fM2vsZ[iBconf][i]->Fill(charge*Z,M2);
+           break;
+         }
+       }
+       
+      }//end kTOF available
+    }//end track loop
+  }//end loop on the events
 }
 
 //_____________________________________________________________________________
@@ -796,3 +794,420 @@ void AliAnalysisNucleiMass::Terminate(Option_t *)
   // Terminate loop
   Printf("Terminate()");
 }
+//_____________________________________________________________________________
+void AliAnalysisNucleiMass::MomVertexCorrection(Double_t p, Double_t *pC, Double_t eta, Int_t FlagPid){
+
+  Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
+
+  for(Int_t iS=0;iS<9;iS++) {
+    if(FlagPid & stdFlagPid[iS]) {
+      if(iS==5) {
+       if(kPvtxCorr==1) pC[iS]=pC[iS]*fPvtxTrueVsReco[0]->Eval(pC[iS],TMath::Abs(eta));//for (bar)d
+       prPvtxTrueVsReco[iBconf][0]->Fill(p,pC[iS]/p);
+      }
+      else if(iS==6) {
+       if(kPvtxCorr==1) pC[iS]=pC[iS]*fPvtxTrueVsReco[1]->Eval(pC[iS],TMath::Abs(eta));//for (bar)t
+       prPvtxTrueVsReco[iBconf][1]->Fill(p,pC[iS]/p);
+      }
+      else if(iS==7) {
+       if(kPvtxCorr==1) pC[iS]=pC[iS]*fPvtxTrueVsReco[2]->Eval(pC[iS],TMath::Abs(eta));//for (bar)He3
+       prPvtxTrueVsReco[iBconf][2]->Fill(p,pC[iS]/p);
+      }
+      else if(iS==8) {
+       if(kPvtxCorr==1) pC[iS]=pC[iS]*fPvtxTrueVsReco[3]->Eval(pC[iS],TMath::Abs(eta));//for (bar)He3
+       prPvtxTrueVsReco[iBconf][3]->Fill(p,pC[iS]/p);
+      }
+    }
+  }
+  
+  return;
+  
+}
+//__________________________________________________________________________________________________
+void AliAnalysisNucleiMass::FillDCAdist(Double_t DCAxy, Double_t DCAz, Double_t charge, Int_t FlagPid, Int_t stdFlagPid[9], Double_t *pC){
+
+  Double_t binP[nbin+1];
+  for(Int_t i=0;i<nbin+1;i++) {
+    binP[i]=0.4+i*0.1;
+  }
+
+  if(charge>0) {
+    for(Int_t iS=0;iS<9;iS++) {
+      if(FlagPid & stdFlagPid[iS]) {
+       if(pC[iS]<1.5) {
+         h2DCAap[iBconf][iS]->Fill(DCAxy,DCAz);
+         h2DCAap[iBconf][iS]->Fill(-DCAxy,-DCAz);
+       }
+       for(Int_t j=0;j<nbin;j++) {
+         if(pC[iS]>binP[j] && pC[iS]<binP[j+1]) {
+           if(TMath::Abs(DCAz)<DCAzCut) {
+             hDCAxy[iBconf][iS][j]->Fill(DCAxy);
+             hDCAxy[iBconf][iS][j]->Fill(-DCAxy);
+           }
+           if(TMath::Abs(DCAxy)<DCAxyCut) {
+             hDCAz[iBconf][iS][j]->Fill(DCAz);
+             hDCAz[iBconf][iS][j]->Fill(-DCAz);
+           }
+           //h2DCA[iBconf][iS][j]->Fill(DCAxy,DCAz);
+           //h2DCA[iBconf][iS][j]->Fill(-DCAxy,-DCAz);
+           break;
+         }
+       }//end loop on the p bins (j)
+      }
+    }//end loop on the particle species (iS)
+  }
+  else {//charge<0
+    for(Int_t iS=0;iS<9;iS++) {
+      if(FlagPid & stdFlagPid[iS]) {
+       if(pC[iS]<1.5) {
+         h2DCAap[iBconf][iS+nPart]->Fill(DCAxy,DCAz);
+         h2DCAap[iBconf][iS+nPart]->Fill(-DCAxy,-DCAz);
+       }
+       for(Int_t j=0;j<nbin;j++) {
+         if(pC[iS]>binP[j] && pC[iS]<binP[j+1]) {
+           if(TMath::Abs(DCAz)<DCAzCut) {
+             hDCAxy[iBconf][iS+nPart][j]->Fill(DCAxy);
+             hDCAxy[iBconf][iS+nPart][j]->Fill(-DCAxy);
+           }
+           if(TMath::Abs(DCAxy)<DCAxyCut) {
+             hDCAz[iBconf][iS+nPart][j]->Fill(DCAz);
+             hDCAz[iBconf][iS+nPart][j]->Fill(-DCAz);
+           }
+           //h2DCA[iBconf][iS+nPart][j]->Fill(DCAxy,DCAz);
+           //h2DCA[iBconf][iS+nPart][j]->Fill(-DCAxy,-DCAz);
+           break;
+         }
+       }//end loop on the p bins (j)
+      }
+    }//end loop on the particle species (iS)
+  } 
+  
+  return;
+}
+//_____________________________________________________________________________
+void AliAnalysisNucleiMass::GetMassFromPvertex(Double_t beta, Double_t p, Double_t &M2) {
+  
+  M2 = p*p*(1-beta*beta)/(beta*beta);
+
+  return;
+  
+}
+//_________________________________________________________________________________________________________________________
+void AliAnalysisNucleiMass::GetZTpc(Double_t dedx, Double_t pTPC, Double_t M2, Double_t &Z2) {
+
+  //z^2_tpc = dedx^{Tpc} / dedx^{exp,Tof}_{z=1}
+  
+  Z2=999.9;
+  
+  Double_t M=999.9;
+  Double_t pTPC_pr=999.9;//rescaling of the pTPC for the proton
+  Double_t expdedx_Tof=999.9;
+  
+  if(M2>0) {
+    M=TMath::Sqrt(M2);
+    pTPC_pr=pTPC*0.938272/M;
+    expdedx_Tof=fPIDResponse->GetTPCResponse().GetExpectedSignal(pTPC_pr,AliPID::kProton);
+    if((dedx/expdedx_Tof)<0) return;
+    Z2=TMath::Power(dedx/expdedx_Tof,0.862);
+  }
+  
+  return;
+}
+//_________________________________________________________________________________________________________________________
+void AliAnalysisNucleiMass::GetMassFromPvertexCorrected(Double_t beta, Double_t *pC, Double_t *Mass2) {
+  
+  for(Int_t iS=0;iS<9;iS++) Mass2[iS] = pC[iS]*pC[iS]*(1-beta*beta)/(beta*beta);
+
+  return;
+}  
+//____________________________________________________________________________________________________________
+void AliAnalysisNucleiMass::GetMassFromExpTimes(Double_t beta, Double_t *IntTimes, Double_t *Mass2) {
+  // m = p_exp/beta/gamma where p_exp = mPDG*beta_exp*gamma_exp; beta_exp = L/t_exp/c = t_e/t_exp ; beta=L/tof/c = t_e/tof
+  // In this way m_tof = mPDG only if tof=t_exp
+  
+  Double_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.875612859,2.808921005,1.404195741,1.863689620};
+    
+  Double_t beta2Exp[9];
+  Double_t p2Exp[9];
+  
+  //Double_t pExp[9];
+  
+  for(Int_t iS=0;iS<9;iS++) {
+    beta2Exp[iS]=IntTimes[0]/IntTimes[iS];//beta = L/tof*c = t_e/tof
+    beta2Exp[iS]=beta2Exp[iS]*beta2Exp[iS];
+    if((1-beta2Exp[iS])==0) {
+      Mass2[iS]=999.9;
+      continue;
+    }
+    p2Exp[iS]=massOverZ[iS]*massOverZ[iS]*beta2Exp[iS]/(1-beta2Exp[iS]);
+    
+    //--------------------for MC corrections
+    if(p2Exp[iS]<0) {
+      Mass2[iS]=999.9;
+      continue;
+    }
+    //pExp[iS]=TMath::Sqrt(p2Exp[iS]);
+    
+    //------------
+    Mass2[iS]=p2Exp[iS]*(1-beta*beta)/(beta*beta);
+  }//end loop on the particle species
+  
+  return;
+}
+//____________________________________________________________________________________________________________
+void AliAnalysisNucleiMass::GetPmeanVsBetaGamma(Double_t *IntTimes, Double_t *pVtx, Int_t FlagPid, Int_t FlagPidTof, Double_t charge) {
+  // m = p_exp/beta/gamma where p_exp = mPDG*beta_exp*gamma_exp; beta_exp = L/t_exp/c = t_e/t_exp ; beta=L/tof/c = t_e/tof
+  // In this way m_tof = mPDG only if tof=t_exp
+  
+  Double_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.875612859,2.808921005,1.404195741,1.863689620};
+    
+  Double_t beta2Exp[9];
+  Double_t p2Exp[9];
+  
+  Double_t pExp[9];
+  
+  Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
+
+  for(Int_t iS=0;iS<9;iS++) {
+    beta2Exp[iS]=IntTimes[0]/IntTimes[iS];//beta = L/tof*c = t_e/tof
+    beta2Exp[iS]=beta2Exp[iS]*beta2Exp[iS];
+    if((1-beta2Exp[iS])==0) {
+      continue;
+    }
+    p2Exp[iS]=massOverZ[iS]*massOverZ[iS]*beta2Exp[iS]/(1-beta2Exp[iS]);
+    
+    if(p2Exp[iS]<0) {
+      continue;
+    }
+    pExp[iS]=TMath::Sqrt(p2Exp[iS]);
+       
+    if((FlagPid & stdFlagPid[iS]) && (FlagPidTof & stdFlagPid[iS])) {
+      if(charge>0){
+       fPmeanVsBetaGamma[iBconf][iS]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
+       prPmeanVsBetaGamma[iBconf][iS]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
+      }
+      else {
+       fPmeanVsBetaGamma[iBconf][iS+nPart]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
+       prPmeanVsBetaGamma[iBconf][iS+nPart]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
+      }
+    }
+  }//end loop on the particle species
+  
+  return;
+  
+}
+//____________________________________________________________________________________________________________
+void AliAnalysisNucleiMass::GetMassFromMeanMom(Double_t beta, Double_t *IntTimes, Double_t *pVtx, Double_t eta, Double_t charge, Double_t *Mass2, Int_t FlagPid, Int_t FlagPidTof) {//Double_t *Mass2, Int_t iCorr
+  // m = p_exp/beta/gamma where p_exp = mPDG*beta_exp*gamma_exp; beta_exp = L/t_exp/c = t_e/t_exp ; beta=L/tof/c = t_e/tof
+  // In this way m_tof = mPDG only if tof=t_exp
+  
+  Double_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.875612859,2.808921005,1.404195741,1.863689620};
+  
+  Double_t beta2Exp[9];
+  Double_t p2Exp[9];
+  
+  Double_t pExp[9];
+  
+  Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
+  
+  for(Int_t iS=0;iS<9;iS++) {
+    if(iS>1) {
+      p2Exp[iS]=pVtx[iS]*fPmeanVsBGcorr[iS-2]->Eval(pVtx[iS]/massOverZ[iS],TMath::Abs(eta));
+      p2Exp[iS]*=p2Exp[iS];
+    }
+    else {
+      beta2Exp[iS]=IntTimes[0]/IntTimes[iS];//beta = L/tof*c = t_e/tof
+      beta2Exp[iS]=beta2Exp[iS]*beta2Exp[iS];
+      if((1-beta2Exp[iS])==0) {
+       Mass2[iS]=999.9;
+       continue;
+      }
+      p2Exp[iS]=massOverZ[iS]*massOverZ[iS]*beta2Exp[iS]/(1-beta2Exp[iS]);
+    }
+    
+    if(p2Exp[iS]<0) {
+      Mass2[iS]=999.9;
+      continue;
+    }
+    pExp[iS]=TMath::Sqrt(p2Exp[iS]);
+    
+    //------------
+    Mass2[iS]=p2Exp[iS]*(1-beta*beta)/(beta*beta);
+    
+    //-----------
+    
+    if(iS>1) {
+      if((FlagPid & stdFlagPid[iS]) && (FlagPidTof & stdFlagPid[iS])) {
+       if(charge>0) {
+         prPmeanVsBGcorr[iBconf][iS-2]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
+       }
+       else if(charge<0) {
+         prPmeanVsBGcorr[iBconf][iS-2+7]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
+       }
+      }
+    }
+  }//end loop on the particle species
+    
+  return;
+  
+}
+//________________________________________________________________________________________
+void AliAnalysisNucleiMass::SetPvtxCorrections(){
+  //for (bar)d
+  fPvtxTrueVsReco[0]=new TF2("fcorr_d","([0]*TMath::Power(x,[1])+[2])*(TMath::Power((TMath::Exp([3]*x)+[4]),[5]*TMath::Power(y,[6])));p_{reco};|#eta|;p_{true}/p_{reco}",0.0001,100,0,1);//for (bar)d
+  fPvtxTrueVsReco[0]->SetParameter(0,0.031263);
+  fPvtxTrueVsReco[0]->SetParameter(1,-3.276770);
+  fPvtxTrueVsReco[0]->SetParameter(2,1.000113);
+  fPvtxTrueVsReco[0]->SetParameter(3,-5.195875);
+  fPvtxTrueVsReco[0]->SetParameter(4,1.000674);
+  fPvtxTrueVsReco[0]->SetParameter(5,2.870503);
+  fPvtxTrueVsReco[0]->SetParameter(6,3.777729);
+  
+  //for (bar)t
+  fPvtxTrueVsReco[1]=new TF2("fcorr_t","([0]*TMath::Power(x,[1])+[2])+[3]*y;p_{reco};|#eta|;p_{true}/p_{reco}",0.0001,100,0,1);//for (bar)He3
+  fPvtxTrueVsReco[1]->SetParameter(0,8.79761e-02);
+  fPvtxTrueVsReco[1]->SetParameter(1,-3.23189e+00);
+  fPvtxTrueVsReco[1]->SetParameter(2,9.99578e-01);
+  fPvtxTrueVsReco[1]->SetParameter(3,0.0);
+  
+  //for (bar)He3
+  fPvtxTrueVsReco[2]=new TF2("fcorr_He","([0]*TMath::Power(x,[1])+[2])*(TMath::Power((TMath::Exp([3]*x)+[4]),[5]*TMath::Power(y,[6])));p_{reco};|#eta|;p_{true}/p_{reco}",0.0001,100,0,1);//for (bar)He3
+  fPvtxTrueVsReco[2]->SetParameter(0,0.037986);
+  fPvtxTrueVsReco[2]->SetParameter(1,-2.707620);
+  fPvtxTrueVsReco[2]->SetParameter(2,1.000742);
+  fPvtxTrueVsReco[2]->SetParameter(3,-4.934743);
+  fPvtxTrueVsReco[2]->SetParameter(4,1.001640);
+  fPvtxTrueVsReco[2]->SetParameter(5,2.744372);
+  fPvtxTrueVsReco[2]->SetParameter(6,3.528561);
+  
+  //for (bar)He4
+  fPvtxTrueVsReco[3]=new TF2("fcorr_He4","([0]*TMath::Power(x,[1])+[2])+[3]*y;p_{reco};|#eta|;p_{true}/p_{reco}",0.0001,100,0,1);//for (bar)He3
+  fPvtxTrueVsReco[3]->SetParameter(0,7.08785e-02);
+  fPvtxTrueVsReco[3]->SetParameter(1,-2.87201e+00);
+  fPvtxTrueVsReco[3]->SetParameter(2,1.00070e+00);
+  fPvtxTrueVsReco[3]->SetParameter(3,0.0);
+  
+  for(Int_t i=0;i<4;i++) {
+    fPvtxTrueVsReco[i]->SetNpx(fPvtxTrueVsReco[i]->GetNpx()*10.0);
+  }
+}
+//________________________________________________________________________________________
+void AliAnalysisNucleiMass::SetPmeanCorrections(){
+  
+  Char_t nameTemp[14][200];
+  snprintf(nameTemp[0],200,"#pi^{+}");
+  snprintf(nameTemp[1],200,"K^{+}");
+  snprintf(nameTemp[2],200,"p");
+  snprintf(nameTemp[3],200,"d");
+  snprintf(nameTemp[4],200,"t");
+  snprintf(nameTemp[5],200,"He3");
+  snprintf(nameTemp[6],200,"He4");
+  snprintf(nameTemp[7],200,"#pi^{-}");
+  snprintf(nameTemp[8],200,"K^{-}");
+  snprintf(nameTemp[9],200,"#bar{p}");
+  snprintf(nameTemp[10],200,"#bar{d}");
+  snprintf(nameTemp[11],200,"#bar{t}");
+  snprintf(nameTemp[12],200,"#bar{He3}");
+  snprintf(nameTemp[13],200,"#bar{He4}");
+  
+  Char_t name_fPmeanVsBGcorr[14][200];
+  for(Int_t i=0;i<14;i++) {
+    snprintf(name_fPmeanVsBGcorr[i],200,"fPmeanVsBGcorr_%s",nameTemp[i]);
+  }
+
+  //Pions
+  fPmeanVsBGcorr[0]=new TF2(name_fPmeanVsBGcorr[0],"(x>[5])*([2]-[0]*TMath::Power(x,[1]))*([3]+[4]*y*y)+(x<=[5])*[6]",0.0001,100,0,0.8);
+  fPmeanVsBGcorr[0]->SetParameter(0,-0.179607);
+  fPmeanVsBGcorr[0]->SetParameter(1,-0.384809);
+  fPmeanVsBGcorr[0]->SetParameter(2,0.885534);
+  fPmeanVsBGcorr[0]->SetParameter(3,0.992710);
+  fPmeanVsBGcorr[0]->SetParameter(4,0.011390);
+  fPmeanVsBGcorr[0]->SetParameter(5,3.231000);
+  fPmeanVsBGcorr[0]->SetParameter(6,0.999900);
+  //Kaons
+  fPmeanVsBGcorr[1]=new TF2(name_fPmeanVsBGcorr[1],"(x>[8])*([2]-[0]*TMath::Power(x,[1]))*TMath::Power([3]+[4]*TMath::Exp([5]*x),[6]+[7]*y*y)+(x<=[8])*[9]",0.0001,20,0,0.8);
+  fPmeanVsBGcorr[1]->SetParameter(0,0.033500);
+  fPmeanVsBGcorr[1]->SetParameter(1,-2.461673);
+  fPmeanVsBGcorr[1]->SetParameter(2,0.996501);
+  fPmeanVsBGcorr[1]->SetParameter(3,1.000000);
+  fPmeanVsBGcorr[1]->SetParameter(4,0.089715);
+  fPmeanVsBGcorr[1]->SetParameter(5,-2.473531);
+  fPmeanVsBGcorr[1]->SetParameter(6,1.000000);
+  fPmeanVsBGcorr[1]->SetParameter(7,-1.562500);
+  fPmeanVsBGcorr[1]->SetParameter(8,0.253000);
+  fPmeanVsBGcorr[1]->SetParameter(9,0.009387);
+  
+  //Protons
+  fPmeanVsBGcorr[2]=new TF2(name_fPmeanVsBGcorr[2],"(x>[8])*([2]-[0]*TMath::Power(x,[1]))*TMath::Power([3]+[4]*TMath::Exp([5]*x),[6]+[7]*y*y)+(x<=[8])*[9]",0.0001,20,0,0.8);
+  fPmeanVsBGcorr[2]->SetParameter(0,0.015081);
+  fPmeanVsBGcorr[2]->SetParameter(1,-2.927557);
+  fPmeanVsBGcorr[2]->SetParameter(2,0.997904);
+  fPmeanVsBGcorr[2]->SetParameter(3,1.000000);
+  fPmeanVsBGcorr[2]->SetParameter(4,0.102697);
+  fPmeanVsBGcorr[2]->SetParameter(5,-3.399528);
+  fPmeanVsBGcorr[2]->SetParameter(6,1.000000);
+  fPmeanVsBGcorr[2]->SetParameter(7,-1.562500);
+  fPmeanVsBGcorr[2]->SetParameter(8,0.239000);
+  fPmeanVsBGcorr[2]->SetParameter(9,0.002054);
+
+  //Deuterons
+  fPmeanVsBGcorr[3]=new TF2(name_fPmeanVsBGcorr[3],"(x>[8])*([2]-[0]*TMath::Power(x,[1]))*TMath::Power([3]+[4]*TMath::Exp([5]*x),[6]+[7]*y*y)+(x<=[8])*[9]",0.0001,20,0,0.8);
+  fPmeanVsBGcorr[3]->SetParameter(0,0.008672);
+  fPmeanVsBGcorr[3]->SetParameter(1,-2.712343);
+  fPmeanVsBGcorr[3]->SetParameter(2,0.997639);
+  fPmeanVsBGcorr[3]->SetParameter(3,1.000000);
+  fPmeanVsBGcorr[3]->SetParameter(4,0.039627);
+  fPmeanVsBGcorr[3]->SetParameter(5,-2.768122);
+  fPmeanVsBGcorr[3]->SetParameter(6,1.000000);
+  fPmeanVsBGcorr[3]->SetParameter(7,-1.562500);
+  fPmeanVsBGcorr[3]->SetParameter(8,0.174000);
+  fPmeanVsBGcorr[3]->SetParameter(9,0.002189);
+
+  //Triton
+  fPmeanVsBGcorr[4]=new TF2(name_fPmeanVsBGcorr[4],"(x>[4])*([2]-[0]*TMath::Power(x,[1])+[3]*y)+(x<=[4])*[5]",0.0001,20,0,0.8);
+  fPmeanVsBGcorr[4]->SetParameter(0,6.79641e-03);
+  fPmeanVsBGcorr[4]->SetParameter(1,-1.92801e+00);
+  fPmeanVsBGcorr[4]->SetParameter(2,1.000000);
+  fPmeanVsBGcorr[4]->SetParameter(3,0.0);
+  fPmeanVsBGcorr[4]->SetParameter(4,0.076);
+  fPmeanVsBGcorr[4]->SetParameter(5,2.25779e-02);
+
+  //Helium-3
+  fPmeanVsBGcorr[5]=new TF2(name_fPmeanVsBGcorr[5],"(x>[8])*([2]-[0]*TMath::Power(x,[1]))*TMath::Power([3]+[4]*TMath::Exp([5]*x),[6]+[7]*y*y)+(x<=[8])*[9]",0.0001,20,0,0.8);
+  fPmeanVsBGcorr[5]->SetParameter(0,0.024339);
+  fPmeanVsBGcorr[5]->SetParameter(1,-2.922613);
+  fPmeanVsBGcorr[5]->SetParameter(2,0.993761);
+  fPmeanVsBGcorr[5]->SetParameter(3,1.000000);
+  fPmeanVsBGcorr[5]->SetParameter(4,1.087549);
+  fPmeanVsBGcorr[5]->SetParameter(5,-6.216154);
+  fPmeanVsBGcorr[5]->SetParameter(6,1.000000);
+  fPmeanVsBGcorr[5]->SetParameter(7,-1.562500);
+  fPmeanVsBGcorr[5]->SetParameter(8,0.282000);
+  fPmeanVsBGcorr[5]->SetParameter(9,0.009711);
+
+  //Helium-4
+  fPmeanVsBGcorr[6]=new TF2(name_fPmeanVsBGcorr[6],"(x>[4])*([2]-[0]*TMath::Power(x,[1])+[3]*y)+(x<=[4])*[5]",0.0001,20,0,0.8);
+  fPmeanVsBGcorr[6]->SetParameter(0,2.34185e-02);
+  fPmeanVsBGcorr[6]->SetParameter(1,-2.31200e+00);
+  fPmeanVsBGcorr[6]->SetParameter(2,1.000000);
+  fPmeanVsBGcorr[6]->SetParameter(3,0.0);
+  fPmeanVsBGcorr[6]->SetParameter(4,0.198);
+  fPmeanVsBGcorr[6]->SetParameter(5,9.9226e-03);
+
+  for(Int_t i=7;i<14;i++) {
+    fPmeanVsBGcorr[i]=(TF2 *)fPmeanVsBGcorr[i-7]->Clone();
+    fPmeanVsBGcorr[i]->SetName(name_fPmeanVsBGcorr[i]);
+  }
+    
+  for(Int_t i=0;i<14;i++) {
+    fPmeanVsBGcorr[i]->SetNpx(fPmeanVsBGcorr[i]->GetNpx()*100.0);
+  }
+
+  return;
+  
+}
+