PWGJE/AliPWG4HighPtTrackQA.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-----------------------------------------------------------------------
  17 // This class stores QA variables as function of pT for different type
  18 // of tracks and track selection criteria
  19 // Output: Histograms for different set of cuts
  20 //-----------------------------------------------------------------------
  21 // Author : M. Verweij - UU
  22 //-----------------------------------------------------------------------
  23
  24 #ifndef ALIPWG4HIGHPTTRACKQA_CXX
  25 #define ALIPWG4HIGHPTTRACKQA_CXX
  26
  27 #include "AliPWG4HighPtTrackQA.h"
  28
  29 #include "TH1.h"
  30 #include "TH2.h"
  31 #include "TH3.h"
  32 #include "TProfile.h"
  33 #include "TList.h"
  34 #include "TFile.h"
  35 #include "TChain.h"
  36 #include "TH3F.h"
  37 #include "TKey.h"
  38 #include "TSystem.h"
  39 #include "TBits.h"
  40
  41 #include "AliAnalysisManager.h"
  42 #include "AliESDInputHandler.h"
  43 #include "AliMCEvent.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliStack.h"
  46 #include "AliESDtrack.h"
  47 #include "AliESDtrackCuts.h"
  48 #include "AliExternalTrackParam.h"
  49 #include "AliLog.h"
  50 #include "AliGenPythiaEventHeader.h"
  51 #include "AliGenCocktailEventHeader.h"
  52 #include "AliCentrality.h"
  53 #include "AliAODVertex.h"
  54 #include "AliAODEvent.h"
  55
  56 using namespace std; //required for resolving the 'cout' symbol
  57
  58 ClassImp(AliPWG4HighPtTrackQA)
  59
  60 AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA()
  61 : AliAnalysisTaskSE(),
  62   fDataType(kESD),
  63   fEvent(0x0),
  64   fESD(0x0),
  65   fVtx(0x0),
  66   fVtxAOD(0x0),
  67   fTrackCuts(0x0),
  68   fTrackCutsITSLoose(0x0),
  69   fTrackCutsTPConly(0x0),
  70   fTrackType(0),
  71   fFilterMask(0),
  72   fIncludeNoITS(kFALSE),
  73   fSigmaConstrainedMax(-1.),
  74   fPtMax(100.),
  75   fIsPbPb(0),
  76   fCentClass(10),
  77   fNVariables(26),
  78   fVariables(0x0),
  79   fITSClusterMap(0),
  80   fAvgTrials(1),
  81   fNEventAll(0),
  82   fNEventSel(0),
  83   fNEventReject(0),
  84   fh1Centrality(0x0),
  85   fh1Xsec(0),
  86   fh1Trials(0),
  87   fh1PtHard(0),
  88   fh1PtHardTrials(0),
  89   fh1NTracksAll(0x0),
  90   fh1NTracksReject(0x0),
  91   fh1NTracksSel(0x0),
  92   fPtAll(0),
  93   fPtSel(0),
  94   fPtPhi(0x0),
  95   fPtEta(0x0),
  96   fPtEtaPhi(0x0),
  97   fPtDCA2D(0x0),
  98   fPtDCAZ(0x0),
  99   fPtNClustersTPC(0x0),
 100   fPtNClustersTPCIter1(0x0),
 101   fPtNClustersTPCIter1Phi(0x0),
 102   fPtNClustersTPCShared(0x0),
 103   fPtNClustersTPCSharedFrac(0x0),
 104   fPtNPointITS(0x0),
 105   fPtNPointITSPhi(0x0),
 106   fPtChi2C(0x0),
 107   fPtNSigmaToVertex(0x0),
 108   fPtRelUncertainty1Pt(0x0),
 109   fPtRelUncertainty1PtNClus(0x0),
 110   fPtRelUncertainty1PtNClusIter1(0x0),
 111   fPtRelUncertainty1PtNPointITS(0x0),
 112   fPtRelUncertainty1PtITSClusterMap(0x0),
 113   fPtRelUncertainty1PtChi2(0x0),
 114   fPtRelUncertainty1PtChi2Iter1(0x0),
 115   fPtRelUncertainty1PtPhi(0x0),
 116   fPtChi2PerClusterTPC(0x0),
 117   fPtChi2PerClusterTPCIter1(0x0),
 118   fPtNCrossedRows(0x0),
 119   fPtNCrossedRowsPhi(0x0),
 120   fPtNCrossedRowsNClusFPhi(0x0),
 121   fPtNCrRNCrRNClusF(0x0),
 122   fPtNCrossedRowsFit(0x0),
 123   fPtNCrossedRowsFitPhi(0x0),
 124   fPtNCrossedRowsNClusFFitPhi(0x0),
 125   fNCrossedRowsNCrossedRowsFit(0x0),
 126   fNClustersNCrossedRows(0x0),
 127   fNClustersNCrossedRowsFit(0x0),
 128   fPtNClustersNClustersFitMap(0x0),
 129   fPtRelUncertainty1PtNCrossedRows(0x0),
 130   fPtRelUncertainty1PtNCrossedRowsFit(0x0),
 131   fPtChi2Gold(0x0),
 132   fPtChi2GGC(0x0),
 133   fPtChi2GoldPhi(0x0),
 134   fPtChi2GGCPhi(0x0),
 135   fChi2GoldChi2GGC(0x0),
 136   fPtChi2ITSPhi(0x0),
 137   fPtSigmaY2(0x0),
 138   fPtSigmaZ2(0x0),
 139   fPtSigmaSnp2(0x0),
 140   fPtSigmaTgl2(0x0),
 141   fPtSigma1Pt2(0x0),
 142   fProfPtSigmaY2(0x0),
 143   fProfPtSigmaZ2(0x0),
 144   fProfPtSigmaSnp2(0x0),
 145   fProfPtSigmaTgl2(0x0),
 146   fProfPtSigma1Pt2(0x0),
 147   fProfPtSigma1Pt(0x0),
 148   fProfPtPtSigma1Pt(0x0),
 149   fHistList(0)
 150 {
 151   //
 152   // Constructor
 153   //
 154   SetNVariables(26);
 155
 156   fPtBinEdges[0][0] = 10.;
 157   fPtBinEdges[0][1] = 1.;
 158   fPtBinEdges[1][0] = 20.;
 159   fPtBinEdges[1][1] = 2.;
 160   fPtBinEdges[2][0] = 100.;
 161   fPtBinEdges[2][1] = 5.;
 162
 163 }
 164 //________________________________________________________________________
 165 AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA(const char *name):
 166   AliAnalysisTaskSE(name),
 167   fDataType(kESD),
 168   fEvent(0x0),
 169   fESD(0x0),
 170   fVtx(0x0),
 171   fVtxAOD(0x0),
 172   fTrackCuts(0x0),
 173   fTrackCutsITSLoose(0x0),
 174   fTrackCutsTPConly(0x0),
 175   fTrackType(0),
 176   fFilterMask(0),
 177   fIncludeNoITS(kFALSE),
 178   fSigmaConstrainedMax(-1.),
 179   fPtMax(100.),
 180   fIsPbPb(0),
 181   fCentClass(10),
 182   fNVariables(26),
 183   fVariables(0x0),
 184   fITSClusterMap(0),
 185   fAvgTrials(1),
 186   fNEventAll(0),
 187   fNEventSel(0),
 188   fNEventReject(0),
 189   fh1Centrality(0x0),
 190   fh1Xsec(0),
 191   fh1Trials(0),
 192   fh1PtHard(0),
 193   fh1PtHardTrials(0),
 194   fh1NTracksAll(0x0),
 195   fh1NTracksReject(0x0),
 196   fh1NTracksSel(0x0),
 197   fPtAll(0),
 198   fPtSel(0),
 199   fPtPhi(0x0),
 200   fPtEta(0x0),
 201   fPtEtaPhi(0x0),
 202   fPtDCA2D(0x0),
 203   fPtDCAZ(0x0),
 204   fPtNClustersTPC(0x0),
 205   fPtNClustersTPCIter1(0x0),
 206   fPtNClustersTPCIter1Phi(0x0),
 207   fPtNClustersTPCShared(0x0),
 208   fPtNClustersTPCSharedFrac(0x0),
 209   fPtNPointITS(0x0),
 210   fPtNPointITSPhi(0x0),
 211   fPtChi2C(0x0),
 212   fPtNSigmaToVertex(0x0),
 213   fPtRelUncertainty1Pt(0x0),
 214   fPtRelUncertainty1PtNClus(0x0),
 215   fPtRelUncertainty1PtNClusIter1(0x0),
 216   fPtRelUncertainty1PtNPointITS(0x0),
 217   fPtRelUncertainty1PtITSClusterMap(0x0),
 218   fPtRelUncertainty1PtChi2(0x0),
 219   fPtRelUncertainty1PtChi2Iter1(0x0),
 220   fPtRelUncertainty1PtPhi(0x0),
 221   fPtChi2PerClusterTPC(0x0),
 222   fPtChi2PerClusterTPCIter1(0x0),
 223   fPtNCrossedRows(0x0),
 224   fPtNCrossedRowsPhi(0x0),
 225   fPtNCrossedRowsNClusFPhi(0x0),
 226   fPtNCrRNCrRNClusF(0x0),
 227   fPtNCrossedRowsFit(0x0),
 228   fPtNCrossedRowsFitPhi(0x0),
 229   fPtNCrossedRowsNClusFFitPhi(0x0),
 230   fNCrossedRowsNCrossedRowsFit(0x0),
 231   fNClustersNCrossedRows(0x0),
 232   fNClustersNCrossedRowsFit(0x0),
 233   fPtNClustersNClustersFitMap(0x0),
 234   fPtRelUncertainty1PtNCrossedRows(0x0),
 235   fPtRelUncertainty1PtNCrossedRowsFit(0x0),
 236   fPtChi2Gold(0x0),
 237   fPtChi2GGC(0x0),
 238   fPtChi2GoldPhi(0x0),
 239   fPtChi2GGCPhi(0x0),
 240   fChi2GoldChi2GGC(0x0),
 241   fPtChi2ITSPhi(0x0),
 242   fPtSigmaY2(0x0),
 243   fPtSigmaZ2(0x0),
 244   fPtSigmaSnp2(0x0),
 245   fPtSigmaTgl2(0x0),
 246   fPtSigma1Pt2(0x0),
 247   fProfPtSigmaY2(0x0),
 248   fProfPtSigmaZ2(0x0),
 249   fProfPtSigmaSnp2(0x0),
 250   fProfPtSigmaTgl2(0x0),
 251   fProfPtSigma1Pt2(0x0),
 252   fProfPtSigma1Pt(0x0),
 253   fProfPtPtSigma1Pt(0x0),
 254   fHistList(0)
 255 {
 256   //
 257   // Constructor. Initialization of Inputs and Outputs
 258   //
 259   AliDebug(2,Form("AliPWG4HighPtTrackQA Calling Constructor"));
 260
 261   SetNVariables(26);
 262
 263   fPtBinEdges[0][0] = 10.;
 264   fPtBinEdges[0][1] = 1.;
 265   fPtBinEdges[1][0] = 20.;
 266   fPtBinEdges[1][1] = 2.;
 267   fPtBinEdges[2][0] = 100.;
 268   fPtBinEdges[2][1] = 5.;
 269
 270   // Input slot #0 works with a TChain ESD
 271   DefineInput(0, TChain::Class());
 272   // Output slot #1 write into a TList
 273   DefineOutput(1, TList::Class());
 274 }
 275
 276 //________________________________________________________________________
 277 void AliPWG4HighPtTrackQA::SetPtBinEdges(Int_t region, Double_t ptmax, Double_t ptBinWidth) {
 278   //
 279   // Set variable bin sizes for pT axis in histos
 280   //
 281
 282   if(region<3) {
 283     fPtBinEdges[region][0] = ptmax;
 284     fPtBinEdges[region][1] = ptBinWidth;
 285   }
 286   else {
 287     AliError("Only 3 regions alowed. Use region 0/1/2\n");
 288     return;
 289   }
 290
 291 }
 292
 293 //________________________________________________________________________
 294 void AliPWG4HighPtTrackQA::UserCreateOutputObjects() {
 295   //Create output objects
 296   AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserCreateOutputObjects \n"));
 297
 298   Bool_t oldStatus = TH1::AddDirectoryStatus();
 299   TH1::AddDirectory(kFALSE);
 300
 301   OpenFile(1);
 302   fHistList = new TList();
 303   fHistList->SetOwner(kTRUE);
 304
 305   Float_t fgkPtMin = 0.;
 306   //  Float_t fgkPtMax = fPtMax;
 307
 308   //fPtBinEdges[region][0] = ptmax of region ; fPtBinEdges[region][1] = binWidth of region
 309   const Float_t ptmin1 =  fgkPtMin;
 310   const Float_t ptmax1 =  fPtBinEdges[0][0];
 311   const Float_t ptmin2 =  ptmax1 ;
 312   const Float_t ptmax2 =  fPtBinEdges[1][0];
 313   const Float_t ptmin3 =  ptmax2 ;
 314   const Float_t ptmax3 =  fPtBinEdges[2][0];//fgkPtMax;
 315   const Int_t nbin11 = (int)((ptmax1-ptmin1)/fPtBinEdges[0][1]);
 316   const Int_t nbin12 = (int)((ptmax2-ptmin2)/fPtBinEdges[1][1])+nbin11;
 317   const Int_t nbin13 = (int)((ptmax3-ptmin3)/fPtBinEdges[2][1])+nbin12;
 318   Int_t fgkNPtBins=nbin13;
 319   //Create array with low edges of each bin
 320   Double_t *binsPt=new Double_t[fgkNPtBins+1];
 321   for(Int_t i=0; i<=fgkNPtBins; i++) {
 322     if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
 323     if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
 324     if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
 325   }
 326
 327   Int_t fgkNPhiBins = 18*6;
 328   Float_t kMinPhi   = 0.;
 329   Float_t kMaxPhi   = 2.*TMath::Pi();
 330   Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
 331   for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
 332
 333   Int_t fgkNEtaBins=20;
 334   Float_t fgkEtaMin = -1.;
 335   Float_t fgkEtaMax = 1.;
 336   Double_t *binsEta=new Double_t[fgkNEtaBins+1];
 337   for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
 338
 339   Int_t fgkNNClustersTPCBins=80;
 340   Float_t fgkNClustersTPCMin = 0.5;
 341   Float_t fgkNClustersTPCMax = 160.5;
 342   Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
 343   for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
 344
 345   Int_t fgkNDCA2DBins=80;
 346   Float_t fgkDCA2DMin = -0.2;
 347   Float_t fgkDCA2DMax = 0.2;
 348   if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==7) {
 349     fgkDCA2DMin = -2.;
 350     fgkDCA2DMax = 2.;
 351   }
 352   Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
 353   for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
 354
 355   Int_t fgkNDCAZBins=80;
 356   Float_t fgkDCAZMin = -2.;
 357   Float_t fgkDCAZMax = 2.;
 358   if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
 359     fgkDCAZMin = -5.;
 360     fgkDCAZMax = 5.;
 361   }
 362   Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
 363   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
 364
 365   Int_t fgkNNPointITSBins=9;
 366   Float_t fgkNPointITSMin = -0.5;
 367   Float_t fgkNPointITSMax = 8.5;
 368   Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
 369   for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
 370
 371   Int_t fgkNITSClusterMapBins=65;
 372   Float_t fgkITSClusterMapMin = -0.5;
 373   Float_t fgkITSClusterMapMax = 64.5;
 374   Double_t *binsITSClusterMap=new Double_t[fgkNITSClusterMapBins+1];
 375   for(Int_t i=0; i<=fgkNITSClusterMapBins; i++) binsITSClusterMap[i]=(Double_t)fgkITSClusterMapMin + (fgkITSClusterMapMax-fgkITSClusterMapMin)/fgkNITSClusterMapBins*(Double_t)i ;
 376
 377
 378   Int_t fgkNNSigmaToVertexBins=9;
 379   Float_t fgkNSigmaToVertexMin = 0.;
 380   Float_t fgkNSigmaToVertexMax = 9.;
 381   Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
 382   for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
 383
 384   Int_t fgkNChi2CBins=10;
 385   //  Float_t fgkChi2CMin = 0.;
 386   //  Float_t fgkChi2CMax = 100.; //10 sigma
 387   Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
 388   for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i] = (Double_t)i * (Double_t)i;
 389
 390   Float_t fgkRel1PtUncertaintyMin = 0.;
 391   Float_t fgkRel1PtUncertaintyMax = 1.;
 392   Float_t binEdgeRel1PtUncertainty1= 0.3;
 393   Int_t fgkNRel1PtUncertaintyBins1 = 45;
 394   Float_t binWidthRel1PtUncertainty1 = (binEdgeRel1PtUncertainty1-fgkRel1PtUncertaintyMin)/((Float_t)fgkNRel1PtUncertaintyBins1);
 395   Int_t fgkNRel1PtUncertaintyBins2 = 35;
 396   Float_t binWidthRel1PtUncertainty2 = (fgkRel1PtUncertaintyMax-binEdgeRel1PtUncertainty1)/((Float_t)fgkNRel1PtUncertaintyBins2);
 397   Int_t fgkNRel1PtUncertaintyBins = fgkNRel1PtUncertaintyBins1 + fgkNRel1PtUncertaintyBins2;
 398
 399   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
 400   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) {
 401     if(i<=fgkNRel1PtUncertaintyBins1)
 402       binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (Double_t)binWidthRel1PtUncertainty1*(Double_t)i ;
 403     if(i<=fgkNRel1PtUncertaintyBins && i>fgkNRel1PtUncertaintyBins1)
 404       binsRel1PtUncertainty[i]=(Double_t)binEdgeRel1PtUncertainty1 + (Double_t)binWidthRel1PtUncertainty2*(Double_t)(i-fgkNRel1PtUncertaintyBins1);
 405   }
 406
 407   Int_t fgkNUncertainty1PtBins = 30;
 408   Float_t fgkUncertainty1PtMin = 0.;
 409   Float_t fgkUncertainty1PtMax = 0.1;
 410   if(fTrackType==1 || fTrackType==2 || fTrackType==4)
 411     fgkUncertainty1PtMax = 0.2;
 412   Double_t *binsUncertainty1Pt=new Double_t[fgkNUncertainty1PtBins+1];
 413   for(Int_t i=0; i<=fgkNUncertainty1PtBins; i++) binsUncertainty1Pt[i]=(Double_t)fgkUncertainty1PtMin + (fgkUncertainty1PtMax-fgkUncertainty1PtMin)/fgkNUncertainty1PtBins*(Double_t)i ;
 414
 415   Float_t fgkChi2PerClusMin = 0.;
 416   Float_t fgkChi2PerClusMax = 4.;
 417   Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
 418   Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
 419   for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
 420
 421   Int_t fgkNCrossedRowsNClusFBins  = 45;
 422   Float_t fgkNCrossedRowsNClusFMin = 0.;
 423   Float_t fgkNCrossedRowsNClusFMax = 1.5;
 424   Double_t *binsNCrossedRowsNClusF=new Double_t[fgkNCrossedRowsNClusFBins+1];
 425   for(Int_t i=0; i<=fgkNCrossedRowsNClusFBins; i++) binsNCrossedRowsNClusF[i]=(Double_t)fgkNCrossedRowsNClusFMin + (fgkNCrossedRowsNClusFMax-fgkNCrossedRowsNClusFMin)/fgkNCrossedRowsNClusFBins*(Double_t)i ;
 426
 427   Float_t fgk1PtMin = 0.;
 428   Float_t fgk1PtMax = 6.;
 429   Float_t binEdge1Pt1 = 1.;
 430   Float_t binWidth1Pt1 = 0.05;
 431   Int_t fgkN1PtBins1 = (int)((binEdge1Pt1-fgk1PtMin)/binWidth1Pt1);
 432   Float_t binWidth1Pt2 = 0.1;
 433   Int_t fgkN1PtBins2 = (int)((fgk1PtMax-binEdge1Pt1)/binWidth1Pt2);
 434   Int_t fgkN1PtBins = fgkN1PtBins1+fgkN1PtBins2;
 435   Double_t *bins1Pt=new Double_t[fgkN1PtBins+1];
 436
 437   for(Int_t i=0; i<=fgkN1PtBins; i++) {
 438     if(i<=fgkN1PtBins1)
 439       bins1Pt[i]=(Double_t)fgk1PtMin + (Double_t)(binEdge1Pt1-fgk1PtMin)/(Double_t)fgkN1PtBins1*(Double_t)i;
 440     if(i<=fgkN1PtBins && i>fgkN1PtBins1)
 441       bins1Pt[i]=(Double_t)binEdge1Pt1 + (Double_t)(fgk1PtMax-binEdge1Pt1)/(Double_t)fgkN1PtBins2*(Double_t)(i-fgkN1PtBins1);
 442   }
 443
 444   Int_t fgkNSigmaY2Bins = 50;
 445   Float_t fgkSigmaY2Min = 0.;
 446   Float_t fgkSigmaY2Max = 1.;
 447   if(fTrackType==1) fgkSigmaY2Max = 4.;
 448   if(fTrackType==2 || fTrackType==4) fgkSigmaY2Max = 0.1;
 449   Double_t *binsSigmaY2=new Double_t[fgkNSigmaY2Bins+1];
 450   for(Int_t i=0; i<=fgkNSigmaY2Bins; i++) binsSigmaY2[i]=(Double_t)fgkSigmaY2Min + (fgkSigmaY2Max-fgkSigmaY2Min)/fgkNSigmaY2Bins*(Double_t)i ;
 451
 452   Int_t fgkNSigmaZ2Bins = 50;
 453   Float_t fgkSigmaZ2Min = 0.;
 454   Float_t fgkSigmaZ2Max = 0.4;
 455   Double_t *binsSigmaZ2=new Double_t[fgkNSigmaZ2Bins+1];
 456   for(Int_t i=0; i<=fgkNSigmaZ2Bins; i++) binsSigmaZ2[i]=(Double_t)fgkSigmaZ2Min + (fgkSigmaZ2Max-fgkSigmaZ2Min)/fgkNSigmaZ2Bins*(Double_t)i ;
 457
 458   Int_t fgkNSigmaSnp2Bins = 50;
 459   Float_t fgkSigmaSnp2Min = 0.;
 460   Float_t fgkSigmaSnp2Max = 0.05;
 461   if(fTrackType==1) fgkSigmaSnp2Max = 0.2;
 462   if(fTrackType==2 || fTrackType==4) fgkSigmaSnp2Max = 0.1;
 463   Double_t *binsSigmaSnp2=new Double_t[fgkNSigmaSnp2Bins+1];
 464   for(Int_t i=0; i<=fgkNSigmaSnp2Bins; i++) binsSigmaSnp2[i]=(Double_t)fgkSigmaSnp2Min + (fgkSigmaSnp2Max-fgkSigmaSnp2Min)/fgkNSigmaSnp2Bins*(Double_t)i ;
 465
 466   Int_t fgkNSigmaTgl2Bins = 50;
 467   Float_t fgkSigmaTgl2Min = 0.;
 468   Float_t fgkSigmaTgl2Max = 0.1;
 469   if(fTrackType==1) fgkSigmaTgl2Max = 0.2;
 470   if(fTrackType==2 || fTrackType==4) fgkSigmaTgl2Max = 0.1;
 471   Double_t *binsSigmaTgl2=new Double_t[fgkNSigmaTgl2Bins+1];
 472   for(Int_t i=0; i<=fgkNSigmaTgl2Bins; i++) binsSigmaTgl2[i]=(Double_t)fgkSigmaTgl2Min + (fgkSigmaTgl2Max-fgkSigmaTgl2Min)/fgkNSigmaTgl2Bins*(Double_t)i ;
 473
 474   Int_t fgkNSigma1Pt2Bins = 50;
 475   Float_t fgkSigma1Pt2Min = 0.;
 476   Float_t fgkSigma1Pt2Max = 1.;
 477   Double_t *binsSigma1Pt2=new Double_t[fgkNSigma1Pt2Bins+1];
 478   for(Int_t i=0; i<=fgkNSigma1Pt2Bins; i++) binsSigma1Pt2[i]=(Double_t)fgkSigma1Pt2Min + (fgkSigma1Pt2Max-fgkSigma1Pt2Min)/fgkNSigma1Pt2Bins*(Double_t)i ;
 479
 480
 481   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 482   fHistList->Add(fNEventAll);
 483   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 484   fHistList->Add(fNEventSel);
 485   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
 486   //Set labels
 487   fNEventReject->Fill("noESD",0);
 488   fNEventReject->Fill("Trigger",0);
 489   fNEventReject->Fill("NTracks<2",0);
 490   fNEventReject->Fill("noVTX",0);
 491   fNEventReject->Fill("VtxStatus",0);
 492   fNEventReject->Fill("NCont<2",0);
 493   fNEventReject->Fill("ZVTX>10",0);
 494   fNEventReject->Fill("cent",0);
 495   fNEventReject->Fill("cent>90",0);
 496   fHistList->Add(fNEventReject);
 497
 498   fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
 499   fHistList->Add(fh1Centrality);
 500
 501   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 502   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 503   fHistList->Add(fh1Xsec);
 504
 505   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 506   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 507   fHistList->Add(fh1Trials);
 508
 509   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 510   fHistList->Add(fh1PtHard);
 511   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 512   fHistList->Add(fh1PtHardTrials);
 513
 514   fh1NTracksAll = new TH1F("fh1NTracksAll","fh1NTracksAll",1,-0.5,0.5);
 515   fHistList->Add(fh1NTracksAll);
 516
 517   fh1NTracksReject = new TH1F("fh1NTracksReject","fh1NTracksReject",1,-0.5,0.5);
 518   fh1NTracksReject->Fill("noHybridTrack",0);
 519   fh1NTracksReject->Fill("noITSrefit",0);
 520   fh1NTracksReject->Fill("noESDtrack",0);
 521   fh1NTracksReject->Fill("noTPCInner",0);
 522   fh1NTracksReject->Fill("FillTPC",0);
 523   fh1NTracksReject->Fill("noTPConly",0);
 524   fh1NTracksReject->Fill("relate",0);
 525   fh1NTracksReject->Fill("trackCuts",0);
 526   fh1NTracksReject->Fill("laser",0);
 527   fh1NTracksReject->Fill("chi2",0);
 528   fHistList->Add(fh1NTracksReject);
 529
 530   fh1NTracksSel = new TH1F("fh1NTracksSel","fh1NTracksSel",1,-0.5,0.5);
 531   fHistList->Add(fh1NTracksSel);
 532
 533   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
 534   fHistList->Add(fPtAll);
 535   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
 536   fHistList->Add(fPtSel);
 537
 538   fPtPhi = new TH2F("fPtPhi","fPtPhi",fgkNPtBins,binsPt,fgkNPhiBins,binsPhi);
 539   fHistList->Add(fPtPhi);
 540
 541   fPtEta = new TH2F("fPtEta","fPtEta",fgkNPtBins,binsPt,fgkNEtaBins,binsEta);
 542   fHistList->Add(fPtEta);
 543
 544   fPtEtaPhi = new TH3F("fPtEtaPhi","fPtEtaPhi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 545   fHistList->Add(fPtEtaPhi);
 546
 547   fPtDCA2D = new TH2F("fPtDCA2D","fPtDCA2D",fgkNPtBins,binsPt,fgkNDCA2DBins,binsDCA2D);
 548   fHistList->Add(fPtDCA2D);
 549
 550   fPtDCAZ = new TH2F("fPtDCAZ","fPtDCAZ",fgkNPtBins,binsPt,fgkNDCAZBins,binsDCAZ);
 551   fHistList->Add(fPtDCAZ);
 552
 553   fPtNClustersTPC = new TH2F("fPtNClustersTPC","fPtNClustersTPC",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 554   fHistList->Add(fPtNClustersTPC);
 555
 556   fPtNClustersTPCIter1 = new TH2F("fPtNClustersTPCIter1","fPtNClustersTPCIter1",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 557   fHistList->Add(fPtNClustersTPCIter1);
 558
 559   fPtNClustersTPCIter1Phi = new TH3F("fPtNClustersTPCIter1Phi","fPtNClustersTPCIter1Phi",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNPhiBins,binsPhi);
 560   fHistList->Add(fPtNClustersTPCIter1Phi);
 561
 562   fPtNClustersTPCShared = new TH2F("fPtNClustersTPCShared","fPtNClustersTPCShared",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 563   fHistList->Add(fPtNClustersTPCShared);
 564
 565   fPtNClustersTPCSharedFrac = new TH2F("fPtNClustersTPCSharedFrac","fPtNClustersTPCSharedFrac",fgkNPtBins,binsPt,fgkNSigma1Pt2Bins,binsSigma1Pt2);
 566   fHistList->Add(fPtNClustersTPCSharedFrac);
 567
 568   fPtNPointITS = new TH2F("fPtNPointITS","fPtNPointITS",fgkNPtBins,binsPt,fgkNNPointITSBins,binsNPointITS);
 569   fHistList->Add(fPtNPointITS);
 570
 571   fPtNPointITSPhi = new TH3F("fPtNPointITSPhi","fPtNPointITSPhi",fgkNPtBins,binsPt,fgkNNPointITSBins,binsNPointITS,fgkNPhiBins,binsPhi);
 572   fHistList->Add(fPtNPointITSPhi);
 573
 574   fPtChi2C = new TH2F("fPtChi2C","fPtChi2C",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
 575   fHistList->Add(fPtChi2C);
 576
 577   fPtNSigmaToVertex = new TH2F("fPtNSigmaToVertex","fPtNSigmaToVertex",fgkNPtBins,binsPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
 578   fHistList->Add(fPtNSigmaToVertex);
 579
 580   fPtRelUncertainty1Pt = new TH2F("fPtRelUncertainty1Pt","fPtRelUncertainty1Pt",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
 581   fHistList->Add(fPtRelUncertainty1Pt);
 582
 583   fPtRelUncertainty1PtNClus = new TH3F("fPtRelUncertainty1PtNClus","fPtRelUncertainty1PtNClus",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
 584   fHistList->Add(fPtRelUncertainty1PtNClus);
 585
 586   fPtRelUncertainty1PtNClusIter1 = new TH3F("fPtRelUncertainty1PtNClusIter1","fPtRelUncertainty1PtNClusIter1",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
 587   fHistList->Add(fPtRelUncertainty1PtNClusIter1);
 588
 589   fPtRelUncertainty1PtNPointITS = new TH3F("fPtRelUncertainty1PtNPointITS","fPtRelUncertainty1PtNPointITS",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNPointITSBins,binsNPointITS);
 590   fHistList->Add(fPtRelUncertainty1PtNPointITS);
 591
 592   fPtRelUncertainty1PtITSClusterMap = new TH3F("fPtRelUncertainty1PtITSClusterMap","fPtRelUncertainty1PtITSClusterMap",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNITSClusterMapBins,binsITSClusterMap);
 593   fHistList->Add(fPtRelUncertainty1PtITSClusterMap);
 594
 595   fPtRelUncertainty1PtChi2 = new TH3F("fPtRelUncertainty1PtChi2","fPtRelUncertainty1PtChi2",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNChi2PerClusBins,binsChi2PerClus);
 596   fHistList->Add(fPtRelUncertainty1PtChi2);
 597
 598   fPtRelUncertainty1PtChi2Iter1 = new TH3F("fPtRelUncertainty1PtChi2Iter1","fPtRelUncertainty1PtChi2Iter1",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNChi2PerClusBins,binsChi2PerClus);
 599   fHistList->Add(fPtRelUncertainty1PtChi2Iter1);
 600
 601   fPtRelUncertainty1PtPhi = new TH3F("fPtRelUncertainty1PtPhi","fPtRelUncertainty1PtPhi",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNPhiBins,binsPhi);
 602   fHistList->Add(fPtRelUncertainty1PtPhi);
 603
 604   fPtChi2PerClusterTPC = new TH2F("fPtChi2PerClusterTPC","fPtChi2PerClusterTPC",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
 605   fHistList->Add(fPtChi2PerClusterTPC);
 606
 607   fPtChi2PerClusterTPCIter1 = new TH2F("fPtChi2PerClusterTPCIter1","fPtChi2PerClusterTPCIter1",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
 608   fHistList->Add(fPtChi2PerClusterTPCIter1);
 609
 610   fPtNCrossedRows = new TH2F("fPtNCrossedRows","fPtNCrossedRows",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 611   fHistList->Add(fPtNCrossedRows);
 612
 613   fPtNCrossedRowsPhi = new TH3F("fPtNCrossedRowsPhi","fPtNCrossedRowsPhi",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNPhiBins,binsPhi);
 614   fHistList->Add(fPtNCrossedRowsPhi);
 615
 616   fPtNCrossedRowsNClusFPhi = new TH3F("fPtNCrossedRowsNClusFPhi","fPtNCrossedRowsNClusFPhi",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF,fgkNPhiBins,binsPhi);
 617   fHistList->Add(fPtNCrossedRowsNClusFPhi);
 618
 619   fPtNCrRNCrRNClusF = new TH3F("fPtNCrRNCrRNClusF","fPtNCrRNCrRNClusF",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
 620   fHistList->Add(fPtNCrRNCrRNClusF);
 621
 622   fPtNCrossedRowsFit = new TH2F("fPtNCrossedRowsFit","fPtNCrossedRowsFit",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 623   fHistList->Add(fPtNCrossedRowsFit);
 624
 625   fPtNCrossedRowsFitPhi = new TH3F("fPtNCrossedRowsFitPhi","fPtNCrossedRowsFitPhi",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNPhiBins,binsPhi);
 626   fHistList->Add(fPtNCrossedRowsFitPhi);
 627
 628   fPtNCrossedRowsNClusFFitPhi = new TH3F("fPtNCrossedRowsNClusFFitPhi","fPtNCrossedRowsNClusFFitPhi",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF,fgkNPhiBins,binsPhi);
 629   fHistList->Add(fPtNCrossedRowsNClusFFitPhi);
 630
 631   fNCrossedRowsNCrossedRowsFit = new TH2F("fNCrossedRowsNCrossedRowsFit","fNCrossedRowsNCrossedRowsFit",fgkNNClustersTPCBins,binsNClustersTPC,fgkNNClustersTPCBins,binsNClustersTPC);
 632   fHistList->Add(fNCrossedRowsNCrossedRowsFit);
 633
 634   fNClustersNCrossedRows = new TH2F("fNClustersNCrossedRows","fNClustersNCrossedRows",fgkNNClustersTPCBins,binsNClustersTPC,fgkNNClustersTPCBins,binsNClustersTPC);
 635   fHistList->Add(fNClustersNCrossedRows);
 636
 637   fNClustersNCrossedRowsFit = new TH2F("fNClustersNCrossedRowsFit","fNClustersNCrossedRowsFit",fgkNNClustersTPCBins,binsNClustersTPC,fgkNNClustersTPCBins,binsNClustersTPC);
 638   fHistList->Add(fNClustersNCrossedRowsFit);
 639
 640   fPtNClustersNClustersFitMap = new TH3F("fPtNClustersNClustersFitMap","fPtNClustersNClustersFitMap;p_{T};N_{cls};N_{cls}^{fit map}",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNNClustersTPCBins,binsNClustersTPC);
 641   fHistList->Add(fPtNClustersNClustersFitMap);
 642
 643   fPtRelUncertainty1PtNCrossedRows = new TH3F("fPtRelUncertainty1PtNCrossedRows","fPtRelUncertainty1PtNCrossedRows",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
 644   fHistList->Add(fPtRelUncertainty1PtNCrossedRows);
 645
 646   fPtRelUncertainty1PtNCrossedRowsFit = new TH3F("fPtRelUncertainty1PtNCrossedRowsFit","fPtRelUncertainty1PtNCrossedRowsFit",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
 647   fHistList->Add(fPtRelUncertainty1PtNCrossedRowsFit);
 648
 649   fPtChi2Gold = new TH2F("fPtChi2Gold","fPtChi2Gold",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
 650   fHistList->Add(fPtChi2Gold);
 651
 652   fPtChi2GGC = new TH2F("fPtChi2GGC","fPtChi2GGC",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
 653   fHistList->Add(fPtChi2GGC);
 654
 655   fPtChi2GoldPhi = new TH3F("fPtChi2GoldPhi","fPtChi2GoldPhi",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C,fgkNPhiBins,binsPhi);
 656   fHistList->Add(fPtChi2GoldPhi);
 657
 658   fPtChi2GGCPhi = new TH3F("fPtChi2GGCPhi","fPtChi2GGCPhi",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C,fgkNPhiBins,binsPhi);
 659   fHistList->Add(fPtChi2GGCPhi);
 660
 661   fChi2GoldChi2GGC = new TH2F("fChi2GoldChi2GGC","fChi2GoldChi2GGC;#chi^{2}_{gold};#chi^{2}_{ggc}",fgkNChi2CBins,binsChi2C,fgkNChi2CBins,binsChi2C);
 662   fHistList->Add(fChi2GoldChi2GGC);
 663
 664   fPtChi2ITSPhi = new TH3F("fPtChi2ITSPhi","fPtChi2ITSPhi;p_{T};#chi^{2}_{ITS};#varphi",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C,fgkNPhiBins,binsPhi);
 665   fHistList->Add(fPtChi2ITSPhi);
 666
 667   fPtSigmaY2 = new TH2F("fPtSigmaY2","fPtSigmaY2",fgkN1PtBins,bins1Pt,fgkNSigmaY2Bins,binsSigmaY2);
 668   fHistList->Add(fPtSigmaY2);
 669
 670   fPtSigmaZ2 = new TH2F("fPtSigmaZ2","fPtSigmaZ2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigmaZ2);
 671   fHistList->Add(fPtSigmaZ2);
 672
 673   fPtSigmaSnp2 = new TH2F("fPtSigmaSnp2","fPtSigmaSnp2",fgkN1PtBins,bins1Pt,fgkNSigmaSnp2Bins,binsSigmaSnp2);
 674   fHistList->Add(fPtSigmaSnp2);
 675
 676   fPtSigmaTgl2 = new TH2F("fPtSigmaTgl2","fPtSigmaTgl2",fgkN1PtBins,bins1Pt,fgkNSigmaTgl2Bins,binsSigmaTgl2);
 677   fHistList->Add(fPtSigmaTgl2);
 678
 679   fPtSigma1Pt2 = new TH2F("fPtSigma1Pt2","fPtSigma1Pt2",fgkN1PtBins,bins1Pt,fgkNSigma1Pt2Bins,binsSigma1Pt2);
 680   fHistList->Add(fPtSigma1Pt2);
 681
 682   fProfPtSigmaY2 = new TProfile("fProfPtSigmaY2","fProfPtSigmaY2",fgkN1PtBins,bins1Pt);
 683   fHistList->Add(fProfPtSigmaY2);
 684
 685   fProfPtSigmaZ2 = new TProfile("fProfPtSigmaZ2","fProfPtSigmaZ2",fgkN1PtBins,bins1Pt);
 686   fHistList->Add(fProfPtSigmaZ2);
 687
 688   fProfPtSigmaSnp2 = new TProfile("fProfPtSigmaSnp2","fProfPtSigmaSnp2",fgkN1PtBins,bins1Pt);
 689   fHistList->Add(fProfPtSigmaSnp2);
 690
 691   fProfPtSigmaTgl2 = new TProfile("fProfPtSigmaTgl2","fProfPtSigmaTgl2",fgkN1PtBins,bins1Pt);
 692   fHistList->Add(fProfPtSigmaTgl2);
 693
 694   fProfPtSigma1Pt2 = new TProfile("fProfPtSigma1Pt2","fProfPtSigma1Pt2",fgkN1PtBins,bins1Pt);
 695   fHistList->Add(fProfPtSigma1Pt2);
 696
 697   fProfPtSigma1Pt = new TProfile("fProfPtSigma1Pt","fProfPtSigma1Pt;p_{T};#sigma(1/p_{T})",fgkNPtBins,binsPt);
 698   fHistList->Add(fProfPtSigma1Pt);
 699
 700   fProfPtPtSigma1Pt = new TProfile("fProfPtPtSigma1Pt","fProfPtPtSigma1Pt;p_{T};p_{T}#sigma(1/p_{T})",fgkNPtBins,binsPt);
 701   fHistList->Add(fProfPtPtSigma1Pt);
 702
 703   TH1::AddDirectory(oldStatus);
 704
 705   PostData(1, fHistList);
 706
 707   if(binsPhi)               delete [] binsPhi;
 708   if(binsPt)                delete [] binsPt;
 709   if(binsNClustersTPC)      delete [] binsNClustersTPC;
 710   if(binsDCA2D)             delete [] binsDCA2D;
 711   if(binsDCAZ)              delete [] binsDCAZ;
 712   if(binsNPointITS)         delete [] binsNPointITS;
 713   if(binsITSClusterMap)     delete [] binsITSClusterMap;
 714   if(binsNSigmaToVertex)    delete [] binsNSigmaToVertex;
 715   if(binsChi2C)             delete [] binsChi2C;
 716   if(binsEta)               delete [] binsEta;
 717   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
 718   if(binsUncertainty1Pt)    delete [] binsUncertainty1Pt;
 719   if(binsChi2PerClus)       delete [] binsChi2PerClus;
 720   if(binsChi2PerClus)       delete [] binsNCrossedRowsNClusF;
 721   if(bins1Pt)               delete [] bins1Pt;
 722   if(binsSigmaY2)           delete [] binsSigmaY2;
 723   if(binsSigmaZ2)           delete [] binsSigmaZ2;
 724   if(binsSigmaSnp2)         delete [] binsSigmaSnp2;
 725   if(binsSigmaTgl2)         delete [] binsSigmaTgl2;
 726   if(binsSigma1Pt2)         delete [] binsSigma1Pt2;
 727 }
 728
 729 //________________________________________________________________________
 730 Bool_t AliPWG4HighPtTrackQA::SelectEvent() {
 731   //
 732   // Decide if event should be selected for analysis
 733   //
 734
 735   // Checks following requirements:
 736   // - fEvent available
 737   // - trigger info from AliPhysicsSelection
 738   // - MCevent available
 739   // - number of reconstructed tracks > 1
 740   // - primary vertex reconstructed
 741   // - z-vertex < 10 cm
 742   // - centrality in case of PbPb
 743
 744   Bool_t selectEvent = kTRUE;
 745
 746   //fEvent object available?
 747   if (!fEvent) {
 748     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
 749     fNEventReject->Fill("noAliVEvent",1);
 750     selectEvent = kFALSE;
 751     return selectEvent;
 752   }
 753
 754   //Check if number of reconstructed tracks is larger than 1
 755   if(!fEvent->GetNumberOfTracks() || fEvent->GetNumberOfTracks()<2)  {
 756     fNEventReject->Fill("NTracks<2",1);
 757     selectEvent = kFALSE;
 758     return selectEvent;
 759   }
 760
 761   //Check if vertex is reconstructed
 762   if(fDataType==kESD&&dynamic_cast<AliESDEvent*>(fEvent)) {
 763     fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexTracks();
 764
 765     if (!fVtx || !fVtx->GetStatus())
 766       fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexSPD();
 767
 768     if(!fVtx) {
 769       fNEventReject->Fill("noVTX",1);
 770       selectEvent = kFALSE;
 771       return selectEvent;
 772     }
 773
 774     if(!fVtx->GetStatus()) {
 775       fNEventReject->Fill("VtxStatus",1);
 776       selectEvent = kFALSE;
 777       return selectEvent;
 778     }
 779
 780     // Need vertex cut
 781     if(fVtx->GetNContributors()<2) {
 782       fNEventReject->Fill("NCont<2",1);
 783       selectEvent = kFALSE;
 784       return selectEvent;
 785     }
 786
 787     //Check if z-vertex < 10 cm
 788     double primVtx[3];
 789     fVtx->GetXYZ(primVtx);
 790     if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 791       fNEventReject->Fill("ZVTX>10",1);
 792       selectEvent = kFALSE;
 793       return selectEvent;
 794     }
 795   }
 796   else if(fDataType==kAOD&&dynamic_cast<AliAODEvent*>(fEvent)) {
 797     fVtxAOD = ((AliAODEvent*)fEvent)->GetPrimaryVertex();
 798     if(!fVtxAOD) {
 799       fNEventReject->Fill("noVTX",1);
 800       selectEvent = kFALSE;
 801       return selectEvent;
 802     }
 803
 804     // Need vertex cut
 805     if(fVtxAOD->GetNContributors()<2) {
 806       fNEventReject->Fill("NCont<2",1);
 807       selectEvent = kFALSE;
 808       return selectEvent;
 809     }
 810
 811     //Check if z-vertex < 10 cm
 812     double primVtx[3];
 813     fVtxAOD->GetXYZ(primVtx);
 814     if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 815       fNEventReject->Fill("ZVTX>10",1);
 816       selectEvent = kFALSE;
 817       return selectEvent;
 818     }
 819
 820   }
 821
 822   //Centrality selection should only be done in case of PbPb
 823   if(IsPbPb()) {
 824     Float_t cent = 0.;
 825     if(fCentClass!=CalculateCentrality(fEvent) && fCentClass!=10) {
 826       fNEventReject->Fill("cent",1);
 827       selectEvent = kFALSE;
 828       return selectEvent;
 829     }
 830     else {
 831       if(fDataType==kESD) {
 832         if(dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()) {
 833           cent = dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()->GetCentralityPercentile("V0M");
 834         }
 835       }
 836       else if(fDataType==kAOD) {
 837         if(dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality())
 838           cent = dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality();
 839        }
 840       if(cent>90.) {
 841         fNEventReject->Fill("cent>90",1);
 842         selectEvent = kFALSE;
 843         return selectEvent;
 844       }
 845       fh1Centrality->Fill(cent);
 846     }
 847   }
 848
 849   return selectEvent;
 850
 851 }
 852
 853 //________________________________________________________________________
 854 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliVEvent *ev){
 855   //
 856   // Get centrality from ESD or AOD
 857   //
 858
 859   if(fDataType==kESD)
 860     return CalculateCentrality(dynamic_cast<AliESDEvent*>(ev));
 861   else if(fDataType==kAOD)
 862     return CalculateCentrality(dynamic_cast<AliAODEvent*>(ev));
 863   else
 864     return 5;
 865 }
 866
 867 //________________________________________________________________________
 868 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliESDEvent *esd){
 869   //
 870   // Get centrality from ESD
 871   //
 872
 873   Float_t cent = -1;
 874
 875   if(esd){
 876     if(esd->GetCentrality()){
 877       cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
 878       if(fDebug>3) printf("centrality: %f\n",cent);
 879     }
 880   }
 881
 882   return GetCentralityClass(cent);
 883
 884 }
 885
 886 //________________________________________________________________________
 887 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(const AliAODEvent *aod){
 888   //
 889   // Get centrality from AOD
 890   //
 891
 892   if(!aod) return 5;
 893   Float_t cent = aod->GetHeader()->GetCentrality();
 894   if(fDebug>3) printf("centrality: %f\n",cent);
 895
 896   return GetCentralityClass(cent);
 897
 898 }
 899
 900 //________________________________________________________________________
 901 Int_t AliPWG4HighPtTrackQA::GetCentralityClass(Float_t cent) const {
 902   //
 903   // Get centrality class
 904   //
 905
 906   if(cent<0)  return 5; // OB - cent sometimes negative
 907   if(cent>80) return 4;
 908   if(cent>50) return 3;
 909   if(cent>30) return 2;
 910   if(cent>10) return 1;
 911   return 0;
 912
 913 }
 914
 915 //________________________________________________________________________
 916 void AliPWG4HighPtTrackQA::UserExec(Option_t *) {
 917   // Main loop
 918   // Called for each event
 919   AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserExec \n"));
 920
 921   fEvent = InputEvent();
 922   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 923
 924   // All events without selection
 925   fNEventAll->Fill(0.);
 926
 927   if(!SelectEvent()) {
 928     // Post output data
 929     PostData(1, fHistList);
 930     return;
 931   }
 932
 933
 934   //Need to keep track of selected events
 935   fNEventSel->Fill(0.);
 936
 937   fVariables = new TArrayF(fNVariables);
 938
 939   if(fDataType==kESD) DoAnalysisESD();
 940   if(fDataType==kAOD) DoAnalysisAOD();
 941
 942   //Delete old fVariables
 943   if(fVariables) delete fVariables;
 944
 945   // Post output data
 946   PostData(1, fHistList);
 947
 948 }
 949
 950 //________________________________________________________________________
 951 void AliPWG4HighPtTrackQA::DoAnalysisESD() {
 952   //
 953   // Run analysis on ESD
 954   //
 955
 956   if(!fESD) {
 957     PostData(1, fHistList);
 958     return;
 959   }
 960
 961   // ---- Get MC Header information (for MC productions in pThard bins) ----
 962   Double_t ptHard = 0.;
 963   Double_t nTrials = 1; // trials for MC trigger weight for real data
 964
 965   AliMCEventHandler *eventHandlerMC = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 966   if (eventHandlerMC) {
 967
 968     if(eventHandlerMC->MCEvent()){
 969       AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(eventHandlerMC->MCEvent());
 970       if(pythiaGenHeader){
 971         nTrials = pythiaGenHeader->Trials();
 972         ptHard  = pythiaGenHeader->GetPtHard();
 973
 974         fh1PtHard->Fill(ptHard);
 975         fh1PtHardTrials->Fill(ptHard,nTrials);
 976
 977         fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 978       }
 979     }
 980   }
 981
 982   Int_t nTracks = fESD->GetNumberOfTracks();
 983   AliDebug(2,Form("nTracks ESD%d", nTracks));
 984
 985   /*
 986     Variables to be put in fVariables
 987     0: pt
 988     1: phi
 989     2: eta
 990     3: dca2D
 991     4: dcaZ
 992     5: nClustersTPC
 993     6: nPointITS
 994     7: chi2C
 995     8: nSigmaToVertex
 996     9: trackLengthTPC
 997     10: chi2PerClusterTPC
 998     11: #crossed rows
 999     12: (#crossed rows)/(#findable clusters)
1000     13: SigmaY2
1001     14: SigmaZ2
1002     15: SigmaSnp2
1003     16: SigmaTgl2
1004     17: Sigma1Pt2
1005     18: NClustersTPCIter1
1006     19: Chi2TPCIter1
1007     20: nClustersTPCShared
1008     21: Golden Chi2 - global vs TPC constrained
1009     22: Chi2 between global and global constrained
1010     23: #crossed rows from fit map
1011     24: (#crossed rows)/(#findable clusters) from fit map
1012     25: chi2ITS
1013   */
1014
1015   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
1016     fh1NTracksAll->Fill(0.);
1017
1018     //Get track for analysis
1019     AliESDtrack *track = 0x0;
1020     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
1021     if(!esdtrack) {
1022       fh1NTracksReject->Fill("noESDtrack",1);
1023       continue;
1024     }
1025     AliESDtrack *origtrack = new AliESDtrack(*esdtrack);
1026     if(!origtrack)
1027       continue;
1028
1029     if(fTrackType==4) {
1030       if (!(fTrackCuts->AcceptTrack(esdtrack))) {
1031         fh1NTracksReject->Fill("trackCuts",1);
1032         if(origtrack) delete origtrack;
1033         continue;
1034       }
1035     }
1036
1037     if(fTrackType==1)
1038       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1039     else if(fTrackType==2 || fTrackType==4) {
1040       track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(fESD),esdtrack->GetID());
1041       if(!track) {
1042         fh1NTracksReject->Fill("noTPConly",1);
1043         if(origtrack) delete origtrack;
1044         continue;
1045       }
1046       AliExternalTrackParam exParam;
1047       Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1048       if( !relate ) {
1049         fh1NTracksReject->Fill("relate",1);
1050         if(track) delete track;
1051         if(origtrack) delete origtrack;
1052         continue;
1053       }
1054       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1055     }
1056     else if(fTrackType==5 || fTrackType==6) {
1057       if(fTrackCuts->AcceptTrack(esdtrack)) {
1058         if(origtrack) delete origtrack;
1059         continue;
1060       }
1061       else {
1062         if( !(fTrackCutsITSLoose->AcceptTrack(esdtrack)) && fTrackCutsTPConly->AcceptTrack(esdtrack) ) {
1063
1064           if(fTrackType==5) {
1065             //use TPConly constrained track
1066             track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1067             if(!track) {
1068               fh1NTracksReject->Fill("noTPConly",1);
1069               if(origtrack) delete origtrack;
1070               continue;
1071             }
1072             AliExternalTrackParam exParam;
1073             Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1074             if( !relate ) {
1075               fh1NTracksReject->Fill("relate",1);
1076               if(track) delete track;
1077               if(origtrack) delete origtrack;
1078               continue;
1079             }
1080             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1081           }
1082           else if(fTrackType==6) {
1083             //use global constrained track
1084             track = new AliESDtrack(*esdtrack);
1085             track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
1086
1087           }
1088         }
1089       }
1090     }
1091     else if(fTrackType==7) {
1092       //use global constrained track
1093       track = new AliESDtrack(*esdtrack);
1094     }
1095     else
1096       track = esdtrack;
1097
1098     if(!track) {
1099       if(origtrack) delete origtrack;
1100       continue;
1101     }
1102
1103     if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
1104       //Cut on chi2 of constrained fit
1105       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
1106         fh1NTracksReject->Fill("chi2",1);
1107         if(track) delete track;
1108         if(origtrack) delete origtrack;
1109         continue;
1110       }
1111     }
1112
1113     fPtAll->Fill(track->Pt());
1114
1115     if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
1116       fh1NTracksReject->Fill("trackCuts",1);
1117       if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
1118         if(track) delete track;
1119       }
1120       if(origtrack) delete origtrack;
1121       continue;
1122     }
1123
1124     if(fTrackType==7) {
1125       if(fTrackCutsITSLoose ) {
1126         if(fTrackCutsITSLoose->AcceptTrack(track) ) {
1127           if(track) delete track;
1128           if(origtrack) delete origtrack;
1129           continue;
1130         }
1131       }
1132
1133       if(esdtrack->GetConstrainedParam())
1134         track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
1135     }
1136
1137     if(!track) {
1138       if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
1139         if(track) delete track;
1140       }
1141       if(origtrack) delete origtrack;
1142       continue;
1143     }
1144
1145     fh1NTracksSel->Fill(0.);
1146
1147     fVariables->Reset(0.);
1148
1149     fVariables->SetAt(track->Pt(),0);
1150     fVariables->SetAt(track->Phi(),1);
1151     fVariables->SetAt(track->Eta(),2);
1152
1153     Float_t dca2D = 0.;
1154     Float_t dcaz  = 0.;
1155
1156     if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
1157       track->GetImpactParametersTPC(dca2D,dcaz); //TPConly
1158     }
1159     else
1160       track->GetImpactParameters(dca2D,dcaz);    //Global
1161
1162     fVariables->SetAt(dca2D,3);
1163     fVariables->SetAt(dcaz,4);
1164
1165     fVariables->SetAt((float)track->GetTPCNcls(),5);
1166
1167     Int_t nPointITS = 0;
1168     fITSClusterMap = track->GetITSClusterMap();
1169     UChar_t itsMap = track->GetITSClusterMap();
1170     for (Int_t i=0; i < 6; i++) {
1171       if (itsMap & (1 << i))
1172         nPointITS ++;
1173     }
1174     fVariables->SetAt((float)nPointITS,6);
1175     Float_t chi2C = (float)track->GetConstrainedChi2();
1176     if(fTrackType==1 || fTrackType==2 || fTrackType==4)
1177       chi2C = (float)track->GetConstrainedChi2TPC();
1178     fVariables->SetAt(chi2C,7);
1179     fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
1180
1181     fVariables->SetAt(GetTrackLengthTPC(track),9);
1182
1183     if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
1184
1185     //fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
1186     fVariables->SetAt(track->GetTPCCrossedRows(),11); //#crossed rows
1187
1188     Float_t crossedRowsTPCNClsF = 1.;//track->GetTPCClusterInfo(2,0);
1189     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
1190     fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
1191     fVariables->SetAt(track->GetSigmaY2(),13);
1192     fVariables->SetAt(track->GetSigmaZ2(),14);
1193     fVariables->SetAt(track->GetSigmaSnp2(),15);
1194     fVariables->SetAt(track->GetSigmaTgl2(),16);
1195     fVariables->SetAt(track->GetSigma1Pt2(),17);
1196
1197     fVariables->SetAt(track->GetTPCNclsIter1(),18);
1198     fVariables->SetAt(track->GetTPCchi2Iter1(),19);
1199
1200     fVariables->SetAt(track->GetTPCnclsS(),20);
1201
1202     Float_t chi2Gold = origtrack->GetChi2TPCConstrainedVsGlobal(fVtx);//GetGoldenChi2(origtrack);
1203     Float_t chi2GGC  = GetGGCChi2(origtrack);
1204
1205     fVariables->SetAt(chi2Gold,21);
1206     fVariables->SetAt(chi2GGC,22);
1207
1208     fVariables->SetAt(GetTPCClusterInfoFitMap(track,2,1),23);
1209     Float_t crossedRowsTPCNClsFFit = 1.;
1210     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsFFit = fVariables->At(23)/track->GetTPCNclsF();
1211     fVariables->SetAt(crossedRowsTPCNClsFFit,24);
1212
1213     fVariables->SetAt(track->GetITSchi2(),25);
1214
1215     TBits fitmap = track->GetTPCFitMap();
1216     fPtNClustersNClustersFitMap->Fill(track->Pt(),track->GetTPCNcls(),(float)fitmap.CountBits());
1217
1218     FillHistograms();
1219
1220     //      int mult = fTrackCuts->CountAcceptedTracks(fESD);
1221
1222     if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
1223       if(track) delete track;
1224     }
1225     if(origtrack) delete origtrack;
1226
1227   }//track loop
1228
1229 }
1230
1231 //________________________________________________________________________
1232 void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
1233   //
1234   // Do QA on AOD input
1235   //
1236   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1237   if(!aod) return;
1238   AliExternalTrackParam exParam;
1239   for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
1240
1241     AliAODTrack *aodtrack = aod->GetTrack(iTrack);
1242     if( !aodtrack->TestFilterMask(fFilterMask) ) {
1243       fh1NTracksReject->Fill("noHybridTrack",1);
1244       continue;
1245     }
1246
1247     if(!fIncludeNoITS) {
1248       if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
1249         fh1NTracksReject->Fill("noITSrefit",1);
1250         continue;
1251       }
1252     }
1253
1254     fVariables->Reset(0.);
1255
1256     fVariables->SetAt(aodtrack->Pt(),0);
1257     fVariables->SetAt(aodtrack->Phi(),1);
1258     fVariables->SetAt(aodtrack->Eta(),2);
1259
1260     Double_t dca[2] = {0.,0.};
1261     if(aodtrack->IsGlobalConstrained()) {
1262       dca[0] = aodtrack->DCA();
1263       dca[1] = aodtrack->ZAtDCA();
1264     } else {
1265       Double_t v[3]   = {0};
1266       Double_t pos[3] = {0};
1267       fVtxAOD->GetXYZ(v);
1268       aodtrack->GetXYZ(pos);
1269       dca[0] = pos[0] - v[0];
1270       dca[1] = pos[1] - v[1];
1271     }
1272     fVariables->SetAt(dca[0],3);
1273     fVariables->SetAt(dca[1],4);
1274     fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
1275     fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
1276     fVariables->SetAt(aodtrack->Chi2perNDF(),7);
1277     fVariables->SetAt(0.,8);
1278     fVariables->SetAt(GetTrackLengthTPC(aodtrack),9);
1279     fVariables->SetAt(aodtrack->Chi2perNDF(),10);
1280     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2,1,0,159,kFALSE),11);
1281     Float_t crossedRowsTPCNClsF = 0.;
1282     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
1283     fVariables->SetAt(crossedRowsTPCNClsF,12);
1284
1285     //get covariance matrix
1286     Double_t cov[21] = {0,};
1287     aodtrack->GetCovMatrix(cov);
1288     Double_t pxpypz[3] = {0,};
1289     aodtrack->PxPyPz(pxpypz);
1290     Double_t xyz[3] = {0,};
1291     aodtrack->GetXYZ(xyz);
1292     Short_t sign = aodtrack->Charge();
1293     exParam.Set(xyz,pxpypz,cov,sign);
1294
1295     fVariables->SetAt(exParam.GetSigmaY2(),13);
1296     fVariables->SetAt(exParam.GetSigmaZ2(),14);
1297     fVariables->SetAt(exParam.GetSigmaSnp2(),15);
1298     fVariables->SetAt(exParam.GetSigmaTgl2(),16);
1299     fVariables->SetAt(exParam.GetSigma1Pt2(),17);
1300
1301     fVariables->SetAt(0.,18); //NClustersTPCIter1
1302     fVariables->SetAt(0.,19); //Chi2TPCIter1
1303
1304     TBits sharedClusterMap = aodtrack->GetTPCSharedMap();
1305     fVariables->SetAt(sharedClusterMap.CountBits(),20);
1306
1307     fVariables->SetAt(0.,21); //not available in AOD golden chi2
1308     fVariables->SetAt(0.,22); //not available in AOD  Chi2 between global and global constrained
1309
1310     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2,1,0,159,kTRUE),23); //not available in AOD #crossed rows from fit map
1311     Float_t crossedRowsTPCNClsFFit = 0.;
1312     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsFFit = fVariables->At(23)/aodtrack->GetTPCNclsF();
1313     fVariables->SetAt(crossedRowsTPCNClsFFit,24); //(#crossed rows)/(#findable clusters) from fit map
1314
1315     fVariables->SetAt(0.,25);
1316
1317     fPtAll->Fill(fVariables->At(0));
1318
1319     FillHistograms();
1320   }
1321 }
1322
1323 //________________________________________________________________________
1324 void AliPWG4HighPtTrackQA::FillHistograms() {
1325   //
1326   // Fill all QA histograms
1327   //
1328
1329   fPtSel->Fill(fVariables->At(0));
1330   fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1331   fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1332   fPtEtaPhi->Fill(fVariables->At(0),fVariables->At(2),fVariables->At(1));
1333   fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1334   fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1335   fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1336   fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1337   fPtNPointITSPhi->Fill(fVariables->At(0),fVariables->At(6),fVariables->At(1));
1338
1339   fPtNClustersTPCIter1->Fill(fVariables->At(0),fVariables->At(18));
1340   fPtNClustersTPCIter1Phi->Fill(fVariables->At(0),fVariables->At(18),fVariables->At(1));
1341   fPtNClustersTPCShared->Fill(fVariables->At(0),fVariables->At(20));
1342   if(fVariables->At(5)>0.)
1343     fPtNClustersTPCSharedFrac->Fill(fVariables->At(0),fVariables->At(20)/fVariables->At(5));
1344
1345   if(fVariables->At(18)>0.)
1346     fPtChi2PerClusterTPCIter1->Fill(fVariables->At(0),fVariables->At(19)/fVariables->At(18));
1347
1348   fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1349   fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1350   fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1351   fPtRelUncertainty1PtNClus->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(5));
1352   fPtRelUncertainty1PtNClusIter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(18));
1353   fPtRelUncertainty1PtNPointITS->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(6));
1354
1355   fPtRelUncertainty1PtITSClusterMap->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),(int)fITSClusterMap);
1356
1357   fPtRelUncertainty1PtChi2->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(10));
1358   if(fVariables->At(18)>0.)
1359     fPtRelUncertainty1PtChi2Iter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(19)/fVariables->At(18));
1360   fPtRelUncertainty1PtPhi->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(1));
1361
1362   fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1363   fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1364   fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1365   fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1366   fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1367
1368   fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1369   fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1370   fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1371   fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1372   fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1373   fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1374   fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1375
1376   fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1377   fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1378   fPtNCrossedRowsPhi->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(1));
1379   fPtNCrossedRowsNClusFPhi->Fill(fVariables->At(0),fVariables->At(12),fVariables->At(1));
1380   fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1381
1382   fPtChi2Gold->Fill(fVariables->At(0),fVariables->At(21));
1383   fPtChi2GGC->Fill(fVariables->At(0),fVariables->At(22));
1384
1385   fPtChi2GoldPhi->Fill(fVariables->At(0),fVariables->At(21),fVariables->At(1));
1386   fPtChi2GGCPhi->Fill(fVariables->At(0),fVariables->At(22),fVariables->At(1));
1387
1388   fChi2GoldChi2GGC->Fill(fVariables->At(21),fVariables->At(22));
1389
1390   fPtNCrossedRowsFit->Fill(fVariables->At(0),fVariables->At(23));
1391   fPtNCrossedRowsFitPhi->Fill(fVariables->At(0),fVariables->At(23),fVariables->At(1));
1392   fPtNCrossedRowsNClusFFitPhi->Fill(fVariables->At(0),fVariables->At(24),fVariables->At(1));
1393   fNCrossedRowsNCrossedRowsFit->Fill(fVariables->At(11),fVariables->At(23));
1394
1395   fNClustersNCrossedRows->Fill(fVariables->At(5),fVariables->At(11));
1396   fNClustersNCrossedRowsFit->Fill(fVariables->At(5),fVariables->At(23));
1397
1398   fPtRelUncertainty1PtNCrossedRows->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(11));
1399   fPtRelUncertainty1PtNCrossedRowsFit->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(23));
1400
1401   if(fVariables->At(6)>0.)
1402     fPtChi2ITSPhi->Fill(fVariables->At(0),fVariables->At(25)/fVariables->At(6),fVariables->At(1));
1403
1404 }
1405
1406 //________________________________________________________________________
1407 Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1408   //
1409   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1410   // This is to called in Notify and should provide the path to the AOD/ESD file
1411   // Copied from AliAnalysisTaskJetSpectrum2
1412   //
1413
1414   TString file(currFile);
1415   fXsec = 0;
1416   fTrials = 1;
1417
1418   if(file.Contains("root_archive.zip#")){
1419     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1420     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1421     file.Replace(pos+1,20,"");
1422   }
1423   else {
1424     // not an archive take the basename....
1425     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1426   }
1427
1428   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
1429   if(!fxsec){
1430     // next trial fetch the histgram file
1431     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1432     if(!fxsec){
1433         // not a severe condition but inciate that we have no information
1434       return kFALSE;
1435     }
1436     else{
1437       // find the tlist we want to be independtent of the name so use the Tkey
1438       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1439       if(!key){
1440         fxsec->Close();
1441         return kFALSE;
1442       }
1443       TList *list = dynamic_cast<TList*>(key->ReadObj());
1444       if(!list){
1445         fxsec->Close();
1446         return kFALSE;
1447       }
1448       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1449       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1450       fxsec->Close();
1451     }
1452   } // no tree pyxsec.root
1453   else {
1454     TTree *xtree = (TTree*)fxsec->Get("Xsection");
1455     if(!xtree){
1456       fxsec->Close();
1457       return kFALSE;
1458     }
1459     UInt_t   ntrials  = 0;
1460     Double_t  xsection  = 0;
1461     xtree->SetBranchAddress("xsection",&xsection);
1462     xtree->SetBranchAddress("ntrials",&ntrials);
1463     xtree->GetEntry(0);
1464     fTrials = ntrials;
1465     fXsec = xsection;
1466     fxsec->Close();
1467   }
1468   return kTRUE;
1469 }
1470 //________________________________________________________________________
1471 Bool_t AliPWG4HighPtTrackQA::Notify()
1472 {
1473   //
1474   // Implemented Notify() to read the cross sections
1475   // and number of trials from pyxsec.root
1476   // Copied from AliAnalysisTaskJetSpectrum2
1477   //
1478
1479   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1480   Float_t xsection = 0;
1481   Float_t ftrials  = 1;
1482
1483   fAvgTrials = 1;
1484   if(tree){
1485     TFile *curfile = tree->GetCurrentFile();
1486     if (!curfile) {
1487       Error("Notify","No current file");
1488       return kFALSE;
1489     }
1490     if(!fh1Xsec||!fh1Trials){
1491       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1492       return kFALSE;
1493     }
1494     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1495     fh1Xsec->Fill("<#sigma>",xsection);
1496     // construct a poor man average trials
1497     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1498     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1499   }
1500   return kTRUE;
1501 }
1502
1503 //________________________________________________________________________
1504 AliGenPythiaEventHeader*  AliPWG4HighPtTrackQA::GetPythiaEventHeader(const AliMCEvent *mcEvent){
1505
1506   if(!mcEvent)return 0;
1507   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1508   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1509   if(!pythiaGenHeader){
1510     // cocktail ??
1511     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1512
1513     if (!genCocktailHeader) {
1514       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1515       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1516       return 0;
1517     }
1518     TList* headerList = genCocktailHeader->GetHeaders();
1519     for (Int_t i=0; i<headerList->GetEntries(); i++) {
1520       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1521       if (pythiaGenHeader)
1522         break;
1523     }
1524     if(!pythiaGenHeader){
1525       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1526       return 0;
1527     }
1528   }
1529   return pythiaGenHeader;
1530
1531 }
1532
1533 //_______________________________________________________________________
1534 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(const AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1, Bool_t useFitMap) const
1535 {
1536   //MV: copied from AliESDtrack since method is not available in AliAODTrack
1537
1538   //
1539   // TPC cluster information
1540   // type 0: get fraction of found/findable clusters with neighbourhood definition
1541   //      1: findable clusters with neighbourhood definition
1542   //      2: found clusters
1543   //
1544   // definition of findable clusters:
1545   //            a cluster is defined as findable if there is another cluster
1546   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1547   //           effects with a very simple algorithm.
1548   //
1549
1550   TBits fTPCClusterMap = 0;
1551   if(useFitMap)
1552     fTPCClusterMap = tr->GetTPCFitMap();
1553   else
1554     fTPCClusterMap = tr->GetTPCClusterMap();
1555
1556   if (type==2) return fTPCClusterMap.CountBits();
1557
1558   Int_t found=0;
1559   Int_t findable=0;
1560   Int_t last=-nNeighbours;
1561
1562   for (Int_t i=row0; i<row1; ++i){
1563     //look to current row
1564     if (fTPCClusterMap[i]) {
1565       last=i;
1566       ++found;
1567       ++findable;
1568       continue;
1569     }
1570     //look to nNeighbours before
1571     if ((i-last)<=nNeighbours) {
1572       ++findable;
1573       continue;
1574     }
1575     //look to nNeighbours after
1576     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1577       if (fTPCClusterMap[j]){
1578         ++findable;
1579         break;
1580       }
1581     }
1582   }
1583   if (type==1) return findable;
1584
1585   if (type==0){
1586     Float_t fraction=0;
1587     if (findable>0)
1588       fraction=(Float_t)found/(Float_t)findable;
1589     else
1590       fraction=0;
1591     return fraction;
1592   }
1593   return 0;  // undefined type - default value
1594 }
1595
1596 //_______________________________________________________________________
1597 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfoFitMap(const AliESDtrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1598 {
1599   //
1600   // TPC cluster information from fit map
1601   // type 0: get fraction of found/findable clusters with neighbourhood definition
1602   //      1: findable clusters with neighbourhood definition
1603   //      2: found clusters
1604   //
1605   // definition of findable clusters:
1606   //            a cluster is defined as findable if there is another cluster
1607   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1608   //           effects with a very simple algorithm.
1609   //
1610
1611   TBits fTPCFitMap = tr->GetTPCFitMap();
1612   if (type==2) return fTPCFitMap.CountBits();
1613
1614   Int_t found=0;
1615   Int_t findable=0;
1616   Int_t last=-nNeighbours;
1617
1618   for (Int_t i=row0; i<row1; ++i){
1619     //look to current row
1620     if (fTPCFitMap[i]) {
1621       last=i;
1622       ++found;
1623       ++findable;
1624       continue;
1625     }
1626     //look to nNeighbours before
1627     if ((i-last)<=nNeighbours) {
1628       ++findable;
1629       continue;
1630     }
1631     //look to nNeighbours after
1632     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1633       if (fTPCFitMap[j]){
1634         ++findable;
1635         break;
1636       }
1637     }
1638   }
1639   if (type==1) return findable;
1640
1641   if (type==0){
1642     Float_t fraction=0;
1643     if (findable>0)
1644       fraction=(Float_t)found/(Float_t)findable;
1645     else
1646       fraction=0;
1647     return fraction;
1648   }
1649   return 0;  // undefined type - default value
1650 }
1651
1652 //_______________________________________________________________________
1653 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(const AliESDtrack *track) const {
1654   //
1655   // returns distance between 1st and last hit in TPC
1656   // distance given in number of padrows
1657   //
1658
1659   TBits fTPCClusterMap = track->GetTPCClusterMap();
1660   int firstHit = 0;
1661   int lastHit = 0;
1662
1663   for(int i=0; i<=159; i++) {
1664     if(fTPCClusterMap[i]>0) firstHit = i;
1665   }
1666   for(int i=159; i>=0; i--) {
1667     if(fTPCClusterMap[i]>0) lastHit = i;
1668   }
1669
1670   Int_t trackLength = lastHit - firstHit;
1671
1672   return trackLength;
1673 }
1674
1675 //_______________________________________________________________________
1676 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(const AliAODTrack *track) const {
1677   //
1678   // returns distance between 1st and last hit in TPC
1679   // distance given in number of padrows
1680   //
1681
1682   TBits fTPCClusterMap = track->GetTPCClusterMap();
1683   int firstHit = 0;
1684   int lastHit = 0;
1685
1686   for(int i=0; i<=159; i++) {
1687     if(fTPCClusterMap[i]>0) firstHit = i;
1688   }
1689   for(int i=159; i>=0; i--) {
1690     if(fTPCClusterMap[i]>0) lastHit = i;
1691   }
1692
1693   Int_t trackLength = lastHit - firstHit;
1694
1695   return trackLength;
1696 }
1697
1698 //_______________________________________________________________________
1699 Float_t AliPWG4HighPtTrackQA::GetGoldenChi2(AliESDtrack *origtrack) {
1700   //
1701   // Return chi2 between global and TPC constrained track
1702   // track should be the global unconstrained track
1703   //
1704
1705   Float_t chi2Gold = 0.;
1706
1707   AliESDtrack *tpcTrack = 0x0;
1708   tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,origtrack->GetID());
1709   if(tpcTrack) {
1710     AliExternalTrackParam exParam;
1711     Bool_t relate = tpcTrack->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1712     if( relate ) {
1713       tpcTrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1714       //          Double_t pTPC[2],covTPC[3];     tpcTrack->PropagateToDCA(fVtx, fESD->GetMagneticField(), 10000,  pTPC, covTPC);
1715     }
1716
1717     tpcTrack->Propagate(origtrack->GetAlpha(), origtrack->GetX(), fESD->GetMagneticField());
1718     chi2Gold = (Float_t)origtrack->GetPredictedChi2(tpcTrack);
1719   }
1720
1721   if(tpcTrack) delete tpcTrack;
1722
1723   return chi2Gold;
1724
1725 }
1726
1727 //_______________________________________________________________________
1728 Float_t AliPWG4HighPtTrackQA::GetGGCChi2(AliESDtrack *origtrack) {
1729   //
1730   // Return chi2 between global and global constrained track
1731   // track should be the global unconstrained track
1732   //
1733
1734   Float_t chi2GGC = 0.;
1735
1736   AliESDtrack *esdtrackC = new AliESDtrack(*origtrack);
1737   if(esdtrackC) {
1738     if(origtrack->GetConstrainedParam()) {
1739       esdtrackC->Set(origtrack->GetConstrainedParam()->GetX(),origtrack->GetConstrainedParam()->GetAlpha(),origtrack->GetConstrainedParam()->GetParameter(),origtrack->GetConstrainedParam()->GetCovariance());
1740       chi2GGC = (Float_t)origtrack->GetPredictedChi2(esdtrackC);
1741     }
1742     delete esdtrackC;
1743   }
1744
1745   return chi2GGC;
1746
1747 }
1748
1749 //________________________________________________________________________
1750 void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1751 {
1752   // The Terminate() function is the last function to be called during
1753   // a query. It always runs on the client, it can be used to present
1754   // the results graphically or save the results to file.
1755
1756 }
1757
1758 #endif