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(0.,7); //ConstrainedChi2TPC -> not available in AOD
1292     fVariables->SetAt(0.,8);
1293     fVariables->SetAt(GetTrackLengthTPC(aodtrack),9);
1294     Float_t chi2pndf = aodtrack->Chi2perNDF();
1295     //if(fVariables->At(5)>0.) chi2pndf = aodtrack->GetTPCchi2()/fVariables->At(5);
1296     fVariables->SetAt(chi2pndf,10);
1297     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2,1,0,159,kFALSE),11);
1298     Float_t crossedRowsTPCNClsF = 0.;
1299     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
1300     fVariables->SetAt(crossedRowsTPCNClsF,12);
1301
1302     //get covariance matrix
1303     Double_t cov[21] = {0,};
1304     aodtrack->GetCovMatrix(cov);
1305     Double_t pxpypz[3] = {0,};
1306     aodtrack->PxPyPz(pxpypz);
1307     Double_t xyz[3] = {0,};
1308     aodtrack->GetXYZ(xyz);
1309     Short_t sign = aodtrack->Charge();
1310     exParam.Set(xyz,pxpypz,cov,sign);
1311
1312     fVariables->SetAt(exParam.GetSigmaY2(),13);
1313     fVariables->SetAt(exParam.GetSigmaZ2(),14);
1314     fVariables->SetAt(exParam.GetSigmaSnp2(),15);
1315     fVariables->SetAt(exParam.GetSigmaTgl2(),16);
1316     fVariables->SetAt(exParam.GetSigma1Pt2(),17);
1317
1318     fVariables->SetAt(0.,18); //NClustersTPCIter1
1319     fVariables->SetAt(0.,19); //Chi2TPCIter1
1320
1321     TBits sharedClusterMap = aodtrack->GetTPCSharedMap();
1322     fVariables->SetAt(sharedClusterMap.CountBits(),20);
1323
1324     fVariables->SetAt(0.,21); //not available in AOD golden chi2
1325     fVariables->SetAt(0.,22); //not available in AOD  Chi2 between global and global constrained
1326
1327     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2,1,0,159,kTRUE),23); //not available in AOD #crossed rows from fit map
1328     Float_t crossedRowsTPCNClsFFit = 0.;
1329     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsFFit = fVariables->At(23)/aodtrack->GetTPCNclsF();
1330     fVariables->SetAt(crossedRowsTPCNClsFFit,24); //(#crossed rows)/(#findable clusters) from fit map
1331
1332     fVariables->SetAt(0.,25);
1333
1334     fPtAll->Fill(fVariables->At(0));
1335
1336     FillHistograms();
1337   }
1338 }
1339
1340 //________________________________________________________________________
1341 void AliPWG4HighPtTrackQA::FillHistograms()
1342 {
1343   //
1344   // Fill all QA histograms
1345   //
1346
1347   fPtSel->Fill(fVariables->At(0));
1348   fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1349   fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1350   fPtEtaPhi->Fill(fVariables->At(0),fVariables->At(2),fVariables->At(1));
1351   fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1352   fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1353   fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1354   fPtNClustersTPCPhi->Fill(fVariables->At(1),fVariables->At(5));
1355   fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1356   fPtNPointITSPhi->Fill(fVariables->At(0),fVariables->At(6),fVariables->At(1));
1357
1358   fPtNClustersTPCIter1->Fill(fVariables->At(0),fVariables->At(18));
1359   fPtNClustersTPCIter1Phi->Fill(fVariables->At(0),fVariables->At(18),fVariables->At(1));
1360   fPtNClustersTPCShared->Fill(fVariables->At(0),fVariables->At(20));
1361   if(fVariables->At(5)>0.)
1362     fPtNClustersTPCSharedFrac->Fill(fVariables->At(0),fVariables->At(20)/fVariables->At(5));
1363
1364   if(fVariables->At(18)>0.)
1365     fPtChi2PerClusterTPCIter1->Fill(fVariables->At(0),fVariables->At(19)/fVariables->At(18));
1366
1367   fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1368   fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1369   fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1370   fPtRelUncertainty1PtNClus->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(5));
1371   fPtRelUncertainty1PtNClusIter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(18));
1372   fPtRelUncertainty1PtNPointITS->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(6));
1373
1374   fPtRelUncertainty1PtITSClusterMap->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),(int)fITSClusterMap);
1375
1376   fPtRelUncertainty1PtChi2->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(10));
1377   if(fVariables->At(18)>0.)
1378     fPtRelUncertainty1PtChi2Iter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(19)/fVariables->At(18));
1379   fPtRelUncertainty1PtPhi->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(1));
1380
1381   fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1382   fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1383   fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1384   fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1385   fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1386
1387   fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1388   fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1389   fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1390   fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1391   fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1392   fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1393   fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1394
1395   fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1396   fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1397   fPtNCrossedRowsPhi->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(1));
1398   fPtNCrossedRowsNClusFPhi->Fill(fVariables->At(0),fVariables->At(12),fVariables->At(1));
1399   fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1400
1401   fPtChi2Gold->Fill(fVariables->At(0),fVariables->At(21));
1402   fPtChi2GGC->Fill(fVariables->At(0),fVariables->At(22));
1403
1404   fPtChi2GoldPhi->Fill(fVariables->At(0),fVariables->At(21),fVariables->At(1));
1405   fPtChi2GGCPhi->Fill(fVariables->At(0),fVariables->At(22),fVariables->At(1));
1406
1407   fChi2GoldChi2GGC->Fill(fVariables->At(21),fVariables->At(22));
1408
1409   fPtNCrossedRowsFit->Fill(fVariables->At(0),fVariables->At(23));
1410   fPtNCrossedRowsFitPhi->Fill(fVariables->At(0),fVariables->At(23),fVariables->At(1));
1411   fPtNCrossedRowsNClusFFitPhi->Fill(fVariables->At(0),fVariables->At(24),fVariables->At(1));
1412   fNCrossedRowsNCrossedRowsFit->Fill(fVariables->At(11),fVariables->At(23));
1413
1414   fNClustersNCrossedRows->Fill(fVariables->At(5),fVariables->At(11));
1415   fNClustersNCrossedRowsFit->Fill(fVariables->At(5),fVariables->At(23));
1416
1417   fPtRelUncertainty1PtNCrossedRows->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(11));
1418   fPtRelUncertainty1PtNCrossedRowsFit->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(23));
1419
1420   if(fVariables->At(6)>0.)
1421     fPtChi2ITSPhi->Fill(fVariables->At(0),fVariables->At(25)/fVariables->At(6),fVariables->At(1));
1422
1423 }
1424
1425 //________________________________________________________________________
1426 Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials)
1427 {
1428   //
1429   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1430   // This is to called in Notify and should provide the path to the AOD/ESD file
1431   // Copied from AliAnalysisTaskJetSpectrum2
1432   //
1433
1434   TString file(currFile);
1435   fXsec = 0;
1436   fTrials = 1;
1437
1438   if(file.Contains("root_archive.zip#")){
1439     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1440     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1441     file.Replace(pos+1,20,"");
1442   }
1443   else {
1444     // not an archive take the basename....
1445     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1446   }
1447
1448   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
1449   if(!fxsec){
1450     // next trial fetch the histgram file
1451     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1452     if(!fxsec){
1453         // not a severe condition but inciate that we have no information
1454       return kFALSE;
1455     }
1456     else{
1457       // find the tlist we want to be independtent of the name so use the Tkey
1458       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1459       if(!key){
1460         fxsec->Close();
1461         return kFALSE;
1462       }
1463       TList *list = dynamic_cast<TList*>(key->ReadObj());
1464       if(!list){
1465         fxsec->Close();
1466         return kFALSE;
1467       }
1468       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1469       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1470       fxsec->Close();
1471     }
1472   } // no tree pyxsec.root
1473   else {
1474     TTree *xtree = (TTree*)fxsec->Get("Xsection");
1475     if(!xtree){
1476       fxsec->Close();
1477       return kFALSE;
1478     }
1479     UInt_t   ntrials  = 0;
1480     Double_t  xsection  = 0;
1481     xtree->SetBranchAddress("xsection",&xsection);
1482     xtree->SetBranchAddress("ntrials",&ntrials);
1483     xtree->GetEntry(0);
1484     fTrials = ntrials;
1485     fXsec = xsection;
1486     fxsec->Close();
1487   }
1488   return kTRUE;
1489 }
1490
1491 //________________________________________________________________________
1492 Bool_t AliPWG4HighPtTrackQA::Notify()
1493 {
1494   //
1495   // Implemented Notify() to read the cross sections
1496   // and number of trials from pyxsec.root
1497   // Copied from AliAnalysisTaskJetSpectrum2
1498   //
1499
1500   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1501   Float_t xsection = 0;
1502   Float_t ftrials  = 1;
1503
1504   fAvgTrials = 1;
1505   if(tree){
1506     TFile *curfile = tree->GetCurrentFile();
1507     if (!curfile) {
1508       Error("Notify","No current file");
1509       return kFALSE;
1510     }
1511     if(!fh1Xsec||!fh1Trials){
1512       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1513       return kFALSE;
1514     }
1515     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1516     fh1Xsec->Fill("<#sigma>",xsection);
1517     // construct a poor man average trials
1518     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1519     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1520   }
1521   return kTRUE;
1522 }
1523
1524 //________________________________________________________________________
1525 AliGenPythiaEventHeader*  AliPWG4HighPtTrackQA::GetPythiaEventHeader(const AliMCEvent *mcEvent)
1526 {
1527
1528   if(!mcEvent)return 0;
1529   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1530   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1531   if(!pythiaGenHeader){
1532     // cocktail ??
1533     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1534
1535     if (!genCocktailHeader) {
1536       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1537       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1538       return 0;
1539     }
1540     TList* headerList = genCocktailHeader->GetHeaders();
1541     for (Int_t i=0; i<headerList->GetEntries(); i++) {
1542       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1543       if (pythiaGenHeader)
1544         break;
1545     }
1546     if(!pythiaGenHeader){
1547       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1548       return 0;
1549     }
1550   }
1551   return pythiaGenHeader;
1552
1553 }
1554
1555 //_______________________________________________________________________
1556 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
1557 {
1558   //MV: copied from AliESDtrack since method is not available in AliAODTrack
1559
1560   //
1561   // TPC cluster information
1562   // type 0: get fraction of found/findable clusters with neighbourhood definition
1563   //      1: findable clusters with neighbourhood definition
1564   //      2: found clusters
1565   //
1566   // definition of findable clusters:
1567   //            a cluster is defined as findable if there is another cluster
1568   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1569   //           effects with a very simple algorithm.
1570   //
1571
1572   TBits fTPCClusterMap = 0;
1573   if(useFitMap)
1574     fTPCClusterMap = tr->GetTPCFitMap();
1575   else
1576     fTPCClusterMap = tr->GetTPCClusterMap();
1577
1578   if (type==2) return fTPCClusterMap.CountBits();
1579
1580   Int_t found=0;
1581   Int_t findable=0;
1582   Int_t last=-nNeighbours;
1583
1584   for (Int_t i=row0; i<row1; ++i){
1585     //look to current row
1586     if (fTPCClusterMap[i]) {
1587       last=i;
1588       ++found;
1589       ++findable;
1590       continue;
1591     }
1592     //look to nNeighbours before
1593     if ((i-last)<=nNeighbours) {
1594       ++findable;
1595       continue;
1596     }
1597     //look to nNeighbours after
1598     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1599       if (fTPCClusterMap[j]){
1600         ++findable;
1601         break;
1602       }
1603     }
1604   }
1605   if (type==1) return findable;
1606
1607   if (type==0){
1608     Float_t fraction=0;
1609     if (findable>0)
1610       fraction=(Float_t)found/(Float_t)findable;
1611     else
1612       fraction=0;
1613     return fraction;
1614   }
1615   return 0;  // undefined type - default value
1616 }
1617
1618 //_______________________________________________________________________
1619 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfoFitMap(const AliESDtrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1620 {
1621   //
1622   // TPC cluster information from fit map
1623   // type 0: get fraction of found/findable clusters with neighbourhood definition
1624   //      1: findable clusters with neighbourhood definition
1625   //      2: found clusters
1626   //
1627   // definition of findable clusters:
1628   //            a cluster is defined as findable if there is another cluster
1629   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1630   //           effects with a very simple algorithm.
1631   //
1632
1633   TBits fTPCFitMap = tr->GetTPCFitMap();
1634   if (type==2) return fTPCFitMap.CountBits();
1635
1636   Int_t found=0;
1637   Int_t findable=0;
1638   Int_t last=-nNeighbours;
1639
1640   for (Int_t i=row0; i<row1; ++i){
1641     //look to current row
1642     if (fTPCFitMap[i]) {
1643       last=i;
1644       ++found;
1645       ++findable;
1646       continue;
1647     }
1648     //look to nNeighbours before
1649     if ((i-last)<=nNeighbours) {
1650       ++findable;
1651       continue;
1652     }
1653     //look to nNeighbours after
1654     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1655       if (fTPCFitMap[j]){
1656         ++findable;
1657         break;
1658       }
1659     }
1660   }
1661   if (type==1) return findable;
1662
1663   if (type==0){
1664     Float_t fraction=0;
1665     if (findable>0)
1666       fraction=(Float_t)found/(Float_t)findable;
1667     else
1668       fraction=0;
1669     return fraction;
1670   }
1671   return 0;  // undefined type - default value
1672 }
1673
1674 //_______________________________________________________________________
1675 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(const AliESDtrack *track) const
1676 {
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 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(const AliAODTrack *track) const
1700 {
1701   //
1702   // returns distance between 1st and last hit in TPC
1703   // distance given in number of padrows
1704   //
1705
1706   TBits fTPCClusterMap = track->GetTPCClusterMap();
1707   int firstHit = 0;
1708   int lastHit = 0;
1709
1710   for(int i=0; i<=159; i++) {
1711     if(fTPCClusterMap[i]>0) firstHit = i;
1712   }
1713   for(int i=159; i>=0; i--) {
1714     if(fTPCClusterMap[i]>0) lastHit = i;
1715   }
1716
1717   Int_t trackLength = lastHit - firstHit;
1718
1719   return trackLength;
1720 }
1721
1722 //_______________________________________________________________________
1723 Float_t AliPWG4HighPtTrackQA::GetGoldenChi2(AliESDtrack *origtrack)
1724 {
1725   //
1726   // Return chi2 between global and TPC constrained track
1727   // track should be the global unconstrained track
1728   //
1729
1730   Float_t chi2Gold = 0.;
1731
1732   AliESDtrack *tpcTrack = 0x0;
1733   tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,origtrack->GetID());
1734   if(tpcTrack) {
1735     AliExternalTrackParam exParam;
1736     Bool_t relate = tpcTrack->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1737     if( relate ) {
1738       tpcTrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1739       //          Double_t pTPC[2],covTPC[3];     tpcTrack->PropagateToDCA(fVtx, fESD->GetMagneticField(), 10000,  pTPC, covTPC);
1740     }
1741
1742     tpcTrack->Propagate(origtrack->GetAlpha(), origtrack->GetX(), fESD->GetMagneticField());
1743     chi2Gold = (Float_t)origtrack->GetPredictedChi2(tpcTrack);
1744   }
1745
1746   if(tpcTrack) delete tpcTrack;
1747
1748   return chi2Gold;
1749
1750 }
1751
1752 //_______________________________________________________________________
1753 Float_t AliPWG4HighPtTrackQA::GetGGCChi2(AliESDtrack *origtrack)
1754 {
1755   //
1756   // Return chi2 between global and global constrained track
1757   // track should be the global unconstrained track
1758   //
1759
1760   Float_t chi2GGC = 0.;
1761
1762   AliESDtrack *esdtrackC = new AliESDtrack(*origtrack);
1763   if(esdtrackC) {
1764     if(origtrack->GetConstrainedParam()) {
1765       esdtrackC->Set(origtrack->GetConstrainedParam()->GetX(),origtrack->GetConstrainedParam()->GetAlpha(),origtrack->GetConstrainedParam()->GetParameter(),origtrack->GetConstrainedParam()->GetCovariance());
1766       chi2GGC = (Float_t)origtrack->GetPredictedChi2(esdtrackC);
1767     }
1768     delete esdtrackC;
1769   }
1770
1771   return chi2GGC;
1772
1773 }
1774
1775 //________________________________________________________________________
1776 void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1777 {
1778   // The Terminate() function is the last function to be called during
1779   // a query. It always runs on the client, it can be used to present
1780   // the results graphically or save the results to file.
1781
1782 }
1783
1784 #endif