PWGLF/SPECTRA/Nuclei/masses/AliAnalysisNucleiMass.cxx

   1 #include "AliAnalysisNucleiMass.h"
   2
   3 // ROOT includes
   4 #include <TMath.h>
   5 #include "TChain.h"
   6
   7 // AliRoot includes
   8 #include "AliInputEventHandler.h"
   9 #include "AliAODEvent.h"
  10 #include "AliESDEvent.h"
  11 #include "AliVEvent.h"
  12 #include "AliAODTrack.h"
  13 #include "AliAODPid.h"
  14 #include "AliCentrality.h"
  15 #include "TH2F.h"
  16 #include "TH2D.h"
  17 #include "TH1F.h"
  18 #include "TF1.h"
  19 #include "TF2.h"
  20 #include "TGraph.h"
  21 #include "TProfile.h"
  22 #include "AliESDtrackCuts.h"
  23 #include "AliAnalysisManager.h"
  24 #include "TFile.h"
  25
  26 ClassImp(AliAnalysisNucleiMass)
  27
  28 //_____________________________________________________________________________
  29 AliAnalysisNucleiMass::AliAnalysisNucleiMass():
  30   AliAnalysisTaskSE(),
  31   Centrality(),
  32   FilterBit(16),
  33   EtaLimit(),
  34   DCAxyCut(0.1),
  35   DCAzCut(1000.),
  36   NsigmaTpcCut(2.0),
  37   NminTpcCluster(0),
  38   iTrdCut(0),
  39   kSignalCheck(1),
  40   iMtof(1),
  41   kPvtxCorr(1),
  42   iBconf(0),
  43   kTOF(0),
  44 //iTriggerSel(-99),
  45   fAOD(NULL),
  46   fESD(NULL),
  47   fEvent(NULL),
  48   fPIDResponse(NULL)
  49 {
  50   Centrality[0]=0.0;
  51   Centrality[1]=100.0;
  52
  53   EtaLimit[0]=-99.0;
  54   EtaLimit[1]=99.0;
  55
  56   fList[0]=new TList();
  57   fList[0]->SetName("results");
  58
  59   fList[1]=new TList();
  60   fList[1]->SetName("results2");
  61 }
  62 //______________________________________________________________________________
  63 AliAnalysisNucleiMass::AliAnalysisNucleiMass(const char *name):
  64   AliAnalysisTaskSE(name),
  65   Centrality(),
  66   FilterBit(16),
  67   EtaLimit(),
  68   DCAxyCut(0.1),
  69   DCAzCut(1000.),
  70   NsigmaTpcCut(2.0),
  71   NminTpcCluster(0),
  72   iTrdCut(0),
  73   kSignalCheck(1),
  74   iMtof(1),
  75   kPvtxCorr(1),
  76   iBconf(0),
  77   kTOF(0),
  78   //iTriggerSel(-99),
  79   fAOD(NULL),
  80   fESD(NULL),
  81   fEvent(NULL),
  82   fPIDResponse(NULL)
  83 {
  84
  85   Centrality[0]=0.0;
  86   Centrality[1]=100.0;
  87
  88   EtaLimit[0]=-99.0;
  89   EtaLimit[1]=99.0;
  90
  91   fList[0]=new TList();
  92   DefineOutput(1, TList::Class());
  93   fList[0]->SetName("results");
  94
  95   fList[1]=new TList();
  96   DefineOutput(2, TList::Class());
  97   fList[1]->SetName("results2");
  98 }
  99 //_____________________________________________________________________________
 100 AliAnalysisNucleiMass::~AliAnalysisNucleiMass()
 101 {
 102   if(fList[0]) delete fList[0];
 103   if(fList[1]) delete fList[1];
 104 }
 105 //______________________________________________________________________________
 106 void AliAnalysisNucleiMass::UserCreateOutputObjects()
 107 {
 108   Char_t namePart[nPart][30];
 109   snprintf(namePart[0],30,"e");
 110   snprintf(namePart[1],30,"#mu");
 111   snprintf(namePart[2],30,"#pi");
 112   snprintf(namePart[3],30,"K");
 113   snprintf(namePart[4],30,"p");
 114   snprintf(namePart[5],30,"d");
 115   snprintf(namePart[6],30,"t");
 116   snprintf(namePart[7],30,"He3");
 117   snprintf(namePart[8],30,"He4");
 118
 119   Char_t name[nSpec][30];
 120   snprintf(name[0],20,"e^{+}");
 121   snprintf(name[1],20,"#mu^{+}");
 122   snprintf(name[2],20,"#pi^{+}");
 123   snprintf(name[3],20,"K^{+}");
 124   snprintf(name[4],20,"p");
 125   snprintf(name[5],20,"d");
 126   snprintf(name[6],20,"t");
 127   snprintf(name[7],20,"He3");
 128   snprintf(name[8],20,"He4");
 129   snprintf(name[9],20,"e^{-}");
 130   snprintf(name[10],20,"#mu^{-}");
 131   snprintf(name[11],20,"#pi^{-}");
 132   snprintf(name[12],20,"K^{-}");
 133   snprintf(name[13],20,"#bar{p}");
 134   snprintf(name[14],20,"#bar{d}");
 135   snprintf(name[15],20,"#bar{t}");
 136   snprintf(name[16],20,"#bar{He3}");
 137   snprintf(name[17],20,"#bar{He4}");
 138
 139   Double_t binP[nbin+1];
 140   for(Int_t i=0;i<nbin+1;i++) {
 141     binP[i]=0.4+0.1*i;
 142   }
 143
 144   Char_t name_nbin[nbin][200];
 145   for(Int_t j=0;j<nbin;j++) {
 146     snprintf(name_nbin[j],200,"%.1f<P<%.1f",binP[j],binP[j+1]);
 147   }
 148
 149   for(Int_t iB=0;iB<nBconf;iB++) {
 150
 151     htemp[iB] = new TH1F("htemp","htemp (avoid the problem with the empty list...);B field",20,-10,10);
 152
 153     //htriggerbits[iB] = new TH1I("htriggerbits","htriggerbits; bits",10,-5,5);
 154     htriggerbits[iB][0] = new TH1I("htriggerbits_0","trigger mask; bits",45,-5,40);
 155     htriggerbits[iB][1] = new TH1I("htriggerbits_1","trigger bits (exclusive); bits",45,-5,40);
 156
 157     hCentrality[iB][0] = new TH1F("hCentrality_Selected","Centrality (selected events);centrality(%)",20,0,100);//20,0,100
 158     hCentrality[iB][1] = new TH1F("hCentrality_Analyzed","Centrality (analyzed events);centrality (%)",20,0,100);//20,0,100
 159
 160     hZvertex[iB][0] = new TH1F("hZvertex_Selected","Vertex distribution of selected events;z vertex (cm)",240,-30,30);
 161     hZvertex[iB][1] = new TH1F("hZvertex_Analyzed","Vertex distribution of analyzed events;z vertex (cm)",240,-30,30);
 162
 163     hEta[iB] = new TH1F("hEta_Analyzed","|#eta| distribution after the track cuts;#eta",200,-1.0,1.0);
 164
 165     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)
 166
 167     Int_t hbins[2];
 168     if(kSignalCheck!=0) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=100; hbins[1]=90;} to reduce RAM consuming (toram)
 169     else {hbins[0]=1; hbins[1]=1;}
 170     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);
 171
 172     hNTpcCluster[iB] = new TH1F("hNTpcCluster","Number of the TPC clusters after the track cuts;n_{cl}^{TPC}",300,0,300);
 173
 174     hNTrdSlices[iB] = new TH1F("hNTrdSlices","Number of the TRD slices after the track cuts;n_{slices}^{TRD}",40,0,40);
 175
 176     if(kSignalCheck==1) {hbins[0]=500; hbins[1]=2000;}
 177     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 178     else if(kSignalCheck==2) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=100; hbins[1]=500;} toram
 179     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);
 180     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);
 181
 182     Char_t name_hDeDxExp[nPart][200];
 183     Char_t title_hDeDxExp[nPart][200];
 184     for(Int_t i=0;i<nPart;i++) {
 185       snprintf(name_hDeDxExp[i],200,"hDeDxExp_%s",namePart[i]);
 186       snprintf(title_hDeDxExp[i],200,"Expected dE/dx of %s in the TPC;p/|z| (GeV/c);dE/dx_{TPC} (a.u.)",namePart[i]);
 187       hDeDxExp[iB][i] = new TProfile(name_hDeDxExp[i],title_hDeDxExp[i],1,0,5,0,1,"");//,500,0,5,0,1000,""); toram
 188     }
 189
 190     Char_t name_fNsigmaTpc[nSpec][200];
 191     Char_t title_fNsigmaTpc[nSpec][200];
 192     if(kSignalCheck==1) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=100; hbins[1]=100;} toram
 193     else {hbins[0]=100; hbins[1]=100;}//temp!
 194     for(Int_t i=0;i<nSpec;i++) {
 195       snprintf(name_fNsigmaTpc[i],200,"NsigmaTpc_%s",name[i]);
 196       snprintf(title_fNsigmaTpc[i],200,"NsigmaTpc_%s;p_{TPC}/|z| (GeV/c);n_{#sigma_{TPC}}^{%s}",name[i],name[i]);
 197       fNsigmaTpc[iB][i] = new TH2F(name_fNsigmaTpc[i],title_fNsigmaTpc[i],hbins[0],0,5,hbins[1],-5,5);
 198     }
 199
 200     if(kSignalCheck>0) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=100; hbins[1]=100;} toram
 201     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 202     Char_t name_fNsigmaTpc_kTOF[nSpec][200];
 203     Char_t title_fNsigmaTpc_kTOF[nSpec][200];
 204     for(Int_t i=0;i<nSpec;i++) {
 205       snprintf(name_fNsigmaTpc_kTOF[i],200,"NsigmaTpc_%s_kTOF",name[i]);
 206       snprintf(title_fNsigmaTpc_kTOF[i],200,"NsigmaTpc_kTOF_%s;p/|z| (GeV/c);n_{#sigma_{TPC}}^{%s}",name[i],name[i]);
 207       fNsigmaTpc_kTOF[iB][i] = new TH2F(name_fNsigmaTpc_kTOF[i],title_fNsigmaTpc_kTOF[i],hbins[0],0,5,hbins[1],-5,5);
 208     }
 209
 210     if(kSignalCheck==1) {hbins[0]=1000; hbins[1]=1300;}
 211     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 212     else if(kSignalCheck==2) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=100; hbins[1]=260;}
 213     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);
 214     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);
 215
 216     Char_t name_hBetaExp[nPart][200];
 217     Char_t title_hBetaExp[nPart][200];
 218     for(Int_t i=0;i<nPart;i++) {
 219       snprintf(name_hBetaExp[i],200,"hBetaTofVsP_Exp_%s",namePart[i]);
 220       snprintf(title_hBetaExp[i],200,"Expected #beta_{TOF} vs p/|z| of %s;p/|z| (GeV/c); #beta_{TOF}",namePart[i]);
 221       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
 222     }
 223
 224     Char_t name_fNsigmaTof[nPart][200];
 225     Char_t title_fNsigmaTof[nPart][200];
 226     if(kSignalCheck==1) {hbins[0]=100; hbins[1]=100;}
 227     else {hbins[0]=1; hbins[1]=1;}
 228     for(Int_t i=0;i<nPart;i++) {
 229       snprintf(name_fNsigmaTof[i],200,"NsigmaTof_%s",namePart[i]);
 230       snprintf(title_fNsigmaTof[i],200,"NsigmaTof_%s;p_{T}/|z| (GeV/c);n_{#sigma_{TOF}}^{%s}",namePart[i],namePart[i]);
 231       fNsigmaTof[iB][i] = new TH2F(name_fNsigmaTof[i],title_fNsigmaTof[i],hbins[0],0,5,hbins[1],-5,5);
 232     }
 233
 234     if(kSignalCheck==1) {hbins[0]=8000; hbins[1]=100;}
 235     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 236     else if(kSignalCheck==2) {hbins[0]=1; hbins[1]=1;}// {hbins[0]=1000; hbins[1]=100;} toram
 237     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);
 238     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);
 239
 240     Char_t name_fM2vsP[1][18][300];
 241     Char_t title_fM2vsP[1][18][300];
 242
 243     for(Int_t i=0;i<nSpec;i++) {
 244       snprintf(name_fM2vsP[0][i],300,"fM2vsPc_%s",name[i]);
 245       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]);
 246
 247       if(kSignalCheck==1) {hbins[0]=8000; hbins[1]=100;}
 248       else {hbins[0]=1; hbins[1]=1;}
 249       fM2vsP[iB][0][i] = new TH2F(name_fM2vsP[0][i],title_fM2vsP[0][i],hbins[0],0,10,hbins[1],0,5);
 250     }
 251
 252     if(kSignalCheck==1) {hbins[0]=4000; hbins[1]=1000;}
 253     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 254     else if(kSignalCheck==2) {hbins[0]=1; hbins[1]=1;}//{hbins[0]=1000 oppure 500; hbins[1]=100;} toram
 255     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);
 256     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);
 257     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);
 258     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);
 259     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);
 260     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);
 261     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);
 262     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);
 263     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);
 264     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);
 265
 266     Char_t name_h2DCAap[18][200];
 267     Char_t title_h2DCAap[18][200];
 268
 269     for(Int_t iS=0;iS<nSpec;iS++) {
 270       snprintf(name_h2DCAap[iS],200,"h2DCAap_%s",name[iS]);
 271       snprintf(title_h2DCAap[iS],200,"h2DCA_%s in for p/z<1.5GeV;DCA_{xy} (cm);DCA_{z} (cm)",name[iS]);
 272       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
 273       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
 274     }
 275
 276     Char_t name_hDCAxy[18][nbin][200];
 277     Char_t title_hDCAxy[18][nbin][200];
 278     Char_t name_hDCAz[18][nbin][200];
 279     Char_t title_hDCAz[18][nbin][200];
 280
 281     //Char_t name_h2DCA[18][nbin][200];
 282     //Char_t title_h2DCA[18][nbin][200];
 283
 284     for(Int_t iS=0;iS<nSpec;iS++) {
 285       for(Int_t j=0;j<nbin;j++) {
 286         snprintf(name_hDCAxy[iS][j],200,"hDCAxy_%s_%s",name[iS],name_nbin[j]);
 287         snprintf(title_hDCAxy[iS][j],200,"hDCAxy_%s_%s in DCAzCut;DCA_{xy} (cm)",name[iS],name_nbin[j]);
 288         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);
 289         else hDCAxy[iB][iS][j] = new TH1D(name_hDCAxy[iS][j],title_hDCAxy[iS][j],1,-3.5,3.5);
 290
 291         snprintf(name_hDCAz[iS][j],200,"hDCAz_%s_%s",name[iS],name_nbin[j]);
 292         snprintf(title_hDCAz[iS][j],200,"hDCAz_%s_%s in DCAxyCut;DCA_{z} (cm)",name[iS],name_nbin[j]);
 293         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);
 294         else hDCAz[iB][iS][j] = new TH1D(name_hDCAz[iS][j],title_hDCAz[iS][j],1,-3.5,3.5);
 295
 296         //snprintf(name_h2DCA[iS][j],200,"h2DCA_%s_%s",name[iS],name_nbin[j]);
 297         //snprintf(title_h2DCA[iS][j],200,"h2DCA_%s_%s;DCA_{xy} (cm);DCA_{z} (cm)",name[iS],name_nbin[j]);
 298         //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);
 299         //else h2DCA[iB][iS][j] = new TH2F(name_h2DCA[iS][j],title_h2DCA[iS][j],1,-4,4,1,-4,4);
 300       }
 301     }
 302
 303     Char_t name_hM2CutDCAxy[18][nbin][200];
 304     Char_t title_hM2CutDCAxy[18][nbin][200];
 305     for(Int_t iS=0;iS<nSpec;iS++) {
 306       for(Int_t j=0;j<nbin;j++) {
 307         snprintf(name_hM2CutDCAxy[iS][j],200,"hM2_CutDCAxy_%s_%s",name[iS],name_nbin[j]);
 308         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]);
 309       }
 310     }
 311
 312     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
 313     const Double_t RangeM2min[nPart]={0.0,0.0,-0.1,0.0,0.0,0.0,0.0,0.0,0.0};
 314     const Double_t RangeM2max[nPart]={1.0,1.0,0.5,2.0,4.0,6.0,12.0,4.0,6.0};
 315
 316     for(Int_t iS=0;iS<nPart;iS++) {
 317       for(Int_t j=0;j<nbin;j++) {
 318         hM2CutDCAxy[iB][iS][j] = new TH1D(name_hM2CutDCAxy[iS][j],title_hM2CutDCAxy[iS][j],BinM2pT[iS],RangeM2min[iS],RangeM2max[iS]);
 319         hM2CutDCAxy[iB][iS+nPart][j] = new TH1D(name_hM2CutDCAxy[iS+nPart][j],title_hM2CutDCAxy[iS+nPart][j],BinM2pT[iS],RangeM2min[iS],RangeM2max[iS]);
 320       }
 321     }
 322
 323     Char_t name_fPmeanVsBetaGamma[18][200];
 324     Char_t title_fPmeanVsBetaGamma[18][200];
 325
 326     if(iMtof==2) {hbins[0]=1; hbins[1]=1;}//if(iMtof==2) {hbins[0]=200; hbins[1]=200;}
 327     else {hbins[0]=1; hbins[1]=1;}
 328     for(Int_t iS=0;iS<nSpec;iS++) {
 329       snprintf(name_fPmeanVsBetaGamma[iS],200,"fPmeanVsPvtx_%s",name[iS]);
 330       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]);
 331       fPmeanVsBetaGamma[iB][iS]=new TH2F(name_fPmeanVsBetaGamma[iS],title_fPmeanVsBetaGamma[iS],hbins[0],0,10,hbins[1],0.8,1.2);
 332     }
 333
 334     Char_t name_prPmeanVsBetaGamma[18][200];
 335     Char_t title_prPmeanVsBetaGamma[18][200];
 336
 337     if(iMtof==2) {hbins[0]=1; hbins[1]=1;}//if(iMtof==2) {hbins[0]=200; hbins[1]=200;}
 338     else {hbins[0]=1; hbins[1]=1;}
 339     for(Int_t iS=0;iS<nSpec;iS++) {
 340       snprintf(name_prPmeanVsBetaGamma[iS],200,"prPmeanVsPvtx_%s",name[iS]);
 341       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]);
 342       prPmeanVsBetaGamma[iB][iS]=new TProfile(name_prPmeanVsBetaGamma[iS],title_prPmeanVsBetaGamma[iS],hbins[0],0,10,0.8,1.2,"");
 343     }
 344
 345     SetPvtxCorrections();
 346
 347     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
 348     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
 349     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
 350     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
 351
 352     SetPmeanCorrections();
 353
 354     Char_t nameTemp[14][200];
 355     snprintf(nameTemp[0],200,"#pi^{+}");
 356     snprintf(nameTemp[1],200,"K^{+}");
 357     snprintf(nameTemp[2],200,"p");
 358     snprintf(nameTemp[3],200,"d");
 359     snprintf(nameTemp[4],200,"t");
 360     snprintf(nameTemp[5],200,"He3");
 361     snprintf(nameTemp[6],200,"He4");
 362     snprintf(nameTemp[7],200,"#pi^{-}");
 363     snprintf(nameTemp[8],200,"K^{-}");
 364     snprintf(nameTemp[9],200,"#bar{p}");
 365     snprintf(nameTemp[10],200,"#bar{d}");
 366     snprintf(nameTemp[11],200,"#bar{t}");
 367     snprintf(nameTemp[12],200,"#bar{He3}");
 368     snprintf(nameTemp[13],200,"#bar{He4}");
 369     Char_t name_prPmeanVsBGcorr[14][200];
 370     Char_t title_prPmeanVsBGcorr[14][200];
 371
 372     hbins[0]=200;
 373     for(Int_t iS=0;iS<14;iS++) {
 374       snprintf(name_prPmeanVsBGcorr[iS],200,"prPmeanVsBGcorr_%s",nameTemp[iS]);
 375       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]);
 376       prPmeanVsBGcorr[iB][iS]=new TProfile(name_prPmeanVsBGcorr[iS],title_prPmeanVsBGcorr[iS],hbins[0],0,20,0.8,1.2,"");
 377     }
 378
 379     fList[iB]->Add(htemp[iB]);
 380     for(Int_t i=0;i<2;i++) fList[iB]->Add(htriggerbits[iB][i]);
 381     for(Int_t i=0;i<2;i++) fList[iB]->Add(hCentrality[iB][i]);
 382     for(Int_t i=0;i<2;i++) fList[iB]->Add(hZvertex[iB][i]);
 383     fList[iB]->Add(hEta[iB]);
 384     fList[iB]->Add(hPhi[iB]);
 385     //fList[iB]->Add(fEtaPhi[iB]);
 386     fList[iB]->Add(hNTpcCluster[iB]);
 387     fList[iB]->Add(hNTrdSlices[iB]);
 388     //for(Int_t i=0;i<2;i++) fList[iB]->Add(fdEdxVSp[iB][i]);
 389     //for(Int_t i=0;i<nPart;i++) fList[iB]->Add(hDeDxExp[iB][i]);
 390     for(Int_t i=0;i<nSpec;i++) fList[iB]->Add(fNsigmaTpc[iB][i]);
 391     for(Int_t i=0;i<nPart;i++) {
 392       if(kSignalCheck!=1)
 393         if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
 394       //fList[iB]->Add(fNsigmaTpc_kTOF[iB][i]);
 395       //fList[iB]->Add(fNsigmaTpc_kTOF[iB][i+nPart]);
 396     }
 397     //for(Int_t i=0;i<2;i++) fList[iB]->Add(fBetaTofVSp[iB][i]);
 398     //for(Int_t i=0;i<nPart;i++) fList[iB]->Add(hBetaExp[iB][i]);
 399     //for(Int_t i=0;i<nPart;i++) fList[iB]->Add(fNsigmaTof[iB][i]);
 400     //for(Int_t i=0;i<2;i++) fList[iB]->Add(fM2vsP_NoTpcCut[iB][0][i]);
 401     for(Int_t i=0;i<nPart;i++) {
 402       if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
 403       //fList[iB]->Add(fM2vsP[iB][0][i]);
 404       //fList[iB]->Add(fM2vsP[iB][0][i+nPart]);
 405     }
 406
 407     for(Int_t i=0;i<2;i++){
 408       //fList[iB]->Add(fPvtxTrueVsReco[i]);
 409       //fList[iB]->Add(prPvtxTrueVsReco[iB][i]);
 410     }
 411     if(iMtof==2) {
 412       for(Int_t i=0;i<nPart;i++){
 413         if(i<2) continue;//e,mu excluded
 414         //fList[iB]->Add(fPmeanVsBetaGamma[iB][i]);
 415         //fList[iB]->Add(prPmeanVsBetaGamma[iB][i]);
 416         //fList[iB]->Add(fPmeanVsBetaGamma[iB][i+nPart]);
 417         //fList[iB]->Add(prPmeanVsBetaGamma[iB][i+nPart]);
 418       }
 419     }
 420     if(iMtof>2) {
 421       //for(Int_t i=0;i<14;i++)fList[iB]->Add(fPmeanVsBGcorr[i]);
 422       //for(Int_t i=0;i<14;i++)fList[iB]->Add(prPmeanVsBGcorr[iB][i]);
 423     }
 424     for(Int_t i=0;i<nPart;i++) {
 425       if(i<5 || i==6 || i==8) continue;//e,mu,pi,K,p,t,he4 excluded//i<5 || i==6 || i==8
 426       fList[iB]->Add(h2DCAap[iB][i]);
 427       fList[iB]->Add(h2DCAap[iB][i+nPart]);
 428     }
 429     /*
 430       for(Int_t i=0;i<nPart;i++) {
 431       if(i<5 || i==6 || i==8) continue;//e,mu,pi,K,p,t,he4 excluded//i<5 || i==6 || i==8
 432       for(Int_t j=0;j<nbin;j++){
 433       fList[iB]->Add(h2DCA[iB][i][j]);
 434       fList[iB]->Add(h2DCA[iB][i+nPart][j]);
 435       }
 436       }
 437     */
 438     //for(Int_t i=0;i<10;i++) fList[iB]->Add(fM2vsZ[iB][i]);
 439     for(Int_t i=0;i<nPart;i++){
 440       if(kSignalCheck!=1)
 441         if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
 442       for(Int_t j=0;j<nbin;j++){
 443         fList[iB]->Add(hDCAxy[iB][i][j]);
 444         fList[iB]->Add(hDCAz[iB][i][j]);
 445         fList[iB]->Add(hM2CutDCAxy[iB][i][j]);
 446         fList[iB]->Add(hDCAxy[iB][i+nPart][j]);
 447         fList[iB]->Add(hDCAz[iB][i+nPart][j]);
 448         fList[iB]->Add(hM2CutDCAxy[iB][i+nPart][j]);
 449       }
 450     }
 451
 452     // Post output data.
 453     PostData(1, fList[0]);
 454     PostData(2, fList[1]);
 455
 456   }//end iB loop
 457 }
 458 //______________________________________________________________________________
 459 void AliAnalysisNucleiMass::UserExec(Option_t *)
 460 {
 461   // Main loop
 462   // Called for each event
 463
 464   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 465   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 466   if(!fAOD && !fESD){
 467     Printf("%s:%d AODEvent and ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 468     return;
 469   }
 470
 471   if(fESD) fEvent = fESD;
 472   else fEvent = fAOD;
 473
 474   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 475   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 476   fPIDResponse=inputHandler->GetPIDResponse();
 477
 478   //--------------------------Magnetic field polarity--------------------
 479   Double_t fBfield=fEvent->GetMagneticField();
 480   if(fBfield<0.0) iBconf=0;//B--
 481   else iBconf=1;//B++
 482   for(Int_t i=0;i<nBconf;i++) htemp[i]->Fill(fBfield);
 483
 484   //--------------------------Centrality--------------------------------
 485   Double_t v0Centr  = -10.;
 486   AliCentrality *centrality = fEvent->GetCentrality();
 487   if (centrality){
 488     v0Centr=centrality->GetCentralityPercentile("V0M"); // VZERO
 489   }
 490   hCentrality[iBconf][0]->Fill(v0Centr);
 491
 492   //-------------------------zVertex determination of event----------------
 493   Double_t zvtx = 9999.9;
 494   const AliVVertex* vtxEVENT = fEvent->GetPrimaryVertex();
 495   if(vtxEVENT->GetNContributors()>0) zvtx = vtxEVENT->GetZ();
 496
 497   hZvertex[iBconf][0]->Fill(zvtx);
 498
 499   //---------------------------EVENT CUTS-----------------------------
 500   if(TMath::Abs(zvtx) < 10.0 && v0Centr>Centrality[0] && v0Centr<Centrality[1]){
 501
 502     //TRIGGER SELECTION
 503     Int_t iTrigger=-2;
 504
 505     if(inputHandler->IsEventSelected() & AliVEvent::kMB) iTrigger = 0;
 506     if(inputHandler->IsEventSelected() & AliVEvent::kCentral) iTrigger = 16;
 507     if(inputHandler->IsEventSelected() & AliVEvent::kSemiCentral) iTrigger = 17;
 508     //if((((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()) & AliVEvent::kAny) iTrigger = 35;
 509
 510     if(iTriggerSel!=-99) {//if a dedicated trigger is required
 511       if(iTrigger!=iTriggerSel) return;
 512     }
 513
 514     for(Int_t i=0;i<32;i++) {
 515       Int_t bit=(1<<i);
 516       if(inputHandler->IsEventSelected() & bit) htriggerbits[iBconf][0]->Fill(i);
 517     }
 518     if(inputHandler->IsEventSelected() & AliVEvent::kAny) htriggerbits[iBconf][0]->Fill(35);
 519     if(inputHandler->IsEventSelected() & AliVEvent::kAnyINT) htriggerbits[iBconf][0]->Fill(36);
 520
 521     htriggerbits[iBconf][1]->Fill(iTrigger);
 522
 523     hCentrality[iBconf][1]->Fill(v0Centr);
 524     hZvertex[iBconf][1]->Fill(zvtx);
 525
 526     Int_t nTracks = fEvent->GetNumberOfTracks();
 527
 528     //----------------------loop on the TRACKS-----------------------------
 529     for(Int_t iT = 0; iT < nTracks; iT++) {
 530       AliVTrack* track = (AliVTrack *) fEvent->GetTrack(iT);
 531
 532       if (!track){
 533         continue;
 534       }
 535
 536      //For the geometrical cuts
 537       Double_t eta = track->Eta();
 538
 539       Bool_t trkFlag = 0;
 540       trkFlag = ((AliAODTrack *) track)->TestFilterBit(FilterBit);
 541       //TestFilterBit(16) -- Standard Cuts with very loose DCA: GetStandardITSTPCTrackCuts2011(kFALSE) && SetMaxDCAToVertexXY(2.4) && SetMaxDCAToVertexZ(3.2) && SetDCaToVertex2D(kTRUE)
 542       //TestFilterBit(32) (STARDARD) -- Standard Cuts with very tight DCA cut ( 7sigma^primaries: 7*(0.0015+0.0050/pt^1.1) ) : GetStandardITSTPCTrackCuts2011().
 543
 544       //Cut on the Minumum Number of the TPC clusters
 545       Bool_t isMinTpcCluster=kFALSE;
 546       Int_t nTpcCluster=0;
 547       nTpcCluster=track->GetTPCNcls();
 548       if(nTpcCluster>NminTpcCluster) isMinTpcCluster=kTRUE;
 549
 550       //-------------------------------------start TRACK CUTS (I): for (II) see below--------
 551       if ((track->Pt() < 0.2) || (eta<EtaLimit[0]) || (eta>EtaLimit[1]) || !trkFlag || !isMinTpcCluster)
 552         continue;
 553
 554       //For the Tpc purity cut
 555       Double_t dedx = track->GetTPCsignal();
 556       if(dedx<10) continue;
 557
 558       Int_t nTrdSlices = track->GetNumberOfTRDslices();
 559       if(nTrdSlices<2 && iTrdCut==1) continue;
 560       if(nTrdSlices>0 && iTrdCut==2) continue;
 561
 562       //-------------------------------------end TRACK CUTS (I)----------------------------------
 563
 564       //-------------------------------------Track info--------------------------------------
 565       Double_t phi= track->Phi();
 566       Double_t charge = (Double_t)track->Charge();
 567       Double_t p = track->P();
 568       Double_t pt = track->Pt();
 569       Double_t tof  = track->GetTOFsignal()-fPIDResponse->GetTOFResponse().GetStartTime(p);
 570       Double_t pTPC = track->GetTPCmomentum();
 571       Double_t beta = 0.0;
 572       Double_t M2 = 999.9;
 573       Double_t Z2 = 999.9;
 574
 575        //Vertex determination
 576       Double_t b[2] = {-99., -99.};
 577       Double_t bCov[3] = {-99., -99., -99.};
 578       if (!track->PropagateToDCA(fEvent->GetPrimaryVertex(), fEvent->GetMagneticField(), 100., b, bCov))
 579         continue;
 580
 581       Double_t DCAxy = b[0];
 582       Double_t DCAz = b[1];
 583
 584       kTOF = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
 585
 586       //-----------------------------TPC info------------------------------
 587       Double_t nsigmaTPC[nPart];
 588       Double_t expdedx[nPart];
 589
 590       Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
 591       Int_t FlagPid = 0;
 592
 593       for(Int_t iS=0;iS<9;iS++){
 594         nsigmaTPC[iS] = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType) iS);
 595         //TPC identification:
 596         if(TMath::Abs(nsigmaTPC[iS])<NsigmaTpcCut) {
 597           FlagPid += ((Int_t)TMath::Power(2,iS));
 598         }
 599       }
 600       //Correction of the momentum to the vertex for (anti)nuclei
 601       Double_t pC[9];
 602       for(Int_t iS=0;iS<9;iS++) pC[iS]=p;
 603       this->MomVertexCorrection(p,pC,eta,FlagPid);
 604
 605       this->FillDCAdist(DCAxy,DCAz,charge,FlagPid,stdFlagPid,pC);
 606
 607       //-------------------------------------start TRACK CUTS (II)-------------------------------------
 608       //Cut on the DCAxy
 609       Bool_t isDCAxyCut=kFALSE;
 610       if(TMath::Abs(DCAxy)<DCAxyCut) isDCAxyCut=kTRUE;
 611
 612       //Cut on the DCAz
 613       Bool_t isDCAzCut=kFALSE;
 614       if(TMath::Abs(DCAz)<DCAzCut) isDCAzCut=kTRUE;
 615
 616       if (!isDCAxyCut || !isDCAzCut)
 617         continue;
 618
 619       //-------------------------------------end TRACK CUTS (II)----------------------------------
 620
 621       hEta[iBconf]->Fill(eta);
 622       hPhi[iBconf]->Fill(phi);
 623       fEtaPhi[iBconf]->Fill(eta,phi);
 624       hNTpcCluster[iBconf]->Fill(nTpcCluster);
 625       hNTrdSlices[iBconf]->Fill(nTrdSlices);
 626
 627       //More TPC info:
 628       for(Int_t iS=0;iS<9;iS++){
 629         expdedx[iS] = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, (AliPID::EParticleType) iS, AliTPCPIDResponse::kdEdxDefault, kTRUE);
 630         hDeDxExp[iBconf][iS]->Fill(pTPC,expdedx[iS]);
 631         nsigmaTPC[iS] = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType) iS);
 632         //fNsigmaTpc[iBconf][iS]->Fill(pTPC,nsigmaTPC[iS]);
 633         if(charge>0) {//positive particle
 634           fNsigmaTpc[iBconf][iS]->Fill(pTPC,nsigmaTPC[iS]);
 635           if(kTOF) fNsigmaTpc_kTOF[iBconf][iS]->Fill(p,nsigmaTPC[iS]);
 636         }
 637         else {//negative particle
 638           fNsigmaTpc[iBconf][iS+nPart]->Fill(pTPC,nsigmaTPC[iS]);
 639           if(kTOF) fNsigmaTpc_kTOF[iBconf][iS+nPart]->Fill(p,nsigmaTPC[iS]);
 640         }
 641         /*
 642           if(TMath::Abs(nsigmaTPC[iS])<NsigmaTpcCut) {
 643           FlagPid += ((Int_t)TMath::Power(2,iS));
 644           }*/
 645       }
 646
 647       if(charge>0) fdEdxVSp[iBconf][0]->Fill(pTPC,dedx);
 648       else fdEdxVSp[iBconf][1]->Fill(pTPC,dedx);
 649
 650       //-----------------------------TOF info------------------------------
 651
 652       Double_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.875612859,2.808921005,1.404195741,1.863689620};
 653
 654       //----------------------------------------kTOF available-----------------------------
 655
 656       if(kTOF) {
 657         Double_t exptimes[9];
 658         track->GetIntegratedTimes(exptimes);
 659         //Integrated times of the Nuclei:
 660         for(Int_t iN=5;iN<9;iN++) {
 661           exptimes[iN] = exptimes[4]*exptimes[4]*(massOverZ[iN]*massOverZ[iN]/p/p+1)/(massOverZ[4]*massOverZ[4]/p/p+1);
 662           exptimes[iN] = TMath::Sqrt(exptimes[iN]);
 663         }
 664
 665         beta=exptimes[0];
 666         beta=beta/tof;//beta = L/tof/c = t_e/tof
 667
 668         Int_t FlagPidTof = 0;
 669         Double_t NsigmaTofCut = 2.0;
 670
 671         Double_t nsigmaTOF[9];
 672         for(Int_t iS=0;iS<9;iS++){
 673           nsigmaTOF[iS] = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType) iS);
 674           fNsigmaTof[iBconf][iS]->Fill(pt,nsigmaTOF[iS]);
 675           if(charge>0) {
 676             hBetaExp[iBconf][iS]->Fill(p,exptimes[0]/exptimes[iS]);
 677           }
 678           else {
 679             hBetaExp[iBconf][iS+nPart]->Fill(p,exptimes[0]/exptimes[iS]);
 680           }
 681
 682           //TOF identification:
 683           if(TMath::Abs(nsigmaTOF[iS])<NsigmaTofCut) {
 684             FlagPidTof += ((Int_t)TMath::Power(2,iS));
 685           }
 686         }
 687
 688         if(charge>0) fBetaTofVSp[iBconf][0]->Fill(p,beta);
 689         else fBetaTofVSp[iBconf][1]->Fill(p,beta);
 690
 691         this->GetMassFromPvertex(beta,p,M2);
 692         this->GetZTpc(dedx,pTPC,M2,Z2);
 693
 694         Double_t Mass2[9];
 695         //-----------------------------M2 as a function of momentum to the primary vertex if iMtof==1---------------------------------
 696         if(iMtof==1) this->GetMassFromPvertexCorrected(beta,pC,Mass2);
 697
 698         if(iMtof==2) this->GetPmeanVsBetaGamma(exptimes,pC,FlagPid,FlagPidTof,charge);
 699
 700         //-----------------------------M2 as a function of expected times---------------------------------
 701         if(iMtof==2) this->GetMassFromExpTimes(beta,exptimes,Mass2);
 702
 703         //-----------------------------M2 as a function of mean momentum calculated from expected time and extrapolated to the (anti)nuclei---------------------------------
 704         if(iMtof>2) this->GetMassFromMeanMom(beta,exptimes,pC,eta,charge,Mass2,FlagPid,FlagPidTof);
 705
 706         //-------------------------------Squared Mass TH2 distributions-----------------------
 707         if(charge>0) {
 708           //without TPC
 709           fM2vsP_NoTpcCut[iBconf][0][0]->Fill(M2,p);
 710           //with TPC
 711           for(Int_t iS=0;iS<9;iS++) {
 712             M2=999.9;
 713             M2=Mass2[iS];
 714             //-----------------
 715             if(FlagPid & stdFlagPid[iS]) {
 716               fM2vsP[iBconf][0][iS]->Fill(M2,pC[iS]);
 717             }
 718           }
 719         }
 720         else {//charge<0
 721           //without TPC
 722           fM2vsP_NoTpcCut[iBconf][0][1]->Fill(M2,p);
 723           //with TPC
 724           for(Int_t iS=0;iS<9;iS++) {
 725             M2=999.9;
 726             M2=Mass2[iS];
 727             //-----------------
 728             if(FlagPid & stdFlagPid[iS]) {
 729               fM2vsP[iBconf][0][iS+nPart]->Fill(M2,pC[iS]);
 730             }
 731           }
 732         }
 733
 734         //------------------------------start Squared Mass TH1 distributions-------------------------
 735         Double_t binP[nbin+1];
 736         for(Int_t i=0;i<nbin+1;i++) {
 737           binP[i]=0.4+i*0.1;
 738         }
 739
 740         if(charge>0) {
 741           for(Int_t iS=0;iS<9;iS++) {
 742             M2=999.9;
 743             M2=Mass2[iS];
 744
 745             if(FlagPid & stdFlagPid[iS]) {
 746               for(Int_t j=0;j<nbin;j++) {
 747                 if(pC[iS]>binP[j] && pC[iS]<binP[j+1]) {
 748                   hM2CutDCAxy[iBconf][iS][j]->Fill(M2);
 749                   break;
 750                 }
 751               }//end loop on the p bins (j)
 752             }
 753           }//end loop on the particle species (iS)
 754         }
 755         else {//charge<0
 756           for(Int_t iS=0;iS<9;iS++) {
 757             M2=999.9;
 758             M2=Mass2[iS];
 759
 760             if(FlagPid & stdFlagPid[iS]) {
 761               for(Int_t j=0;j<nbin;j++) {
 762                 if(pC[iS]>binP[j] && pC[iS]<binP[j+1]) {
 763                   hM2CutDCAxy[iBconf][iS+nPart][j]->Fill(M2);
 764                   break;
 765                 }
 766               }//end loop on the p bins (j)
 767             }
 768           }//end loop on the particle species (iS)
 769         }
 770
 771         //-------------------------------------------------M2/Z2 vs Z-------------------------
 772
 773
 774         Double_t binCutPt[10] = {0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0};
 775         Double_t Z=999.9;
 776         if(Z2>0) Z=TMath::Sqrt(Z2);
 777
 778         fM2vsZ[iBconf][0]->Fill(charge*TMath::Sqrt(Z2),M2);
 779         for(Int_t i=1;i<10;i++) {
 780           if(pt>binCutPt[i-1] && pt<binCutPt[i]){
 781             fM2vsZ[iBconf][i]->Fill(charge*Z,M2);
 782             break;
 783           }
 784         }
 785
 786       }//end kTOF available
 787     }//end track loop
 788   }//end loop on the events
 789 }
 790
 791 //_____________________________________________________________________________
 792 void AliAnalysisNucleiMass::Terminate(Option_t *)
 793 {
 794   // Terminate loop
 795   Printf("Terminate()");
 796 }
 797 //_____________________________________________________________________________
 798 void AliAnalysisNucleiMass::MomVertexCorrection(Double_t p, Double_t *pC, Double_t eta, Int_t FlagPid){
 799
 800   Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
 801
 802   for(Int_t iS=0;iS<9;iS++) {
 803     if(FlagPid & stdFlagPid[iS]) {
 804       if(iS==5) {
 805         if(kPvtxCorr==1) pC[iS]=pC[iS]*fPvtxTrueVsReco[0]->Eval(pC[iS],TMath::Abs(eta));//for (bar)d
 806         prPvtxTrueVsReco[iBconf][0]->Fill(p,pC[iS]/p);
 807       }
 808       else if(iS==6) {
 809         if(kPvtxCorr==1) pC[iS]=pC[iS]*fPvtxTrueVsReco[1]->Eval(pC[iS],TMath::Abs(eta));//for (bar)t
 810         prPvtxTrueVsReco[iBconf][1]->Fill(p,pC[iS]/p);
 811       }
 812       else if(iS==7) {
 813         if(kPvtxCorr==1) pC[iS]=pC[iS]*fPvtxTrueVsReco[2]->Eval(pC[iS],TMath::Abs(eta));//for (bar)He3
 814         prPvtxTrueVsReco[iBconf][2]->Fill(p,pC[iS]/p);
 815       }
 816       else if(iS==8) {
 817         if(kPvtxCorr==1) pC[iS]=pC[iS]*fPvtxTrueVsReco[3]->Eval(pC[iS],TMath::Abs(eta));//for (bar)He3
 818         prPvtxTrueVsReco[iBconf][3]->Fill(p,pC[iS]/p);
 819       }
 820     }
 821   }
 822
 823   return;
 824
 825 }
 826 //__________________________________________________________________________________________________
 827 void AliAnalysisNucleiMass::FillDCAdist(Double_t DCAxy, Double_t DCAz, Double_t charge, Int_t FlagPid, Int_t stdFlagPid[9], Double_t *pC){
 828
 829   Double_t binP[nbin+1];
 830   for(Int_t i=0;i<nbin+1;i++) {
 831     binP[i]=0.4+i*0.1;
 832   }
 833
 834   if(charge>0) {
 835     for(Int_t iS=0;iS<9;iS++) {
 836       if(FlagPid & stdFlagPid[iS]) {
 837         if(pC[iS]<1.5) {
 838           h2DCAap[iBconf][iS]->Fill(DCAxy,DCAz);
 839           h2DCAap[iBconf][iS]->Fill(-DCAxy,-DCAz);
 840         }
 841         for(Int_t j=0;j<nbin;j++) {
 842           if(pC[iS]>binP[j] && pC[iS]<binP[j+1]) {
 843             if(TMath::Abs(DCAz)<DCAzCut) {
 844               hDCAxy[iBconf][iS][j]->Fill(DCAxy);
 845               hDCAxy[iBconf][iS][j]->Fill(-DCAxy);
 846             }
 847             if(TMath::Abs(DCAxy)<DCAxyCut) {
 848               hDCAz[iBconf][iS][j]->Fill(DCAz);
 849               hDCAz[iBconf][iS][j]->Fill(-DCAz);
 850             }
 851             //h2DCA[iBconf][iS][j]->Fill(DCAxy,DCAz);
 852             //h2DCA[iBconf][iS][j]->Fill(-DCAxy,-DCAz);
 853             break;
 854           }
 855         }//end loop on the p bins (j)
 856       }
 857     }//end loop on the particle species (iS)
 858   }
 859   else {//charge<0
 860     for(Int_t iS=0;iS<9;iS++) {
 861       if(FlagPid & stdFlagPid[iS]) {
 862         if(pC[iS]<1.5) {
 863           h2DCAap[iBconf][iS+nPart]->Fill(DCAxy,DCAz);
 864           h2DCAap[iBconf][iS+nPart]->Fill(-DCAxy,-DCAz);
 865         }
 866         for(Int_t j=0;j<nbin;j++) {
 867           if(pC[iS]>binP[j] && pC[iS]<binP[j+1]) {
 868             if(TMath::Abs(DCAz)<DCAzCut) {
 869               hDCAxy[iBconf][iS+nPart][j]->Fill(DCAxy);
 870               hDCAxy[iBconf][iS+nPart][j]->Fill(-DCAxy);
 871             }
 872             if(TMath::Abs(DCAxy)<DCAxyCut) {
 873               hDCAz[iBconf][iS+nPart][j]->Fill(DCAz);
 874               hDCAz[iBconf][iS+nPart][j]->Fill(-DCAz);
 875             }
 876             //h2DCA[iBconf][iS+nPart][j]->Fill(DCAxy,DCAz);
 877             //h2DCA[iBconf][iS+nPart][j]->Fill(-DCAxy,-DCAz);
 878             break;
 879           }
 880         }//end loop on the p bins (j)
 881       }
 882     }//end loop on the particle species (iS)
 883   }
 884
 885   return;
 886 }
 887 //_____________________________________________________________________________
 888 void AliAnalysisNucleiMass::GetMassFromPvertex(Double_t beta, Double_t p, Double_t &M2) {
 889
 890   M2 = p*p*(1-beta*beta)/(beta*beta);
 891
 892   return;
 893
 894 }
 895 //_________________________________________________________________________________________________________________________
 896 void AliAnalysisNucleiMass::GetZTpc(Double_t dedx, Double_t pTPC, Double_t M2, Double_t &Z2) {
 897
 898   //z^2_tpc = dedx^{Tpc} / dedx^{exp,Tof}_{z=1}
 899
 900   Z2=999.9;
 901
 902   Double_t M=999.9;
 903   Double_t pTPC_pr=999.9;//rescaling of the pTPC for the proton
 904   Double_t expdedx_Tof=999.9;
 905
 906   if(M2>0) {
 907     M=TMath::Sqrt(M2);
 908     pTPC_pr=pTPC*0.938272/M;
 909     expdedx_Tof=fPIDResponse->GetTPCResponse().GetExpectedSignal(pTPC_pr,AliPID::kProton);
 910     if((dedx/expdedx_Tof)<0) return;
 911     Z2=TMath::Power(dedx/expdedx_Tof,0.862);
 912   }
 913
 914   return;
 915 }
 916 //_________________________________________________________________________________________________________________________
 917 void AliAnalysisNucleiMass::GetMassFromPvertexCorrected(Double_t beta, Double_t *pC, Double_t *Mass2) {
 918
 919   for(Int_t iS=0;iS<9;iS++) Mass2[iS] = pC[iS]*pC[iS]*(1-beta*beta)/(beta*beta);
 920
 921   return;
 922 }
 923 //____________________________________________________________________________________________________________
 924 void AliAnalysisNucleiMass::GetMassFromExpTimes(Double_t beta, Double_t *IntTimes, Double_t *Mass2) {
 925
 926   // 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
 927   // In this way m_tof = mPDG only if tof=t_exp
 928
 929   Double_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.875612859,2.808921005,1.404195741,1.863689620};
 930
 931   Double_t beta2Exp[9];
 932   Double_t p2Exp[9];
 933
 934   //Double_t pExp[9];
 935
 936   for(Int_t iS=0;iS<9;iS++) {
 937     beta2Exp[iS]=IntTimes[0]/IntTimes[iS];//beta = L/tof*c = t_e/tof
 938     beta2Exp[iS]=beta2Exp[iS]*beta2Exp[iS];
 939     if((1-beta2Exp[iS])==0) {
 940       Mass2[iS]=999.9;
 941       continue;
 942     }
 943     p2Exp[iS]=massOverZ[iS]*massOverZ[iS]*beta2Exp[iS]/(1-beta2Exp[iS]);
 944
 945     //--------------------for MC corrections
 946     if(p2Exp[iS]<0) {
 947       Mass2[iS]=999.9;
 948       continue;
 949     }
 950     //pExp[iS]=TMath::Sqrt(p2Exp[iS]);
 951
 952     //------------
 953     Mass2[iS]=p2Exp[iS]*(1-beta*beta)/(beta*beta);
 954   }//end loop on the particle species
 955
 956   return;
 957 }
 958 //____________________________________________________________________________________________________________
 959 void AliAnalysisNucleiMass::GetPmeanVsBetaGamma(Double_t *IntTimes, Double_t *pVtx, Int_t FlagPid, Int_t FlagPidTof, Double_t charge) {
 960
 961   // 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
 962   // In this way m_tof = mPDG only if tof=t_exp
 963
 964   Double_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.875612859,2.808921005,1.404195741,1.863689620};
 965
 966   Double_t beta2Exp[9];
 967   Double_t p2Exp[9];
 968
 969   Double_t pExp[9];
 970
 971   Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
 972
 973   for(Int_t iS=0;iS<9;iS++) {
 974     beta2Exp[iS]=IntTimes[0]/IntTimes[iS];//beta = L/tof*c = t_e/tof
 975     beta2Exp[iS]=beta2Exp[iS]*beta2Exp[iS];
 976     if((1-beta2Exp[iS])==0) {
 977       continue;
 978     }
 979     p2Exp[iS]=massOverZ[iS]*massOverZ[iS]*beta2Exp[iS]/(1-beta2Exp[iS]);
 980
 981     if(p2Exp[iS]<0) {
 982       continue;
 983     }
 984     pExp[iS]=TMath::Sqrt(p2Exp[iS]);
 985
 986     if((FlagPid & stdFlagPid[iS]) && (FlagPidTof & stdFlagPid[iS])) {
 987       if(charge>0){
 988         fPmeanVsBetaGamma[iBconf][iS]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
 989         prPmeanVsBetaGamma[iBconf][iS]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
 990       }
 991       else {
 992         fPmeanVsBetaGamma[iBconf][iS+nPart]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
 993         prPmeanVsBetaGamma[iBconf][iS+nPart]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
 994       }
 995     }
 996   }//end loop on the particle species
 997
 998   return;
 999
1000 }
1001 //____________________________________________________________________________________________________________
1002 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
1003
1004   // 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
1005   // In this way m_tof = mPDG only if tof=t_exp
1006
1007   Double_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.875612859,2.808921005,1.404195741,1.863689620};
1008
1009   Double_t beta2Exp[9];
1010   Double_t p2Exp[9];
1011
1012   Double_t pExp[9];
1013
1014   Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
1015
1016   for(Int_t iS=0;iS<9;iS++) {
1017     if(iS>1) {
1018       p2Exp[iS]=pVtx[iS]*fPmeanVsBGcorr[iS-2]->Eval(pVtx[iS]/massOverZ[iS],TMath::Abs(eta));
1019       p2Exp[iS]*=p2Exp[iS];
1020     }
1021     else {
1022       beta2Exp[iS]=IntTimes[0]/IntTimes[iS];//beta = L/tof*c = t_e/tof
1023       beta2Exp[iS]=beta2Exp[iS]*beta2Exp[iS];
1024       if((1-beta2Exp[iS])==0) {
1025         Mass2[iS]=999.9;
1026         continue;
1027       }
1028       p2Exp[iS]=massOverZ[iS]*massOverZ[iS]*beta2Exp[iS]/(1-beta2Exp[iS]);
1029     }
1030
1031     if(p2Exp[iS]<0) {
1032       Mass2[iS]=999.9;
1033       continue;
1034     }
1035     pExp[iS]=TMath::Sqrt(p2Exp[iS]);
1036
1037     //------------
1038     Mass2[iS]=p2Exp[iS]*(1-beta*beta)/(beta*beta);
1039
1040     //-----------
1041
1042     if(iS>1) {
1043       if((FlagPid & stdFlagPid[iS]) && (FlagPidTof & stdFlagPid[iS])) {
1044         if(charge>0) {
1045           prPmeanVsBGcorr[iBconf][iS-2]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
1046         }
1047         else if(charge<0) {
1048           prPmeanVsBGcorr[iBconf][iS-2+7]->Fill(pVtx[iS]/massOverZ[iS],pExp[iS]/pVtx[iS]);
1049         }
1050       }
1051     }
1052   }//end loop on the particle species
1053
1054   return;
1055
1056 }
1057 //________________________________________________________________________________________
1058 void AliAnalysisNucleiMass::SetPvtxCorrections(){
1059   //for (bar)d
1060   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
1061   fPvtxTrueVsReco[0]->SetParameter(0,0.031263);
1062   fPvtxTrueVsReco[0]->SetParameter(1,-3.276770);
1063   fPvtxTrueVsReco[0]->SetParameter(2,1.000113);
1064   fPvtxTrueVsReco[0]->SetParameter(3,-5.195875);
1065   fPvtxTrueVsReco[0]->SetParameter(4,1.000674);
1066   fPvtxTrueVsReco[0]->SetParameter(5,2.870503);
1067   fPvtxTrueVsReco[0]->SetParameter(6,3.777729);
1068
1069   //for (bar)t
1070   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
1071   fPvtxTrueVsReco[1]->SetParameter(0,8.79761e-02);
1072   fPvtxTrueVsReco[1]->SetParameter(1,-3.23189e+00);
1073   fPvtxTrueVsReco[1]->SetParameter(2,9.99578e-01);
1074   fPvtxTrueVsReco[1]->SetParameter(3,0.0);
1075
1076   //for (bar)He3
1077   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
1078   fPvtxTrueVsReco[2]->SetParameter(0,0.037986);
1079   fPvtxTrueVsReco[2]->SetParameter(1,-2.707620);
1080   fPvtxTrueVsReco[2]->SetParameter(2,1.000742);
1081   fPvtxTrueVsReco[2]->SetParameter(3,-4.934743);
1082   fPvtxTrueVsReco[2]->SetParameter(4,1.001640);
1083   fPvtxTrueVsReco[2]->SetParameter(5,2.744372);
1084   fPvtxTrueVsReco[2]->SetParameter(6,3.528561);
1085
1086   //for (bar)He4
1087   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
1088   fPvtxTrueVsReco[3]->SetParameter(0,7.08785e-02);
1089   fPvtxTrueVsReco[3]->SetParameter(1,-2.87201e+00);
1090   fPvtxTrueVsReco[3]->SetParameter(2,1.00070e+00);
1091   fPvtxTrueVsReco[3]->SetParameter(3,0.0);
1092
1093   for(Int_t i=0;i<4;i++) {
1094     fPvtxTrueVsReco[i]->SetNpx(fPvtxTrueVsReco[i]->GetNpx()*10.0);
1095   }
1096 }
1097 //________________________________________________________________________________________
1098 void AliAnalysisNucleiMass::SetPmeanCorrections(){
1099
1100   Char_t nameTemp[14][200];
1101   snprintf(nameTemp[0],200,"#pi^{+}");
1102   snprintf(nameTemp[1],200,"K^{+}");
1103   snprintf(nameTemp[2],200,"p");
1104   snprintf(nameTemp[3],200,"d");
1105   snprintf(nameTemp[4],200,"t");
1106   snprintf(nameTemp[5],200,"He3");
1107   snprintf(nameTemp[6],200,"He4");
1108   snprintf(nameTemp[7],200,"#pi^{-}");
1109   snprintf(nameTemp[8],200,"K^{-}");
1110   snprintf(nameTemp[9],200,"#bar{p}");
1111   snprintf(nameTemp[10],200,"#bar{d}");
1112   snprintf(nameTemp[11],200,"#bar{t}");
1113   snprintf(nameTemp[12],200,"#bar{He3}");
1114   snprintf(nameTemp[13],200,"#bar{He4}");
1115
1116   Char_t name_fPmeanVsBGcorr[14][200];
1117   for(Int_t i=0;i<14;i++) {
1118     snprintf(name_fPmeanVsBGcorr[i],200,"fPmeanVsBGcorr_%s",nameTemp[i]);
1119   }
1120
1121   //Pions
1122   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);
1123   fPmeanVsBGcorr[0]->SetParameter(0,-0.179607);
1124   fPmeanVsBGcorr[0]->SetParameter(1,-0.384809);
1125   fPmeanVsBGcorr[0]->SetParameter(2,0.885534);
1126   fPmeanVsBGcorr[0]->SetParameter(3,0.992710);
1127   fPmeanVsBGcorr[0]->SetParameter(4,0.011390);
1128   fPmeanVsBGcorr[0]->SetParameter(5,3.231000);
1129   fPmeanVsBGcorr[0]->SetParameter(6,0.999900);
1130
1131   //Kaons
1132   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);
1133   fPmeanVsBGcorr[1]->SetParameter(0,0.033500);
1134   fPmeanVsBGcorr[1]->SetParameter(1,-2.461673);
1135   fPmeanVsBGcorr[1]->SetParameter(2,0.996501);
1136   fPmeanVsBGcorr[1]->SetParameter(3,1.000000);
1137   fPmeanVsBGcorr[1]->SetParameter(4,0.089715);
1138   fPmeanVsBGcorr[1]->SetParameter(5,-2.473531);
1139   fPmeanVsBGcorr[1]->SetParameter(6,1.000000);
1140   fPmeanVsBGcorr[1]->SetParameter(7,-1.562500);
1141   fPmeanVsBGcorr[1]->SetParameter(8,0.253000);
1142   fPmeanVsBGcorr[1]->SetParameter(9,0.009387);
1143
1144   //Protons
1145   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);
1146   fPmeanVsBGcorr[2]->SetParameter(0,0.015081);
1147   fPmeanVsBGcorr[2]->SetParameter(1,-2.927557);
1148   fPmeanVsBGcorr[2]->SetParameter(2,0.997904);
1149   fPmeanVsBGcorr[2]->SetParameter(3,1.000000);
1150   fPmeanVsBGcorr[2]->SetParameter(4,0.102697);
1151   fPmeanVsBGcorr[2]->SetParameter(5,-3.399528);
1152   fPmeanVsBGcorr[2]->SetParameter(6,1.000000);
1153   fPmeanVsBGcorr[2]->SetParameter(7,-1.562500);
1154   fPmeanVsBGcorr[2]->SetParameter(8,0.239000);
1155   fPmeanVsBGcorr[2]->SetParameter(9,0.002054);
1156
1157   //Deuterons
1158   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);
1159   fPmeanVsBGcorr[3]->SetParameter(0,0.008672);
1160   fPmeanVsBGcorr[3]->SetParameter(1,-2.712343);
1161   fPmeanVsBGcorr[3]->SetParameter(2,0.997639);
1162   fPmeanVsBGcorr[3]->SetParameter(3,1.000000);
1163   fPmeanVsBGcorr[3]->SetParameter(4,0.039627);
1164   fPmeanVsBGcorr[3]->SetParameter(5,-2.768122);
1165   fPmeanVsBGcorr[3]->SetParameter(6,1.000000);
1166   fPmeanVsBGcorr[3]->SetParameter(7,-1.562500);
1167   fPmeanVsBGcorr[3]->SetParameter(8,0.174000);
1168   fPmeanVsBGcorr[3]->SetParameter(9,0.002189);
1169
1170   //Triton
1171   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);
1172   fPmeanVsBGcorr[4]->SetParameter(0,6.79641e-03);
1173   fPmeanVsBGcorr[4]->SetParameter(1,-1.92801e+00);
1174   fPmeanVsBGcorr[4]->SetParameter(2,1.000000);
1175   fPmeanVsBGcorr[4]->SetParameter(3,0.0);
1176   fPmeanVsBGcorr[4]->SetParameter(4,0.076);
1177   fPmeanVsBGcorr[4]->SetParameter(5,2.25779e-02);
1178
1179   //Helium-3
1180   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);
1181   fPmeanVsBGcorr[5]->SetParameter(0,0.024339);
1182   fPmeanVsBGcorr[5]->SetParameter(1,-2.922613);
1183   fPmeanVsBGcorr[5]->SetParameter(2,0.993761);
1184   fPmeanVsBGcorr[5]->SetParameter(3,1.000000);
1185   fPmeanVsBGcorr[5]->SetParameter(4,1.087549);
1186   fPmeanVsBGcorr[5]->SetParameter(5,-6.216154);
1187   fPmeanVsBGcorr[5]->SetParameter(6,1.000000);
1188   fPmeanVsBGcorr[5]->SetParameter(7,-1.562500);
1189   fPmeanVsBGcorr[5]->SetParameter(8,0.282000);
1190   fPmeanVsBGcorr[5]->SetParameter(9,0.009711);
1191
1192   //Helium-4
1193   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);
1194   fPmeanVsBGcorr[6]->SetParameter(0,2.34185e-02);
1195   fPmeanVsBGcorr[6]->SetParameter(1,-2.31200e+00);
1196   fPmeanVsBGcorr[6]->SetParameter(2,1.000000);
1197   fPmeanVsBGcorr[6]->SetParameter(3,0.0);
1198   fPmeanVsBGcorr[6]->SetParameter(4,0.198);
1199   fPmeanVsBGcorr[6]->SetParameter(5,9.9226e-03);
1200
1201   for(Int_t i=7;i<14;i++) {
1202     fPmeanVsBGcorr[i]=(TF2 *)fPmeanVsBGcorr[i-7]->Clone();
1203     fPmeanVsBGcorr[i]->SetName(name_fPmeanVsBGcorr[i]);
1204   }
1205
1206   for(Int_t i=0;i<14;i++) {
1207     fPmeanVsBGcorr[i]->SetNpx(fPmeanVsBGcorr[i]->GetNpx()*100.0);
1208   }
1209
1210   return;
1211
1212 }
1213