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