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