PWG1/ITS/AliAnalysisTaskITSsaTracks.cxx

   1 #include "AliAnalysisTaskSE.h"
   2 #include "AliTrackPointArray.h"
   3 #include "AliAnalysisManager.h"
   4 #include "AliAnalysisDataContainer.h"
   5 #include "AliESDEvent.h"
   6 #include "AliESDfriend.h"
   7 #include "AliStack.h"
   8 #include "AliPID.h"
   9 #include "AliITSPIDResponse.h"
  10 #include "AliMCEventHandler.h"
  11 #include "AliMCEvent.h"
  12 #include <TParticle.h>
  13 #include <TSystem.h>
  14 #include <TTree.h>
  15 #include <TH1F.h>
  16 #include <TH2F.h>
  17 #include <TChain.h>
  18 #include "AliESDInputHandlerRP.h"
  19 #include "AliAnalysisTaskITSsaTracks.h"
  20
  21 /**************************************************************************
  22  * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
  23  *                                                                        *
  24  * Author: The ALICE Off-line Project.                                    *
  25  * Contributors are mentioned in the code where appropriate.              *
  26  *                                                                        *
  27  * Permission to use, copy, modify and distribute this software and its   *
  28  * documentation strictly for non-commercial purposes is hereby granted   *
  29  * without fee, provided that the above copyright notice appears in all   *
  30  * copies and that both the copyright notice and this permission notice   *
  31  * appear in the supporting documentation. The authors make no claims     *
  32  * about the suitability of this software for any purpose. It is          *
  33  * provided "as is" without express or implied warranty.                  *
  34  **************************************************************************/
  35
  36 //*************************************************************************
  37 // Implementation of class AliAnalysisTaskITSsaTracks
  38 // AliAnalysisTaskSE to extract QA and performance histos for ITS standalone tracks
  39 //
  40 //
  41 // Authors: L. Milano, milano@to.infn.it
  42 //          F. Prino, prino@to.infn.it
  43 //
  44 //*************************************************************************
  45
  46 ClassImp(AliAnalysisTaskITSsaTracks)
  47 //______________________________________________________________________________
  48 AliAnalysisTaskITSsaTracks::AliAnalysisTaskITSsaTracks() : AliAnalysisTaskSE("ITSsa resolution"),
  49   fOutput(0),
  50   fHistNEvents(0),
  51   fHistPtTPCITSAll(0),
  52   fHistPtTPCITSGood(0),
  53   fHistPtTPCITSFake(0),
  54   fHistPtITSsaAll(0),
  55   fHistPtITSsaGood(0),
  56   fHistPtITSsaFake(0),
  57   fHistPtITSpureSAAll(0),
  58   fHistPtITSpureSAGood(0),
  59   fHistPtITSpureSAFake(0),
  60   fHistdedxvsPtITSpureSA3cls(0),
  61   fHistdedxvsPITSpureSA3cls(0),
  62   fHistdedxvsPtITSpureSA4cls(0),
  63   fHistdedxvsPITSpureSA4cls(0),
  64   fNPtBins(kMaxPtBins),
  65   fMinITSpts(4),
  66   fMinSPDpts(1),
  67   fMinITSptsForPid(3),
  68   fMinITSptsForMatch(6),
  69   fMinTPCpts(50),
  70   fReadMC(kFALSE),
  71   fUseMCId(kFALSE)
  72 {
  73   //
  74   for(Int_t ilay=0; ilay<6;ilay++) fRequirePoint[ilay]=kFALSE;
  75   DefineInput(0, TChain::Class());
  76   DefineOutput(1, TList::Class());
  77   const Int_t nbins = 29;
  78   Double_t xbins[nbins+1]={0.06,0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,
  79                            0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,
  80                            0.85,0.90,0.95,1.00,1.20,1.40,1.60,1.80,1.90,2.00};
  81   SetPtBins(nbins,xbins);
  82 }
  83
  84
  85 //___________________________________________________________________________
  86 AliAnalysisTaskITSsaTracks::~AliAnalysisTaskITSsaTracks(){
  87   //
  88   if (fOutput) {
  89     delete fOutput;
  90     fOutput = 0;
  91   }
  92 }
  93
  94 //________________________________________________________________________
  95 void AliAnalysisTaskITSsaTracks::SetPtBins(Int_t n, Double_t* lim){
  96   // define pt bins for analysis
  97   if(n>kMaxPtBins){
  98     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
  99     return;
 100   }else{
 101     fNPtBins=n;
 102     for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
 103     for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
 104   }
 105 }
 106 //___________________________________________________________________________
 107 void AliAnalysisTaskITSsaTracks::UserCreateOutputObjects() {
 108   // create output histos
 109
 110   fOutput = new TList();
 111   fOutput->SetOwner();
 112   fOutput->SetName("OutputHistos");
 113
 114   fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
 115   fHistNEvents->Sumw2();
 116   fHistNEvents->SetMinimum(0);
 117   fOutput->Add(fHistNEvents);
 118
 119   //binning for the dedx histogram
 120   const Int_t hnbinsdedx=400;
 121   Double_t hxmindedx = 0.01;
 122   Double_t hxmaxdedx = 10;
 123   Double_t hlogxmindedx = TMath::Log10(hxmindedx);
 124   Double_t hlogxmaxdedx = TMath::Log10(hxmaxdedx);
 125   Double_t hbinwidthdedx = (hlogxmaxdedx-hlogxmindedx)/hnbinsdedx;
 126   Double_t hxbinsdedx[hnbinsdedx+1];
 127   hxbinsdedx[0] = 0.01;
 128   for (Int_t i=1;i<=hnbinsdedx;i++) {
 129     hxbinsdedx[i] = hxmindedx + TMath::Power(10,hlogxmindedx+i*hbinwidthdedx);
 130   }
 131
 132   fHistdedxvsPtITSpureSA3cls = new TH2F("hdedxvsPtITSpureSA3cls","",hnbinsdedx,hxbinsdedx,900,0,1000);
 133   fHistdedxvsPtITSpureSA3cls->Sumw2();
 134   fOutput->Add(fHistdedxvsPtITSpureSA3cls);
 135
 136   fHistdedxvsPITSpureSA3cls = new TH2F("hdedxvsPITSpureSA3cls","",hnbinsdedx,hxbinsdedx,900,0,1000);
 137   fHistdedxvsPITSpureSA3cls->Sumw2();
 138   fOutput->Add(fHistdedxvsPITSpureSA3cls);
 139
 140   fHistdedxvsPtITSpureSA4cls = new TH2F("hdedxvsPtITSpureSA4cls","",hnbinsdedx,hxbinsdedx,900,0,1000);
 141   fHistdedxvsPtITSpureSA4cls->Sumw2();
 142   fOutput->Add(fHistdedxvsPtITSpureSA4cls);
 143
 144   fHistdedxvsPITSpureSA4cls = new TH2F("hdedxvsPITSpureSA4cls","",hnbinsdedx,hxbinsdedx,900,0,1000);
 145   fHistdedxvsPITSpureSA4cls->Sumw2();
 146   fOutput->Add(fHistdedxvsPITSpureSA4cls);
 147
 148   TString spname[3]={"Pion","Kaon","Proton"};
 149   TString hisname;
 150   const Int_t nbins = fNPtBins;
 151   Double_t xbins[nbins+1];
 152   for(Int_t ibin=0; ibin<=nbins; ibin++) xbins[ibin]=fPtLimits[ibin];
 153
 154
 155   fHistPtTPCITSAll = new TH1F("hPtTPCITSAll","",100,0.,2.);
 156   fHistPtTPCITSAll->Sumw2();
 157   fOutput->Add(fHistPtTPCITSAll);
 158   fHistPtTPCITSGood = new TH1F("hPtTPCITSGood","",100,0.,2.);
 159   fHistPtTPCITSGood->Sumw2();
 160   fOutput->Add(fHistPtTPCITSGood);
 161   fHistPtTPCITSFake = new TH1F("hPtTPCITSFake","",100,0.,2.);
 162   fHistPtTPCITSFake->Sumw2();
 163   fOutput->Add(fHistPtTPCITSFake);
 164
 165   fHistPtITSsaAll  = new TH1F("hPtITSsaAll","",100,0.,2.);
 166   fHistPtITSsaAll->Sumw2();
 167   fOutput->Add(fHistPtITSsaAll);
 168   fHistPtITSsaGood  = new TH1F("hPtITSsaGood","",100,0.,2.);
 169   fHistPtITSsaGood->Sumw2();
 170   fOutput->Add(fHistPtITSsaGood);
 171   fHistPtITSsaFake  = new TH1F("hPtITSsaFake","",100,0.,2.);
 172   fHistPtITSsaFake->Sumw2();
 173   fOutput->Add(fHistPtITSsaFake);
 174
 175   fHistPtITSpureSAAll = new TH1F("hPtITSpureSAAll","",100,0.,2.);
 176   fHistPtITSpureSAAll->Sumw2();
 177   fOutput->Add(fHistPtITSpureSAAll);
 178   fHistPtITSpureSAGood = new TH1F("hPtITSpureSAGood","",100,0.,2.);
 179   fHistPtITSpureSAGood->Sumw2();
 180   fOutput->Add(fHistPtITSpureSAGood);
 181   fHistPtITSpureSAFake = new TH1F("hPtITSpureSAFake","",100,0.,2.);
 182   fHistPtITSpureSAFake->Sumw2();
 183   fOutput->Add(fHistPtITSpureSAFake);
 184
 185   for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
 186
 187     hisname.Form("hPtTPCITS%s",spname[iSpec].Data());
 188     fHistPtTPCITS[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
 189     fHistPtTPCITS[iSpec]->Sumw2();
 190     fOutput->Add(fHistPtTPCITS[iSpec]);
 191
 192     hisname.Form("hPtITSsa%s",spname[iSpec].Data());
 193     fHistPtITSsa[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
 194     fHistPtITSsa[iSpec]->Sumw2();
 195     fOutput->Add(fHistPtITSsa[iSpec]);
 196
 197     hisname.Form("hPtITSpureSA%s",spname[iSpec].Data());
 198     fHistPtITSpureSA[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
 199     fHistPtITSpureSA[iSpec]->Sumw2();
 200     fOutput->Add(fHistPtITSpureSA[iSpec]);
 201
 202     //---
 203
 204     hisname.Form("hEtaPhiTPCITS%s",spname[iSpec].Data());
 205     fHistEtaPhiTPCITS[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
 206     fHistEtaPhiTPCITS[iSpec]->Sumw2();
 207     fOutput->Add(fHistEtaPhiTPCITS[iSpec]);
 208
 209     hisname.Form("hEtaPhiITSsa%s",spname[iSpec].Data());
 210     fHistEtaPhiITSsa[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
 211     fHistEtaPhiITSsa[iSpec]->Sumw2();
 212     fOutput->Add(fHistEtaPhiITSsa[iSpec]);
 213
 214     hisname.Form("hEtaPhiITSpureSA%s",spname[iSpec].Data());
 215     fHistEtaPhiITSpureSA[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
 216     fHistEtaPhiITSpureSA[iSpec]->Sumw2();
 217     fOutput->Add(fHistEtaPhiITSpureSA[iSpec]);
 218
 219     //---
 220
 221     hisname.Form("hNcluTPCITS%s",spname[iSpec].Data());
 222     fHistNcluTPCITS[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
 223     fHistNcluTPCITS[iSpec]->Sumw2();
 224     fOutput->Add(fHistNcluTPCITS[iSpec]);
 225
 226     hisname.Form("hNcluITSsa%s",spname[iSpec].Data());
 227     fHistNcluITSsa[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
 228     fHistNcluITSsa[iSpec]->Sumw2();
 229     fOutput->Add(fHistNcluITSsa[iSpec]);
 230
 231     hisname.Form("hNcluITSpureSA%s",spname[iSpec].Data());
 232     fHistNcluITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
 233     fHistNcluITSpureSA[iSpec]->Sumw2();
 234     fOutput->Add(fHistNcluITSpureSA[iSpec]);
 235
 236     //---
 237
 238     hisname.Form("hd0rphiITSpureSA%s",spname[iSpec].Data());
 239     fHistd0rphiITSpureSA[iSpec] = new TH2F(hisname.Data(),"",nbins,xbins,2000,-1,1);
 240     fHistd0rphiITSpureSA[iSpec]->Sumw2();
 241     fOutput->Add(fHistd0rphiITSpureSA[iSpec]);
 242
 243     hisname.Form("hd0zITSpureSA%s",spname[iSpec].Data());
 244     fHistd0zITSpureSA[iSpec] = new TH2F(hisname.Data(),"",nbins,xbins,2000,-1,1);
 245     fHistd0zITSpureSA[iSpec]->Sumw2();
 246     fOutput->Add(fHistd0zITSpureSA[iSpec]);
 247
 248     //---
 249
 250     hisname.Form("hCluInLayTPCITS%s",spname[iSpec].Data());
 251     fHistCluInLayTPCITS[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
 252     fHistCluInLayTPCITS[iSpec]->Sumw2();
 253     fOutput->Add(fHistCluInLayTPCITS[iSpec]);
 254
 255     hisname.Form("hCluInLayITSsa%s",spname[iSpec].Data());
 256     fHistCluInLayITSsa[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
 257     fHistCluInLayITSsa[iSpec]->Sumw2();
 258     fOutput->Add(fHistCluInLayITSsa[iSpec]);
 259
 260     hisname.Form("hCluInLayITSpureSA%s",spname[iSpec].Data());
 261     fHistCluInLayITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
 262     fHistCluInLayITSpureSA[iSpec]->Sumw2();
 263     fOutput->Add(fHistCluInLayITSpureSA[iSpec]);
 264
 265     hisname.Form("hOuterLayITSpureSA%s",spname[iSpec].Data());
 266     fHistOuterLayITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
 267     fHistOuterLayITSpureSA[iSpec]->Sumw2();
 268     fOutput->Add(fHistOuterLayITSpureSA[iSpec]);
 269
 270     //---
 271
 272     hisname.Form("hPtResid%s",spname[iSpec].Data());
 273     fHistPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
 274     fHistPtResid[iSpec]->Sumw2();
 275     fOutput->Add(fHistPtResid[iSpec]);
 276     hisname.Form("hPtRelResid%s",spname[iSpec].Data());
 277     fHistPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
 278     fHistPtRelResid[iSpec]->Sumw2();
 279     fOutput->Add(fHistPtRelResid[iSpec]);
 280
 281     hisname.Form("hInvPtResid%s",spname[iSpec].Data());
 282     fHistInvPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
 283     fHistInvPtResid[iSpec]->Sumw2();
 284     fOutput->Add(fHistInvPtResid[iSpec]);
 285     hisname.Form("hInvPtRelResid%s",spname[iSpec].Data());
 286     fHistInvPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
 287     fHistInvPtRelResid[iSpec]->Sumw2();
 288     fOutput->Add(fHistInvPtRelResid[iSpec]);
 289
 290     hisname.Form("hMCPtResid%s",spname[iSpec].Data());
 291     fHistMCPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
 292     fHistMCPtResid[iSpec]->Sumw2();
 293     fOutput->Add(fHistMCPtResid[iSpec]);
 294     hisname.Form("hMCPtRelResid%s",spname[iSpec].Data());
 295     fHistMCPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
 296     fHistMCPtRelResid[iSpec]->Sumw2();
 297     fOutput->Add(fHistMCPtRelResid[iSpec]);
 298
 299     hisname.Form("hMCInvPtResid%s",spname[iSpec].Data());
 300     fHistMCInvPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
 301     fHistMCInvPtResid[iSpec]->Sumw2();
 302     fOutput->Add(fHistMCInvPtResid[iSpec]);
 303     hisname.Form("hMCInvPtRelResid%s",spname[iSpec].Data());
 304     fHistMCInvPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
 305     fHistMCInvPtRelResid[iSpec]->Sumw2();
 306     fOutput->Add(fHistMCInvPtRelResid[iSpec]);
 307
 308   }
 309
 310 }
 311 //______________________________________________________________________________
 312 void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
 313 {
 314   //
 315
 316   AliESDEvent *esd = (AliESDEvent*) (InputEvent());
 317
 318
 319   if(!esd) {
 320     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
 321     return;
 322   }
 323
 324
 325   if(!ESDfriend()) {
 326     printf("AliAnalysisTaskSDDRP::Exec(): bad ESDfriend\n");
 327     return;
 328   }
 329
 330   AliStack* stack=0;
 331
 332   if(fReadMC){
 333     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 334     if (!eventHandler) {
 335       Printf("ERROR: Could not retrieve MC event handler");
 336       return;
 337     }
 338     AliMCEvent* mcEvent = eventHandler->MCEvent();
 339     if (!mcEvent) {
 340       Printf("ERROR: Could not retrieve MC event");
 341       return;
 342     }
 343     stack = mcEvent->Stack();
 344     if (!stack) {
 345       Printf("ERROR: stack not available");
 346       return;
 347     }
 348   }
 349
 350   PostData(1, fOutput);
 351
 352   fHistNEvents->Fill(0);
 353   const AliESDVertex *spdv=esd->GetPrimaryVertexSPD();
 354   if(spdv->GetNContributors()<=0) return;
 355   fHistNEvents->Fill(1);
 356
 357   Int_t ntracks = esd->GetNumberOfTracks();
 358   for (Int_t iTrack=0; iTrack < ntracks; iTrack++) {
 359     AliESDtrack * track = esd->GetTrack(iTrack);
 360     if (!track) continue;
 361     Int_t status=track->GetStatus();
 362     if(!(status&AliESDtrack::kITSrefit)) continue;
 363     Bool_t isSA=kTRUE;
 364     if(status&AliESDtrack::kTPCin) isSA=kFALSE;
 365     Int_t nTPCclus=track->GetNcls(1);
 366     Int_t nITSclus=track->GetNcls(0);
 367     if(nITSclus<fMinITSpts)continue;
 368
 369     UChar_t clumap=track->GetITSClusterMap();
 370     Int_t nSPDPoints=0;
 371     for(Int_t i=0; i<2; i++){
 372       if(clumap&(1<<i)) ++nSPDPoints;
 373     }
 374     if(nSPDPoints<fMinSPDpts) continue;
 375
 376     Int_t nPointsForPid=0;
 377     for(Int_t i=2; i<6; i++){
 378       if(clumap&(1<<i)) ++nPointsForPid;
 379     }
 380
 381     Float_t pttrack=track->Pt();
 382     Float_t ptrack=track->P();
 383     Float_t dedx=track->GetITSsignal();
 384     Int_t hadronSpecie=-1;
 385     if(status&AliESDtrack::kTPCin){
 386       fHistPtTPCITSAll->Fill(pttrack);
 387     }else{
 388       if(status&AliESDtrack::kITSpureSA){
 389         fHistPtITSpureSAAll->Fill(pttrack);
 390       }else{
 391         fHistPtITSsaAll->Fill(pttrack);
 392       }
 393     }
 394
 395     if(fReadMC && fUseMCId){
 396       Int_t trlabel=track->GetLabel();
 397       if(trlabel>=0){
 398         if(status&AliESDtrack::kTPCin){
 399           fHistPtTPCITSGood->Fill(pttrack);
 400         }else{
 401           if(status&AliESDtrack::kITSpureSA){
 402             fHistPtITSpureSAGood->Fill(pttrack);
 403           }else{
 404             fHistPtITSsaGood->Fill(pttrack);
 405           }
 406         }
 407       }else{
 408         if(status&AliESDtrack::kTPCin){
 409           fHistPtTPCITSFake->Fill(pttrack);
 410         }else{
 411           if(status&AliESDtrack::kITSpureSA){
 412             fHistPtITSpureSAFake->Fill(pttrack);
 413           }else{
 414             fHistPtITSsaFake->Fill(pttrack);
 415           }
 416         }
 417         continue; // fake track
 418       }
 419       TParticle* part = stack->Particle(trlabel);
 420       Int_t pdg=TMath::Abs(part->GetPdgCode());
 421       if(pdg==211) hadronSpecie=kPion;
 422       else if(pdg==321) hadronSpecie=kKaon;
 423       else if(pdg==2212) hadronSpecie=kProton;
 424     }else{
 425       if(nPointsForPid<fMinITSptsForPid)continue;
 426       AliITSPIDResponse pidits(fReadMC);
 427       Float_t nSigmaPion=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kPion,nPointsForPid,isSA);
 428       if(nSigmaPion>-2. && nSigmaPion<2.){
 429         hadronSpecie=kPion;
 430       }else{
 431         Float_t nSigmaKaon=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kKaon,nPointsForPid,isSA);
 432         if(nSigmaKaon>-2. && nSigmaKaon<2.){
 433           hadronSpecie=kKaon;
 434         }else{
 435           Float_t nSigmaProton=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kProton,nPointsForPid,isSA);
 436           if(nSigmaProton>-2. && nSigmaProton<2.){
 437             hadronSpecie=kProton;
 438           }
 439         }
 440       }
 441     }
 442     if(hadronSpecie<0) continue;
 443     if(status&AliESDtrack::kTPCin){ // TPC+ITS tracks
 444       fHistPtTPCITS[hadronSpecie]->Fill(pttrack);
 445       fHistEtaPhiTPCITS[hadronSpecie]->Fill(track->Eta(),track->Phi());
 446       fHistNcluTPCITS[hadronSpecie]->Fill(pttrack,nITSclus);
 447       fHistCluInLayTPCITS[hadronSpecie]->Fill(-1.);
 448       for(Int_t iBit=0; iBit<6; iBit++){
 449         if(clumap&(1<<iBit)) fHistCluInLayTPCITS[hadronSpecie]->Fill(pttrack,iBit);
 450       }
 451     }else{ // ITS standalone and pureSA tracks
 452       if(status&AliESDtrack::kITSpureSA){
 453         Float_t impactXY=-999, impactZ=-999;
 454         track->GetImpactParameters(impactXY, impactZ);
 455         fHistPtITSpureSA[hadronSpecie]->Fill(pttrack);
 456         fHistEtaPhiITSpureSA[hadronSpecie]->Fill(track->Eta(),track->Phi());
 457         fHistNcluITSpureSA[hadronSpecie]->Fill(pttrack,nITSclus);
 458         fHistd0rphiITSpureSA[hadronSpecie]->Fill(pttrack,impactXY);
 459         fHistd0zITSpureSA[hadronSpecie]->Fill(pttrack,impactZ);
 460         if(nPointsForPid==3){
 461           fHistdedxvsPtITSpureSA3cls->Fill(pttrack,dedx);
 462           fHistdedxvsPITSpureSA3cls->Fill(ptrack,dedx);
 463         }
 464         if(nPointsForPid==4){
 465           fHistdedxvsPtITSpureSA4cls->Fill(pttrack,dedx);
 466           fHistdedxvsPITSpureSA4cls->Fill(ptrack,dedx);
 467         }
 468         fHistCluInLayITSpureSA[hadronSpecie]->Fill(-1.);
 469         Int_t outerLay=-1;
 470         for(Int_t iBit=0; iBit<6; iBit++){
 471           if(clumap&(1<<iBit)){
 472             fHistCluInLayITSpureSA[hadronSpecie]->Fill(pttrack,iBit);
 473             if(iBit>outerLay) outerLay=iBit;
 474           }
 475         }
 476         fHistOuterLayITSpureSA[hadronSpecie]->Fill(pttrack,outerLay);
 477
 478
 479         if(fReadMC){
 480           Int_t trlabel=track->GetLabel();
 481           if(trlabel<0)continue; // fake track
 482           TParticle* part = stack->Particle(trlabel);
 483           Float_t ptgen=part->Pt();
 484           Float_t invpttrack=track->OneOverPt();
 485           Float_t invptgen=0.;
 486           if(ptgen>0.) invptgen=1./ptgen;
 487           fHistMCPtResid[hadronSpecie]->Fill(pttrack,pttrack-ptgen);
 488           fHistMCPtRelResid[hadronSpecie]->Fill(pttrack,(pttrack-ptgen)/ptgen);
 489           fHistMCInvPtResid[hadronSpecie]->Fill(pttrack,invpttrack-invptgen);
 490           fHistMCInvPtRelResid[hadronSpecie]->Fill(pttrack,(invpttrack-invptgen)/invptgen);
 491         }
 492       }else{
 493         fHistPtITSsa[hadronSpecie]->Fill(pttrack);
 494         fHistEtaPhiITSsa[hadronSpecie]->Fill(track->Eta(),track->Phi());
 495         fHistNcluITSsa[hadronSpecie]->Fill(pttrack,nITSclus);
 496         fHistCluInLayITSsa[hadronSpecie]->Fill(-1.);
 497         for(Int_t iBit=0; iBit<6; iBit++){
 498           if(clumap&(1<<iBit)) fHistCluInLayITSsa[hadronSpecie]->Fill(pttrack,iBit);
 499         }
 500       }
 501     }
 502
 503     if(nITSclus<fMinITSptsForMatch || nTPCclus<fMinTPCpts) continue;
 504     Int_t idxMI[12],idxSA[12];
 505     for(Int_t icl=0; icl<12; icl++){
 506       idxMI[icl]=-1;
 507       idxSA[icl]=-1;
 508     }
 509     Int_t ncls=track->GetClusters(0,idxMI);
 510     if(fMinITSpts<6){
 511       Bool_t accept=kTRUE;
 512       for(Int_t ilay=0; ilay<6; ilay++){
 513         if(fRequirePoint[ilay]){
 514           Int_t mask = 1<<ilay;
 515           if(!(clumap & mask)){
 516             accept=kFALSE;
 517             break;
 518           }
 519         }
 520       }
 521       if(!accept) continue;
 522     }
 523     // Sort
 524     for(Int_t i=0;i<12;i++){
 525       for(Int_t j=i+1;j<12;j++){
 526         if(idxMI[j]>idxMI[i]){
 527           Int_t tmp=idxMI[j];
 528           idxMI[j]=idxMI[i];
 529           idxMI[i]=tmp;
 530         }
 531       }
 532     }
 533     //    for(Int_t i=0; i<12; i++) printf("%d ",idxMI[i]);
 534     //    printf("\n");
 535     if(idxMI[0]<0 && idxMI[0]==idxMI[1]) continue;
 536     Bool_t matched=kFALSE;
 537     for (Int_t iTrack2 = 0; iTrack2 < ntracks; iTrack2++) {
 538       if(matched) break;
 539       if(iTrack2==iTrack) continue;
 540       AliESDtrack* track2 = esd->GetTrack(iTrack2);
 541       Int_t status2=track2->GetStatus();
 542       if(!(status2&AliESDtrack::kITSrefit)) continue;
 543       if(!(status2&AliESDtrack::kITSpureSA)) continue;
 544       Int_t clumap2=track2->GetITSClusterMap();
 545       Int_t nITSclus2=track2->GetNcls(0);
 546       Int_t nTPCclus2=track2->GetNcls(1);
 547       if(nITSclus2<fMinITSpts || nTPCclus2>0) continue;
 548       Int_t ncls2=track2->GetClusters(0,idxSA);
 549       if(ncls2!=ncls) continue;
 550       if(fMinITSpts<6){
 551         Bool_t accept=kTRUE;
 552         for(Int_t ilay=0; ilay<6; ilay++){
 553           if(fRequirePoint[ilay]){
 554             Int_t mask = 1<<ilay;
 555             if(!(clumap2 & mask)){
 556               accept=kFALSE;
 557               break;
 558             }
 559           }
 560         }
 561         if(!accept) continue;
 562       }
 563       // Sort
 564       for(Int_t i=0;i<12;i++){
 565         for(Int_t j=i+1;j<12;j++){
 566           if(idxSA[j]>idxSA[i]){
 567             Int_t tmp=idxSA[j];
 568             idxSA[j]=idxSA[i];
 569             idxSA[i]=tmp;
 570           }
 571         }
 572       }
 573       Int_t match=0;
 574       for(Int_t icl=0; icl<ncls; icl++){
 575         if(idxSA[icl]!=idxMI[icl]){
 576           match=0;
 577           break;
 578         }
 579         else match++;
 580       }
 581       if(match==ncls && match>0){
 582         matched=kTRUE;
 583         Float_t pt1=track->Pt();
 584         Float_t pt2=track2->Pt();
 585         Float_t ptm1=track->OneOverPt();
 586         Float_t ptm2=track2->OneOverPt();
 587         fHistPtResid[hadronSpecie]->Fill(pt1,pt2-pt1);
 588         fHistPtRelResid[hadronSpecie]->Fill(pt1,(pt2-pt1)/pt1);
 589         fHistInvPtResid[hadronSpecie]->Fill(pt1,ptm2-ptm1);
 590         fHistInvPtRelResid[hadronSpecie]->Fill(pt1,(ptm2-ptm1)/ptm1);
 591       }
 592     }
 593   }
 594
 595   PostData(1,fOutput);
 596
 597 }
 598 //______________________________________________________________________________
 599 void AliAnalysisTaskITSsaTracks::Terminate(Option_t */*option*/)
 600 {
 601   // Terminate analysis
 602   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 603   if (!fOutput) {
 604     printf("ERROR: fOutput not available\n");
 605     return;
 606   }
 607   fHistNEvents= dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 608   fHistdedxvsPtITSpureSA3cls= dynamic_cast<TH2F*>(fOutput->FindObject("hdedxvsPtITSpureSA3cls%s"));
 609   fHistdedxvsPITSpureSA3cls= dynamic_cast<TH2F*>(fOutput->FindObject("hdedxvsPITSpureSA3cls%s"));
 610   fHistdedxvsPtITSpureSA4cls= dynamic_cast<TH2F*>(fOutput->FindObject("hdedxvsPtITSpureSA4cls%s"));
 611   fHistdedxvsPITSpureSA4cls= dynamic_cast<TH2F*>(fOutput->FindObject("hdedxvsPITSpureSA4cls%s"));
 612
 613   TString spname[3]={"Pion","Kaon","Proton"};
 614
 615   fHistPtTPCITSAll =dynamic_cast<TH1F*>(fOutput->FindObject("hPtTPCITSAll"));
 616   fHistPtTPCITSGood =dynamic_cast<TH1F*>(fOutput->FindObject("hPtTPCITSGood"));
 617   fHistPtTPCITSFake =dynamic_cast<TH1F*>(fOutput->FindObject("hPtTPCITSFake"));
 618
 619   fHistPtITSsaAll =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSsaAll"));
 620   fHistPtITSsaGood =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSsaGood"));
 621   fHistPtITSsaFake =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSsaFake"));
 622
 623   fHistPtITSpureSAAll =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSpureSAAll"));
 624   fHistPtITSpureSAGood =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSpureSAGood"));
 625   fHistPtITSpureSAFake =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSpureSAFake"));
 626
 627   for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
 628     fHistPtTPCITS[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtTPCITS%s",spname[iSpec].Data())));
 629     fHistPtITSsa[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtITSsa%s",spname[iSpec].Data())));
 630     fHistPtITSpureSA[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtITSpureSA%s",spname[iSpec].Data())));
 631
 632     fHistEtaPhiTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiTPCITS%s",spname[iSpec].Data())));
 633     fHistEtaPhiITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiITSsa%s",spname[iSpec].Data())));
 634     fHistEtaPhiITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiITSpureSA%s",spname[iSpec].Data())));
 635
 636     fHistNcluTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluTPCITS%s",spname[iSpec].Data())));
 637     fHistNcluITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluITSsa%s",spname[iSpec].Data())));
 638     fHistNcluITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluITSpureSA%s",spname[iSpec].Data())));
 639
 640     fHistd0rphiITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hd0rphiITSpureSA%s",spname[iSpec].Data())));
 641     fHistd0zITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hd0zITSpureSA%s",spname[iSpec].Data())));
 642
 643     fHistCluInLayTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayTPCITS%s",spname[iSpec].Data())));
 644     fHistCluInLayITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayITSsa%s",spname[iSpec].Data())));
 645     fHistCluInLayITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayITSpureSA%s",spname[iSpec].Data())));
 646
 647     fHistOuterLayITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hOuterLayITSpureSA%s",spname[iSpec].Data())));
 648
 649     fHistPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hPtResid%s",spname[iSpec].Data())));
 650     fHistPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hPtRelResid%s",spname[iSpec].Data())));
 651     fHistInvPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hInvPtResid%s",spname[iSpec].Data())));
 652     fHistInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hInvPtRelResid%s",spname[iSpec].Data())));
 653     fHistMCPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCPtResid%s",spname[iSpec].Data())));
 654     fHistMCPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCPtRelResid%s",spname[iSpec].Data())));
 655     fHistMCInvPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtResid%s",spname[iSpec].Data())));
 656     fHistMCInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtRelResid%s",spname[iSpec].Data())));
 657
 658   }
 659   return;
 660 }
 661
 662
 663
 664
 665
 666