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 <TNtuple.h>
  16 #include <TH1F.h>
  17 #include <TH2F.h>
  18 #include <TChain.h>
  19 #include "AliESDInputHandlerRP.h"
  20 #include "AliAnalysisTaskITSsaTracks.h"
  21
  22 /**************************************************************************
  23  * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
  24  *                                                                        *
  25  * Author: The ALICE Off-line Project.                                    *
  26  * Contributors are mentioned in the code where appropriate.              *
  27  *                                                                        *
  28  * Permission to use, copy, modify and distribute this software and its   *
  29  * documentation strictly for non-commercial purposes is hereby granted   *
  30  * without fee, provided that the above copyright notice appears in all   *
  31  * copies and that both the copyright notice and this permission notice   *
  32  * appear in the supporting documentation. The authors make no claims     *
  33  * about the suitability of this software for any purpose. It is          *
  34  * provided "as is" without express or implied warranty.                  *
  35  **************************************************************************/
  36
  37 //*************************************************************************
  38 // Implementation of class AliAnalysisTaskITSsaTracks
  39 // AliAnalysisTaskSE to extract QA and performance histos for ITS standalone tracks
  40 //
  41 //
  42 // Authors: L. Milano, milano@to.infn.it
  43 //          F. Prino, prino@to.infn.it
  44 //
  45 //*************************************************************************
  46
  47 ClassImp(AliAnalysisTaskITSsaTracks)
  48 //______________________________________________________________________________
  49 AliAnalysisTaskITSsaTracks::AliAnalysisTaskITSsaTracks() : AliAnalysisTaskSE("ITSsa resolution"),
  50   fOutput(0),
  51   fHistNEvents(0),
  52   fHistMCPhiResid(0),
  53   fHistPhiResid(0),
  54   fNtupleTracks(0),
  55   fNPtBins(kMaxPtBins),
  56   fMinITSpts(4),
  57   fMinSPDpts(1),
  58   fMinPtsforPid(3),
  59   fMinTPCpts(50),
  60   fMaxITSChi2Clu(100.),
  61   fFillNtuple(kFALSE),
  62   fReadMC(kFALSE),
  63   fUseMCId(kFALSE)
  64 {
  65   //
  66   for(Int_t ilay=0; ilay<6;ilay++) fRequirePoint[ilay]=kFALSE;
  67   DefineInput(0, TChain::Class());
  68   DefineOutput(1, TList::Class());
  69   const Int_t nbins = 29;
  70   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,
  71                            0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,
  72                            0.85,0.90,0.95,1.00,1.20,1.40,1.60,1.80,1.90,2.00};
  73   SetPtBins(nbins,xbins);
  74 }
  75
  76
  77 //___________________________________________________________________________
  78 AliAnalysisTaskITSsaTracks::~AliAnalysisTaskITSsaTracks(){
  79   //
  80   delete fHistNEvents;
  81   for(Int_t iType=0; iType<kNtrackTypes; iType++){
  82     delete fHistPt[iType];
  83     delete fHistPtGood[iType];
  84     delete fHistPtFake[iType];
  85     delete fHistEtaPhi[iType];
  86     delete fHistEtaPhiGood[iType];
  87     delete fHistEtaPhiFake[iType];
  88     delete fHistChi2[iType];
  89     delete fHistChi2Good[iType];
  90     delete fHistChi2Fake[iType];
  91     delete fHistNclu[iType];
  92     delete fHistNcluGood[iType];
  93     delete fHistNcluFake[iType];
  94     delete fHistdedxvsPt3cls[iType];
  95     delete fHistdedxvsPt4cls[iType];
  96   }
  97   for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
  98     delete fHistPtTPCITS[iSpec];
  99     delete fHistPtITSsa[iSpec];
 100     delete fHistPtITSpureSA[iSpec];
 101     delete fHistEtaPhiTPCITS[iSpec];
 102     delete fHistEtaPhiITSsa[iSpec];
 103     delete fHistEtaPhiITSpureSA[iSpec];
 104     delete fHistNcluTPCITS[iSpec];
 105     delete fHistNcluITSsa[iSpec];
 106     delete fHistNcluITSpureSA[iSpec];
 107     delete fHistd0rphiITSpureSA[iSpec];
 108     delete fHistd0zITSpureSA[iSpec];
 109     delete fHistCluInLayTPCITS[iSpec];
 110     delete fHistCluInLayITSsa[iSpec];
 111     delete fHistCluInLayITSpureSA[iSpec];
 112     delete fHistOuterLayITSpureSA[iSpec];
 113     delete fHistPtResid[iSpec];
 114     delete fHistPtRelResid[iSpec];
 115     delete fHistInvPtResid[iSpec];
 116     delete fHistInvPtRelResid[iSpec];
 117     delete fHistMCPtResid[iSpec];
 118     delete fHistMCPtRelResid[iSpec];
 119     delete fHistMCInvPtResid[iSpec];
 120     delete fHistMCInvPtRelResid[iSpec];
 121   }
 122
 123   delete fHistMCPhiResid;
 124   delete fHistPhiResid;
 125   delete fNtupleTracks;
 126   if (fOutput) {
 127     delete fOutput;
 128     fOutput = 0;
 129   }
 130 }
 131
 132 //________________________________________________________________________
 133 void AliAnalysisTaskITSsaTracks::SetPtBins(Int_t n, Double_t* lim){
 134   // define pt bins for analysis
 135   if(n>kMaxPtBins){
 136     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
 137     return;
 138   }else{
 139     fNPtBins=n;
 140     for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
 141     for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
 142   }
 143 }
 144 //___________________________________________________________________________
 145 void AliAnalysisTaskITSsaTracks::UserCreateOutputObjects() {
 146   // create output histos
 147
 148   fOutput = new TList();
 149   fOutput->SetOwner();
 150   fOutput->SetName("OutputHistos");
 151
 152   fNtupleTracks= new TNtuple("ntupleTracks","ntupleTracks","pt:px:py:pz:eta:phi:d0xy:d0z:dedx3:dedx4:dedx5:dedx6:truncmean:chi2:status:ITSrefit:TPCin:TPCrefit:ITSpureSA:nITSclu:nTPCclu:clumap:spdin:spdout:sddin:sddout:ssdin:ssdout:label:ptgen:pxgen:pygen:pzgen:etagen:phigen:ntracks");
 153   // kinematics: pt, p eta,phi        ->  6 variables
 154   // impact parameters: d0xy, d0z     ->  2 variables
 155   // dE/dx: 4 Layers + trunc mean     ->  5 variables
 156   // chi2 and track status:           ->  6 variables
 157   // cluster infos:                   ->  9 variables
 158   // MC info:                         ->  7 variables
 159   // Multiplicity                     ->  1 variable
 160   // Total:                              36
 161   fOutput->Add(fNtupleTracks);
 162
 163   fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
 164   fHistNEvents->Sumw2();
 165   fHistNEvents->SetMinimum(0);
 166   fOutput->Add(fHistNEvents);
 167
 168   //binning for the dedx histogram
 169   const Int_t hnbinsdedx=400;
 170   Double_t hxmindedx = 0.01;
 171   Double_t hxmaxdedx = 10;
 172   Double_t hlogxmindedx = TMath::Log10(hxmindedx);
 173   Double_t hlogxmaxdedx = TMath::Log10(hxmaxdedx);
 174   Double_t hbinwidthdedx = (hlogxmaxdedx-hlogxmindedx)/hnbinsdedx;
 175   Double_t hxbinsdedx[hnbinsdedx+1];
 176   hxbinsdedx[0] = 0.01;
 177   for (Int_t i=1;i<=hnbinsdedx;i++) {
 178     hxbinsdedx[i] = hxmindedx + TMath::Power(10,hlogxmindedx+i*hbinwidthdedx);
 179   }
 180
 181   TString tit[kNtrackTypes]={"TPCITS","ITSsa","ITSpureSA"};
 182   for(Int_t iType=0; iType<kNtrackTypes; iType++){
 183     fHistPt[iType]=new TH1F(Form("hPt%s",tit[iType].Data()),"",100,0.,2.);
 184     fHistPt[iType]->Sumw2();
 185     fOutput->Add(fHistPt[iType]);
 186     fHistPtGood[iType]=new TH1F(Form("hPtGood%s",tit[iType].Data()),"",100,0.,2.);
 187     fHistPtGood[iType]->Sumw2();
 188     fOutput->Add(fHistPtGood[iType]);
 189     fHistPtFake[iType]=new TH1F(Form("hPtFake%s",tit[iType].Data()),"",100,0.,2.);
 190     fHistPtFake[iType]->Sumw2();
 191     fOutput->Add(fHistPtFake[iType]);
 192
 193     fHistEtaPhi[iType] = new TH2F(Form("hEtaPhi%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
 194     fHistEtaPhi[iType]->Sumw2();
 195     fOutput->Add(fHistEtaPhi[iType]);
 196     fHistEtaPhiGood[iType] = new TH2F(Form("hEtaPhiGood%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
 197     fHistEtaPhiGood[iType]->Sumw2();
 198     fOutput->Add(fHistEtaPhiGood[iType]);
 199     fHistEtaPhiFake[iType] = new TH2F(Form("hEtaPhiFake%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
 200     fHistEtaPhiFake[iType]->Sumw2();
 201     fOutput->Add(fHistEtaPhiFake[iType]);
 202
 203     fHistChi2[iType]=new TH1F(Form("hChi2%s",tit[iType].Data()),"",100,0.,10.);
 204     fHistChi2[iType]->Sumw2();
 205     fOutput->Add(fHistChi2[iType]);
 206     fHistChi2Good[iType]=new TH1F(Form("hChi2Good%s",tit[iType].Data()),"",100,0.,10.);
 207     fHistChi2Good[iType]->Sumw2();
 208     fOutput->Add(fHistChi2Good[iType]);
 209     fHistChi2Fake[iType]=new TH1F(Form("hChi2Fake%s",tit[iType].Data()),"",100,0.,10.);
 210     fHistChi2Fake[iType]->Sumw2();
 211     fOutput->Add(fHistChi2Fake[iType]);
 212
 213     fHistNclu[iType]=new TH1F(Form("hNclu%s",tit[iType].Data()),"",7,-0.5,6.5);
 214     fHistNclu[iType]->Sumw2();
 215     fOutput->Add(fHistNclu[iType]);
 216     fHistNcluGood[iType]=new TH1F(Form("hNcluGood%s",tit[iType].Data()),"",7,-0.5,6.5);
 217     fHistNcluGood[iType]->Sumw2();
 218     fOutput->Add(fHistNcluGood[iType]);
 219     fHistNcluFake[iType]=new TH1F(Form("hNcluFake%s",tit[iType].Data()),"",7,-0.5,6.5);
 220     fHistNcluFake[iType]->Sumw2();
 221     fOutput->Add(fHistNcluFake[iType]);
 222
 223     fHistdedxvsPt3cls[iType] = new TH2F(Form("hdedxvsPt3cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
 224     fHistdedxvsPt3cls[iType]->Sumw2();
 225     fOutput->Add(fHistdedxvsPt3cls[iType]);
 226
 227     fHistdedxvsPt4cls[iType] = new TH2F(Form("hdedxvsPt4cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
 228     fHistdedxvsPt4cls[iType]->Sumw2();
 229     fOutput->Add(fHistdedxvsPt4cls[iType]);
 230   }
 231
 232   //-----------------------------------------------------------
 233
 234   TString spname[3]={"Pion","Kaon","Proton"};
 235   TString hisname;
 236   const Int_t nbins = fNPtBins;
 237   Double_t xbins[nbins+1];
 238   for(Int_t ibin=0; ibin<=nbins; ibin++) xbins[ibin]=fPtLimits[ibin];
 239
 240   for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
 241
 242     hisname.Form("hPtTPCITS%s",spname[iSpec].Data());
 243     fHistPtTPCITS[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
 244     fHistPtTPCITS[iSpec]->Sumw2();
 245     fOutput->Add(fHistPtTPCITS[iSpec]);
 246
 247     hisname.Form("hPtITSsa%s",spname[iSpec].Data());
 248     fHistPtITSsa[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
 249     fHistPtITSsa[iSpec]->Sumw2();
 250     fOutput->Add(fHistPtITSsa[iSpec]);
 251
 252     hisname.Form("hPtITSpureSA%s",spname[iSpec].Data());
 253     fHistPtITSpureSA[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
 254     fHistPtITSpureSA[iSpec]->Sumw2();
 255     fOutput->Add(fHistPtITSpureSA[iSpec]);
 256
 257     //---
 258
 259     hisname.Form("hEtaPhiTPCITS%s",spname[iSpec].Data());
 260     fHistEtaPhiTPCITS[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
 261     fHistEtaPhiTPCITS[iSpec]->Sumw2();
 262     fOutput->Add(fHistEtaPhiTPCITS[iSpec]);
 263
 264     hisname.Form("hEtaPhiITSsa%s",spname[iSpec].Data());
 265     fHistEtaPhiITSsa[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
 266     fHistEtaPhiITSsa[iSpec]->Sumw2();
 267     fOutput->Add(fHistEtaPhiITSsa[iSpec]);
 268
 269     hisname.Form("hEtaPhiITSpureSA%s",spname[iSpec].Data());
 270     fHistEtaPhiITSpureSA[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
 271     fHistEtaPhiITSpureSA[iSpec]->Sumw2();
 272     fOutput->Add(fHistEtaPhiITSpureSA[iSpec]);
 273
 274     //---
 275
 276     hisname.Form("hNcluTPCITS%s",spname[iSpec].Data());
 277     fHistNcluTPCITS[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
 278     fHistNcluTPCITS[iSpec]->Sumw2();
 279     fOutput->Add(fHistNcluTPCITS[iSpec]);
 280
 281     hisname.Form("hNcluITSsa%s",spname[iSpec].Data());
 282     fHistNcluITSsa[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
 283     fHistNcluITSsa[iSpec]->Sumw2();
 284     fOutput->Add(fHistNcluITSsa[iSpec]);
 285
 286     hisname.Form("hNcluITSpureSA%s",spname[iSpec].Data());
 287     fHistNcluITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
 288     fHistNcluITSpureSA[iSpec]->Sumw2();
 289     fOutput->Add(fHistNcluITSpureSA[iSpec]);
 290
 291     //---
 292
 293     hisname.Form("hd0rphiITSpureSA%s",spname[iSpec].Data());
 294     fHistd0rphiITSpureSA[iSpec] = new TH2F(hisname.Data(),"",nbins,xbins,2000,-1,1);
 295     fHistd0rphiITSpureSA[iSpec]->Sumw2();
 296     fOutput->Add(fHistd0rphiITSpureSA[iSpec]);
 297
 298     hisname.Form("hd0zITSpureSA%s",spname[iSpec].Data());
 299     fHistd0zITSpureSA[iSpec] = new TH2F(hisname.Data(),"",nbins,xbins,2000,-1,1);
 300     fHistd0zITSpureSA[iSpec]->Sumw2();
 301     fOutput->Add(fHistd0zITSpureSA[iSpec]);
 302
 303     //---
 304
 305     hisname.Form("hCluInLayTPCITS%s",spname[iSpec].Data());
 306     fHistCluInLayTPCITS[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
 307     fHistCluInLayTPCITS[iSpec]->Sumw2();
 308     fOutput->Add(fHistCluInLayTPCITS[iSpec]);
 309
 310     hisname.Form("hCluInLayITSsa%s",spname[iSpec].Data());
 311     fHistCluInLayITSsa[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
 312     fHistCluInLayITSsa[iSpec]->Sumw2();
 313     fOutput->Add(fHistCluInLayITSsa[iSpec]);
 314
 315     hisname.Form("hCluInLayITSpureSA%s",spname[iSpec].Data());
 316     fHistCluInLayITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
 317     fHistCluInLayITSpureSA[iSpec]->Sumw2();
 318     fOutput->Add(fHistCluInLayITSpureSA[iSpec]);
 319
 320     hisname.Form("hOuterLayITSpureSA%s",spname[iSpec].Data());
 321     fHistOuterLayITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
 322     fHistOuterLayITSpureSA[iSpec]->Sumw2();
 323     fOutput->Add(fHistOuterLayITSpureSA[iSpec]);
 324
 325     //---
 326
 327     hisname.Form("hPtResid%s",spname[iSpec].Data());
 328     fHistPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
 329     fHistPtResid[iSpec]->Sumw2();
 330     fOutput->Add(fHistPtResid[iSpec]);
 331     hisname.Form("hPtRelResid%s",spname[iSpec].Data());
 332     fHistPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
 333     fHistPtRelResid[iSpec]->Sumw2();
 334     fOutput->Add(fHistPtRelResid[iSpec]);
 335
 336     hisname.Form("hInvPtResid%s",spname[iSpec].Data());
 337     fHistInvPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
 338     fHistInvPtResid[iSpec]->Sumw2();
 339     fOutput->Add(fHistInvPtResid[iSpec]);
 340     hisname.Form("hInvPtRelResid%s",spname[iSpec].Data());
 341     fHistInvPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
 342     fHistInvPtRelResid[iSpec]->Sumw2();
 343     fOutput->Add(fHistInvPtRelResid[iSpec]);
 344
 345     hisname.Form("hMCPtResid%s",spname[iSpec].Data());
 346     fHistMCPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
 347     fHistMCPtResid[iSpec]->Sumw2();
 348     fOutput->Add(fHistMCPtResid[iSpec]);
 349     hisname.Form("hMCPtRelResid%s",spname[iSpec].Data());
 350     fHistMCPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
 351     fHistMCPtRelResid[iSpec]->Sumw2();
 352     fOutput->Add(fHistMCPtRelResid[iSpec]);
 353
 354     hisname.Form("hMCInvPtResid%s",spname[iSpec].Data());
 355     fHistMCInvPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
 356     fHistMCInvPtResid[iSpec]->Sumw2();
 357     fOutput->Add(fHistMCInvPtResid[iSpec]);
 358     hisname.Form("hMCInvPtRelResid%s",spname[iSpec].Data());
 359     fHistMCInvPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
 360     fHistMCInvPtRelResid[iSpec]->Sumw2();
 361     fOutput->Add(fHistMCInvPtRelResid[iSpec]);
 362
 363   }
 364
 365   fHistMCPhiResid=new TH2F("hMCPhiResid","",nbins,xbins,100,-0.5,0.5);
 366   fHistMCPhiResid->Sumw2();
 367   fOutput->Add(fHistMCPhiResid);
 368   fHistPhiResid=new TH2F("hPhiResid","",nbins,xbins,100,-0.5,0.5);
 369   fHistPhiResid->Sumw2();
 370   fOutput->Add(fHistPhiResid);
 371
 372   PostData(1,fOutput);
 373
 374 }
 375 //______________________________________________________________________________
 376 void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
 377 {
 378   //
 379
 380   AliESDEvent *esd = (AliESDEvent*) (InputEvent());
 381
 382
 383   if(!esd) {
 384     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
 385     return;
 386   }
 387
 388
 389   if(!ESDfriend()) {
 390     printf("AliAnalysisTaskSDDRP::Exec(): bad ESDfriend\n");
 391     return;
 392   }
 393
 394   AliStack* stack=0;
 395
 396   if(fReadMC){
 397     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 398     if (!eventHandler) {
 399       Printf("ERROR: Could not retrieve MC event handler");
 400       return;
 401     }
 402     AliMCEvent* mcEvent = eventHandler->MCEvent();
 403     if (!mcEvent) {
 404       Printf("ERROR: Could not retrieve MC event");
 405       return;
 406     }
 407     stack = mcEvent->Stack();
 408     if (!stack) {
 409       Printf("ERROR: stack not available");
 410       return;
 411     }
 412   }
 413
 414
 415   fHistNEvents->Fill(0);
 416   const AliESDVertex *spdv=esd->GetPrimaryVertexSPD();
 417   if(spdv->GetNContributors()<=0) return;
 418   fHistNEvents->Fill(1);
 419
 420   const Int_t ntSize=36;
 421   Float_t xnt[ntSize];
 422
 423   Int_t ntracks = esd->GetNumberOfTracks();
 424   for (Int_t iTrack=0; iTrack < ntracks; iTrack++) {
 425     AliESDtrack * track = esd->GetTrack(iTrack);
 426     if (!track) continue;
 427     Float_t pttrack=track->Pt();
 428     Float_t ptrack=track->P();
 429     Float_t pxtrack=track->Px();
 430     Float_t pytrack=track->Py();
 431     Float_t pztrack=track->Pz();
 432     Float_t impactXY=-999, impactZ=-999;
 433     track->GetImpactParameters(impactXY, impactZ);
 434     Double_t dedxlay[4];
 435     track->GetITSdEdxSamples(dedxlay);
 436     Float_t dedx=track->GetITSsignal();
 437     Float_t chi2=track->GetITSchi2();
 438     Int_t status=track->GetStatus();
 439     Bool_t isITSrefit=kFALSE;
 440     Bool_t isTPCin=kFALSE;
 441     Bool_t isTPCrefit=kFALSE;
 442     Bool_t isPureSA=kFALSE;
 443     if(status&AliESDtrack::kITSrefit) isITSrefit=kTRUE;
 444     if(status&AliESDtrack::kTPCin) isTPCin=kTRUE;
 445     if(status&AliESDtrack::kTPCrefit) isTPCrefit=kTRUE;
 446     if(status&AliESDtrack::kITSpureSA) isPureSA=kTRUE;
 447     Bool_t isSA=kTRUE;
 448     if(status&AliESDtrack::kTPCin) isSA=kFALSE;
 449     Int_t nTPCclus=track->GetNcls(1);
 450     Int_t nITSclus=track->GetNcls(0);
 451     UChar_t clumap=track->GetITSClusterMap();
 452     Int_t idet;
 453     Float_t xloc,zloc;
 454     Int_t statusLay[6];
 455     for(Int_t iLay=0; iLay<6;iLay++) track->GetITSModuleIndexInfo(iLay,idet,statusLay[iLay],xloc,zloc);
 456     Int_t trlabel=track->GetLabel();
 457     Float_t ptgen=-999.;
 458     Float_t pgen=-999.;
 459     Float_t pxgen=-999.;
 460     Float_t pygen=-999.;
 461     Float_t pzgen=-999.;
 462     Float_t etagen=-999.;
 463     Float_t phigen=-999.;
 464     if(fReadMC){
 465       TParticle* part = stack->Particle(TMath::Abs(trlabel));
 466       ptgen=part->Pt();
 467       pgen=part->P();
 468       pxgen=part->Px();
 469       pygen=part->Py();
 470       pzgen=part->Pz();
 471       etagen=part->Eta();
 472       phigen=part->Phi();
 473     }
 474
 475     Int_t indexn=0;
 476     xnt[indexn++]=pttrack;
 477     xnt[indexn++]=pxtrack;
 478     xnt[indexn++]=pytrack;
 479     xnt[indexn++]=pztrack;
 480     xnt[indexn++]=track->Eta();
 481     xnt[indexn++]=track->Phi();
 482     xnt[indexn++]=impactXY;
 483     xnt[indexn++]=impactZ;
 484     for(Int_t iLay=0; iLay<4; iLay++) xnt[indexn++]=dedxlay[iLay];
 485     xnt[indexn++]=dedx;
 486     xnt[indexn++]=chi2;
 487     xnt[indexn++]=status;
 488     xnt[indexn++]=isITSrefit;
 489     xnt[indexn++]=isTPCin;
 490     xnt[indexn++]=isTPCrefit;
 491     xnt[indexn++]=isPureSA;
 492     xnt[indexn++]=nITSclus;
 493     xnt[indexn++]=nTPCclus;
 494     xnt[indexn++]=clumap;
 495     for(Int_t iLay=0; iLay<6; iLay++) xnt[indexn++]=(Float_t)statusLay[iLay];
 496     xnt[indexn++]=trlabel;
 497     xnt[indexn++]=ptgen;
 498     xnt[indexn++]=pxgen;
 499     xnt[indexn++]=pygen;
 500     xnt[indexn++]=pzgen;
 501     xnt[indexn++]=etagen;
 502     xnt[indexn++]=phigen;
 503     xnt[indexn++]=ntracks;
 504
 505     if(indexn>ntSize) printf("AliAnalysisTaskITSsaTracks: ERROR ntuple insexout of range\n");
 506     if(fFillNtuple) fNtupleTracks->Fill(xnt);
 507
 508     if(!(status&AliESDtrack::kITSrefit)) continue;
 509     if(nITSclus<fMinITSpts)continue;
 510     if((chi2/nITSclus) > fMaxITSChi2Clu) continue;
 511
 512     Int_t nSPDPoints=0;
 513     for(Int_t i=0; i<2; i++){
 514       if(clumap&(1<<i)) ++nSPDPoints;
 515     }
 516     if(nSPDPoints<fMinSPDpts) continue;
 517
 518     Int_t nPointsForPid=0;
 519     for(Int_t i=2; i<6; i++){
 520       if(clumap&(1<<i)) ++nPointsForPid;
 521     }
 522
 523     Int_t iTrackType=-1;
 524     if(status&AliESDtrack::kTPCin){
 525       iTrackType=kTypeTPCITS;
 526     }else{
 527       if(status&AliESDtrack::kITSpureSA) iTrackType=kTypeITSpureSA;
 528       else iTrackType=kTypeITSsa;
 529     }
 530     if(iTrackType<0 || iTrackType>=kNtrackTypes) continue;
 531     fHistPt[iTrackType]->Fill(pttrack);
 532     fHistEtaPhi[iTrackType]->Fill(track->Eta(),track->Phi());
 533     fHistChi2[iTrackType]->Fill(chi2/nITSclus);
 534     fHistNclu[iTrackType]->Fill(nITSclus);
 535     if(nPointsForPid==3){
 536       fHistdedxvsPt3cls[iTrackType]->Fill(pttrack,dedx);
 537     }
 538     if(nPointsForPid==4){
 539       fHistdedxvsPt4cls[iTrackType]->Fill(pttrack,dedx);
 540     }
 541     if(fReadMC){
 542       if(trlabel<0){
 543         fHistPtFake[iTrackType]->Fill(pttrack);
 544         fHistEtaPhiFake[iTrackType]->Fill(track->Eta(),track->Phi());
 545         fHistChi2Fake[iTrackType]->Fill(chi2/nITSclus);
 546         fHistNcluFake[iTrackType]->Fill(nITSclus);
 547       }else{
 548         fHistPtGood[iTrackType]->Fill(pttrack);
 549         fHistEtaPhiGood[iTrackType]->Fill(track->Eta(),track->Phi());
 550         fHistChi2Good[iTrackType]->Fill(chi2/nITSclus);
 551         fHistNcluGood[iTrackType]->Fill(nITSclus);
 552       }
 553     }
 554
 555     Int_t hadronSpecie=-1;
 556     if(fReadMC && fUseMCId){
 557       if(trlabel>=0){
 558         TParticle* part = stack->Particle(trlabel);
 559         Int_t pdg=TMath::Abs(part->GetPdgCode());
 560         if(pdg==211) hadronSpecie=kPion;
 561         else if(pdg==321) hadronSpecie=kKaon;
 562         else if(pdg==2212) hadronSpecie=kProton;
 563       }
 564     }else{
 565       if(nPointsForPid<fMinPtsforPid)continue;
 566       AliITSPIDResponse pidits(fReadMC);
 567       Float_t nSigmaPion=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kPion,nPointsForPid,isSA);
 568       if(nSigmaPion>-2. && nSigmaPion<2.){
 569         hadronSpecie=kPion;
 570       }else{
 571         Float_t nSigmaKaon=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kKaon,nPointsForPid,isSA);
 572         if(nSigmaKaon>-2. && nSigmaKaon<2.){
 573           hadronSpecie=kKaon;
 574         }else{
 575           Float_t nSigmaProton=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kProton,nPointsForPid,isSA);
 576           if(nSigmaProton>-2. && nSigmaProton<2.){
 577             hadronSpecie=kProton;
 578           }
 579         }
 580       }
 581     }
 582     if(hadronSpecie<0) continue;
 583     if(iTrackType==kTypeTPCITS){ // TPC+ITS tracks
 584       fHistPtTPCITS[hadronSpecie]->Fill(pttrack);
 585       fHistEtaPhiTPCITS[hadronSpecie]->Fill(track->Eta(),track->Phi());
 586       fHistNcluTPCITS[hadronSpecie]->Fill(pttrack,nITSclus);
 587       fHistCluInLayTPCITS[hadronSpecie]->Fill(-1.);
 588       for(Int_t iBit=0; iBit<6; iBit++){
 589         if(clumap&(1<<iBit)) fHistCluInLayTPCITS[hadronSpecie]->Fill(pttrack,iBit);
 590       }
 591     }else if(iTrackType==kTypeITSpureSA){ // TPC+ITS tracks
 592       fHistPtITSpureSA[hadronSpecie]->Fill(pttrack);
 593       fHistEtaPhiITSpureSA[hadronSpecie]->Fill(track->Eta(),track->Phi());
 594       fHistNcluITSpureSA[hadronSpecie]->Fill(pttrack,nITSclus);
 595       fHistd0rphiITSpureSA[hadronSpecie]->Fill(pttrack,impactXY);
 596       fHistd0zITSpureSA[hadronSpecie]->Fill(pttrack,impactZ);
 597       fHistCluInLayITSpureSA[hadronSpecie]->Fill(-1.);
 598       Int_t outerLay=-1;
 599       for(Int_t iBit=0; iBit<6; iBit++){
 600         if(clumap&(1<<iBit)){
 601           fHistCluInLayITSpureSA[hadronSpecie]->Fill(pttrack,iBit);
 602           if(iBit>outerLay) outerLay=iBit;
 603         }
 604         }
 605       fHistOuterLayITSpureSA[hadronSpecie]->Fill(pttrack,outerLay);
 606
 607
 608       if(fReadMC){
 609         if(trlabel>=0){
 610           TParticle* part = stack->Particle(trlabel);
 611           ptgen=part->Pt();
 612           Float_t invpttrack=track->OneOverPt();
 613           Float_t invptgen=0.;
 614           if(ptgen>0.) invptgen=1./ptgen;
 615           fHistMCPtResid[hadronSpecie]->Fill(pttrack,pttrack-ptgen);
 616           fHistMCPtRelResid[hadronSpecie]->Fill(pttrack,(pttrack-ptgen)/ptgen);
 617           fHistMCInvPtResid[hadronSpecie]->Fill(pttrack,invpttrack-invptgen);
 618           fHistMCInvPtRelResid[hadronSpecie]->Fill(pttrack,(invpttrack-invptgen)/invptgen);
 619           Float_t deltaphi=track->Phi()-part->Phi();
 620           if(deltaphi<-TMath::Pi()) deltaphi+=2*TMath::Pi();
 621           if(deltaphi>TMath::Pi()) deltaphi-=2*TMath::Pi();
 622           fHistMCPhiResid->Fill(pttrack,deltaphi);
 623         }
 624       }
 625     }else if(iTrackType==kTypeITSsa){ // TPC+ITS tracks
 626       fHistPtITSsa[hadronSpecie]->Fill(pttrack);
 627       fHistEtaPhiITSsa[hadronSpecie]->Fill(track->Eta(),track->Phi());
 628       fHistNcluITSsa[hadronSpecie]->Fill(pttrack,nITSclus);
 629       fHistCluInLayITSsa[hadronSpecie]->Fill(-1.);
 630       for(Int_t iBit=0; iBit<6; iBit++){
 631         if(clumap&(1<<iBit)) fHistCluInLayITSsa[hadronSpecie]->Fill(pttrack,iBit);
 632       }
 633     }
 634
 635     // match TPCITS with ITSpureSA
 636     if(nITSclus<fMinITSpts || nTPCclus<fMinTPCpts) continue;
 637     Int_t idxMI[12],idxSA[12];
 638     for(Int_t icl=0; icl<12; icl++){
 639       idxMI[icl]=-1;
 640       idxSA[icl]=-1;
 641     }
 642     Int_t ncls=track->GetClusters(0,idxMI);
 643     if(fMinITSpts<6){
 644       Bool_t accept=kTRUE;
 645       for(Int_t ilay=0; ilay<6; ilay++){
 646         if(fRequirePoint[ilay]){
 647           Int_t mask = 1<<ilay;
 648           if(!(clumap & mask)){
 649             accept=kFALSE;
 650             break;
 651           }
 652         }
 653       }
 654       if(!accept) continue;
 655     }
 656     // Sort
 657     for(Int_t i=0;i<12;i++){
 658       for(Int_t j=i+1;j<12;j++){
 659         if(idxMI[j]>idxMI[i]){
 660           Int_t tmp=idxMI[j];
 661           idxMI[j]=idxMI[i];
 662           idxMI[i]=tmp;
 663         }
 664       }
 665     }
 666     //    for(Int_t i=0; i<12; i++) printf("%d ",idxMI[i]);
 667     //    printf("\n");
 668     if(idxMI[0]<0 && idxMI[0]==idxMI[1]) continue;
 669     Bool_t matched=kFALSE;
 670     for (Int_t iTrack2 = 0; iTrack2 < ntracks; iTrack2++) {
 671       if(matched) break;
 672       if(iTrack2==iTrack) continue;
 673       AliESDtrack* track2 = esd->GetTrack(iTrack2);
 674       Int_t status2=track2->GetStatus();
 675       if(!(status2&AliESDtrack::kITSrefit)) continue;
 676       if(!(status2&AliESDtrack::kITSpureSA)) continue;
 677       Int_t clumap2=track2->GetITSClusterMap();
 678       Int_t nITSclus2=track2->GetNcls(0);
 679       Int_t nTPCclus2=track2->GetNcls(1);
 680       if(nITSclus2<fMinITSpts || nTPCclus2>0) continue;
 681       Int_t ncls2=track2->GetClusters(0,idxSA);
 682       if(ncls2!=ncls) continue;
 683       if(fMinITSpts<6){
 684         Bool_t accept=kTRUE;
 685         for(Int_t ilay=0; ilay<6; ilay++){
 686           if(fRequirePoint[ilay]){
 687             Int_t mask = 1<<ilay;
 688             if(!(clumap2 & mask)){
 689               accept=kFALSE;
 690               break;
 691             }
 692           }
 693         }
 694         if(!accept) continue;
 695       }
 696       // Sort
 697       for(Int_t i=0;i<12;i++){
 698         for(Int_t j=i+1;j<12;j++){
 699           if(idxSA[j]>idxSA[i]){
 700             Int_t tmp=idxSA[j];
 701             idxSA[j]=idxSA[i];
 702             idxSA[i]=tmp;
 703           }
 704         }
 705       }
 706       Int_t match=0;
 707       for(Int_t icl=0; icl<ncls; icl++){
 708         if(idxSA[icl]!=idxMI[icl]){
 709           match=0;
 710           break;
 711         }
 712         else match++;
 713       }
 714       if(match==ncls && match>0){
 715         matched=kTRUE;
 716         Float_t pt1=track->Pt();
 717         Float_t pt2=track2->Pt();
 718         Float_t ptm1=track->OneOverPt();
 719         Float_t ptm2=track2->OneOverPt();
 720         fHistPtResid[hadronSpecie]->Fill(pt1,pt2-pt1);
 721         fHistPtRelResid[hadronSpecie]->Fill(pt1,(pt2-pt1)/pt1);
 722         fHistInvPtResid[hadronSpecie]->Fill(pt1,ptm2-ptm1);
 723         fHistInvPtRelResid[hadronSpecie]->Fill(pt1,(ptm2-ptm1)/ptm1);
 724         Float_t deltaphi=track->Phi()-track2->Phi();
 725         if(deltaphi<-TMath::Pi()) deltaphi+=2*TMath::Pi();
 726         if(deltaphi>TMath::Pi()) deltaphi-=2*TMath::Pi();
 727         fHistPhiResid->Fill(pt1,deltaphi);
 728       }
 729     }
 730   }
 731
 732   PostData(1,fOutput);
 733
 734 }
 735 //______________________________________________________________________________
 736 void AliAnalysisTaskITSsaTracks::Terminate(Option_t */*option*/)
 737 {
 738   // Terminate analysis
 739   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 740   if (!fOutput) {
 741     printf("ERROR: fOutput not available\n");
 742     return;
 743   }
 744   fNtupleTracks = dynamic_cast<TNtuple*>(fOutput->FindObject("ntupleTracks"));
 745   fHistNEvents= dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 746
 747   TString spname[3]={"Pion","Kaon","Proton"};
 748
 749   for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
 750     fHistPtTPCITS[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtTPCITS%s",spname[iSpec].Data())));
 751     fHistPtITSsa[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtITSsa%s",spname[iSpec].Data())));
 752     fHistPtITSpureSA[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtITSpureSA%s",spname[iSpec].Data())));
 753
 754     fHistEtaPhiTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiTPCITS%s",spname[iSpec].Data())));
 755     fHistEtaPhiITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiITSsa%s",spname[iSpec].Data())));
 756     fHistEtaPhiITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiITSpureSA%s",spname[iSpec].Data())));
 757
 758     fHistNcluTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluTPCITS%s",spname[iSpec].Data())));
 759     fHistNcluITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluITSsa%s",spname[iSpec].Data())));
 760     fHistNcluITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluITSpureSA%s",spname[iSpec].Data())));
 761
 762     fHistd0rphiITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hd0rphiITSpureSA%s",spname[iSpec].Data())));
 763     fHistd0zITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hd0zITSpureSA%s",spname[iSpec].Data())));
 764
 765     fHistCluInLayTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayTPCITS%s",spname[iSpec].Data())));
 766     fHistCluInLayITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayITSsa%s",spname[iSpec].Data())));
 767     fHistCluInLayITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayITSpureSA%s",spname[iSpec].Data())));
 768
 769     fHistOuterLayITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hOuterLayITSpureSA%s",spname[iSpec].Data())));
 770
 771     fHistPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hPtResid%s",spname[iSpec].Data())));
 772     fHistPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hPtRelResid%s",spname[iSpec].Data())));
 773     fHistInvPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hInvPtResid%s",spname[iSpec].Data())));
 774     fHistInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hInvPtRelResid%s",spname[iSpec].Data())));
 775     fHistMCPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCPtResid%s",spname[iSpec].Data())));
 776     fHistMCPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCPtRelResid%s",spname[iSpec].Data())));
 777     fHistMCInvPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtResid%s",spname[iSpec].Data())));
 778     fHistMCInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtRelResid%s",spname[iSpec].Data())));
 779
 780   }
 781   fHistMCPhiResid= dynamic_cast<TH2F*>(fOutput->FindObject("hMCPhiResid"));
 782   fHistPhiResid= dynamic_cast<TH2F*>(fOutput->FindObject("hPhiResid"));
 783   return;
 784 }
 785
 786
 787
 788
 789