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