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