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