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