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