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