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