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 <TStyle.h>
  28 #include <TSystem.h>
  29 #include <TTree.h>
  30 #include <TH1F.h>
  31 #include <TRandom3.h>
  32 #include <TH2F.h>
  33 #include <TChain.h>
  34 #include <TNtuple.h>
  35 #include <TParticle.h>
  36 #include <Rtypes.h>
  37 #include "AliAnalysisTaskSE.h"
  38 #include "AliAnalysisManager.h"
  39 #include "AliAnalysisDataContainer.h"
  40 #include "AliESDEvent.h"
  41 #include "AliESDInputHandler.h"
  42 #include "AliESDtrack.h"
  43 #include "AliStack.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliMCEvent.h"
  46 #include "AliPhysicsSelection.h"
  47 #include "AliAnalysisTaskSEITSsaSpectra.h"
  48
  49 ClassImp(AliAnalysisTaskSEITSsaSpectra)
  50
  51 //________________________________________________________________________
  52 AliAnalysisTaskSEITSsaSpectra::AliAnalysisTaskSEITSsaSpectra():
  53 AliAnalysisTaskSE("Task CFits"),
  54   fESD(0),
  55   fOutput(0),
  56   fHistNEvents(0),
  57   fHistNTracks(0),
  58   fHistDEDX(0),
  59   fHistDEDXdouble(0),
  60   fHistBeforeEvSel(0),
  61   fHistAfterEvSel(0),
  62   fMinSPDPts(1),
  63   fMinNdEdxSamples(3),
  64   fMindEdx(0.),
  65   fMinNSigma(3.),
  66   fMaxY(0.5),
  67   fMaxChi2Clu(1.),
  68   fNSigmaDCAxy(7.),
  69   fNSigmaDCAz(7.),
  70   fMC(kFALSE),
  71   fSmearMC(kFALSE),
  72   fSmearP(0.),
  73   fSmeardEdx(0.),
  74   fRandGener(0),
  75   fFillNtuple(kFALSE),
  76   fNtupleNSigma(0),
  77   fNtupleMC(0)
  78 {
  79   // Constructor
  80   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};
  81   for(Int_t iBin=0; iBin<kNbins+1; iBin++) fPtBinLimits[iBin]=xbins[iBin];
  82   fRandGener=new TRandom3(0);
  83
  84   DefineInput(0, TChain::Class());
  85   DefineOutput(1, TList::Class());
  86   Printf("end of AliAnalysisTaskSEITSsaSpectra");
  87 }
  88
  89 //___________________________________________________________________________
  90 AliAnalysisTaskSEITSsaSpectra::~AliAnalysisTaskSEITSsaSpectra(){
  91   // Destructor
  92   if (fOutput) {
  93     delete fOutput;
  94     fOutput = 0;
  95   }
  96   if(fRandGener) delete fRandGener;
  97 }
  98
  99 //________________________________________________________________________
 100 Double_t AliAnalysisTaskSEITSsaSpectra::CookdEdx(Double_t *s){
 101   // truncated mean for the dEdx
 102   Int_t nc=0;
 103   Double_t dedx[4]={0.,0.,0.,0.};
 104   for (Int_t il=0; il<4; il++) { // count good (>0) dE/dx values
 105     if(s[il]>fMindEdx){
 106       dedx[nc]= s[il];
 107       nc++;
 108     }
 109   }
 110   if(nc<fMinNdEdxSamples) return -1.;
 111
 112   Double_t tmp;
 113   Int_t swap; // sort in ascending order
 114   do {
 115     swap=0;
 116     for (Int_t i=0; i<nc-1; i++) {
 117       if (dedx[i]<=dedx[i+1]) continue;
 118       tmp=dedx[i];
 119       dedx[i]=dedx[i+1];
 120       dedx[i+1]=tmp;
 121       swap++;
 122     }
 123   } while (swap);
 124
 125   Double_t sumamp=0,sumweight=0;
 126   Double_t weight[4]={1.,1.,0.,0.};
 127   if(nc==3) weight[1]=0.5;
 128   else if(nc<3) weight[1]=0.;
 129   for (Int_t i=0; i<nc; i++) {
 130     sumamp+= dedx[i]*weight[i];
 131     sumweight+=weight[i];
 132   }
 133   return sumamp/sumweight;
 134 }
 135
 136
 137 //________________________________________________________________________
 138 Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt, Bool_t optMC){
 139   // cut on transverse impact parameter updaated on 20-5-2010
 140   // from the study of L. Milano, F. Prino on the ITS standalone tracks
 141   // using the common binning of the TPC tracks
 142   Double_t xyP[3];
 143   Double_t zP[3];
 144   if(optMC){
 145     xyP[0]=88.63;//SIMpass4a12
 146     xyP[1]=19.57;
 147     xyP[2]=1.65;
 148     zP[0]=140.98;
 149     zP[1]=62.33;
 150     zP[2]=1.15;
 151   }
 152   else{
 153     xyP[0]=85.28;//DATApass6
 154     xyP[1]=25.78;
 155     xyP[2]=1.55;
 156     zP[0]=146.80;
 157     zP[1]=70.07;
 158     zP[2]=1.11;
 159   }
 160   Double_t xySigma = xyP[0] + xyP[1]/TMath::Power(TMath::Abs(pt),xyP[2]);
 161   Double_t xyMax = fNSigmaDCAxy*xySigma; //in micron
 162   if((TMath::Abs(impactXY)*10000)>xyMax) return kFALSE;
 163
 164   Double_t zSigma = zP[0] + zP[1]/TMath::Power(TMath::Abs(pt),zP[2]);
 165   Double_t zMax = fNSigmaDCAz*zSigma; //in micron
 166   if((TMath::Abs(impactZ)*10000)>zMax) return kFALSE;
 167
 168   return kTRUE;
 169 }
 170
 171 //________________________________________________________________________
 172 Double_t AliAnalysisTaskSEITSsaSpectra::Eta2y(Double_t pt, Double_t m, Double_t eta){
 173   Double_t mt = TMath::Sqrt(m*m + pt*pt);
 174   return TMath::ASinH(pt/mt*TMath::SinH(eta));
 175 }
 176
 177
 178 //________________________________________________________________________
 179 Double_t AliAnalysisTaskSEITSsaSpectra::BetheBloch(Double_t bg,Bool_t optMC) {
 180   // BB PHOBOS parametrization tuned by G. Ortona on 900 GeV pp data of 2009
 181   Double_t par[5];
 182   if(optMC){//Double_t par[5]={139.1,23.36,0.06052,0.2043,-0.0004999};
 183     par[0]=139.1;
 184     par[1]=23.36;
 185     par[2]=0.06052;
 186     par[3]=0.2043;
 187     par[4]=-0.0004999;
 188   }else {
 189     //Double_t par[5]={5.33458e+04,1.65303e+01,2.60065e-03,3.59533e-04,7.51168e-05};}
 190   par[0]=5.33458e+04;
 191   par[1]=1.65303e+01;
 192   par[2]=2.60065e-03;
 193   par[3]=3.59533e-04;
 194   par[4]=7.51168e-05;
 195   }
 196   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 197   Double_t gamma=bg/beta;
 198   Double_t eff=1.0;
 199   if(bg<par[2]) eff=(bg-par[3])*(bg-par[3])+par[4];
 200   else eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 201   return (par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 202 }
 203
 204
 205 //________________________________________________________________________
 206 void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){
 207   fOutput = new TList();
 208   fOutput->SetOwner();
 209   fOutput->SetName("Spiderman");
 210
 211   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events",6,0.5,6.5);
 212   fHistNEvents->Sumw2();
 213   fHistNEvents->SetMinimum(0);
 214   fOutput->Add(fHistNEvents);
 215
 216   fHistNTracks = new TH1F("fHistNTracks", "Number of ITSsa tracks",20,0.5,20.5);
 217   fHistNTracks->Sumw2();
 218   fHistNTracks->SetMinimum(0);
 219   fOutput->Add(fHistNTracks);
 220
 221   //binning for the histogram
 222   const Int_t hnbins=400;
 223   Double_t hxmin = 0.01;
 224   Double_t hxmax = 10;
 225   Double_t hlogxmin = TMath::Log10(hxmin);
 226   Double_t hlogxmax = TMath::Log10(hxmax);
 227   Double_t hbinwidth = (hlogxmax-hlogxmin)/hnbins;
 228   Double_t hxbins[hnbins+1];
 229   hxbins[0] = 0.01;
 230   for (Int_t i=1;i<=hnbins;i++) {
 231     hxbins[i] = hxmin + TMath::Power(10,hlogxmin+i*hbinwidth);
 232   }
 233
 234   fHistDEDX = new TH2F("fHistDEDX","",hnbins,hxbins,900,0,1000);
 235   fOutput->Add(fHistDEDX);
 236
 237   fHistDEDXdouble = new TH2F("fHistDEDXdouble","",500,-5,5,900,0,1000);
 238   fOutput->Add(fHistDEDXdouble);
 239
 240
 241   fHistBeforeEvSel = new TH1F("fHistBeforeEvSel","fHistBeforeEvSel",kNbins,fPtBinLimits);
 242   fHistAfterEvSel = new TH1F("fHistAfterEvSel","fHistAfterEvSel",kNbins,fPtBinLimits);
 243   fOutput->Add(fHistBeforeEvSel);
 244   fOutput->Add(fHistAfterEvSel);
 245
 246
 247
 248   for(Int_t j=0;j<3;j++){
 249     fHistMCpos[j] = new TH1F(Form("fHistMCpos%d",j),Form("fHistMCpos%d",j),kNbins,fPtBinLimits);
 250     fHistMCneg[j] = new TH1F(Form("fHistMCneg%d",j),Form("fHistMCneg%d",j),kNbins,fPtBinLimits);
 251     fOutput->Add(fHistMCneg[j]);
 252                 fOutput->Add(fHistMCpos[j]);
 253   }
 254
 255
 256   for(Int_t j=0;j<3;j++){
 257     fHistMCposBefEvSel[j] = new TH1F(Form("fHistMCposBefEvSel%d",j),Form("fHistMCposBefEvSel%d",j),kNbins,fPtBinLimits);
 258     fHistMCnegBefEvSel[j] = new TH1F(Form("fHistMCnegBefEvSel%d",j),Form("fHistMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
 259     fOutput->Add(fHistMCnegBefEvSel[j]);
 260     fOutput->Add(fHistMCposBefEvSel[j]);
 261   }
 262
 263   for(Int_t i=0; i<4; i++){
 264     fHistCharge[i] = new TH1F(Form("fHistChargeLay%d",i),Form("fHistChargeLay%d",i),100,0,300);
 265     fOutput->Add(fHistCharge[i]);
 266   }
 267
 268   for(Int_t i=0; i<kNbins; i++){
 269     fHistPosPi[i] = new TH1F(Form("fHistPosPi%d",i),Form("fHistPosPi%d",i),175,-3.5,3.5);
 270     fHistPosK[i]  = new TH1F(Form("fHistPosK%d",i),Form("fHistPosK%d",i),175,-3.5,3.5);
 271     fHistPosP[i]  = new TH1F(Form("fHistPosP%d",i),Form("fHistPosP%d",i),175,-3.5,3.5);
 272     fHistNegPi[i] = new TH1F(Form("fHistNegPi%d",i),Form("fHistNegPi%d",i),175,-3.5,3.5);
 273     fHistNegK[i]  = new TH1F(Form("fHistNegK%d",i),Form("fHistNegK%d",i),175,-3.5,3.5);
 274     fHistNegP[i]  = new TH1F(Form("fHistNegP%d",i),Form("fHistNegP%d",i),175,-3.5,3.5);
 275
 276     fHistDCAPosPi[i] = new TH1F(Form("fHistDCAPosPi%d",i),Form("fHistDCAPosPi%d",i),2000,-1,1);  //DCA distr.
 277     fHistDCAPosK[i]  = new TH1F(Form("fHistDCAPosK%d",i),Form("fHistDCAPosK%d",i),2000,-1,1);
 278     fHistDCAPosP[i]  = new TH1F(Form("fHistDCAPosP%d",i),Form("fHistDCAPosP%d",i),2000,-1,1);
 279     fHistDCANegPi[i] = new TH1F(Form("fHistDCANegPi%d",i),Form("fHistDCANegPi%d",i),2000,-1,1);
 280     fHistDCANegK[i]  = new TH1F(Form("fHistDCANegK%d",i),Form("fHistDCANegK%d",i),2000,-1,1);
 281     fHistDCANegP[i]  = new TH1F(Form("fHistDCANegP%d",i),Form("fHistDCANegP%d",i),2000,-1,1);
 282
 283     fHistMCPosPi[i] = new TH1F(Form("fHistMCPosPi%d",i),Form("fHistMCPosPi%d",i),175,-3.5,3.5); //MC truth
 284     fHistMCPosK[i]  = new TH1F(Form("fHistMCPosK%d",i),Form("fHistMCPosK%d",i),175,-3.5,3.5);
 285     fHistMCPosP[i]  = new TH1F(Form("fHistMCPosP%d",i),Form("fHistMCPosP%d",i),175,-3.5,3.5);
 286     fHistMCNegPi[i] = new TH1F(Form("fHistMCNegPi%d",i),Form("fHistMCNegPi%d",i),175,-3.5,3.5);
 287     fHistMCNegK[i]  = new TH1F(Form("fHistMCNegK%d",i),Form("fHistMCNegK%d",i),175,-3.5,3.5);
 288     fHistMCNegP[i]  = new TH1F(Form("fHistMCNegP%d",i),Form("fHistMCNegP%d",i),175,-3.5,3.5);
 289     fOutput->Add(fHistPosPi[i]);
 290     fOutput->Add(fHistPosK[i]);
 291     fOutput->Add(fHistPosP[i]);
 292     fOutput->Add(fHistNegPi[i]);
 293     fOutput->Add(fHistNegK[i]);
 294     fOutput->Add(fHistNegP[i]);
 295
 296     fOutput->Add(fHistDCAPosPi[i]);//DCA distr.
 297     fOutput->Add(fHistDCAPosK[i]);
 298     fOutput->Add(fHistDCAPosP[i]);
 299     fOutput->Add(fHistDCANegPi[i]);
 300     fOutput->Add(fHistDCANegK[i]);
 301     fOutput->Add(fHistDCANegP[i]);
 302
 303     fOutput->Add(fHistMCPosPi[i]);//MC truth
 304     fOutput->Add(fHistMCPosK[i]);
 305     fOutput->Add(fHistMCPosP[i]);
 306     fOutput->Add(fHistMCNegPi[i]);
 307     fOutput->Add(fHistMCNegK[i]);
 308     fOutput->Add(fHistMCNegP[i]);
 309   }
 310
 311   //NSigma Histos
 312   for(Int_t j=0;j<3;j++){
 313     fHistPosNSigma[j] = new TH1F(Form("hHistPosNSigma%d",j),Form("hHistPosNSigma%d",j),kNbins,fPtBinLimits);
 314     fHistNegNSigma[j] = new TH1F(Form("hHistNegNSigma%d",j),Form("hHistNegNSigma%d",j),kNbins,fPtBinLimits);
 315     fHistPosNSigmaPrim[j] = new TH1F(Form("hHistPosNSigmaPrim%d",j),Form("hHistPosNSigmaPrim%d",j),kNbins,fPtBinLimits);
 316     fHistNegNSigmaPrim[j] = new TH1F(Form("hHistNegNSigmaPrim%d",j),Form("hHistNegNSigmaPrim%d",j),kNbins,fPtBinLimits);
 317     fHistPosNSigmaPrimMC[j] = new TH1F(Form("hHistPosNSigmaPrimMC%d",j),Form("hHistPosNSigmaPrimMC%d",j),kNbins,fPtBinLimits);
 318     fHistNegNSigmaPrimMC[j] = new TH1F(Form("hHistNegNSigmaPrimMC%d",j),Form("hHistNegNSigmaPrimMC%d",j),kNbins,fPtBinLimits);
 319     fOutput->Add(fHistPosNSigma[j]);
 320     fOutput->Add(fHistNegNSigma[j]);
 321     fOutput->Add(fHistPosNSigmaPrim[j]);
 322     fOutput->Add(fHistNegNSigmaPrim[j]);
 323     fOutput->Add(fHistPosNSigmaPrimMC[j]);
 324     fOutput->Add(fHistNegNSigmaPrimMC[j]);
 325   }
 326
 327   fNtupleNSigma = new TNtuple("fNtupleNSigma","fNtupleNSigma","p:pt:dedx:ncls:sign:run:eta:impactXY:impactZ:isph:pdgcode:mfl");
 328   fOutput->Add(fNtupleNSigma);
 329   fNtupleMC = new TNtuple("fNtupleMC","fNtupleMC","ptMC:pdgcode:signMC:etaMC:yMC:isph:evSel:run");
 330   fOutput->Add(fNtupleMC);
 331
 332   Printf("end of CreateOutputObjects");
 333 }
 334
 335 //________________________________________________________________________
 336 void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
 337
 338   fESD=(AliESDEvent*)InputEvent();
 339   if(!fESD) {
 340     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
 341     return;
 342   }
 343
 344   //binning for the dEdx distributions
 345
 346   //variables
 347   Float_t pdgmass[3]={0.13957,0.493677,0.938272}; //mass for pi, K, P (Gev/c^2)
 348   Int_t listcode[3]={211,321,2212};//code for pi, K, P (Gev/c^2)
 349   Double_t s[4];
 350   Float_t ptMC=-999,yMC=-999,code=-999, signMC=-999,isph=-999,mfl=-999.;
 351   Float_t impactXY=-999, impactZ=-999;
 352   Int_t evSel=1;
 353   AliESDtrack* track;
 354   UInt_t status;
 355   AliStack* stack=0;
 356   TParticle *part=0;
 357   TParticlePDG *pdgPart=0;
 358
 359   //Nsigma Method
 360   Float_t resodedx[4];
 361   if(fMC){
 362     resodedx[0]=0.13;
 363     resodedx[1]=0.13;
 364     resodedx[2]=0.134;
 365     resodedx[3]=0.127;
 366   }else{
 367     resodedx[0]=0.23;
 368     resodedx[1]=0.18;
 369     resodedx[2]=0.16;
 370     resodedx[3]=0.14;
 371   }
 372   /////////////////////
 373   if(fMC){
 374     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 375     if (!eventHandler) {
 376       Printf("ERROR: Could not retrieve MC event handler");
 377       return;
 378     }
 379     AliMCEvent* mcEvent = eventHandler->MCEvent();
 380     if (!mcEvent) {
 381       Printf("ERROR: Could not retrieve MC event");
 382       return;
 383     }
 384     stack = mcEvent->Stack();
 385     if (!stack) {
 386       printf("ERROR: stack not available\n");
 387       return;
 388     }
 389   }
 390   //flags for MC
 391   Int_t nTrackMC=0;
 392   if(stack) nTrackMC = stack->GetNtrack();
 393   const AliESDVertex *vtx =  fESD->GetPrimaryVertexSPD();
 394
 395   //event selection
 396   fHistNEvents->Fill(1);
 397   if(!vtx)evSel=0;
 398   else{
 399     fHistNEvents->Fill(2);
 400     if(vtx->GetNContributors()<0) evSel=0;
 401     else{
 402       fHistNEvents->Fill(3);
 403       if(TMath::Abs(vtx->GetZv())>10) evSel=0;
 404       else{
 405         fHistNEvents->Fill(4);
 406         if(vtx->GetZRes()>0.5) evSel=0;
 407         else{
 408           fHistNEvents->Fill(5);
 409           if(vtx->IsFromVertexerZ() && vtx->GetDispersion()>0.03) evSel=0;
 410           else fHistNEvents->Fill(6);
 411         }
 412       }
 413     }
 414   }
 415
 416   /////first loop on stack, before event selection, filling MC ntuple
 417
 418   for(Int_t imc=0; imc<nTrackMC; imc++){
 419     part = stack->Particle(imc);
 420     if(!stack->IsPhysicalPrimary(imc))continue;//no secondary in the MC sample
 421     isph=1.;
 422     pdgPart = part->GetPDG();
 423     if(pdgPart->Charge()==0) continue; //no neutral particles
 424     if(TMath::Abs(part->Eta()) > 0.9) continue; //pseudorapidity-acceptance cut
 425     if(part->Energy() != TMath::Abs(part->Pz())) yMC = 0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
 426     if(TMath::Abs(yMC) > fMaxY) continue; //rapidity cut
 427
 428     if(pdgPart->Charge()>0) signMC=1;
 429     else signMC=-1;
 430     ptMC=part->Pt();
 431     code=pdgPart->PdgCode();
 432
 433
 434     //filling MC ntuple
 435     if(TMath::Abs(code)==211 || TMath::Abs(code)==321 || TMath::Abs(code)==2212){
 436       Float_t xntMC[8];
 437       Int_t indexMC=0;
 438       xntMC[indexMC++]=(Float_t)ptMC;
 439       xntMC[indexMC++]=(Float_t)code;
 440       xntMC[indexMC++]=(Float_t)signMC;
 441       xntMC[indexMC++]=(Float_t)part->Eta();
 442       xntMC[indexMC++]=(Float_t)yMC;
 443       xntMC[indexMC++]=(Float_t)isph;
 444       xntMC[indexMC++]=(Float_t)evSel;
 445       xntMC[indexMC++]=(Float_t)fESD->GetRunNumber();
 446
 447       if(fFillNtuple) fNtupleMC->Fill(xntMC);
 448     }
 449
 450     for(Int_t j=0; j<3; j++){
 451       if(TMath::Abs(code)==listcode[j]){
 452         if(signMC>0) fHistMCposBefEvSel[j]->Fill(TMath::Abs(ptMC));
 453         else  fHistMCnegBefEvSel[j]->Fill(TMath::Abs(ptMC));
 454       }
 455     }
 456     if(evSel==1){
 457       for(Int_t j=0; j<3; j++){
 458         if(TMath::Abs(code)==listcode[j]){
 459           if(signMC>0) fHistMCpos[j]->Fill(TMath::Abs(ptMC));
 460           else  fHistMCneg[j]->Fill(TMath::Abs(ptMC));
 461         }
 462       }
 463     }
 464   }
 465
 466   if(evSel==0)return;
 467
 468   //loop on tracks
 469   for (Int_t iTrack=0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
 470     isph=-999.;
 471     code=-999;
 472     mfl=-999;
 473
 474     track = (AliESDtrack*)fESD->GetTrack(iTrack);
 475     if (!track) continue;
 476
 477     //track selection
 478     fHistNTracks->Fill(1);
 479     status=track->GetStatus();
 480     if((status&AliESDtrack::kITSpureSA)==0) continue; //its standalone
 481     fHistNTracks->Fill(2);
 482     if((status&AliESDtrack::kITSrefit)==0) continue; //its refit
 483     fHistNTracks->Fill(3);
 484     if(track->GetSign()==0.) continue; //no neutral particles
 485     fHistNTracks->Fill(4);
 486
 487
 488     //cluster in ITS
 489     UInt_t clumap = track->GetITSClusterMap();
 490     Int_t nSPD=0;
 491     for(Int_t il=0; il<2; il++) if(TESTBIT(clumap,il)) nSPD++;
 492     if(nSPD<fMinSPDPts) continue;
 493     fHistNTracks->Fill(5);
 494     Int_t count=0;
 495     for(Int_t j=2;j<6;j++) if(TESTBIT(clumap,j)) count++;
 496     if(count<fMinNdEdxSamples) continue; //at least 3 points on SSD/SDD
 497     fHistNTracks->Fill(6);
 498     //chisquare/nclusters
 499     Int_t nclu=nSPD+count;
 500     if(track->GetITSchi2()/nclu > fMaxChi2Clu) continue;
 501     fHistNTracks->Fill(7);
 502     //pseudorapidity and rapidity
 503     if(TMath::Abs(track->Eta()) > 0.9) continue;
 504     fHistNTracks->Fill(8);
 505     //truncated mean
 506     //if(fMC) for(Int_t j=0;j<2;j++) s[j]*=3.34/5.43;//correction for SDD miscalibration of the MCpass4
 507     track->GetITSdEdxSamples(s);
 508     Double_t dedx = CookdEdx(s);
 509     if(dedx<0) continue;
 510     fHistNTracks->Fill(9);
 511
 512
 513     Float_t pt = track->Pt();
 514     Int_t theBin=-1;
 515     for(Int_t m=0; m<kNbins; m++){
 516       if(TMath::Abs(pt) > fPtBinLimits[m] && TMath::Abs(pt) < fPtBinLimits[m+1]){
 517         theBin=m;
 518         break;
 519       }
 520     }
 521     track->GetImpactParameters(impactXY, impactZ);
 522
 523     //Filling Ntuple
 524     //information from the MC kinematics
 525     if(fMC){
 526       if(track->GetLabel()<0)isph=-1.;
 527       if(track->GetLabel()>=0){
 528         part = (TParticle*)stack->Particle(track->GetLabel());
 529         pdgPart = part->GetPDG();
 530         code = pdgPart->PdgCode();
 531         if(stack->IsPhysicalPrimary(track->GetLabel()))isph=1.;
 532         else{
 533           isph=0.;
 534           TParticle* moth = stack->Particle(part->GetFirstMother());
 535           Float_t codemoth = TMath::Abs(moth->GetPdgCode());
 536           mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
 537         }
 538       }
 539     }
 540     Float_t xnt[12];
 541     Int_t index=0;
 542     xnt[index++]=(Float_t)track->GetP();
 543     xnt[index++]=(Float_t)track->Pt();
 544     xnt[index++]=(Float_t)dedx;
 545     xnt[index++]=(Float_t)count;
 546     xnt[index++]=(Float_t)track->GetSign();
 547     xnt[index++]=(Float_t)fESD->GetRunNumber();
 548     xnt[index++]=(Float_t)track->Eta();
 549     xnt[index++]=(Float_t)impactXY;
 550     xnt[index++]=(Float_t)impactZ;
 551     xnt[index++]=(Float_t)isph;
 552     xnt[index++]=(Float_t)code;
 553     xnt[index]=(Float_t)mfl;
 554
 555     if(fFillNtuple) fNtupleNSigma->Fill(xnt);
 556
 557
 558
 559     //Compute y and bb
 560     Double_t y[3],bbtheo[3],logdiff[3];
 561     Float_t p=track->GetP();
 562     if(fMC && fSmearMC){
 563       dedx=fRandGener->Gaus(dedx,fSmeardEdx*dedx);
 564       p=fRandGener->Gaus(p,fSmearP*p);
 565     }
 566
 567     for(Int_t i=0;i<3;i++){
 568       y[i] = Eta2y(pt,pdgmass[i],track->Eta());
 569       bbtheo[i]=BetheBloch(p/pdgmass[i],fMC);
 570       logdiff[i]=TMath::Log(dedx) - TMath::Log(bbtheo[i]);
 571     }
 572
 573     //NSigma Method
 574     Int_t resocls=(Int_t)count-1;
 575
 576     Double_t nsigmas[3];
 577     Double_t min=999999.;
 578     Int_t minPos=-1;
 579     for(Int_t isp=0; isp<3; isp++){
 580       Double_t bb=bbtheo[isp];
 581       nsigmas[isp]=TMath::Abs((dedx-bb)/(resodedx[resocls]*bb));
 582       if(nsigmas[isp]<min){
 583         min=nsigmas[isp];
 584         minPos=isp;
 585       }
 586     }
 587     Double_t yPart=y[minPos];
 588
 589     if(min<fMinNSigma && yPart<fMaxY){
 590       //DCA distributions, before the DCA cuts
 591       if(theBin>=0 && theBin<kNbins){
 592         if(track->GetSign()>0){
 593           if(minPos==0) fHistDCAPosPi[theBin]->Fill(impactXY);
 594           else if(minPos==1) fHistDCAPosK[theBin]->Fill(impactXY);
 595           else if(minPos==2) fHistDCAPosP[theBin]->Fill(impactXY);
 596         }else{
 597           if(minPos==0) fHistDCANegPi[theBin]->Fill(impactXY);
 598           else if(minPos==1) fHistDCANegK[theBin]->Fill(impactXY);
 599           else if(minPos==2) fHistDCANegP[theBin]->Fill(impactXY);
 600         }
 601       }
 602     }
 603
 604     //DCA cut on xy and z
 605     if(!DCAcut(impactXY,impactZ,pt,fMC)) continue;
 606     fHistNTracks->Fill(10);
 607
 608     if(min<fMinNSigma && yPart<fMaxY){
 609       if(track->GetSign()>0) fHistPosNSigma[minPos]->Fill(pt);
 610       else fHistNegNSigma[minPos]->Fill(pt);
 611       if(fMC){
 612         if(isph==1){
 613           if(track->GetSign()>0) fHistPosNSigmaPrim[minPos]->Fill(pt);
 614           else fHistNegNSigmaPrim[minPos]->Fill(pt);
 615           if(TMath::Abs(code)==listcode[minPos]){
 616             if(track->GetSign()>0) fHistPosNSigmaPrimMC[minPos]->Fill(pt);
 617             else fHistNegNSigmaPrimMC[minPos]->Fill(pt);
 618           }
 619         }
 620       }
 621     }
 622
 623     if(theBin>=0 && theBin<kNbins){
 624       if(track->GetSign()>0){
 625         if(TMath::Abs(y[0]) < fMaxY)fHistPosPi[theBin]->Fill(logdiff[0]);
 626         if(TMath::Abs(y[1]) < fMaxY)fHistPosK[theBin]->Fill(logdiff[1]);
 627         if(TMath::Abs(y[2]) < fMaxY)fHistPosP[theBin]->Fill(logdiff[2]);
 628         if(fMC){
 629           if((TMath::Abs(y[0])<fMaxY) && (TMath::Abs(code)==211))  fHistMCPosPi[theBin]->Fill(logdiff[0]);
 630           if((TMath::Abs(y[1])<fMaxY) && (TMath::Abs(code)==321))  fHistMCPosK[theBin]->Fill(logdiff[1]);
 631           if((TMath::Abs(y[2])<fMaxY) && (TMath::Abs(code)==2212)) fHistMCPosP[theBin]->Fill(logdiff[2]);
 632         }
 633       }else{
 634         if(TMath::Abs(y[0]) < fMaxY)fHistNegPi[theBin]->Fill(logdiff[0]);
 635         if(TMath::Abs(y[1]) < fMaxY)fHistNegK[theBin]->Fill(logdiff[1]);
 636         if(TMath::Abs(y[2]) < fMaxY)fHistNegP[theBin]->Fill(logdiff[2]);
 637         if(fMC){
 638           if((TMath::Abs(y[0])<fMaxY) && (TMath::Abs(code)==211))  fHistMCNegPi[theBin]->Fill(logdiff[0]);
 639           if((TMath::Abs(y[1])<fMaxY) && (TMath::Abs(code)==321))  fHistMCNegK[theBin]->Fill(logdiff[1]);
 640           if((TMath::Abs(y[2])<fMaxY) && (TMath::Abs(code)==2212)) fHistMCNegP[theBin]->Fill(logdiff[2]);
 641         }
 642       }
 643     }
 644
 645
 646     //fill propaganda plot with dedx
 647     fHistDEDX->Fill(track->GetP(),dedx);
 648     fHistDEDXdouble->Fill(track->GetP()*track->GetSign(),dedx);
 649
 650     //fill charge distribution histo to check the calibration
 651     for(Int_t j=0;j<4;j++){
 652       if(s[j]<5) continue;
 653       fHistCharge[j]->Fill(s[j]);
 654     }
 655   }
 656
 657   // Post output data.
 658   PostData(1,fOutput);
 659   Printf("............. end of Exec");
 660 }
 661
 662 //________________________________________________________________________
 663 void AliAnalysisTaskSEITSsaSpectra::Terminate(Option_t *) {
 664
 665   fOutput = dynamic_cast<TList*>(GetOutputData(1));
 666   if (!fOutput) {
 667     printf("ERROR: fOutput not available\n");
 668     return;
 669   }
 670   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 671   fHistNTracks = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNTracks"));
 672   fHistDEDX = dynamic_cast<TH2F*>(fOutput->FindObject("fHistDEDX"));
 673   fHistDEDXdouble = dynamic_cast<TH2F*>(fOutput->FindObject("fHistDEDXdouble"));
 674
 675   fHistBeforeEvSel = dynamic_cast<TH1F*>(fOutput->FindObject("fHistBeforeEvSel"));
 676   fHistAfterEvSel = dynamic_cast<TH1F*>(fOutput->FindObject("fHistAfterEvSel"));
 677
 678
 679   for(Int_t j=0;j<3;j++){
 680     fHistMCpos[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCpos%d",j)));
 681     fHistMCneg[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCneg%d",j)));
 682   }
 683
 684
 685   for(Int_t j=0;j<3;j++){
 686     fHistMCposBefEvSel[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCposBefEvSel%d",j)));
 687     fHistMCnegBefEvSel[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCnegBefEvSel%d",j)));
 688   }
 689
 690   for(Int_t i=0; i<4; i++){
 691     fHistCharge[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistChargeLay%d",i)));
 692   }
 693
 694   for(Int_t i=0; i<kNbins; i++){
 695     fHistPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPosPi%d",i)));
 696     fHistPosK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPosK%d",i)));
 697     fHistPosP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPosP%d",i)));
 698     fHistNegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistNegPi%d",i)));
 699     fHistNegK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistNegK%d",i)));
 700     fHistNegP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistNegP%d",i)));
 701
 702     fHistDCAPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCAPosPi%d",i)));
 703     fHistDCAPosK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCAPosK%d",i)));
 704     fHistDCAPosP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCAPosP%d",i)));
 705     fHistDCANegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCANegPi%d",i)));
 706     fHistDCANegK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCANegK%d",i)));
 707     fHistDCANegP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCANegP%d",i)));
 708
 709
 710     if(fMC){
 711       fHistMCPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosPi%d",i)));
 712       fHistMCPosK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosK%d",i)));
 713       fHistMCPosP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosP%d",i)));
 714       fHistMCNegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegPi%d",i)));
 715       fHistMCNegK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegK%d",i)));
 716       fHistMCNegP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegP%d",i)));
 717     }
 718   }
 719
 720   for(Int_t j=0;j<3;j++){
 721     fHistPosNSigma[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigma%d",j)));
 722     fHistNegNSigma[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigma%d",j)));
 723     fHistPosNSigmaPrim[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaPrim%d",j)));
 724     fHistNegNSigmaPrim[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaPrim%d",j)));
 725     fHistPosNSigmaPrimMC[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaPrimMC%d",j)));
 726     fHistNegNSigmaPrimMC[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaPrimMC%d",j)));
 727   }
 728
 729   fNtupleNSigma = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleNSigma"));
 730   fNtupleMC = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleMC"));
 731
 732   Printf("end of Terminate");
 733   return;
 734 }