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 "TGraph.h"
  20 #include "TProfile.h"
  21 #include "AliESDtrackCuts.h"
  22 #include "AliAnalysisManager.h"
  23 #include "TFile.h"
  24
  25 ClassImp(AliAnalysisNucleiMass)//...
  26
  27 //_____________________________________________________________________________
  28 AliAnalysisNucleiMass::AliAnalysisNucleiMass():
  29   AliAnalysisTaskSE(),
  30 //Centrality(NULL),
  31   FilterBit(16),
  32 //  EtaLimit(NULL),
  33   DCAxyCut(0.1),
  34   DCAzCut(1000.),
  35   NsigmaTpcCut(2.0),
  36   NminTpcCluster(0),
  37   iTrdCut(0),
  38   kSignalCheck(1),
  39   iMtof(1),
  40   iBconf(0),
  41   kTOF(0),
  42   fAOD(NULL),
  43   fESD(NULL),
  44   fEvent(NULL),
  45   fPIDResponse(NULL)
  46 {
  47   fList[0]=new TList();
  48   fList[0]->SetName("results");
  49
  50   fList[1]=new TList();
  51   fList[1]->SetName("results2");
  52 }
  53 //______________________________________________________________________________
  54 AliAnalysisNucleiMass::AliAnalysisNucleiMass(const char *name):
  55   AliAnalysisTaskSE(name),
  56   //Centrality(NULL),
  57   FilterBit(16),
  58   //EtaLimit(NULL),
  59   DCAxyCut(0.1),
  60   DCAzCut(1000.),
  61   NsigmaTpcCut(2.0),
  62   NminTpcCluster(0),
  63   iTrdCut(0),
  64   kSignalCheck(1),
  65   iMtof(1),
  66   iBconf(0),
  67   kTOF(0),
  68   fAOD(NULL),
  69   fESD(NULL),
  70   fEvent(NULL),
  71   fPIDResponse(NULL)
  72 {
  73   fList[0]=new TList();
  74   DefineOutput(1, TList::Class());
  75   fList[0]->SetName("results");
  76
  77   fList[1]=new TList();
  78   DefineOutput(2, TList::Class());
  79   fList[1]->SetName("results2");
  80 }
  81 //_____________________________________________________________________________
  82 AliAnalysisNucleiMass::~AliAnalysisNucleiMass()
  83 {
  84   if(fList[0]) delete fList[0];
  85   if(fList[1]) delete fList[1];
  86 }
  87 //______________________________________________________________________________
  88 void AliAnalysisNucleiMass::UserCreateOutputObjects()
  89 {
  90   Char_t namePart[nPart][30];
  91   snprintf(namePart[0],30,"e");
  92   snprintf(namePart[1],30,"#mu");
  93   snprintf(namePart[2],30,"#pi");
  94   snprintf(namePart[3],30,"K");
  95   snprintf(namePart[4],30,"p");
  96   snprintf(namePart[5],30,"d");
  97   snprintf(namePart[6],30,"t");
  98   snprintf(namePart[7],30,"He3");
  99   snprintf(namePart[8],30,"He4");
 100
 101   Char_t name[nSpec][30];
 102   snprintf(name[0],20,"e^{+}");
 103   snprintf(name[1],20,"#mu^{+}");
 104   snprintf(name[2],20,"#pi^{+}");
 105   snprintf(name[3],20,"K^{+}");
 106   snprintf(name[4],20,"p");
 107   snprintf(name[5],20,"d");
 108   snprintf(name[6],20,"t");
 109   snprintf(name[7],20,"He3");
 110   snprintf(name[8],20,"He4");
 111   snprintf(name[9],20,"e^{-}");
 112   snprintf(name[10],20,"#mu^{-}");
 113   snprintf(name[11],20,"#pi^{-}");
 114   snprintf(name[12],20,"K^{-}");
 115   snprintf(name[13],20,"#bar{p}");
 116   snprintf(name[14],20,"#bar{d}");
 117   snprintf(name[15],20,"#bar{t}");
 118   snprintf(name[16],20,"#bar{He3}");
 119   snprintf(name[17],20,"#bar{He4}");
 120
 121   Double_t binPt[nbin+1];
 122   for(Int_t i=0;i<nbin+1;i++) {
 123     binPt[i]=0.4+0.1*i;
 124   }
 125
 126   Char_t name_nbin[nbin][200];
 127   for(Int_t j=0;j<nbin;j++) {
 128     snprintf(name_nbin[j],200,"%.1f<Pt<%.1f",binPt[j],binPt[j+1]);
 129   }
 130
 131   for(Int_t iB=0;iB<nBconf;iB++) {
 132
 133     htemp[iB] = new TH1F("htemp","htemp (avoid the problem with the empty list...);B field",20,-10,10);
 134
 135     hCentrality[iB][0] = new TH1F("hCentrality_Selected","Centrality (selected events);centrality(%)",20,0,100);
 136     hCentrality[iB][1] = new TH1F("hCentrality_Analyzed","Centrality (analyzed events);centrality (%)",20,0,100);
 137
 138     hZvertex[iB][0] = new TH1F("hZvertex_Selected","Vertex distribution of selected events;z vertex (cm)",240,-30,30);
 139     hZvertex[iB][1] = new TH1F("hZvertex_Analyzed","Vertex distribution of analyzed events;z vertex (cm)",240,-30,30);
 140
 141     hEta[iB] = new TH1F("hEta_Analyzed","|#eta| distribution after the track cuts;|#eta|",100,0.0,1.0);
 142
 143     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)
 144
 145     Int_t hbins[2];
 146     if(kSignalCheck>1) {hbins[0]=100; hbins[1]=90;}
 147     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 148     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);
 149
 150     hNTpcCluster[iB] = new TH1F("hNTpcCluster","Number of the TPC clusters after the track cuts;n_{cl}^{TPC}",300,0,300);
 151
 152     hNTrdSlices[iB] = new TH1F("hNTrdSlices","Number of the TRD slices after the track cuts;n_{slices}^{TRD}",40,0,40);
 153
 154     if(kSignalCheck==1) {hbins[0]=500; hbins[1]=2000;}
 155     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 156     else if(kSignalCheck==2) {hbins[0]=100; hbins[1]=500;}
 157     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);
 158     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);
 159
 160     Char_t name_hDeDxExp[nPart][200];
 161     Char_t title_hDeDxExp[nPart][200];
 162     for(Int_t i=0;i<nPart;i++) {
 163       snprintf(name_hDeDxExp[i],200,"hDeDxExp_%s",namePart[i]);
 164       snprintf(title_hDeDxExp[i],200,"Expected dE/dx of %s in the TPC;p/|z| (GeV/c);dE/dx_{TPC} (a.u.)",namePart[i]);
 165       hDeDxExp[iB][i] = new TProfile(name_hDeDxExp[i],title_hDeDxExp[i],500,0,5,0,1000,"");
 166     }
 167
 168     Char_t name_fNsigmaTpc[nPart][200];
 169     Char_t title_fNsigmaTpc[nPart][200];
 170     if(kSignalCheck==1) {hbins[0]=100; hbins[1]=100;}
 171     else {hbins[0]=1; hbins[1]=1;}
 172     for(Int_t i=0;i<nPart;i++) {
 173       snprintf(name_fNsigmaTpc[i],200,"NsigmaTpc_%s",namePart[i]);
 174       snprintf(title_fNsigmaTpc[i],200,"NsigmaTpc_%s;p_{TPC}/|z| (GeV/c);n_{#sigma_{TPC}}^{%s}",namePart[i],namePart[i]);
 175       fNsigmaTpc[iB][i] = new TH2F(name_fNsigmaTpc[i],title_fNsigmaTpc[i],hbins[0],0,5,hbins[1],-5,5);
 176     }
 177
 178     if(kSignalCheck>1) {hbins[0]=100; hbins[1]=100;}
 179     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 180     Char_t name_fNsigmaTpc_kTOF[nSpec][200];
 181     Char_t title_fNsigmaTpc_kTOF[nSpec][200];
 182     for(Int_t i=0;i<nSpec;i++) {
 183       snprintf(name_fNsigmaTpc_kTOF[i],200,"NsigmaTpc_%s_kTOF",name[i]);
 184       snprintf(title_fNsigmaTpc_kTOF[i],200,"NsigmaTpc_kTOF_%s in DCAxyCut;p_{T}/|z| (GeV/c);n_{#sigma_{TPC}}^{%s}",name[i],name[i]);
 185       fNsigmaTpc_kTOF[iB][i] = new TH2F(name_fNsigmaTpc_kTOF[i],title_fNsigmaTpc_kTOF[i],hbins[0],0,5,hbins[1],-5,5);
 186     }
 187
 188     if(kSignalCheck==1) {hbins[0]=1000; hbins[1]=1300;}
 189     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 190     else if(kSignalCheck==2) {hbins[0]=100; hbins[1]=260;}
 191     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);
 192     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);
 193
 194     Char_t name_hBetaExp[nPart][200];
 195     Char_t title_hBetaExp[nPart][200];
 196     for(Int_t i=0;i<nPart;i++) {
 197       snprintf(name_hBetaExp[i],200,"hBetaTofVsP_Exp_%s",namePart[i]);
 198       snprintf(title_hBetaExp[i],200,"Expected #beta_{TOF} vs p/|z| of %s;p/|z| (GeV/c); #beta_{TOF}",namePart[i]);
 199       hBetaExp[iB][i] = new TProfile(name_hBetaExp[i],title_hBetaExp[i],400,0,5,0.4,1.05,"");
 200     }
 201
 202     Char_t name_fNsigmaTof[nPart][200];
 203     Char_t title_fNsigmaTof[nPart][200];
 204     if(kSignalCheck==1) {hbins[0]=100; hbins[1]=100;}
 205     else {hbins[0]=1; hbins[1]=1;}
 206     for(Int_t i=0;i<nPart;i++) {
 207       snprintf(name_fNsigmaTof[i],200,"NsigmaTof_%s",namePart[i]);
 208       snprintf(title_fNsigmaTof[i],200,"NsigmaTof_%s;p_{T}/|z| (GeV/c);n_{#sigma_{TOF}}^{%s}",namePart[i],namePart[i]);
 209       fNsigmaTof[iB][i] = new TH2F(name_fNsigmaTof[i],title_fNsigmaTof[i],hbins[0],0,5,hbins[1],-5,5);
 210     }
 211
 212     Char_t name_fNsigmaTof_DcaCut[nSpec][200];
 213     Char_t title_fNsigmaTof_DcaCut[nSpec][200];
 214     if(kSignalCheck>1) {hbins[0]=100; hbins[1]=100;}
 215     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 216     for(Int_t i=0;i<nSpec;i++) {
 217       snprintf(name_fNsigmaTof_DcaCut[i],200,"NsigmaTof_DcaCut_%s",name[i]);
 218       snprintf(title_fNsigmaTof_DcaCut[i],200,"NsigmaTof_%s with DCAxyCut;p_{T}/|z| (GeV/c);n_{#sigma_{TOF}}^{%s}",name[i],name[i]);
 219       fNsigmaTof_DcaCut[iB][i] = new TH2F(name_fNsigmaTof_DcaCut[i],title_fNsigmaTof_DcaCut[i],hbins[0],0,5,hbins[1],-5,5);
 220     }
 221
 222     if(kSignalCheck==1) {hbins[0]=8000; hbins[1]=100;}
 223     else {hbins[0]=1; hbins[1]=1;}
 224     fM2vsPt_NoTpcCut[iB][0][0] = new TH2F("fM2vsPt_NoTpcCut_pos","m^{2}/z^{2}_{TOF} vs p_{T}/|z| (positive charge);m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4});p_{T}/|z| (GeV/c)",hbins[0],0,10,hbins[1],0,5);
 225     fM2vsPt_NoTpcCut[iB][0][1] = new TH2F("fM2vsPt_NoTpcCut_neg","m^{2}/z^{2}_{TOF} vs p_{T}/|z| (negative charge);m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4});p_{T}/|z| (GeV/c)",hbins[0],0,10,hbins[1],0,5);
 226
 227     if(kSignalCheck==1) {hbins[0]=8000; hbins[1]=100;}
 228     else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 229     else if(kSignalCheck==2) {hbins[0]=1000; hbins[1]=100;}
 230     fM2vsPt_NoTpcCut[iB][1][0] = new TH2F("fM2vsPt_NoTpcCut_DCAxyCut_pos","m^{2}/z^{2}_{TOF} vs p_{T}/|z| (positive charge) with DCAxy cut;m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4});p_{T}/|z| (GeV/c)",hbins[0],0,10,hbins[1],0,5);
 231     fM2vsPt_NoTpcCut[iB][1][1] = new TH2F("fM2vsPt_NoTpcCut_DCAxyCut_neg","m^{2}/z^{2}_{TOF} vs p_{T}/|z| (negative charge) with DCAxy cut;m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4});p_{T}/|z| (GeV/c)",hbins[0],0,10,hbins[1],0,5);
 232
 233     Char_t name_fM2vsPt[2][18][300];
 234     Char_t title_fM2vsPt[2][18][300];
 235
 236     for(Int_t i=0;i<nSpec;i++) {
 237       snprintf(name_fM2vsPt[0][i],300,"fM2vsPt_%s",name[i]);
 238       snprintf(title_fM2vsPt[0][i],300,"m^{2}/z^{2}_{TOF} vs p_{T}/|z| of %s with a NsigmaTpcCut;m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4});p_{T}/|z| (GeV/c)",name[i]);
 239
 240       snprintf(name_fM2vsPt[1][i],300,"fM2vsPt_%s_DCAxyCut",name[i]);
 241       snprintf(title_fM2vsPt[1][i],300,"m^{2}/z^{2}_{TOF} vs p_{T}/|z| of %s with a NsigmaTpcCut and with the DCAxy cut;m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4});p_{T}/|z| (GeV/c)",name[i]);
 242
 243       if(kSignalCheck==1) {hbins[0]=8000; hbins[1]=100;}
 244       else {hbins[0]=1; hbins[1]=1;}
 245       fM2vsPt[iB][0][i] = new TH2F(name_fM2vsPt[0][i],title_fM2vsPt[0][i],hbins[0],0,10,hbins[1],0,5);
 246
 247       if(kSignalCheck==1) {hbins[0]=8000; hbins[1]=100;}
 248       else if(kSignalCheck==0) {hbins[0]=1; hbins[1]=1;}
 249       else if(kSignalCheck==2) {hbins[0]=1000; hbins[1]=100;}
 250       fM2vsPt[iB][1][i] = new TH2F(name_fM2vsPt[1][i],title_fM2vsPt[1][i],hbins[0],0,10,hbins[1],0,5);
 251     }
 252
 253     if(kSignalCheck==1) {hbins[0]=4000; hbins[1]=1000;}
 254     else {hbins[0]=1; hbins[1]=1;}
 255     fM2vsZ[iB][0] = new TH2F("fM2vsZ","m^{2}/z^{2}_{TOF} vs z_{TPC} Integrated p_{T};z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
 256     fM2vsZ[iB][1] = new TH2F("fM2vsZ_0.5pT1.0","m^{2}/z^{2}_{TOF} vs z_{TPC} 0.5<pT<1.0;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
 257     fM2vsZ[iB][2] = new TH2F("fM2vsZ_1.0pT1.5","m^{2}/z^{2}_{TOF} vs z_{TPC} 1.0<pT<1.5;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
 258     fM2vsZ[iB][3] = new TH2F("fM2vsZ_1.5pT2.0","m^{2}/z^{2}_{TOF} vs z_{TPC} 1.5<pT<2.0;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
 259     fM2vsZ[iB][4] = new TH2F("fM2vsZ_2.0pT2.5","m^{2}/z^{2}_{TOF} vs z_{TPC} 2.0<pT<2.5;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
 260     fM2vsZ[iB][5] = new TH2F("fM2vsZ_2.5pT3.0","m^{2}/z^{2}_{TOF} vs z_{TPC} 2.5<pT<3.0;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
 261     fM2vsZ[iB][6] = new TH2F("fM2vsZ_3.0pT3.5","m^{2}/z^{2}_{TOF} vs z_{TPC} 3.0<pT<3.5;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
 262     fM2vsZ[iB][7] = new TH2F("fM2vsZ_3.5pT4.0","m^{2}/z^{2}_{TOF} vs z_{TPC} 3.5<pT<4.0;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
 263     fM2vsZ[iB][8] = new TH2F("fM2vsZ_4.0pT4.5","m^{2}/z^{2}_{TOF} vs z_{TPC} 4.0<pT<4.5;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
 264     fM2vsZ[iB][9] = new TH2F("fM2vsZ_4.5pT5.0","m^{2}/z^{2}_{TOF} vs z_{TPC} 2.0<pT<2.5;z_{TPC};m^{2}/z^{2}_{TOF} (GeV^{2}/c^{4})",hbins[0],-4,4,hbins[1],0,10);
 265
 266     Char_t name_hDCAxy[18][nbin][200];
 267     Char_t title_hDCAxy[18][nbin][200];
 268     Char_t name_hDCAz[18][nbin][200];
 269     Char_t title_hDCAz[18][nbin][200];
 270     for(Int_t iS=0;iS<nSpec;iS++) {
 271       for(Int_t j=0;j<nbin;j++) {
 272         snprintf(name_hDCAxy[iS][j],200,"hDCAxy_%s_%s",name[iS],name_nbin[j]);
 273         snprintf(title_hDCAxy[iS][j],200,"hDCAxy_%s_%s;DCA_{xy} (cm)",name[iS],name_nbin[j]);
 274         hDCAxy[iB][iS][j] = new TH1D(name_hDCAxy[iS][j],title_hDCAxy[iS][j],875,-3.5,3.5);
 275
 276         snprintf(name_hDCAz[iS][j],200,"hDCAz_%s_%s",name[iS],name_nbin[j]);
 277         snprintf(title_hDCAz[iS][j],200,"hDCAz_%s_%s;DCA_{z} (cm)",name[iS],name_nbin[j]);
 278         hDCAz[iB][iS][j] = new TH1D(name_hDCAz[iS][j],title_hDCAz[iS][j],875,-3.5,3.5);
 279       }
 280     }
 281
 282     Char_t name_hM2CutDCAxy[18][nbin][200];
 283     Char_t title_hM2CutDCAxy[18][nbin][200];
 284     Char_t name_hM2CutGroundDCAxy[18][nbin][200];
 285     Char_t title_hM2CutGroundDCAxy[18][nbin][200];
 286     for(Int_t iS=0;iS<nSpec;iS++) {
 287       for(Int_t j=0;j<nbin;j++) {
 288         snprintf(name_hM2CutDCAxy[iS][j],200,"hM2_CutDCAxy_%s_%s",name[iS],name_nbin[j]);
 289         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]);
 290         snprintf(name_hM2CutGroundDCAxy[iS][j],200,"hM2_GroundCatDCAxy_%s_%s",name[iS],name_nbin[j]);
 291         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]);
 292       }
 293     }
 294
 295     const Int_t BinM2pT[nPart]={1,1,600,250,500,500,1000,400,600};
 296     const Double_t RangeM2min[nPart]={0.0,0.0,-0.1,0.0,0.0,0.0,0.0,0.0,0.0};
 297     const Double_t RangeM2max[nPart]={1.0,1.0,0.5,2.0,4.0,6.0,12.0,4.0,6.0};
 298
 299     for(Int_t iS=0;iS<nPart;iS++) {
 300       for(Int_t j=0;j<nbin;j++) {
 301
 302         hM2CutDCAxy[iB][iS][j] = new TH1D(name_hM2CutDCAxy[iS][j],title_hM2CutDCAxy[iS][j],BinM2pT[iS],RangeM2min[iS],RangeM2max[iS]);
 303         hM2CutGroundDCAxy[iB][iS][j] = new TH1D(name_hM2CutGroundDCAxy[iS][j],title_hM2CutGroundDCAxy[iS][j],BinM2pT[iS],RangeM2min[iS],RangeM2max[iS]);
 304
 305         hM2CutDCAxy[iB][iS+nPart][j] = new TH1D(name_hM2CutDCAxy[iS+nPart][j],title_hM2CutDCAxy[iS+nPart][j],BinM2pT[iS],RangeM2min[iS],RangeM2max[iS]);
 306         hM2CutGroundDCAxy[iB][iS+nPart][j] = new TH1D(name_hM2CutGroundDCAxy[iS+nPart][j],title_hM2CutGroundDCAxy[iS+nPart][j],BinM2pT[iS],RangeM2min[iS],RangeM2max[iS]);
 307       }
 308     }
 309
 310     //Parameterizations:
 311     fPmeanVsPexp[0]=new TF1("fPmeanVsPexp_p","[2]-[0]*TMath::Exp(-(TMath::Max(x,[3])*[1]))",0,20);
 312     fPmeanVsPexp[1]=new TF1("fPmeanVsPexp_d","[2]-[0]*TMath::Exp(-(TMath::Max(x,[3])*[1]))",0,20);
 313     fPmeanVsPexp[2]=new TF1("fPmeanVsPexp_He3","[2]-[0]*TMath::Exp(-(TMath::Max(x,[3])*[1]))",0,20);
 314
 315     Double_t fpars_p[4]={5.14500484596484148e-03,9.74729863202270397e-01,0.0,1.00607413672776569e+00};
 316     Double_t fpars_d[4]={3.16023942908439243e-02,1.24005027514358490e+00,-1.50000000000000003e-03,1.40607413672776560e+00};
 317     Double_t fpars_He3[4]={2.73329079591698026e-02,1.53005942367188852e+00,-4.10231310888738848e-03,1.20607413672776564e+00};
 318     fPmeanVsPexp[0]->SetParameters(fpars_p);
 319     fPmeanVsPexp[1]->SetParameters(fpars_d);
 320     fPmeanVsPexp[2]->SetParameters(fpars_He3);
 321
 322     /*Char_t title_Xaxis[3][200];
 323     Char_t title_Yaxis[3][200];
 324     snprintf(title_Xaxis[0],200,"p(t_{exp}^{%s})",namePart[4]);
 325     snprintf(title_Yaxis[0],200,"p(t_{TOF})-p(t_{exp}^{%s})/p(t_{exp}^{%s})",namePart[4],namePart[4]);
 326     snprintf(title_Xaxis[1],200,"p(t_{exp}^{%s})",namePart[5]);
 327     snprintf(title_Yaxis[1],200,"p(t_{TOF})-p(t_{exp}^{%s})/p(t_{exp}^{%s})",namePart[5],namePart[5]);
 328     snprintf(title_Xaxis[2],200,"p(t_{exp}^{%s})",namePart[7]);
 329     snprintf(title_Yaxis[2],200,"p(t_{TOF})-p(t_{exp}^{%s})/p(t_{exp}^{%s})",namePart[7],namePart[7]);
 330     for(Int_t i=0;i<3;i++){
 331       fPmeanVsPexp[i]->GetXaxis()->SetTitle(title_Xaxis[i]);
 332       fPmeanVsPexp[i]->GetYaxis()->SetTitle(title_Yaxis[i]);
 333       fPmeanVsPexp[i]->SetTitle("Parameterization calculated with Monte Carlo (LHC13d15)");
 334       }*/
 335     //end parameterizations
 336
 337     Char_t name_fPmeanVsBetaGamma[18][200];
 338     Char_t title_fPmeanVsBetaGamma[18][200];
 339
 340     hbins[0]=200; hbins[1]=200;
 341     for(Int_t iS=0;iS<nSpec;iS++) {
 342       snprintf(name_fPmeanVsBetaGamma[iS],200,"fPmeanVsPvtx_%s",name[iS]);
 343       snprintf(title_fPmeanVsBetaGamma[iS],200,"<p>/p_{vtx} vs #beta#gamma of %s;p_{vtx}/m_{%s};<p>_{%s}/p_{vtx}",name[iS],name[iS],name[iS]);
 344       fPmeanVsBetaGamma[iB][iS]=new TH2F(name_fPmeanVsBetaGamma[iS],title_fPmeanVsBetaGamma[iS],hbins[0],0,10,hbins[1],0.8,1.2);
 345     }
 346
 347     Char_t name_prPmeanVsBetaGamma[18][200];
 348     Char_t title_prPmeanVsBetaGamma[18][200];
 349
 350     for(Int_t iS=0;iS<nSpec;iS++) {
 351       snprintf(name_prPmeanVsBetaGamma[iS],200,"prPmeanVsPvtx_%s",name[iS]);
 352       snprintf(title_prPmeanVsBetaGamma[iS],200,"<p>/p_{vtx} vs #beta#gamma of %s;p_{vtx}/m_{%s};<p>_{%s}/p_{vtx}",name[iS],name[iS],name[iS]);
 353       prPmeanVsBetaGamma[iB][iS]=new TProfile(name_prPmeanVsBetaGamma[iS],title_prPmeanVsBetaGamma[iS],hbins[0],0,10,0.8,1.2,"");
 354     }
 355
 356     fList[iB]->Add(htemp[iB]);
 357     for(Int_t i=0;i<2;i++) fList[iB]->Add(hCentrality[iB][i]);
 358     for(Int_t i=0;i<2;i++) fList[iB]->Add(hZvertex[iB][i]);
 359     fList[iB]->Add(hEta[iB]);
 360     fList[iB]->Add(hPhi[iB]);
 361     fList[iB]->Add(fEtaPhi[iB]);
 362     fList[iB]->Add(hNTpcCluster[iB]);
 363     fList[iB]->Add(hNTrdSlices[iB]);
 364     for(Int_t i=0;i<2;i++) fList[iB]->Add(fdEdxVSp[iB][i]);
 365     for(Int_t i=0;i<nPart;i++) fList[iB]->Add(hDeDxExp[iB][i]);
 366     for(Int_t i=0;i<nPart;i++) fList[iB]->Add(fNsigmaTpc[iB][i]);
 367     for(Int_t i=0;i<nPart;i++) {
 368       if(kSignalCheck!=1)
 369         if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
 370       fList[iB]->Add(fNsigmaTpc_kTOF[iB][i]);
 371       fList[iB]->Add(fNsigmaTpc_kTOF[iB][i+nPart]);
 372     }
 373     for(Int_t i=0;i<2;i++) fList[iB]->Add(fBetaTofVSp[iB][i]);
 374     for(Int_t i=0;i<nPart;i++) fList[iB]->Add(hBetaExp[iB][i]);
 375     for(Int_t i=0;i<nPart;i++) fList[iB]->Add(fNsigmaTof[iB][i]);
 376     for(Int_t i=0;i<nPart;i++) {
 377       if(kSignalCheck!=1)
 378         if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
 379       fList[iB]->Add(fNsigmaTof_DcaCut[iB][i]);
 380       fList[iB]->Add(fNsigmaTof_DcaCut[iB][i+nPart]);
 381     }
 382     for(Int_t i=0;i<2;i++) fList[iB]->Add(fM2vsPt_NoTpcCut[iB][0][i]);
 383     for(Int_t i=0;i<2;i++) fList[iB]->Add(fM2vsPt_NoTpcCut[iB][1][i]);
 384     for(Int_t i=0;i<nPart;i++) {
 385       if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
 386       fList[iB]->Add(fM2vsPt[iB][0][i]);
 387       fList[iB]->Add(fM2vsPt[iB][0][i+nPart]);
 388     }
 389     for(Int_t i=0;i<nPart;i++){
 390       if(kSignalCheck!=1)
 391         if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
 392       fList[iB]->Add(fM2vsPt[iB][1][i]);
 393       fList[iB]->Add(fM2vsPt[iB][1][i+nPart]);
 394     }
 395     if(iMtof!=1) {
 396       for(Int_t i=0;i<nPart;i++){
 397         if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
 398         fList[iB]->Add(fPmeanVsBetaGamma[iB][i]);
 399         fList[iB]->Add(prPmeanVsBetaGamma[iB][i]);
 400         fList[iB]->Add(fPmeanVsBetaGamma[iB][i+nPart]);
 401         fList[iB]->Add(prPmeanVsBetaGamma[iB][i+nPart]);
 402       }
 403     }
 404     if(iMtof==8) for(Int_t i=0;i<3;i++) fList[iB]->Add(fPmeanVsPexp[i]);
 405     else if(iMtof==4) for(Int_t i=1;i<3;i++) fList[iB]->Add(fPmeanVsPexp[i]);
 406     for(Int_t i=0;i<10;i++) fList[iB]->Add(fM2vsZ[iB][i]);
 407     for(Int_t i=0;i<nPart;i++){
 408       if(kSignalCheck!=1)
 409         if(i<3 || i==6 || i==8) continue;//e,mu,pi,t,he4 excluded
 410       for(Int_t j=0;j<nbin;j++){
 411         fList[iB]->Add(hDCAxy[iB][i][j]);
 412         fList[iB]->Add(hDCAz[iB][i][j]);
 413         fList[iB]->Add(hM2CutDCAxy[iB][i][j]);
 414         fList[iB]->Add(hM2CutGroundDCAxy[iB][i][j]);
 415         fList[iB]->Add(hDCAxy[iB][i+nPart][j]);
 416         fList[iB]->Add(hDCAz[iB][i+nPart][j]);
 417         fList[iB]->Add(hM2CutDCAxy[iB][i+nPart][j]);
 418         fList[iB]->Add(hM2CutGroundDCAxy[iB][i+nPart][j]);
 419       }
 420     }
 421
 422     // Post output data.
 423     PostData(1, fList[0]);
 424     PostData(2, fList[1]);
 425
 426   }//end iB loop
 427 }
 428 //______________________________________________________________________________
 429 void AliAnalysisNucleiMass::UserExec(Option_t *)
 430 {
 431   // Main loop
 432   // Called for each event
 433
 434   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 435   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 436   if(!fAOD && !fESD){
 437     Printf("%s:%d AODEvent and ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 438     return;
 439   }
 440
 441   if(fESD) fEvent = fESD;
 442   else fEvent = fAOD;
 443
 444   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 445   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 446   fPIDResponse=inputHandler->GetPIDResponse();
 447
 448   //--------------------------Magnetic field polarity--------------------
 449   Double_t fBfield=fEvent->GetMagneticField();
 450   if(fBfield<0.0) iBconf=0;//B--
 451   else iBconf=1;//B++
 452   for(Int_t i=0;i<nBconf;i++) htemp[i]->Fill(fBfield);
 453
 454   //--------------------------Centrality--------------------------------
 455   Double_t v0Centr  = -10.;
 456   AliCentrality *centrality = fEvent->GetCentrality();
 457   if (centrality){
 458     v0Centr=centrality->GetCentralityPercentile("V0M"); // VZERO
 459   }
 460   hCentrality[iBconf][0]->Fill(v0Centr);
 461
 462   //-------------------------zVertex determination of event----------------
 463   Double_t zvtx = 9999.9;
 464   const AliVVertex* vtxEVENT = fEvent->GetPrimaryVertex();
 465   if(vtxEVENT->GetNContributors()>0) zvtx = vtxEVENT->GetZ();
 466
 467   hZvertex[iBconf][0]->Fill(zvtx);
 468
 469   //---------------------------EVENT CUTS-----------------------------
 470   if(TMath::Abs(zvtx) < 10.0 && v0Centr>Centrality[0] && v0Centr<Centrality[1]){
 471
 472     hCentrality[iBconf][1]->Fill(v0Centr);
 473     hZvertex[iBconf][1]->Fill(zvtx);
 474
 475     Int_t nTracks = fEvent->GetNumberOfTracks();
 476
 477     //----------------------loop on the TRACKS-----------------------------
 478     for(Int_t iT = 0; iT < nTracks; iT++) {
 479       AliVTrack* track = (AliVTrack *) fEvent->GetTrack(iT);
 480
 481       if (!track){
 482         continue;
 483       }
 484
 485      //For the geometrical cuts
 486       Double_t etaAbs = TMath::Abs(track->Eta());
 487
 488       Bool_t trkFlag = 0;
 489       trkFlag = ((AliAODTrack *) track)->TestFilterBit(FilterBit);
 490       //TestFilterBit(16) -- Standard Cuts with very loose DCA: GetStandardITSTPCTrackCuts2011(kFALSE) && SetMaxDCAToVertexXY(2.4) && SetMaxDCAToVertexZ(3.2) && SetDCaToVertex2D(kTRUE)
 491       //TestFilterBit(32) (STARDARD) -- Standard Cuts with very tight DCA cut ( 7sigma^primaries: 7*(0.0015+0.0050/pt^1.1) ) : GetStandardITSTPCTrackCuts2011().
 492
 493       //Cut on the Minumum Number of the TPC clusters
 494       Bool_t isMinTpcCluster=kFALSE;
 495       Int_t nTpcCluster=0;
 496       nTpcCluster=track->GetTPCNcls();
 497       if(nTpcCluster>NminTpcCluster) isMinTpcCluster=kTRUE;
 498
 499       //-------------------------------------start TRACK CUTS----------------------------------
 500       if ((track->Pt() < 0.2) || (etaAbs<EtaLimit[0]) || (etaAbs>EtaLimit[1]) || !trkFlag || !isMinTpcCluster)
 501         continue;
 502
 503       //Vertex determination
 504       Double_t b[2] = {-99., -99.};
 505       Double_t bCov[3] = {-99., -99., -99.};
 506       if (!track->PropagateToDCA(fEvent->GetPrimaryVertex(), fEvent->GetMagneticField(), 100., b, bCov))
 507         continue;
 508
 509       Double_t DCAxy = b[0];
 510       Double_t DCAz = b[1];
 511
 512       //Cut on the DCAz
 513       Bool_t isDCAzCut=kFALSE;
 514       if(DCAz<DCAzCut) isDCAzCut=kTRUE;
 515
 516       if(!isDCAzCut)
 517         continue;
 518
 519       //For the Tpc purity cut
 520       Double_t dedx = track->GetTPCsignal();
 521       if(dedx<10) continue;
 522
 523       Int_t nTrdSlices = track->GetNumberOfTRDslices();
 524       if(nTrdSlices<2 && iTrdCut==1) continue;
 525       if(nTrdSlices>0 && iTrdCut==2) continue;
 526
 527       //-------------------------------------end TRACK CUTS----------------------------------
 528
 529       //Track info:
 530       Double_t phi= track->Phi();
 531
 532       hEta[iBconf]->Fill(etaAbs);
 533       hPhi[iBconf]->Fill(phi);
 534       fEtaPhi[iBconf]->Fill(etaAbs,phi);
 535       hNTpcCluster[iBconf]->Fill(nTpcCluster);
 536       hNTrdSlices[iBconf]->Fill(nTrdSlices);
 537
 538       Double_t charge = (Double_t)track->Charge();
 539       Double_t p = track->P();
 540       Double_t pt = track->Pt();
 541       Double_t tof  = track->GetTOFsignal()-fPIDResponse->GetTOFResponse().GetStartTime(p);
 542       Double_t pTPC = track->GetTPCmomentum();
 543       Double_t beta = 0.0;
 544       Double_t M2 = 999.9;
 545       Double_t Z2 = 999.9;
 546
 547       kTOF = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
 548
 549       //-----------------------------TPC info------------------------------
 550       Double_t nsigmaTPC[nPart];
 551       Double_t expdedx[nPart];
 552
 553       Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
 554       Int_t FlagPid = 0;
 555
 556       for(Int_t iS=0;iS<9;iS++){
 557         expdedx[iS] = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, (AliPID::EParticleType) iS, AliTPCPIDResponse::kdEdxDefault, kTRUE);
 558         hDeDxExp[iBconf][iS]->Fill(pTPC,expdedx[iS]);
 559         nsigmaTPC[iS] = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType) iS);
 560         fNsigmaTpc[iBconf][iS]->Fill(pTPC,nsigmaTPC[iS]);
 561         if(charge>0) {//positive particle
 562           if(kTOF && (TMath::Abs(DCAxy)<DCAxyCut)) fNsigmaTpc_kTOF[iBconf][iS]->Fill(pt,nsigmaTPC[iS]);
 563         }
 564         else {//negative particle
 565           if(kTOF && (TMath::Abs(DCAxy)<DCAxyCut)) fNsigmaTpc_kTOF[iBconf][iS+nPart]->Fill(pt,nsigmaTPC[iS]);
 566         }
 567
 568         //TPC identification:
 569         if(TMath::Abs(nsigmaTPC[iS])<NsigmaTpcCut) {
 570           FlagPid += ((Int_t)TMath::Power(2,iS));
 571         }
 572       }
 573
 574       if(charge>0) fdEdxVSp[iBconf][0]->Fill(pTPC,dedx);
 575       else fdEdxVSp[iBconf][1]->Fill(pTPC,dedx);
 576
 577       //-----------------------------TOF info------------------------------
 578
 579       Double_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.875612859,2.808921005,1.404195741,1.863689620};
 580
 581       //----------------------------------------kTOF available-----------------------------
 582
 583       if(kTOF) {
 584         Double_t exptimes[9];
 585         track->GetIntegratedTimes(exptimes);
 586         //Integrated times of the Nuclei:
 587         for(Int_t iN=5;iN<9;iN++) {
 588           exptimes[iN] = exptimes[4]*exptimes[4]*(massOverZ[iN]*massOverZ[iN]/p/p+1)/(massOverZ[4]*massOverZ[4]/p/p+1);
 589           exptimes[iN] = TMath::Sqrt(exptimes[iN]);
 590         }
 591
 592         beta=exptimes[0];
 593         beta=beta/tof;//beta = L/tof/c = t_e/tof
 594
 595         Double_t nsigmaTOF[9];
 596         for(Int_t iS=0;iS<9;iS++){
 597           nsigmaTOF[iS] = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType) iS);
 598           fNsigmaTof[iBconf][iS]->Fill(pt,nsigmaTOF[iS]);
 599           if(charge>0) {
 600             hBetaExp[iBconf][iS]->Fill(p,exptimes[0]/exptimes[iS]);
 601             if(TMath::Abs(DCAxy)<DCAxyCut) fNsigmaTof_DcaCut[iBconf][iS]->Fill(pt,nsigmaTOF[iS]);
 602           }
 603           else {
 604             hBetaExp[iBconf][iS+nPart]->Fill(p,exptimes[0]/exptimes[iS]);
 605             if(TMath::Abs(DCAxy)<DCAxyCut) fNsigmaTof_DcaCut[iBconf][iS+nPart]->Fill(pt,nsigmaTOF[iS]);
 606           }
 607         }
 608         if(charge>0) fBetaTofVSp[iBconf][0]->Fill(p,beta);
 609         else fBetaTofVSp[iBconf][1]->Fill(p,beta);
 610
 611         this->GetMassFromPvertex(beta,p,M2);
 612         this->GetZTpc(dedx,pTPC,M2,Z2);
 613
 614         //-----------------------------M2 as a function of expected times, if iMtof>1---------------------------------
 615         Double_t Mass2[9];
 616         if(iMtof>1) this->GetMassFromExpTimes(beta,exptimes,Mass2,iMtof,p,FlagPid,charge);
 617
 618         //-------------------------------Squared Mass TH2 distributions-----------------------
 619         if(charge>0) {
 620           //without TPC
 621           fM2vsPt_NoTpcCut[iBconf][0][0]->Fill(M2,pt);
 622           if(TMath::Abs(DCAxy)<DCAxyCut) fM2vsPt_NoTpcCut[iBconf][1][0]->Fill(M2,pt);
 623           //with TPC
 624           for(Int_t iS=0;iS<9;iS++) {
 625             //-----------------------------M2 as a function of expected times, if iMtof>1---------------------------------
 626             if(iMtof>1) {
 627               M2=999.9;
 628               M2=Mass2[iS];
 629             }
 630             //-----------------
 631             if(FlagPid & stdFlagPid[iS]) {
 632               fM2vsPt[iBconf][0][iS]->Fill(M2,pt);
 633               if(TMath::Abs(DCAxy)<DCAxyCut) fM2vsPt[iBconf][1][iS]->Fill(M2,pt);
 634             }
 635           }
 636         }
 637         else {//charge<0
 638           //without TPC
 639           fM2vsPt_NoTpcCut[iBconf][0][1]->Fill(M2,pt);
 640           if(TMath::Abs(DCAxy)<DCAxyCut) fM2vsPt_NoTpcCut[iBconf][1][1]->Fill(M2,pt);
 641            //with TPC
 642           for(Int_t iS=0;iS<9;iS++) {
 643             //-----------------------------M2 as a function of expected times, if iMtof>1---------------------------------
 644             if(iMtof>1) {
 645               M2=999.9;
 646               M2=Mass2[iS];
 647             }
 648             //-----------------
 649             if(FlagPid & stdFlagPid[iS]) {
 650               fM2vsPt[iBconf][0][iS+nPart]->Fill(M2,pt);
 651               if(TMath::Abs(DCAxy)<DCAxyCut) fM2vsPt[iBconf][1][iS+nPart]->Fill(M2,pt);
 652             }
 653           }
 654         }
 655
 656         //------------------------------start DCA and Squared Mass TH1 distributions-------------------------
 657         Double_t binPt[nbin+1];
 658         for(Int_t i=0;i<nbin+1;i++) {
 659           binPt[i]=0.4+i*0.1;
 660         }
 661
 662         if(charge>0) {
 663           for(Int_t iS=0;iS<9;iS++) {
 664             //-----------------------------M2 as a function of expected times, if iMtof>1---------------------------------
 665             if(iMtof>1) {
 666               M2=999.9;
 667               M2=Mass2[iS];
 668             }
 669             //-----------------
 670             if(FlagPid & stdFlagPid[iS]) {
 671               for(Int_t j=0;j<nbin;j++) {
 672                 if(pt>binPt[j] && pt<binPt[j+1]) {
 673                   hDCAxy[iBconf][iS][j]->Fill(DCAxy);
 674                   hDCAxy[iBconf][iS][j]->Fill(-DCAxy);
 675                   hDCAz[iBconf][iS][j]->Fill(DCAz);
 676                   hDCAz[iBconf][iS][j]->Fill(-DCAz);
 677                   if(TMath::Abs(DCAxy)<DCAxyCut) {
 678                     hM2CutDCAxy[iBconf][iS][j]->Fill(M2);
 679                   }
 680                   if(TMath::Abs(DCAxy+0.5)<DCAxyCut) {
 681                     hM2CutGroundDCAxy[iBconf][iS][j]->Fill(M2);
 682                   }
 683                   break;
 684                 }
 685               }//end loop on the pT bins (j)
 686             }
 687           }//end loop on the particle species (iS)
 688         }
 689         else {//charge<0
 690           for(Int_t iS=0;iS<9;iS++) {
 691             //-----------------------------M2 as a function of expected times, if iMtof>1---------------------------------
 692             if(iMtof>1) {
 693               M2=999.9;
 694               M2=Mass2[iS];
 695             }
 696             //-----------------
 697             if(FlagPid & stdFlagPid[iS]) {
 698               for(Int_t j=0;j<nbin;j++) {
 699                 if(pt>binPt[j] && pt<binPt[j+1]) {
 700                   hDCAxy[iBconf][iS+nPart][j]->Fill(DCAxy);
 701                   hDCAxy[iBconf][iS+nPart][j]->Fill(-DCAxy);
 702                   hDCAz[iBconf][iS+nPart][j]->Fill(DCAz);
 703                   hDCAz[iBconf][iS+nPart][j]->Fill(-DCAz);
 704                   if(TMath::Abs(DCAxy)<DCAxyCut) {
 705                     hM2CutDCAxy[iBconf][iS+nPart][j]->Fill(M2);
 706                   }
 707                   if(TMath::Abs(DCAxy+0.5)<DCAxyCut) {
 708                     hM2CutGroundDCAxy[iBconf][iS+nPart][j]->Fill(M2);
 709                   }
 710                   break;
 711                 }
 712               }//end loop on the pT bins (j)
 713             }
 714           }//end loop on the particle species (iS)
 715         }
 716
 717         //-------------------------------------------------M2/Z2 vs Z-------------------------
 718
 719
 720         Double_t binCutPt[10] = {0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0};
 721         Double_t Z=999.9;
 722         if(Z2>0) Z=TMath::Sqrt(Z2);
 723
 724         fM2vsZ[iBconf][0]->Fill(charge*TMath::Sqrt(Z2),M2);
 725         for(Int_t i=1;i<10;i++) {
 726           if(pt>binCutPt[i-1] && pt<binCutPt[i]){
 727             fM2vsZ[iBconf][i]->Fill(charge*Z,M2);
 728             break;
 729           }
 730         }
 731
 732
 733
 734       }//end kTOF available
 735     }//end track loop
 736   }//end loop on the events
 737 }
 738
 739 //_____________________________________________________________________________
 740 void AliAnalysisNucleiMass::Terminate(Option_t *)
 741 {
 742   // Terminate loop
 743   Printf("Terminate()");
 744 }
 745 //_____________________________________________________________________________
 746 void AliAnalysisNucleiMass::GetMassFromPvertex(Double_t beta, Double_t p, Double_t &M2) {
 747
 748   M2 = p*p*(1-beta*beta)/(beta*beta);
 749
 750   return;
 751
 752 }
 753 //____________________________________________________________________________________________________________
 754 void AliAnalysisNucleiMass::GetMassFromExpTimes(Double_t beta, Double_t *IntTimes, Double_t *Mass2, Int_t iCorr, Double_t pVtx, Int_t FlagPid, Double_t charge) {
 755
 756   // 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
 757   // In this way m_tof = mPDG only if tof=t_exp
 758
 759   Double_t massOverZ[9] = {0.000511,0.105658,0.139570,0.493677,0.938272,1.875612859,2.808921005,1.404195741,1.863689620};
 760
 761   Double_t beta2Exp[9];
 762   Double_t p2Exp[9];
 763
 764   Double_t pExp[9];
 765   Double_t CorrFactor=0.0;
 766
 767   Int_t stdFlagPid[9] = {1,2,4,8,16,32,64,128,256};//e,#mu,#pi,K,p,d,t,3He,4He
 768
 769   for(Int_t iS=0;iS<9;iS++) {
 770     beta2Exp[iS]=IntTimes[0]/IntTimes[iS];//beta = L/tof*c = t_e/tof
 771     beta2Exp[iS]=beta2Exp[iS]*beta2Exp[iS];
 772     if((1-beta2Exp[iS])==0) {
 773       Mass2[iS]=999.9;
 774       continue;
 775     }
 776     p2Exp[iS]=massOverZ[iS]*massOverZ[iS]*beta2Exp[iS]/(1-beta2Exp[iS]);
 777
 778     //--------------------for MC corrections
 779     if(p2Exp[iS]<0) {
 780       Mass2[iS]=999.9;
 781       continue;
 782     }
 783     pExp[iS]=TMath::Sqrt(p2Exp[iS]);
 784
 785     CorrFactor=0.0;
 786     if(iCorr & 12) {//iCorr==4 || iCorr==8
 787       if(iCorr==8 && iS==4) CorrFactor=fPmeanVsPexp[0]->Eval(pExp[iS]);
 788
 789       if(iS==5) CorrFactor=fPmeanVsPexp[1]->Eval(pExp[iS]);
 790       else if(iS==7) CorrFactor=fPmeanVsPexp[2]->Eval(pExp[iS]);
 791       CorrFactor=pExp[iS]*CorrFactor;
 792       pExp[iS]=pExp[iS]+CorrFactor;//CorrFactor is negative so pExp(Corrected)<pExp
 793     }
 794     p2Exp[iS]=pExp[iS]*pExp[iS];
 795     //------------
 796     Mass2[iS]=p2Exp[iS]*(1-beta*beta)/(beta*beta);
 797
 798     //------------
 799     if(FlagPid & stdFlagPid[iS]) {
 800       if(charge>0){
 801         fPmeanVsBetaGamma[iBconf][iS]->Fill(pVtx/massOverZ[iS],pExp[iS]/pVtx);
 802         prPmeanVsBetaGamma[iBconf][iS]->Fill(pVtx/massOverZ[iS],pExp[iS]/pVtx);
 803       }
 804       else {
 805         fPmeanVsBetaGamma[iBconf][iS+nPart]->Fill(pVtx/massOverZ[iS],pExp[iS]/pVtx);
 806         prPmeanVsBetaGamma[iBconf][iS+nPart]->Fill(pVtx/massOverZ[iS],pExp[iS]/pVtx);
 807       }
 808     }
 809   }//end loop on the particle species
 810
 811   return;
 812
 813 }
 814 //_________________________________________________________________________________________________________________________
 815 void AliAnalysisNucleiMass::GetZTpc(Double_t dedx, Double_t pTPC, Double_t M2, Double_t &Z2) {
 816
 817   //z^2_tpc = dedx^{Tpc} / dedx^{exp,Tof}_{z=1}
 818
 819   Z2=999.9;
 820
 821   Double_t M=999.9;
 822   Double_t pTPC_pr=999.9;//rescaling of the pTPC for the proton
 823   Double_t expdedx_Tof=999.9;
 824
 825   if(M2>0) {
 826     M=TMath::Sqrt(M2);
 827     pTPC_pr=pTPC*0.938272/M;
 828     expdedx_Tof=fPIDResponse->GetTPCResponse().GetExpectedSignal(pTPC_pr,AliPID::kProton);
 829     if((dedx/expdedx_Tof)<0) return;
 830     Z2=TMath::Power(dedx/expdedx_Tof,0.862);
 831   }
 832
 833   return;
 834 }