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