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