PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////
  17 // AliAnalysisTaskSE for the extraction of the various histograms to
  18 // study the pt spectra of identified hadrons:
  19 // - log(dEdx)-log(dEdxBB) distributions for pions, kaons and protons in pt bins
  20 // - Pt distributions of pions, kaons and protons with nSigma PID
  21 // Authors:
  22 // E. Biolcati, biolcati@to.infn.it
  23 // L. Milano, milano@to.infn.it
  24 // F. Prino, prino@to.infn.it
  25 ///////////////////////////////////////////////////////////////////////////
  26
  27 #include <TH1F.h>
  28 #include <TF1.h>
  29 #include <TRandom3.h>
  30 #include <TH2F.h>
  31 #include <TChain.h>
  32 #include <TNtuple.h>
  33 #include <TParticle.h>
  34 #include <Rtypes.h>
  35 #include "AliAnalysisTaskSE.h"
  36 #include "AliAnalysisManager.h"
  37 #include "AliAnalysisDataContainer.h"
  38 #include "AliESDEvent.h"
  39 #include "AliESDInputHandler.h"
  40 #include "AliESDtrack.h"
  41 #include "AliStack.h"
  42 #include "AliMCEventHandler.h"
  43 #include "AliMCEvent.h"
  44 #include "AliPhysicsSelection.h"
  45 #include "AliAnalysisTaskSEITSsaSpectra.h"
  46 #include "AliESDtrackCuts.h"
  47 #include "AliCentrality.h"
  48 #include "AliMultiplicity.h"
  49 #include "AliESDUtils.h"
  50 #include "AliITSPIDResponse.h"
  51
  52 ClassImp(AliAnalysisTaskSEITSsaSpectra)
  53 /* $Id$ */
  54
  55 //________________________________________________________________________
  56 AliAnalysisTaskSEITSsaSpectra::AliAnalysisTaskSEITSsaSpectra():
  57 AliAnalysisTaskSE("taskITSsaSpectra"),
  58   fESD(0),
  59   fOutput(0),
  60   fListCuts(0),
  61   fHistNEvents(0),
  62   fHistMult(0),
  63   fHistCen(0),
  64   fHistNTracks(0),
  65   fHistNTracksPos(0),
  66   fHistNTracksNeg(0),
  67   fHistDEDX(0),
  68   fHistDEDXdouble(0),
  69   fHistBeforeEvSel(0),
  70   fHistAfterEvSel(0),
  71   fDCAxyCutFunc(0x0),
  72   fDCAzCutFunc(0x0),
  73   fITSPIDResponse(0),
  74   fMinSPDPts(1),
  75   fMinNdEdxSamples(3),
  76   fMindEdx(0.),
  77   fMinNSigma(1.5),
  78   fMaxY(0.5),
  79   fMaxChi2Clu(2.5),
  80   fNSigmaDCAxy(7.),
  81   fNSigmaDCAz(7.),
  82   fEtaRange(0.8),
  83   fLowMult(-1),
  84   fUpMult(-1),
  85   fLowCentrality(-1.0),
  86   fUpCentrality(-1.0),
  87   fMultEstimator(0),
  88   fHImode(0),
  89   fYear(2010),
  90   fMC(kFALSE),
  91   fSmearMC(kFALSE),
  92   fSmearP(0.),
  93   fSmeardEdx(0.),
  94   fRandGener(0),
  95   fFillNtuple(kFALSE),
  96   fLowEnergypp(kFALSE),
  97   fNtupleNSigma(0),
  98   fNtupleMC(0)
  99 {
 100   // Constructor
 101   Double_t xbins[kNbins+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};
 102   for(Int_t iBin=0; iBin<kNbins+1; iBin++) fPtBinLimits[iBin]=xbins[iBin];
 103   fRandGener=new TRandom3(0);
 104   DefineInput(0, TChain::Class());
 105   DefineOutput(1, TList::Class());
 106   DefineOutput(2,TList::Class());
 107   Printf("end of AliAnalysisTaskSEITSsaSpectra");
 108 }
 109
 110 //___________________________________________________________________________
 111 AliAnalysisTaskSEITSsaSpectra::~AliAnalysisTaskSEITSsaSpectra(){
 112   // Destructor
 113   if (fOutput) {
 114     delete fOutput;
 115     fOutput = 0;
 116   }
 117   if (fListCuts) {
 118     delete fListCuts;
 119     fListCuts = 0;
 120   }
 121   if(fRandGener) delete fRandGener;
 122   if(fITSPIDResponse) delete fITSPIDResponse;
 123   delete fDCAxyCutFunc;
 124   delete fDCAzCutFunc;
 125 }
 126
 127 //________________________________________________________________________
 128 Double_t AliAnalysisTaskSEITSsaSpectra::CookdEdx(Double_t *s) const {
 129   // truncated mean for the dEdx
 130   Int_t nc=0;
 131   Double_t dedx[4]={0.,0.,0.,0.};
 132   for (Int_t il=0; il<4; il++) { // count good (>0) dE/dx values
 133     if(s[il]>fMindEdx){
 134       dedx[nc]= s[il];
 135       nc++;
 136     }
 137   }
 138   if(nc<fMinNdEdxSamples) return -1.;
 139
 140   Double_t tmp;
 141   Int_t swap; // sort in ascending order
 142   do {
 143     swap=0;
 144     for (Int_t i=0; i<nc-1; i++) {
 145       if (dedx[i]<=dedx[i+1]) continue;
 146       tmp=dedx[i];
 147       dedx[i]=dedx[i+1];
 148       dedx[i+1]=tmp;
 149       swap++;
 150     }
 151   } while (swap);
 152
 153   Double_t sumamp=0,sumweight=0;
 154   Double_t weight[4]={1.,1.,0.,0.};
 155   if(nc==3) weight[1]=0.5;
 156   else if(nc<3) weight[1]=0.;
 157   for (Int_t i=0; i<nc; i++) {
 158     sumamp+= dedx[i]*weight[i];
 159     sumweight+=weight[i];
 160   }
 161   return sumamp/sumweight;
 162 }
 163
 164
 165 //________________________________________________________________________
 166 Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt) const {
 167   // cut on transverse impact parameter updaated on 20-5-2010
 168   // from the study of L. Milano, F. Prino on the ITS standalone tracks
 169   // using the common binning of the TPC tracks
 170
 171   Double_t xyMax = fDCAxyCutFunc->Eval(pt); //in micron
 172   if((TMath::Abs(impactXY)*10000)>xyMax) return kFALSE;
 173   Double_t zMax = fDCAzCutFunc->Eval(pt); //in micron
 174   if((TMath::Abs(impactZ)*10000)>zMax) return kFALSE;
 175
 176   return kTRUE;
 177 }
 178
 179 //________________________________________________________________________
 180 Double_t AliAnalysisTaskSEITSsaSpectra::Eta2y(Double_t pt, Double_t m, Double_t eta) const {
 181   // convert eta to y
 182   Double_t mt = TMath::Sqrt(m*m + pt*pt);
 183   return TMath::ASinH(pt/mt*TMath::SinH(eta));
 184 }
 185
 186
 187 //________________________________________________________________________
 188 void AliAnalysisTaskSEITSsaSpectra::Init(){
 189   //
 190   // Initialization
 191   //
 192   fListCuts=new TList();
 193   fListCuts->SetOwner();
 194   Double_t xyP[3];
 195   Double_t zP[3];
 196   if(fMC){
 197     if(fYear==2009){
 198       xyP[0]=88.63; //MC LHC10a12
 199       xyP[1]=19.57;
 200       xyP[2]=1.65;
 201       zP[0]=140.98;
 202       zP[1]=62.33;
 203       zP[2]=1.15;
 204     }else{
 205       xyP[0]=36.; //MC LHC10d1
 206       xyP[1]=43.9;
 207       xyP[2]=1.3;
 208       zP[0]=111.9;
 209       zP[1]=59.8;
 210       zP[2]=1.2;
 211     }
 212   }else{
 213     if(fYear==2009){
 214       xyP[0]=85.28;//DATA 900 GeV pass6
 215       xyP[1]=25.78;
 216       xyP[2]=1.55;
 217       zP[0]=146.80;
 218       zP[1]=70.07;
 219       zP[2]=1.11;
 220     }else{
 221       xyP[0]=32.7;//DATA 7 TeV pass2
 222       xyP[1]=44.8;
 223       xyP[2]=1.3;
 224       zP[0]=117.3;
 225       zP[1]=66.8;
 226       zP[2]=1.2;
 227     }
 228   }
 229   fDCAxyCutFunc = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
 230   for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutFunc->SetParameter(ipar,xyP[ipar]);
 231   fDCAxyCutFunc->SetParameter(3,fNSigmaDCAxy);
 232
 233
 234
 235   fDCAzCutFunc = new TF1("fDCAzCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
 236   for(Int_t ipar=0; ipar<3; ipar++) fDCAzCutFunc->SetParameter(ipar,zP[ipar]);
 237   fDCAzCutFunc->SetParameter(3,fNSigmaDCAz);
 238
 239   fListCuts->Add(fDCAxyCutFunc);
 240   fListCuts->Add(fDCAzCutFunc);
 241
 242   PostData(2,fListCuts);
 243 }
 244
 245 //________________________________________________________________________
 246 void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){
 247   // Create a TList with histograms and a TNtuple
 248   // Called once
 249
 250   fOutput = new TList();
 251   fOutput->SetOwner();
 252   fOutput->SetName("Spiderman");
 253
 254   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events",9,-1.5,7.5);
 255   fHistNEvents->Sumw2();
 256   fHistNEvents->SetMinimum(0);
 257   fHistNEvents->GetXaxis()->SetBinLabel(1,"Read from ESD");
 258   fHistNEvents->GetXaxis()->SetBinLabel(2,"Pass Phys. Sel. + Trig");
 259   fHistNEvents->GetXaxis()->SetBinLabel(3,"SDD read out");
 260   fHistNEvents->GetXaxis()->SetBinLabel(4,"In mult. range");
 261   fHistNEvents->GetXaxis()->SetBinLabel(5,"With SPD vertex");
 262   fHistNEvents->GetXaxis()->SetBinLabel(6,"Vertex contributors >0");
 263   fHistNEvents->GetXaxis()->SetBinLabel(7,"|zVertex|<10");
 264   fHistNEvents->GetXaxis()->SetBinLabel(8,"Error on zVertex<0.5");
 265   fHistNEvents->GetXaxis()->SetBinLabel(9,"Good Z vertex");
 266   fOutput->Add(fHistNEvents);
 267
 268   fHistMult = new TH1F("fHistMult", "Event Multiplicity",3000,-0.5,2999.5);
 269   fHistMult->Sumw2();
 270   fHistMult->SetMinimum(0);
 271   if(fMultEstimator==0) fHistMult->GetXaxis()->SetTitle("Multiplicity |#eta|<0.8");
 272   else if(fMultEstimator==1) fHistMult->GetXaxis()->SetTitle("Tracklets |#eta|<0.8");
 273   else if(fMultEstimator==2) fHistMult->GetXaxis()->SetTitle("Clusters on SPD1");
 274   fOutput->Add(fHistMult);
 275
 276   fHistCen = new TH1F("fHistCen", "Event Centrality",101,-0.5,100.5);
 277   fHistCen->Sumw2();
 278   fHistCen->SetMinimum(0);
 279   fOutput->Add(fHistCen);
 280
 281   fHistNTracks = new TH1F("fHistNTracks", "Number of ITSsa tracks",20,0.5,20.5);
 282   fHistNTracks->Sumw2();
 283   fHistNTracks->SetMinimum(0);
 284   fOutput->Add(fHistNTracks);
 285
 286   fHistNTracksPos = new TH1F("fHistNTracksPos", "Number of positive ITSsa tracks",20,0.5,20.5);
 287   fHistNTracksPos->Sumw2();
 288   fHistNTracksPos->SetMinimum(0);
 289   fOutput->Add(fHistNTracksPos);
 290
 291   fHistNTracksNeg = new TH1F("fHistNTracksNeg", "Number of negative ITSsa tracks",20,0.5,20.5);
 292   fHistNTracksNeg->Sumw2();
 293   fHistNTracksNeg->SetMinimum(0);
 294   fOutput->Add(fHistNTracksNeg);
 295
 296   //binning for the histogram
 297   const Int_t hnbins=400;
 298   Double_t hxmin = 0.01;
 299   Double_t hxmax = 10;
 300   Double_t hlogxmin = TMath::Log10(hxmin);
 301   Double_t hlogxmax = TMath::Log10(hxmax);
 302   Double_t hbinwidth = (hlogxmax-hlogxmin)/hnbins;
 303   Double_t hxbins[hnbins+1];
 304   hxbins[0] = 0.01;
 305   for (Int_t i=1;i<=hnbins;i++) {
 306     hxbins[i] = hxmin + TMath::Power(10,hlogxmin+i*hbinwidth);
 307   }
 308
 309   fHistDEDX = new TH2F("fHistDEDX","",hnbins,hxbins,900,0,1000);
 310   fOutput->Add(fHistDEDX);
 311
 312   fHistDEDXdouble = new TH2F("fHistDEDXdouble","",500,-5,5,900,0,1000);
 313   fOutput->Add(fHistDEDXdouble);
 314
 315   fHistBeforeEvSel = new TH1F("fHistBeforeEvSel","fHistBeforeEvSel",kNbins,fPtBinLimits);
 316   fHistAfterEvSel = new TH1F("fHistAfterEvSel","fHistAfterEvSel",kNbins,fPtBinLimits);
 317   fOutput->Add(fHistBeforeEvSel);
 318   fOutput->Add(fHistAfterEvSel);
 319
 320
 321
 322   for(Int_t j=0;j<3;j++){
 323
 324     fHistPosNSigmaSep[j] = new TH2F(Form("fHistPosNSigmaSep%d",j),"",hnbins,hxbins,1000,-10,10);
 325     fOutput->Add(fHistPosNSigmaSep[j]);
 326     fHistNegNSigmaSep[j] = new TH2F(Form("fHistNegNSigmaSep%d",j),"",hnbins,hxbins,1000,-10,10);
 327     fOutput->Add(fHistNegNSigmaSep[j]);
 328
 329     fHistPrimMCpos[j] = new TH1F(Form("fHistPrimMCpos%d",j),Form("fHistPrimMCpos%d",j),kNbins,fPtBinLimits);
 330     fHistPrimMCneg[j] = new TH1F(Form("fHistPrimMCneg%d",j),Form("fHistPrimMCneg%d",j),kNbins,fPtBinLimits);
 331     fOutput->Add(fHistPrimMCneg[j]);
 332     fOutput->Add(fHistPrimMCpos[j]);
 333     fHistSecStrMCpos[j] = new TH1F(Form("fHistSecStrMCpos%d",j),Form("fHistSecStrMCpos%d",j),kNbins,fPtBinLimits);
 334     fHistSecStrMCneg[j] = new TH1F(Form("fHistSecStrMCneg%d",j),Form("fHistSecStrMCneg%d",j),kNbins,fPtBinLimits);
 335     fOutput->Add(fHistSecStrMCneg[j]);
 336     fOutput->Add(fHistSecStrMCpos[j]);
 337     fHistSecMatMCpos[j] = new TH1F(Form("fHistSecMatMCpos%d",j),Form("fHistSecMatMCpos%d",j),kNbins,fPtBinLimits);
 338     fHistSecMatMCneg[j] = new TH1F(Form("fHistSecMatMCneg%d",j),Form("fHistSecMatMCneg%d",j),kNbins,fPtBinLimits);
 339     fOutput->Add(fHistSecMatMCneg[j]);
 340     fOutput->Add(fHistSecMatMCpos[j]);
 341     //
 342     fHistPrimMCposBefEvSel[j] = new TH1F(Form("fHistPrimMCposBefEvSel%d",j),Form("fHistPrimMCposBefEvSel%d",j),kNbins,fPtBinLimits);
 343     fHistPrimMCposBefEvSelEta[j] = new TH1F(Form("fHistPrimMCposBefEvSelEta%d",j),Form("fHistPrimMCposBefEvSelEta%d",j),kNbins,fPtBinLimits);
 344     fHistPrimMCposBefEvSelEtaY[j] = new TH1F(Form("fHistPrimMCposBefEvSelEtaY%d",j),Form("fHistPrimMCposBefEvSelEtaY%d",j),kNbins,fPtBinLimits);
 345     fHistPrimMCnegBefEvSel[j] = new TH1F(Form("fHistPrimMCnegBefEvSel%d",j),Form("fHistPrimMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
 346     fHistPrimMCnegBefEvSelEta[j] = new TH1F(Form("fHistPrimMCnegBefEvSelEta%d",j),Form("fHistPrimMCnegBefEvSelEta%d",j),kNbins,fPtBinLimits);
 347     fHistPrimMCnegBefEvSelEtaY[j] = new TH1F(Form("fHistPrimMCnegBefEvSelEtaY%d",j),Form("fHistPrimMCnegBefEvSelEtaY%d",j),kNbins,fPtBinLimits);
 348     fOutput->Add(fHistPrimMCposBefEvSel[j]);
 349     fOutput->Add(fHistPrimMCposBefEvSelEta[j]);
 350     fOutput->Add(fHistPrimMCposBefEvSelEtaY[j]);
 351     fOutput->Add(fHistPrimMCnegBefEvSel[j]);
 352     fOutput->Add(fHistPrimMCnegBefEvSelEta[j]);
 353     fOutput->Add(fHistPrimMCnegBefEvSelEtaY[j]);
 354     fHistSecStrMCposBefEvSel[j] = new TH1F(Form("fHistSecStrMCposBefEvSel%d",j),Form("fHistSecStrMCposBefEvSel%d",j),kNbins,fPtBinLimits);
 355     fHistSecStrMCnegBefEvSel[j] = new TH1F(Form("fHistSecStrMCnegBefEvSel%d",j),Form("fHistSecStrMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
 356     fOutput->Add(fHistSecStrMCnegBefEvSel[j]);
 357     fOutput->Add(fHistSecStrMCposBefEvSel[j]);
 358     fHistSecMatMCposBefEvSel[j] = new TH1F(Form("fHistSecMatMCposBefEvSel%d",j),Form("fHistSecMatMCposBefEvSel%d",j),kNbins,fPtBinLimits);
 359     fHistSecMatMCnegBefEvSel[j] = new TH1F(Form("fHistSecMatMCnegBefEvSel%d",j),Form("fHistSecMatMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
 360     fOutput->Add(fHistSecMatMCnegBefEvSel[j]);
 361     fOutput->Add(fHistSecMatMCposBefEvSel[j]);
 362     //
 363     fHistPrimMCposReco[j] = new TH1F(Form("fHistPrimMCposReco%d",j),Form("fHistPrimMCposReco%d",j),kNbins,fPtBinLimits);
 364     fHistPrimMCposRecoEtaY[j] = new TH1F(Form("fHistPrimMCposRecoEtaY%d",j),Form("fHistPrimMCposRecoEtaY%d",j),kNbins,fPtBinLimits);
 365     fHistPrimMCnegReco[j] = new TH1F(Form("fHistPrimMCnegReco%d",j),Form("fHistPrimMCnegReco%d",j),kNbins,fPtBinLimits);
 366     fHistPrimMCnegRecoEtaY[j] = new TH1F(Form("fHistPrimMCnegRecoEtaY%d",j),Form("fHistPrimMCnegRecoEtaY%d",j),kNbins,fPtBinLimits);
 367     fOutput->Add(fHistPrimMCposReco[j]);
 368     fOutput->Add(fHistPrimMCposRecoEtaY[j]);
 369     fOutput->Add(fHistPrimMCnegReco[j]);
 370     fOutput->Add(fHistPrimMCnegRecoEtaY[j]);
 371     fHistSecStrMCposReco[j] = new TH1F(Form("fHistSecStrMCposReco%d",j),Form("fHistSecStrMCposReco%d",j),kNbins,fPtBinLimits);
 372     fHistSecStrMCnegReco[j] = new TH1F(Form("fHistSecStrMCnegReco%d",j),Form("fHistSecStrMCnegReco%d",j),kNbins,fPtBinLimits);
 373     fOutput->Add(fHistSecStrMCnegReco[j]);
 374     fOutput->Add(fHistSecStrMCposReco[j]);
 375     fHistSecMatMCposReco[j] = new TH1F(Form("fHistSecMatMCposReco%d",j),Form("fHistSecMatMCposReco%d",j),kNbins,fPtBinLimits);
 376     fHistSecMatMCnegReco[j] = new TH1F(Form("fHistSecMatMCnegReco%d",j),Form("fHistSecMatMCnegReco%d",j),kNbins,fPtBinLimits);
 377     fOutput->Add(fHistSecMatMCnegReco[j]);
 378     fOutput->Add(fHistSecMatMCposReco[j]);
 379
 380   }
 381
 382   for(Int_t i=0; i<4; i++){
 383     fHistCharge[i] = new TH1F(Form("fHistChargeLay%d",i),Form("fHistChargeLay%d",i),100,0,300);
 384     fOutput->Add(fHistCharge[i]);
 385   }
 386
 387   for(Int_t i=0; i<kNbins; i++){
 388     fHistPosPi[i] = new TH1F(Form("fHistPosPi%d",i),Form("fHistPosPi%d",i),175,-3.5,3.5);
 389     fHistPosK[i]  = new TH1F(Form("fHistPosK%d",i),Form("fHistPosK%d",i),175,-3.5,3.5);
 390     fHistPosP[i]  = new TH1F(Form("fHistPosP%d",i),Form("fHistPosP%d",i),175,-3.5,3.5);
 391     fHistNegPi[i] = new TH1F(Form("fHistNegPi%d",i),Form("fHistNegPi%d",i),175,-3.5,3.5);
 392     fHistNegK[i]  = new TH1F(Form("fHistNegK%d",i),Form("fHistNegK%d",i),175,-3.5,3.5);
 393     fHistNegP[i]  = new TH1F(Form("fHistNegP%d",i),Form("fHistNegP%d",i),175,-3.5,3.5);
 394
 395     fHistDCAPosPi[i] = new TH1F(Form("fHistDCAPosPi%d",i),Form("fHistDCAPosPi%d",i),2000,-1,1);  //DCA distr. with NSigma PID
 396     fHistDCAPosK[i]  = new TH1F(Form("fHistDCAPosK%d",i),Form("fHistDCAPosK%d",i),2000,-1,1);
 397     fHistDCAPosP[i]  = new TH1F(Form("fHistDCAPosP%d",i),Form("fHistDCAPosP%d",i),2000,-1,1);
 398     fHistDCANegPi[i] = new TH1F(Form("fHistDCANegPi%d",i),Form("fHistDCANegPi%d",i),2000,-1,1);
 399     fHistDCANegK[i]  = new TH1F(Form("fHistDCANegK%d",i),Form("fHistDCANegK%d",i),2000,-1,1);
 400     fHistDCANegP[i]  = new TH1F(Form("fHistDCANegP%d",i),Form("fHistDCANegP%d",i),2000,-1,1);
 401
 402     fHistMCPrimDCAPosPi[i] = new TH1F(Form("fHistMCPrimDCAPosPi%d",i),Form("fHistMCPrimDCAPosPi%d",i),2000,-1,1);  //DCA distr. with MC truth
 403     fHistMCPrimDCAPosK[i]  = new TH1F(Form("fHistMCPrimDCAPosK%d",i),Form("fHistMCPrimDCAPosK%d",i),2000,-1,1);
 404     fHistMCPrimDCAPosP[i]  = new TH1F(Form("fHistMCPrimDCAPosP%d",i),Form("fHistMCPrimDCAPosP%d",i),2000,-1,1);
 405     fHistMCPrimDCANegPi[i] = new TH1F(Form("fHistMCPrimDCANegPi%d",i),Form("fHistMCPrimDCANegPi%d",i),2000,-1,1);
 406     fHistMCPrimDCANegK[i]  = new TH1F(Form("fHistMCPrimDCANegK%d",i),Form("fHistMCPrimDCANegK%d",i),2000,-1,1);
 407     fHistMCPrimDCANegP[i]  = new TH1F(Form("fHistMCPrimDCANegP%d",i),Form("fHistMCPrimDCANegP%d",i),2000,-1,1);
 408
 409     fHistMCSecStDCAPosPi[i] = new TH1F(Form("fHistMCSecStDCAPosPi%d",i),Form("fHistMCSecStDCAPosPi%d",i),2000,-1,1);  //DCA distr. with MC truth
 410     fHistMCSecStDCAPosK[i]  = new TH1F(Form("fHistMCSecStDCAPosK%d",i),Form("fHistMCSecStDCAPosK%d",i),2000,-1,1);
 411     fHistMCSecStDCAPosP[i]  = new TH1F(Form("fHistMCSecStDCAPosP%d",i),Form("fHistMCSecStDCAPosP%d",i),2000,-1,1);
 412     fHistMCSecStDCANegPi[i] = new TH1F(Form("fHistMCSecStDCANegPi%d",i),Form("fHistMCSecStDCANegPi%d",i),2000,-1,1);
 413     fHistMCSecStDCANegK[i]  = new TH1F(Form("fHistMCSecStDCANegK%d",i),Form("fHistMCSecStDCANegK%d",i),2000,-1,1);
 414     fHistMCSecStDCANegP[i]  = new TH1F(Form("fHistMCSecStDCANegP%d",i),Form("fHistMCSecStDCANegP%d",i),2000,-1,1);
 415
 416     fHistMCSecMatDCAPosPi[i] = new TH1F(Form("fHistMCSecMatDCAPosPi%d",i),Form("fHistMCSecMatDCAPosPi%d",i),2000,-1,1);  //DCA distr. with MC truth
 417     fHistMCSecMatDCAPosK[i]  = new TH1F(Form("fHistMCSecMatDCAPosK%d",i),Form("fHistMCSecMatDCAPosK%d",i),2000,-1,1);
 418     fHistMCSecMatDCAPosP[i]  = new TH1F(Form("fHistMCSecMatDCAPosP%d",i),Form("fHistMCSecMatDCAPosP%d",i),2000,-1,1);
 419     fHistMCSecMatDCANegPi[i] = new TH1F(Form("fHistMCSecMatDCANegPi%d",i),Form("fHistMCSecMatDCANegPi%d",i),2000,-1,1);
 420     fHistMCSecMatDCANegK[i]  = new TH1F(Form("fHistMCSecMatDCANegK%d",i),Form("fHistMCSecMatDCANegK%d",i),2000,-1,1);
 421     fHistMCSecMatDCANegP[i]  = new TH1F(Form("fHistMCSecMatDCANegP%d",i),Form("fHistMCSecMatDCANegP%d",i),2000,-1,1);
 422
 423     fHistMCPosOtherHypPion[i] = new TH1F(Form("fHistMCPosOtherHypPion%d",i),Form("fHistMCPosOtherHypPion%d",i),175,-3.5,3.5);   //MC truth
 424     fHistMCPosOtherHypKaon[i] = new TH1F(Form("fHistMCPosOtherHypKaon%d",i),Form("fHistMCPosOtherHypKaon%d",i),175,-3.5,3.5);
 425     fHistMCPosOtherHypProton[i] = new TH1F(Form("fHistMCPosOtherHypProton%d",i),Form("fHistMCPosOtherHypProton%d",i),175,-3.5,3.5);
 426     fHistMCPosElHypPion[i] = new TH1F(Form("fHistMCPosElHypPion%d",i),Form("fHistMCPosElHypPion%d",i),175,-3.5,3.5);
 427     fHistMCPosElHypKaon[i] = new TH1F(Form("fHistMCPosElHypKaon%d",i),Form("fHistMCPosElHypKaon%d",i),175,-3.5,3.5);
 428     fHistMCPosElHypProton[i] = new TH1F(Form("fHistMCPosElHypProton%d",i),Form("fHistMCPosElHypProton%d",i),175,-3.5,3.5);
 429     fHistMCPosPiHypPion[i] = new TH1F(Form("fHistMCPosPiHypPion%d",i),Form("fHistMCPosPiHypPion%d",i),175,-3.5,3.5);
 430     fHistMCPosPiHypKaon[i] = new TH1F(Form("fHistMCPosPiHypKaon%d",i),Form("fHistMCPosPiHypKaon%d",i),175,-3.5,3.5);
 431     fHistMCPosPiHypProton[i] = new TH1F(Form("fHistMCPosPiHypProton%d",i),Form("fHistMCPosPiHypProton%d",i),175,-3.5,3.5);
 432     fHistMCPosKHypPion[i] = new TH1F(Form("fHistMCPosKHypPion%d",i),Form("fHistMCPosKHypPion%d",i),175,-3.5,3.5);
 433     fHistMCPosKHypKaon[i]  = new TH1F(Form("fHistMCPosKHypKaon%d",i),Form("fHistMCPosKHypKaon%d",i),175,-3.5,3.5);
 434     fHistMCPosKHypProton[i] = new TH1F(Form("fHistMCPosKHypProton%d",i),Form("fHistMCPosKHypProton%d",i),175,-3.5,3.5);
 435     fHistMCPosPHypPion[i] = new TH1F(Form("fHistMCPosPHypPion%d",i),Form("fHistMCPosPHypPion%d",i),175,-3.5,3.5);
 436     fHistMCPosPHypKaon[i] = new TH1F(Form("fHistMCPosPHypKaon%d",i),Form("fHistMCPosPHypKaon%d",i),175,-3.5,3.5);
 437     fHistMCPosPHypProton[i]  = new TH1F(Form("fHistMCPosPHypProton%d",i),Form("fHistMCPosPHypProton%d",i),175,-3.5,3.5);
 438
 439     fHistMCNegOtherHypPion[i] = new TH1F(Form("fHistMCNegOtherHypPion%d",i),Form("fHistMCNegOtherHypPion%d",i),175,-3.5,3.5);   //MC truth
 440     fHistMCNegOtherHypKaon[i] = new TH1F(Form("fHistMCNegOtherHypKaon%d",i),Form("fHistMCNegOtherHypKaon%d",i),175,-3.5,3.5);
 441     fHistMCNegOtherHypProton[i] = new TH1F(Form("fHistMCNegOtherHypProton%d",i),Form("fHistMCNegOtherHypProton%d",i),175,-3.5,3.5);
 442     fHistMCNegElHypPion[i] = new TH1F(Form("fHistMCNegElHypPion%d",i),Form("fHistMCNegElHypPion%d",i),175,-3.5,3.5);
 443     fHistMCNegElHypKaon[i] = new TH1F(Form("fHistMCNegElHypKaon%d",i),Form("fHistMCNegElHypKaon%d",i),175,-3.5,3.5);
 444     fHistMCNegElHypProton[i] = new TH1F(Form("fHistMCNegElHypProton%d",i),Form("fHistMCNegElHypProton%d",i),175,-3.5,3.5);
 445     fHistMCNegPiHypPion[i] = new TH1F(Form("fHistMCNegPiHypPion%d",i),Form("fHistMCNegPiHypPion%d",i),175,-3.5,3.5);
 446     fHistMCNegPiHypKaon[i] = new TH1F(Form("fHistMCNegPiHypKaon%d",i),Form("fHistMCNegPiHypKaon%d",i),175,-3.5,3.5);
 447     fHistMCNegPiHypProton[i] = new TH1F(Form("fHistMCNegPiHypProton%d",i),Form("fHistMCNegPiHypProton%d",i),175,-3.5,3.5);
 448     fHistMCNegKHypPion[i] = new TH1F(Form("fHistMCNegKHypPion%d",i),Form("fHistMCNegKHypPion%d",i),175,-3.5,3.5);
 449     fHistMCNegKHypKaon[i]  = new TH1F(Form("fHistMCNegKHypKaon%d",i),Form("fHistMCNegKHypKaon%d",i),175,-3.5,3.5);
 450     fHistMCNegKHypProton[i] = new TH1F(Form("fHistMCNegKHypProton%d",i),Form("fHistMCNegKHypProton%d",i),175,-3.5,3.5);
 451     fHistMCNegPHypPion[i] = new TH1F(Form("fHistMCNegPHypPion%d",i),Form("fHistMCNegPHypPion%d",i),175,-3.5,3.5);
 452     fHistMCNegPHypKaon[i] = new TH1F(Form("fHistMCNegPHypKaon%d",i),Form("fHistMCNegPHypKaon%d",i),175,-3.5,3.5);
 453     fHistMCNegPHypProton[i]  = new TH1F(Form("fHistMCNegPHypProton%d",i),Form("fHistMCNegPHypProton%d",i),175,-3.5,3.5);
 454
 455
 456     fOutput->Add(fHistPosPi[i]);
 457     fOutput->Add(fHistPosK[i]);
 458     fOutput->Add(fHistPosP[i]);
 459     fOutput->Add(fHistNegPi[i]);
 460     fOutput->Add(fHistNegK[i]);
 461     fOutput->Add(fHistNegP[i]);
 462
 463     fOutput->Add(fHistDCAPosPi[i]); //DCA distr
 464     fOutput->Add(fHistDCAPosK[i]);
 465     fOutput->Add(fHistDCAPosP[i]);
 466     fOutput->Add(fHistDCANegPi[i]);
 467     fOutput->Add(fHistDCANegK[i]);
 468     fOutput->Add(fHistDCANegP[i]);
 469
 470     fOutput->Add(fHistMCPrimDCAPosPi[i]);//DCA distr.
 471     fOutput->Add(fHistMCPrimDCAPosK[i]);
 472     fOutput->Add(fHistMCPrimDCAPosP[i]);
 473     fOutput->Add(fHistMCPrimDCANegPi[i]);
 474     fOutput->Add(fHistMCPrimDCANegK[i]);
 475     fOutput->Add(fHistMCPrimDCANegP[i]);
 476
 477     fOutput->Add(fHistMCSecStDCAPosPi[i]);//DCA distr.
 478     fOutput->Add(fHistMCSecStDCAPosK[i]);
 479     fOutput->Add(fHistMCSecStDCAPosP[i]);
 480     fOutput->Add(fHistMCSecStDCANegPi[i]);
 481     fOutput->Add(fHistMCSecStDCANegK[i]);
 482     fOutput->Add(fHistMCSecStDCANegP[i]);
 483
 484     fOutput->Add(fHistMCSecMatDCAPosPi[i]);//DCA distr.
 485     fOutput->Add(fHistMCSecMatDCAPosK[i]);
 486     fOutput->Add(fHistMCSecMatDCAPosP[i]);
 487     fOutput->Add(fHistMCSecMatDCANegPi[i]);
 488     fOutput->Add(fHistMCSecMatDCANegK[i]);
 489     fOutput->Add(fHistMCSecMatDCANegP[i]);
 490
 491     fOutput->Add(fHistMCPosOtherHypPion[i]);//MC truth
 492     fOutput->Add(fHistMCPosOtherHypKaon[i]);
 493     fOutput->Add(fHistMCPosOtherHypProton[i]);
 494     fOutput->Add(fHistMCPosElHypPion[i]);
 495     fOutput->Add(fHistMCPosElHypKaon[i]);
 496     fOutput->Add(fHistMCPosElHypProton[i]);
 497     fOutput->Add(fHistMCPosPiHypPion[i]);
 498     fOutput->Add(fHistMCPosPiHypKaon[i]);
 499     fOutput->Add(fHistMCPosPiHypProton[i]);
 500     fOutput->Add(fHistMCPosKHypPion[i]);
 501     fOutput->Add(fHistMCPosKHypKaon[i]);
 502     fOutput->Add(fHistMCPosKHypProton[i]);
 503     fOutput->Add(fHistMCPosPHypPion[i]);
 504     fOutput->Add(fHistMCPosPHypKaon[i]);
 505     fOutput->Add(fHistMCPosPHypProton[i]);
 506
 507     fOutput->Add(fHistMCNegOtherHypPion[i]);//MC truth
 508     fOutput->Add(fHistMCNegOtherHypKaon[i]);
 509     fOutput->Add(fHistMCNegOtherHypProton[i]);
 510     fOutput->Add(fHistMCNegElHypPion[i]);
 511     fOutput->Add(fHistMCNegElHypKaon[i]);
 512     fOutput->Add(fHistMCNegElHypProton[i]);
 513     fOutput->Add(fHistMCNegPiHypPion[i]);
 514     fOutput->Add(fHistMCNegPiHypKaon[i]);
 515     fOutput->Add(fHistMCNegPiHypProton[i]);
 516     fOutput->Add(fHistMCNegKHypPion[i]);
 517     fOutput->Add(fHistMCNegKHypKaon[i]);
 518     fOutput->Add(fHistMCNegKHypProton[i]);
 519     fOutput->Add(fHistMCNegPHypPion[i]);
 520     fOutput->Add(fHistMCNegPHypKaon[i]);
 521     fOutput->Add(fHistMCNegPHypProton[i]);
 522
 523   }
 524
 525   //NSigma Histos
 526   for(Int_t j=0;j<3;j++){
 527
 528     fHistPosNSigmaMean[j] = new TH1F(Form("hHistPosNSigmaMean%d",j),Form("hHistPosNSigmaMean%d",j),kNbins,fPtBinLimits);
 529     fHistNegNSigmaMean[j] = new TH1F(Form("hHistNegNSigmaMean%d",j),Form("hHistNegNSigmaMean%d",j),kNbins,fPtBinLimits);
 530     fHistPosNSigmaMCMean[j] = new TH1F(Form("hHistPosNSigmaMCMean%d",j),Form("hHistPosNSigmaMCMean%d",j),kNbins,fPtBinLimits);
 531     fHistNegNSigmaMCMean[j] = new TH1F(Form("hHistNegNSigmaMCMean%d",j),Form("hHistNegNSigmaMCMean%d",j),kNbins,fPtBinLimits);
 532     fHistPosNSigmaPrimMean[j] = new TH1F(Form("hHistPosNSigmaPrimMean%d",j),Form("hHistPosNSigmaPrimMean%d",j),kNbins,fPtBinLimits);
 533     fHistNegNSigmaPrimMean[j] = new TH1F(Form("hHistNegNSigmaPrimMean%d",j),Form("hHistNegNSigmaPrimMean%d",j),kNbins,fPtBinLimits);
 534     fHistPosNSigmaPrimMCMean[j] = new TH1F(Form("hHistPosNSigmaPrimMCMean%d",j),Form("hHistPosNSigmaPrimMCMean%d",j),kNbins,fPtBinLimits);
 535     fHistNegNSigmaPrimMCMean[j] = new TH1F(Form("hHistNegNSigmaPrimMCMean%d",j),Form("hHistNegNSigmaPrimMCMean%d",j),kNbins,fPtBinLimits);
 536     fOutput->Add(fHistPosNSigmaMean[j]);
 537     fOutput->Add(fHistNegNSigmaMean[j]);
 538     fOutput->Add(fHistPosNSigmaMCMean[j]);
 539     fOutput->Add(fHistNegNSigmaMCMean[j]);
 540     fOutput->Add(fHistPosNSigmaPrimMean[j]);
 541     fOutput->Add(fHistNegNSigmaPrimMean[j]);
 542     fOutput->Add(fHistPosNSigmaPrimMCMean[j]);
 543     fOutput->Add(fHistNegNSigmaPrimMCMean[j]);
 544
 545     fHistPosNSigma[j] = new TH1F(Form("hHistPosNSigma%d",j),Form("hHistPosNSigma%d",j),kNbins,fPtBinLimits);
 546     fHistNegNSigma[j] = new TH1F(Form("hHistNegNSigma%d",j),Form("hHistNegNSigma%d",j),kNbins,fPtBinLimits);
 547     fHistPosNSigmaMC[j] = new TH1F(Form("hHistPosNSigmaMC%d",j),Form("hHistPosNSigmaMC%d",j),kNbins,fPtBinLimits);
 548     fHistNegNSigmaMC[j] = new TH1F(Form("hHistNegNSigmaMC%d",j),Form("hHistNegNSigmaMC%d",j),kNbins,fPtBinLimits);
 549     fHistPosNSigmaPrim[j] = new TH1F(Form("hHistPosNSigmaPrim%d",j),Form("hHistPosNSigmaPrim%d",j),kNbins,fPtBinLimits);
 550     fHistNegNSigmaPrim[j] = new TH1F(Form("hHistNegNSigmaPrim%d",j),Form("hHistNegNSigmaPrim%d",j),kNbins,fPtBinLimits);
 551     fHistPosNSigmaPrimMC[j] = new TH1F(Form("hHistPosNSigmaPrimMC%d",j),Form("hHistPosNSigmaPrimMC%d",j),kNbins,fPtBinLimits);
 552     fHistNegNSigmaPrimMC[j] = new TH1F(Form("hHistNegNSigmaPrimMC%d",j),Form("hHistNegNSigmaPrimMC%d",j),kNbins,fPtBinLimits);
 553     fOutput->Add(fHistPosNSigma[j]);
 554     fOutput->Add(fHistNegNSigma[j]);
 555     fOutput->Add(fHistPosNSigmaMC[j]);
 556     fOutput->Add(fHistNegNSigmaMC[j]);
 557     fOutput->Add(fHistPosNSigmaPrim[j]);
 558     fOutput->Add(fHistNegNSigmaPrim[j]);
 559     fOutput->Add(fHistPosNSigmaPrimMC[j]);
 560     fOutput->Add(fHistNegNSigmaPrimMC[j]);
 561
 562   }
 563
 564   fNtupleNSigma = new TNtuple("fNtupleNSigma","fNtupleNSigma","p:pt:dedx:dedx3:dedx4:dedx5:dedx6:ncls:nclspid:sign:run:eta:impactXY:impactZ:isph:pdgcode:mfl:uniqueID:chi2ncls");
 565   fOutput->Add(fNtupleNSigma);
 566   fNtupleMC = new TNtuple("fNtupleMC","fNtupleMC","ptMC:pdgcode:signMC:etaMC:yMC:isph:evSel:run");
 567   fOutput->Add(fNtupleMC);
 568
 569
 570
 571   // Post output data.
 572   PostData(1,fOutput);
 573
 574   Printf("end of CreateOutputObjects");
 575 }
 576
 577 //________________________________________________________________________
 578 void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
 579   // Main loop
 580   // Called for each event
 581
 582   ///////////////////////////////////////
 583   //variables
 584   Float_t pdgmass[4]={0.13957,0.493677,0.938272,1.8756}; //mass for pi, K, P (Gev/c^2)
 585   Int_t listcode[3]={211,321,2212};//code for pi, K, P (Gev/c^2)
 586   Double_t dedxLay[4];
 587   Float_t ptMC=-999;
 588   Int_t code=-999, signMC=-999,isph=-999,mfl=-999,uniqueID=-999;
 589   Float_t impactXY=-999, impactZ=-999;
 590   Int_t evSel=1;
 591   AliESDtrack* track;
 592   UInt_t status;
 593   AliStack* stack=0;
 594   TParticle *part=0;
 595   TParticlePDG *pdgPart=0;
 596
 597   /////////////////////
 598
 599   fESD=(AliESDEvent*)InputEvent();
 600   if(!fESD) {
 601     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
 602     return;
 603   }
 604   fHistNEvents->Fill(-1);
 605
 606   UInt_t maskPhysSel = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 607   TString firedTriggerClasses=fESD->GetFiredTriggerClasses();
 608   //  if(!firedTriggerClasses.Contains("CINT1B")) return;
 609   if((maskPhysSel & AliVEvent::kMB)==0) return;
 610   fHistNEvents->Fill(0);
 611
 612   if(fLowEnergypp && !fMC){ // remove events without SDD in pp 2.76 TeV
 613     if(!firedTriggerClasses.Contains("ALL")) return;
 614   }
 615   fHistNEvents->Fill(1);
 616
 617
 618   if(fMC){
 619     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 620     if (!eventHandler) {
 621       Printf("ERROR: Could not retrieve MC event handler");
 622       return;
 623     }
 624     AliMCEvent* mcEvent = eventHandler->MCEvent();
 625     if (!mcEvent) {
 626       Printf("ERROR: Could not retrieve MC event");
 627       return;
 628     }
 629     stack = mcEvent->Stack();
 630     if (!stack) {
 631       printf("ERROR: stack not available\n");
 632       return;
 633     }
 634   }
 635   if(!fITSPIDResponse){
 636     fITSPIDResponse=new AliITSPIDResponse(fMC);
 637   }
 638
 639   //flags for MC
 640   Int_t nTrackMC=0;
 641   if(stack) nTrackMC = stack->GetNtrack();
 642   const AliESDVertex *vtx =  fESD->GetPrimaryVertexSPD();
 643
 644   ///////////selection of the centrality or multiplicity bin
 645
 646   //selection on the event centrality
 647   if(fHImode){
 648     if(!(fLowCentrality<0.0)&&fUpCentrality>0.0)
 649       {
 650         AliCentrality *centrality = fESD->GetCentrality();
 651         if(!centrality->IsEventInCentralityClass(fLowCentrality,fUpCentrality,"V0M"))
 652           return;
 653         Printf("Centrality of the event: %.1f",centrality->GetCentralityPercentile("V0M"));
 654         Printf("Centrality cut: %.1f to %.1f",fLowCentrality,fUpCentrality);
 655         fHistCen->Fill(centrality->GetCentralityPercentile("V0M"));
 656       }
 657   }
 658
 659   //selection on the event multiplicity based on global tracks
 660   Int_t multiplicity = -1;
 661   if(fMultEstimator==0){
 662     // tracks+tracklets
 663     multiplicity=AliESDtrackCuts::GetReferenceMultiplicity(fESD, AliESDtrackCuts::kTrackletsITSTPC, 0.8);
 664   }else if(fMultEstimator==1){
 665     // tracklets
 666     multiplicity=AliESDtrackCuts::GetReferenceMultiplicity(fESD, AliESDtrackCuts::kTracklets, 0.8);
 667   }else if(fMultEstimator==2){
 668     // clusters in SPD1
 669     const AliMultiplicity *mult = fESD->GetMultiplicity();
 670     Float_t nClu1 = (Float_t)mult->GetNumberOfITSClusters(1);
 671     multiplicity =(Int_t)(AliESDUtils::GetCorrSPD2(nClu1,vtx->GetZ())+0.5);
 672   }
 673   if(fLowMult>-1 && multiplicity<fLowMult) return;
 674   if(fUpMult>-1 && multiplicity>fUpMult) return;
 675   fHistMult->Fill(multiplicity);
 676   fHistNEvents->Fill(2);
 677
 678   //event selection
 679   if(!vtx)evSel=0;
 680   else{
 681     fHistNEvents->Fill(3);
 682     if(vtx->GetNContributors()<0) evSel=0;
 683     else{
 684       fHistNEvents->Fill(4);
 685       if(TMath::Abs(vtx->GetZv())>10) evSel=0;
 686       else{
 687         fHistNEvents->Fill(5);
 688         if(vtx->GetZRes()>0.5) evSel=0;
 689         else{
 690           fHistNEvents->Fill(6);
 691           if(vtx->IsFromVertexerZ() && vtx->GetDispersion()>0.03) evSel=0;
 692           else fHistNEvents->Fill(7);
 693         }
 694       }
 695     }
 696   }
 697
 698
 699   /////first loop on stack, before event selection, filling MC ntuple
 700
 701   for(Int_t imc=0; imc<nTrackMC; imc++){
 702     part = stack->Particle(imc);
 703     isph=1;
 704     if(!stack->IsPhysicalPrimary(imc)) isph=0;
 705     pdgPart = part->GetPDG();
 706     if(!pdgPart)continue;
 707     if(pdgPart->Charge()==0) continue; //no neutral particles
 708     Float_t yMC=-999.;
 709     if(part->Energy() != TMath::Abs(part->Pz())) yMC = 0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
 710     if(pdgPart->Charge()>0) signMC=1;
 711     else signMC=-1;
 712     ptMC=part->Pt();
 713     code=pdgPart->PdgCode();
 714     Int_t jpart=-1;
 715     for(Int_t j=0; j<3; j++){
 716       if(TMath::Abs(code)==listcode[j]){
 717         jpart=j;
 718         break;
 719       }
 720     }
 721
 722     if(jpart>=0 && stack->IsPhysicalPrimary(imc) && TMath::Abs(part->Eta())< fEtaRange){
 723       if(signMC>0) fHistPrimMCposBefEvSelEta[jpart]->Fill(ptMC);
 724       else  fHistPrimMCnegBefEvSelEta[jpart]->Fill(ptMC);
 725     }
 726
 727     if(TMath::Abs(yMC) > fMaxY) continue; //rapidity cut
 728
 729     Int_t indexMoth=part->GetFirstMother();
 730     if(indexMoth>=0){
 731       TParticle* moth = stack->Particle(indexMoth);
 732       Float_t codemoth = TMath::Abs(moth->GetPdgCode());
 733       mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
 734     }
 735     uniqueID = part->GetUniqueID();
 736
 737     //filling MC ntuple
 738     if(TMath::Abs(code)==211 || TMath::Abs(code)==321 || TMath::Abs(code)==2212){
 739       Float_t xntMC[8];
 740       Int_t indexMC=0;
 741       xntMC[indexMC++]=(Float_t)ptMC;
 742       xntMC[indexMC++]=(Float_t)code;
 743       xntMC[indexMC++]=(Float_t)signMC;
 744       xntMC[indexMC++]=(Float_t)part->Eta();
 745       xntMC[indexMC++]=(Float_t)yMC;
 746       xntMC[indexMC++]=(Float_t)isph;
 747       xntMC[indexMC++]=(Float_t)evSel;
 748       xntMC[indexMC++]=(Float_t)fESD->GetRunNumber();
 749
 750       if(fFillNtuple) fNtupleMC->Fill(xntMC);
 751     }
 752
 753     if(jpart>=0){
 754       if(stack->IsPhysicalPrimary(imc)){
 755         if(signMC>0){
 756           fHistPrimMCposBefEvSel[jpart]->Fill(ptMC);
 757           if(TMath::Abs(part->Eta())<fEtaRange) fHistPrimMCposBefEvSelEtaY[jpart]->Fill(ptMC);
 758         }
 759         else{
 760           fHistPrimMCnegBefEvSel[jpart]->Fill(ptMC);
 761           if(TMath::Abs(part->Eta())<fEtaRange) fHistPrimMCnegBefEvSelEtaY[jpart]->Fill(ptMC);
 762         }
 763         if(evSel==1){
 764           if(signMC>0) fHistPrimMCpos[jpart]->Fill(ptMC);
 765           else  fHistPrimMCneg[jpart]->Fill(ptMC);
 766         }
 767       }else{
 768         if(mfl==3 && uniqueID == kPDecay){ // If a particle is not a physical primary, check if it comes from weak decay
 769           if(signMC>0) fHistSecStrMCposBefEvSel[jpart]->Fill(ptMC);
 770           else  fHistSecStrMCnegBefEvSel[jpart]->Fill(ptMC);
 771           if(evSel==1){
 772             if(signMC>0) fHistSecStrMCpos[jpart]->Fill(ptMC);
 773             else  fHistSecStrMCneg[jpart]->Fill(ptMC);
 774           }
 775         }else{
 776           if(signMC>0) fHistSecMatMCposBefEvSel[jpart]->Fill(ptMC);
 777           else  fHistSecMatMCnegBefEvSel[jpart]->Fill(ptMC);
 778           if(evSel==1){
 779             if(signMC>0) fHistSecMatMCpos[jpart]->Fill(ptMC);
 780             else  fHistSecMatMCneg[jpart]->Fill(ptMC);
 781           }
 782         }
 783       }
 784     }
 785   }
 786
 787
 788   if(evSel==0)return;  //event selection
 789
 790
 791   //loop on tracks
 792   for (Int_t iTrack=0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
 793     isph=-999;
 794     code=-999;
 795     mfl=-999;
 796     uniqueID=-999;
 797
 798     track = (AliESDtrack*)fESD->GetTrack(iTrack);
 799     if (!track) continue;
 800
 801     //track selection
 802     Int_t countBinTrk=1;
 803     TString label;
 804
 805     label="no selection";
 806     fHistNTracks->Fill(countBinTrk);
 807     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
 808     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
 809     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
 810     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
 811     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
 812     countBinTrk++;
 813
 814     status=track->GetStatus();
 815     if((status&AliESDtrack::kITSpureSA)==0) continue; //its standalone
 816
 817     label="ITSsa";
 818     fHistNTracks->Fill(countBinTrk);
 819     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
 820     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
 821     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
 822     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
 823     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
 824     countBinTrk++;
 825
 826     if((status&AliESDtrack::kITSrefit)==0) continue; //its refit
 827
 828     label="ITSrefit";
 829     fHistNTracks->Fill(countBinTrk);
 830     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
 831     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
 832     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
 833     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
 834     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
 835     countBinTrk++;
 836
 837     if(TMath::Abs(track->GetSign())<0.0001) continue; //no neutral particles
 838
 839     label="neutral particle";
 840     fHistNTracks->Fill(countBinTrk);
 841     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
 842     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
 843     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
 844     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
 845     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
 846     countBinTrk++;
 847
 848     //cluster in ITS
 849     UInt_t clumap = track->GetITSClusterMap();
 850     Int_t nSPD=0;
 851     for(Int_t il=0; il<2; il++) if(TESTBIT(clumap,il)) nSPD++;
 852     if(nSPD<fMinSPDPts) continue;
 853
 854     label="SPDcls";
 855     fHistNTracks->Fill(countBinTrk);
 856     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
 857     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
 858     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
 859     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
 860     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
 861     countBinTrk++;
 862
 863     Int_t nPtsForPid=0;
 864     for(Int_t j=2;j<6;j++) if(TESTBIT(clumap,j)) nPtsForPid++;
 865     if(nPtsForPid<fMinNdEdxSamples) continue; //at least 3 points on SSD/SDD
 866
 867     label="SDD+SSD cls";
 868     fHistNTracks->Fill(countBinTrk);
 869     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
 870     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
 871     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
 872     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
 873     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
 874     countBinTrk++;
 875
 876     //chisquare/nclusters
 877     Int_t nclu=nSPD+nPtsForPid;
 878     if(track->GetITSchi2()/nclu > fMaxChi2Clu) continue;
 879
 880     label="chi2/ncls";
 881     fHistNTracks->Fill(countBinTrk);
 882     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
 883     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
 884     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
 885     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
 886     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
 887     countBinTrk++;
 888
 889     //pseudorapidity and rapidity
 890     if(TMath::Abs(track->Eta()) > fEtaRange) continue;
 891
 892     label="eta";
 893     fHistNTracks->Fill(countBinTrk);
 894     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
 895     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
 896     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
 897     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
 898     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
 899     countBinTrk++;
 900
 901     //truncated mean
 902     track->GetITSdEdxSamples(dedxLay);
 903     Double_t dedx = CookdEdx(dedxLay);
 904     if(dedx<0) continue;
 905
 906     label="de/dx<0";
 907     fHistNTracks->Fill(countBinTrk);
 908     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
 909     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
 910     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
 911     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
 912     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
 913     countBinTrk++;
 914
 915     Float_t pt = track->Pt();
 916     Int_t theBin=-1;
 917     for(Int_t m=0; m<kNbins; m++){
 918       if(TMath::Abs(pt) > fPtBinLimits[m] && TMath::Abs(pt) < fPtBinLimits[m+1]){
 919         theBin=m;
 920         break;
 921       }
 922     }
 923     track->GetImpactParameters(impactXY, impactZ);
 924
 925     //Filling Ntuple
 926     //information from the MC kinematics
 927     if(fMC){
 928       if(track->GetLabel()<0)isph=-1;
 929       if(track->GetLabel()>=0){
 930         part = (TParticle*)stack->Particle(track->GetLabel());
 931         pdgPart = part->GetPDG();
 932         code = pdgPart->PdgCode();
 933         if(stack->IsPhysicalPrimary(track->GetLabel())) isph=1;
 934         else{
 935           isph=0;
 936           Int_t indexMoth=part->GetFirstMother();
 937           if(indexMoth>=0){
 938             TParticle* moth = stack->Particle(indexMoth);
 939             Float_t codemoth = TMath::Abs(moth->GetPdgCode());
 940             mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
 941           }
 942           uniqueID = part->GetUniqueID();
 943         }
 944
 945         //Filling DCA distribution with MC truth
 946
 947         if(theBin>=0 && theBin<kNbins){
 948           if(isph==1){//primaries in MC
 949             if(track->GetSign()>0){
 950               if(TMath::Abs(code)==listcode[0]) fHistMCPrimDCAPosPi[theBin]->Fill(impactXY);
 951               if(TMath::Abs(code)==listcode[1]) fHistMCPrimDCAPosK[theBin]->Fill(impactXY);
 952               if(TMath::Abs(code)==listcode[2]) fHistMCPrimDCAPosP[theBin]->Fill(impactXY);
 953             }else{
 954               if(TMath::Abs(code)==listcode[0]) fHistMCPrimDCANegPi[theBin]->Fill(impactXY);
 955               if(TMath::Abs(code)==listcode[1]) fHistMCPrimDCANegK[theBin]->Fill(impactXY);
 956               if(TMath::Abs(code)==listcode[2]) fHistMCPrimDCANegP[theBin]->Fill(impactXY);
 957             }
 958           }
 959
 960           if(isph==0){//primaries in MC
 961             if(mfl==3 && uniqueID == kPDecay){
 962               if(track->GetSign()>0){
 963                 if(TMath::Abs(code)==listcode[0]) fHistMCSecStDCAPosPi[theBin]->Fill(impactXY);
 964                 if(TMath::Abs(code)==listcode[1]) fHistMCSecStDCAPosK[theBin]->Fill(impactXY);
 965                 if(TMath::Abs(code)==listcode[2]) fHistMCSecStDCAPosP[theBin]->Fill(impactXY);
 966               }else{
 967                 if(TMath::Abs(code)==listcode[0]) fHistMCSecStDCANegPi[theBin]->Fill(impactXY);
 968                 if(TMath::Abs(code)==listcode[1]) fHistMCSecStDCANegK[theBin]->Fill(impactXY);
 969                 if(TMath::Abs(code)==listcode[2]) fHistMCSecStDCANegP[theBin]->Fill(impactXY);
 970               }
 971             }else{
 972               if(track->GetSign()>0){
 973                 if(TMath::Abs(code)==listcode[0]) fHistMCSecMatDCAPosPi[theBin]->Fill(impactXY);
 974                 if(TMath::Abs(code)==listcode[1]) fHistMCSecMatDCAPosK[theBin]->Fill(impactXY);
 975                 if(TMath::Abs(code)==listcode[2]) fHistMCSecMatDCAPosP[theBin]->Fill(impactXY);
 976               }else{
 977                 if(TMath::Abs(code)==listcode[0]) fHistMCSecMatDCANegPi[theBin]->Fill(impactXY);
 978                 if(TMath::Abs(code)==listcode[1]) fHistMCSecMatDCANegK[theBin]->Fill(impactXY);
 979                 if(TMath::Abs(code)==listcode[2]) fHistMCSecMatDCANegP[theBin]->Fill(impactXY);
 980               }
 981             }
 982           }
 983         }
 984       }
 985     }
 986     Float_t xnt[19];
 987     Int_t index=0;
 988     xnt[index++]=(Float_t)track->GetP();
 989     xnt[index++]=(Float_t)track->Pt();
 990     xnt[index++]=(Float_t)dedx;
 991     xnt[index++]=(Float_t)dedxLay[0];
 992     xnt[index++]=(Float_t)dedxLay[1];
 993     xnt[index++]=(Float_t)dedxLay[2];
 994     xnt[index++]=(Float_t)dedxLay[3];
 995     xnt[index++]=(Float_t)nclu;
 996     xnt[index++]=(Float_t)nPtsForPid;
 997     xnt[index++]=(Float_t)track->GetSign();
 998     xnt[index++]=(Float_t)fESD->GetRunNumber();
 999     xnt[index++]=(Float_t)track->Eta();
1000     xnt[index++]=(Float_t)impactXY;
1001     xnt[index++]=(Float_t)impactZ;
1002     xnt[index++]=(Float_t)isph;
1003     xnt[index++]=(Float_t)code;
1004     xnt[index++]=(Float_t)mfl;
1005     xnt[index++]=(Float_t)uniqueID;
1006     xnt[index]=(Float_t)track->GetITSchi2()/nclu;
1007
1008     if(fFillNtuple) fNtupleNSigma->Fill(xnt);
1009
1010
1011
1012     //Compute y and bb
1013     Double_t y[4],bbtheo[4],logdiff[4];
1014     Float_t p=track->GetP();
1015     if(fMC && fSmearMC){
1016       dedx=fRandGener->Gaus(dedx,fSmeardEdx*dedx);
1017       p=fRandGener->Gaus(p,fSmearP*p);
1018     }
1019
1020     //Nsigma Method
1021     Float_t resodedx[4];
1022     for(Int_t ires=0;ires<4;ires++){
1023       resodedx[ires]=fITSPIDResponse->GetResolution(1,ires+1,kTRUE);
1024     }
1025
1026     for(Int_t i=0;i<4;i++){
1027       y[i] = Eta2y(pt,pdgmass[i],track->Eta());
1028       //bbtheo[i]=fITSPIDResponse->Bethe(p,pdgmass[i],kTRUE); //Pure PHOBOS BB
1029       bbtheo[i]=fITSPIDResponse->BetheITSsaHybrid(p,pdgmass[i]); //PHOBOS + polinomial at low pt (below beta*gamma = 0.76)
1030       logdiff[i]=TMath::Log(dedx) - TMath::Log(bbtheo[i]);
1031     }
1032
1033     Int_t resocls=(Int_t)nPtsForPid-1;
1034
1035     //NSigma Method, with asymmetric bands
1036     Int_t minPosMean=-1;
1037     for(Int_t isp=0; isp<3; isp++){
1038       if(dedx<bbtheo[0])continue;
1039       Double_t bb=bbtheo[isp];
1040       if(dedx<bb){
1041         Double_t bbdistance=TMath::Abs((bbtheo[isp]-bbtheo[isp-1])/2);
1042         Double_t nsigma=TMath::Abs((dedx-bb)/bbdistance);
1043         if(nsigma<1.)minPosMean=isp;
1044       }
1045       else{
1046         Double_t bbdistance=TMath::Abs((bbtheo[isp]-bbtheo[isp+1])/2);
1047         Double_t nsigma=TMath::Abs((dedx-bb)/bbdistance);
1048         if(nsigma<1.)minPosMean=isp;
1049       }
1050     }
1051     if(dedx<bbtheo[0] && TMath::Abs((dedx-bbtheo[0])/(resodedx[resocls]*bbtheo[0]))<2)minPosMean=0;
1052
1053     //NSigma method with simmetric bands
1054
1055     Double_t nsigmas[3];
1056     Double_t min=999999.;
1057     Int_t minPos=-1;
1058     for(Int_t isp=0; isp<3; isp++){
1059       Double_t bb=bbtheo[isp];
1060       nsigmas[isp]=TMath::Abs((dedx-bb)/(resodedx[resocls]*bb));
1061       if(nsigmas[isp]<min){
1062         min=nsigmas[isp];
1063         minPos=isp;
1064       }
1065       //Filling histos with nsigma separation
1066       if(track->GetSign()>0)fHistPosNSigmaSep[isp]->Fill(track->GetP(),((dedx-bb)/(resodedx[resocls]*bb)));
1067       else fHistNegNSigmaSep[isp]->Fill(track->GetP(),((dedx-bb)/(resodedx[resocls]*bb)));
1068     }
1069
1070     // y calculation
1071     Double_t yPartMean=-999;
1072     Double_t yPart=-999;
1073     if(minPosMean>-1) yPartMean=y[minPosMean];
1074     if(minPos>-1) yPart=y[minPos];
1075
1076     if(TMath::Abs(yPartMean)<fMaxY){
1077       //DCA distributions, before the DCA cuts, based on asymmetrinc nsigma approach
1078       if(theBin>=0 && theBin<kNbins){
1079         if(track->GetSign()>0){
1080           if(minPosMean==0) fHistDCAPosPi[theBin]->Fill(impactXY);
1081           else if(minPosMean==1) fHistDCAPosK[theBin]->Fill(impactXY);
1082           else if(minPosMean==2) fHistDCAPosP[theBin]->Fill(impactXY);
1083         }else{
1084           if(minPosMean==0) fHistDCANegPi[theBin]->Fill(impactXY);
1085           else if(minPosMean==1) fHistDCANegK[theBin]->Fill(impactXY);
1086           else if(minPosMean==2) fHistDCANegP[theBin]->Fill(impactXY);
1087         }
1088       }
1089     }
1090
1091     //DCA cut on xy and z
1092     if(!DCAcut(impactXY,impactZ,pt)) continue;
1093
1094     label="DCA";
1095     fHistNTracks->Fill(countBinTrk);
1096     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
1097     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
1098     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
1099     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
1100     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
1101     countBinTrk++;
1102
1103
1104     Int_t jpart=-1;
1105
1106     //Filling Histos for Reco Efficiency
1107     //information from the MC kinematics
1108     if(fMC){
1109       if(track->GetLabel()<0)isph=-1;
1110       if(track->GetLabel()>=0){
1111         part = (TParticle*)stack->Particle(track->GetLabel());
1112         pdgPart = part->GetPDG();
1113         code = pdgPart->PdgCode();
1114         for(Int_t j=0; j<3; j++){
1115           if(TMath::Abs(code)==listcode[j]){
1116             jpart=j;
1117             break;
1118           }
1119         }
1120         if(jpart>=0){
1121           if(pdgPart->Charge()>0) signMC=1;
1122           else signMC=-1;
1123           ptMC=part->Pt();
1124           if(stack->IsPhysicalPrimary(track->GetLabel())){
1125             if(signMC>0){
1126               fHistPrimMCposReco[jpart]->Fill(ptMC);
1127               if(TMath::Abs(y[jpart])<fMaxY) fHistPrimMCposRecoEtaY[jpart]->Fill(ptMC);
1128             }
1129             else{
1130               fHistPrimMCnegReco[jpart]->Fill(ptMC);
1131               if(TMath::Abs(y[jpart])<fMaxY) fHistPrimMCnegRecoEtaY[jpart]->Fill(ptMC);
1132             }
1133           }else{
1134             Int_t indexMoth=part->GetFirstMother();
1135             if(indexMoth>=0){
1136               TParticle* moth = stack->Particle(indexMoth);
1137               Float_t codemoth = TMath::Abs(moth->GetPdgCode());
1138               mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
1139             }
1140             uniqueID = part->GetUniqueID();
1141             if(mfl==3 && uniqueID == kPDecay){ // strangeness
1142               if(signMC>0) fHistSecStrMCposReco[jpart]->Fill(ptMC);
1143               else  fHistSecStrMCnegReco[jpart]->Fill(ptMC);
1144             }else{
1145               if(signMC>0) fHistSecMatMCposReco[jpart]->Fill(ptMC);
1146               else  fHistSecMatMCnegReco[jpart]->Fill(ptMC);
1147             }
1148           }
1149         }
1150       }
1151     }
1152
1153     //Nsigma histos with MC truth
1154
1155     //asymmetric approach
1156     if(TMath::Abs(yPartMean)<fMaxY && minPosMean>-1){
1157       //nsigma histos
1158       if(track->GetSign()>0) fHistPosNSigmaMean[minPosMean]->Fill(pt);
1159       else fHistNegNSigmaMean[minPosMean]->Fill(pt);
1160       if(fMC){
1161         //nsigma histos with MC truth on PID
1162         if(TMath::Abs(code)==listcode[minPosMean]){
1163           if(track->GetSign()>0) fHistPosNSigmaMCMean[minPosMean]->Fill(pt);
1164           else fHistNegNSigmaMCMean[minPosMean]->Fill(pt);
1165         }
1166         //nsigma histos with MC truth on IsPhysicalPrimary
1167         if(isph==1){
1168           if(track->GetSign()>0) fHistPosNSigmaPrimMean[minPosMean]->Fill(pt);
1169           else fHistNegNSigmaPrimMean[minPosMean]->Fill(pt);
1170           //nsigma histos with MC truth on IsPhysicalPrimary and PID
1171           if(TMath::Abs(code)==listcode[minPosMean]){
1172             if(track->GetSign()>0) fHistPosNSigmaPrimMCMean[minPosMean]->Fill(pt);
1173             else fHistNegNSigmaPrimMCMean[minPosMean]->Fill(pt);
1174           }
1175         }
1176       }
1177     }
1178
1179     //symmetric bands
1180     if(min<fMinNSigma && TMath::Abs(yPart)<fMaxY){
1181       //nsigma histos
1182       if(track->GetSign()>0) fHistPosNSigma[minPos]->Fill(pt);
1183       else fHistNegNSigma[minPos]->Fill(pt);
1184       if(fMC){
1185         //nsigma histos with MC truth on PID
1186         if(TMath::Abs(code)==listcode[minPos]){
1187           if(track->GetSign()>0) fHistPosNSigmaMC[minPos]->Fill(pt);
1188           else fHistNegNSigmaMC[minPos]->Fill(pt);
1189         }
1190         //nsigma histos with MC truth on IsPhysicalPrimary
1191         if(isph==1){
1192           if(track->GetSign()>0) fHistPosNSigmaPrim[minPos]->Fill(pt);
1193           else fHistNegNSigmaPrim[minPos]->Fill(pt);
1194           //nsigma histos with MC truth on IsPhysicalPrimary and PID
1195           if(TMath::Abs(code)==listcode[minPos]){
1196             if(track->GetSign()>0) fHistPosNSigmaPrimMC[minPos]->Fill(pt);
1197             else fHistNegNSigmaPrimMC[minPos]->Fill(pt);
1198           }
1199         }
1200       }
1201     }
1202
1203
1204     //integral approach histograms
1205     if(theBin>=0 && theBin<kNbins){
1206       if(track->GetSign()>0){
1207         if(TMath::Abs(y[0]) < fMaxY)fHistPosPi[theBin]->Fill(logdiff[0]);
1208         if(TMath::Abs(y[1]) < fMaxY)fHistPosK[theBin]->Fill(logdiff[1]);
1209         if(TMath::Abs(y[2]) < fMaxY)fHistPosP[theBin]->Fill(logdiff[2]);
1210         if(fMC){
1211           if(TMath::Abs(y[0])<fMaxY){
1212             if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypPion[theBin]->Fill(logdiff[0]);
1213             if(TMath::Abs(code)==11)fHistMCPosElHypPion[theBin]->Fill(logdiff[0]);
1214             if(TMath::Abs(code)==211)fHistMCPosPiHypPion[theBin]->Fill(logdiff[0]);
1215             if(TMath::Abs(code)==321)fHistMCPosKHypPion[theBin]->Fill(logdiff[0]);
1216             if(TMath::Abs(code)==2212)fHistMCPosPHypPion[theBin]->Fill(logdiff[0]);
1217           }
1218           if(TMath::Abs(y[1])<fMaxY){
1219             if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypKaon[theBin]->Fill(logdiff[1]);
1220             if(TMath::Abs(code)==11)fHistMCPosElHypKaon[theBin]->Fill(logdiff[1]);
1221             if(TMath::Abs(code)==211)fHistMCPosPiHypKaon[theBin]->Fill(logdiff[1]);
1222             if(TMath::Abs(code)==321)fHistMCPosKHypKaon[theBin]->Fill(logdiff[1]);
1223             if(TMath::Abs(code)==2212)fHistMCPosPHypKaon[theBin]->Fill(logdiff[1]);
1224           }
1225           if(TMath::Abs(y[2])<fMaxY){
1226               if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypProton[theBin]->Fill(logdiff[2]);
1227               if(TMath::Abs(code)==11)fHistMCPosElHypProton[theBin]->Fill(logdiff[2]);
1228               if(TMath::Abs(code)==211)fHistMCPosPiHypProton[theBin]->Fill(logdiff[2]);
1229               if(TMath::Abs(code)==321)fHistMCPosKHypProton[theBin]->Fill(logdiff[2]);
1230               if(TMath::Abs(code)==2212)fHistMCPosPHypProton[theBin]->Fill(logdiff[2]);
1231           }
1232         }
1233       }else{
1234         if(TMath::Abs(y[0]) < fMaxY)fHistNegPi[theBin]->Fill(logdiff[0]);
1235         if(TMath::Abs(y[1]) < fMaxY)fHistNegK[theBin]->Fill(logdiff[1]);
1236         if(TMath::Abs(y[2]) < fMaxY)fHistNegP[theBin]->Fill(logdiff[2]);
1237         if(fMC){
1238           if(TMath::Abs(y[0])<fMaxY){
1239             if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypPion[theBin]->Fill(logdiff[0]);
1240             if(TMath::Abs(code)==11)fHistMCNegElHypPion[theBin]->Fill(logdiff[0]);
1241             if(TMath::Abs(code)==211)fHistMCNegPiHypPion[theBin]->Fill(logdiff[0]);
1242             if(TMath::Abs(code)==321)fHistMCNegKHypPion[theBin]->Fill(logdiff[0]);
1243             if(TMath::Abs(code)==2212)fHistMCNegPHypPion[theBin]->Fill(logdiff[0]);
1244           }
1245           if(TMath::Abs(y[1])<fMaxY){
1246             if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypKaon[theBin]->Fill(logdiff[1]);
1247             if(TMath::Abs(code)==11)fHistMCNegElHypKaon[theBin]->Fill(logdiff[1]);
1248             if(TMath::Abs(code)==211)fHistMCNegPiHypKaon[theBin]->Fill(logdiff[1]);
1249             if(TMath::Abs(code)==321)fHistMCNegKHypKaon[theBin]->Fill(logdiff[1]);
1250             if(TMath::Abs(code)==2212)fHistMCNegPHypKaon[theBin]->Fill(logdiff[1]);
1251           }
1252           if(TMath::Abs(y[2])<fMaxY){
1253             if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypProton[theBin]->Fill(logdiff[2]);
1254             if(TMath::Abs(code)==11)fHistMCNegElHypProton[theBin]->Fill(logdiff[2]);
1255             if(TMath::Abs(code)==211)fHistMCNegPiHypProton[theBin]->Fill(logdiff[2]);
1256             if(TMath::Abs(code)==321)fHistMCNegKHypProton[theBin]->Fill(logdiff[2]);
1257             if(TMath::Abs(code)==2212)fHistMCNegPHypProton[theBin]->Fill(logdiff[2]);
1258           }
1259         }
1260       }
1261     }
1262
1263
1264     //fill propaganda plot with dedx
1265     fHistDEDX->Fill(track->GetP(),dedx);
1266     fHistDEDXdouble->Fill(track->GetP()*track->GetSign(),dedx);
1267
1268     //fill charge distribution histo to check the calibration
1269     for(Int_t j=0;j<4;j++){
1270       if(dedxLay[j]<5) continue;
1271       fHistCharge[j]->Fill(dedxLay[j]);
1272     }
1273   }
1274
1275   // Post output data.
1276   PostData(1,fOutput);
1277   PostData(2,fListCuts);
1278   //  Printf("............. end of Exec");
1279 }
1280
1281 //________________________________________________________________________
1282 void AliAnalysisTaskSEITSsaSpectra::Terminate(Option_t *) {
1283   // Merge output
1284   // Called once at the end of the query
1285
1286   fOutput = dynamic_cast<TList*>(GetOutputData(1));
1287   if (!fOutput) {
1288     printf("ERROR: fOutput not available\n");
1289     return;
1290   }
1291   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1292   printf("Number of Analyzed Events = %f\n",fHistNEvents->GetBinContent(1));
1293
1294   Printf("end of Terminate");
1295   return;
1296 }