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("noESDtrack",0);
 510   fh1NTracksReject->Fill("noTPCInner",0);
 511   fh1NTracksReject->Fill("FillTPC",0);
 512   fh1NTracksReject->Fill("noTPConly",0);
 513   fh1NTracksReject->Fill("relate",0);
 514   fh1NTracksReject->Fill("trackCuts",0);
 515   fh1NTracksReject->Fill("laser",0);
 516   fh1NTracksReject->Fill("chi2",0);
 517   fHistList->Add(fh1NTracksReject);
 518
 519   fh1NTracksSel = new TH1F("fh1NTracksSel","fh1NTracksSel",1,-0.5,0.5);
 520   fHistList->Add(fh1NTracksSel);
 521
 522   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
 523   fHistList->Add(fPtAll);
 524   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
 525   fHistList->Add(fPtSel);
 526
 527   fPtPhi = new TH2F("fPtPhi","fPtPhi",fgkNPtBins,binsPt,fgkNPhiBins,binsPhi);
 528   fHistList->Add(fPtPhi);
 529
 530   fPtEta = new TH2F("fPtEta","fPtEta",fgkNPtBins,binsPt,fgkNEtaBins,binsEta);
 531   fHistList->Add(fPtEta);
 532
 533   fPtDCA2D = new TH2F("fPtDCA2D","fPtDCA2D",fgkNPtBins,binsPt,fgkNDCA2DBins,binsDCA2D);
 534   fHistList->Add(fPtDCA2D);
 535
 536   fPtDCAZ = new TH2F("fPtDCAZ","fPtDCAZ",fgkNPtBins,binsPt,fgkNDCAZBins,binsDCAZ);
 537   fHistList->Add(fPtDCAZ);
 538
 539   fPtNClustersTPC = new TH2F("fPtNClustersTPC","fPtNClustersTPC",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 540   fHistList->Add(fPtNClustersTPC);
 541
 542   fPtNClustersTPCIter1 = new TH2F("fPtNClustersTPCIter1","fPtNClustersTPCIter1",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 543   fHistList->Add(fPtNClustersTPCIter1);
 544
 545   fPtNClustersTPCIter1Phi = new TH3F("fPtNClustersTPCIter1Phi","fPtNClustersTPCIter1Phi",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNPhiBins,binsPhi);
 546   fHistList->Add(fPtNClustersTPCIter1Phi);
 547
 548   fPtNClustersTPCShared = new TH2F("fPtNClustersTPCShared","fPtNClustersTPCShared",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 549   fHistList->Add(fPtNClustersTPCShared);
 550
 551   fPtNClustersTPCSharedFrac = new TH2F("fPtNClustersTPCSharedFrac","fPtNClustersTPCSharedFrac",fgkNPtBins,binsPt,fgkNSigma1Pt2Bins,binsSigma1Pt2);
 552   fHistList->Add(fPtNClustersTPCSharedFrac);
 553
 554   fPtNPointITS = new TH2F("fPtNPointITS","fPtNPointITS",fgkNPtBins,binsPt,fgkNNPointITSBins,binsNPointITS);
 555   fHistList->Add(fPtNPointITS);
 556
 557   fPtChi2C = new TH2F("fPtChi2C","fPtChi2C",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
 558   fHistList->Add(fPtChi2C);
 559
 560   fPtNSigmaToVertex = new TH2F("fPtNSigmaToVertex","fPtNSigmaToVertex",fgkNPtBins,binsPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
 561   fHistList->Add(fPtNSigmaToVertex);
 562
 563   fPtRelUncertainty1Pt = new TH2F("fPtRelUncertainty1Pt","fPtRelUncertainty1Pt",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
 564   fHistList->Add(fPtRelUncertainty1Pt);
 565
 566   fPtRelUncertainty1PtNClus = new TH3F("fPtRelUncertainty1PtNClus","fPtRelUncertainty1PtNClus",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
 567   fHistList->Add(fPtRelUncertainty1PtNClus);
 568
 569   fPtRelUncertainty1PtNClusIter1 = new TH3F("fPtRelUncertainty1PtNClusIter1","fPtRelUncertainty1PtNClusIter1",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
 570   fHistList->Add(fPtRelUncertainty1PtNClusIter1);
 571
 572   fPtRelUncertainty1PtNPointITS = new TH3F("fPtRelUncertainty1PtNPointITS","fPtRelUncertainty1PtNPointITS",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNPointITSBins,binsNPointITS);
 573   fHistList->Add(fPtRelUncertainty1PtNPointITS);
 574
 575   fPtRelUncertainty1PtITSClusterMap = new TH3F("fPtRelUncertainty1PtITSClusterMap","fPtRelUncertainty1PtITSClusterMap",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNITSClusterMapBins,binsITSClusterMap);
 576   fHistList->Add(fPtRelUncertainty1PtITSClusterMap);
 577
 578   fPtRelUncertainty1PtChi2 = new TH3F("fPtRelUncertainty1PtChi2","fPtRelUncertainty1PtChi2",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNChi2PerClusBins,binsChi2PerClus);
 579   fHistList->Add(fPtRelUncertainty1PtChi2);
 580
 581   fPtRelUncertainty1PtChi2Iter1 = new TH3F("fPtRelUncertainty1PtChi2Iter1","fPtRelUncertainty1PtChi2Iter1",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNChi2PerClusBins,binsChi2PerClus);
 582   fHistList->Add(fPtRelUncertainty1PtChi2Iter1);
 583
 584   fPtRelUncertainty1PtPhi = new TH3F("fPtRelUncertainty1PtPhi","fPtRelUncertainty1PtPhi",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNPhiBins,binsPhi);
 585   fHistList->Add(fPtRelUncertainty1PtPhi);
 586
 587   fPtUncertainty1Pt = new TH2F("fPtUncertainty1Pt","fPtUncertainty1Pt",fgkNPtBins,binsPt,fgkNUncertainty1PtBins,binsUncertainty1Pt);
 588   fHistList->Add(fPtUncertainty1Pt);
 589
 590   fPtChi2PerClusterTPC = new TH2F("fPtChi2PerClusterTPC","fPtChi2PerClusterTPC",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
 591   fHistList->Add(fPtChi2PerClusterTPC);
 592
 593   fPtChi2PerClusterTPCIter1 = new TH2F("fPtChi2PerClusterTPCIter1","fPtChi2PerClusterTPCIter1",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
 594   fHistList->Add(fPtChi2PerClusterTPCIter1);
 595
 596   fPtNCrossedRows = new TH2F("fPtNCrossedRows","fPtNCrossedRows",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 597   fHistList->Add(fPtNCrossedRows);
 598
 599   fPtNCrossedRowsPhi = new TH3F("fPtNCrossedRowsPhi","fPtNCrossedRowsPhi",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNPhiBins,binsPhi);
 600   fHistList->Add(fPtNCrossedRowsPhi);
 601
 602   fPtNCrossedRowsNClusFPhi = new TH3F("fPtNCrossedRowsNClusFPhi","fPtNCrossedRowsNClusFPhi",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF,fgkNPhiBins,binsPhi);
 603   fHistList->Add(fPtNCrossedRowsNClusFPhi);
 604
 605   fPtNCrRNCrRNClusF = new TH3F("fPtNCrRNCrRNClusF","fPtNCrRNCrRNClusF",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
 606   fHistList->Add(fPtNCrRNCrRNClusF);
 607
 608   fPtNCrossedRowsFit = new TH2F("fPtNCrossedRowsFit","fPtNCrossedRowsFit",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 609   fHistList->Add(fPtNCrossedRowsFit);
 610
 611   fPtNCrossedRowsFitPhi = new TH3F("fPtNCrossedRowsFitPhi","fPtNCrossedRowsFitPhi",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNPhiBins,binsPhi);
 612   fHistList->Add(fPtNCrossedRowsFitPhi);
 613
 614   fPtNCrossedRowsNClusFFitPhi = new TH3F("fPtNCrossedRowsNClusFFitPhi","fPtNCrossedRowsNClusFFitPhi",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF,fgkNPhiBins,binsPhi);
 615   fHistList->Add(fPtNCrossedRowsNClusFFitPhi);
 616
 617   fNCrossedRowsNCrossedRowsFit = new TH2F("fNCrossedRowsNCrossedRowsFit","fNCrossedRowsNCrossedRowsFit",fgkNNClustersTPCBins,binsNClustersTPC,fgkNNClustersTPCBins,binsNClustersTPC);
 618   fHistList->Add(fNCrossedRowsNCrossedRowsFit);
 619
 620   fNClustersNCrossedRows = new TH2F("fNClustersNCrossedRows","fNClustersNCrossedRows",fgkNNClustersTPCBins,binsNClustersTPC,fgkNNClustersTPCBins,binsNClustersTPC);
 621   fHistList->Add(fNClustersNCrossedRows);
 622
 623   fNClustersNCrossedRowsFit = new TH2F("fNClustersNCrossedRowsFit","fNClustersNCrossedRowsFit",fgkNNClustersTPCBins,binsNClustersTPC,fgkNNClustersTPCBins,binsNClustersTPC);
 624   fHistList->Add(fNClustersNCrossedRowsFit);
 625
 626   fPtRelUncertainty1PtNCrossedRows = new TH3F("fPtRelUncertainty1PtNCrossedRows","fPtRelUncertainty1PtNCrossedRows",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
 627   fHistList->Add(fPtRelUncertainty1PtNCrossedRows);
 628
 629   fPtRelUncertainty1PtNCrossedRowsFit = new TH3F("fPtRelUncertainty1PtNCrossedRowsFit","fPtRelUncertainty1PtNCrossedRowsFit",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNNClustersTPCBins,binsNClustersTPC);
 630   fHistList->Add(fPtRelUncertainty1PtNCrossedRowsFit);
 631
 632   fPtChi2Gold = new TH2F("fPtChi2Gold","fPtChi2Gold",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
 633   fHistList->Add(fPtChi2Gold);
 634
 635   fPtChi2GGC = new TH2F("fPtChi2GGC","fPtChi2GGC",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
 636   fHistList->Add(fPtChi2GGC);
 637
 638   fPtChi2GoldPhi = new TH3F("fPtChi2GoldPhi","fPtChi2GoldPhi",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C,fgkNPhiBins,binsPhi);
 639   fHistList->Add(fPtChi2GoldPhi);
 640
 641   fPtChi2GGCPhi = new TH3F("fPtChi2GGCPhi","fPtChi2GGCPhi",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C,fgkNPhiBins,binsPhi);
 642   fHistList->Add(fPtChi2GGCPhi);
 643
 644   fChi2GoldChi2GGC = new TH2F("fChi2GoldChi2GGC","fChi2GoldChi2GGC;#chi^{2}_{gold};#chi^{2}_{ggc}",fgkNChi2CBins,binsChi2C,fgkNChi2CBins,binsChi2C);
 645   fHistList->Add(fChi2GoldChi2GGC);
 646
 647
 648   fPtSigmaY2 = new TH2F("fPtSigmaY2","fPtSigmaY2",fgkN1PtBins,bins1Pt,fgkNSigmaY2Bins,binsSigmaY2);
 649   fHistList->Add(fPtSigmaY2);
 650
 651   fPtSigmaZ2 = new TH2F("fPtSigmaZ2","fPtSigmaZ2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigmaZ2);
 652   fHistList->Add(fPtSigmaZ2);
 653
 654   fPtSigmaSnp2 = new TH2F("fPtSigmaSnp2","fPtSigmaSnp2",fgkN1PtBins,bins1Pt,fgkNSigmaSnp2Bins,binsSigmaSnp2);
 655   fHistList->Add(fPtSigmaSnp2);
 656
 657   fPtSigmaTgl2 = new TH2F("fPtSigmaTgl2","fPtSigmaTgl2",fgkN1PtBins,bins1Pt,fgkNSigmaTgl2Bins,binsSigmaTgl2);
 658   fHistList->Add(fPtSigmaTgl2);
 659
 660   fPtSigma1Pt2 = new TH2F("fPtSigma1Pt2","fPtSigma1Pt2",fgkN1PtBins,bins1Pt,fgkNSigma1Pt2Bins,binsSigma1Pt2);
 661   fHistList->Add(fPtSigma1Pt2);
 662
 663   fProfPtSigmaY2 = new TProfile("fProfPtSigmaY2","fProfPtSigmaY2",fgkN1PtBins,bins1Pt);
 664   fHistList->Add(fProfPtSigmaY2);
 665
 666   fProfPtSigmaZ2 = new TProfile("fProfPtSigmaZ2","fProfPtSigmaZ2",fgkN1PtBins,bins1Pt);
 667   fHistList->Add(fProfPtSigmaZ2);
 668
 669   fProfPtSigmaSnp2 = new TProfile("fProfPtSigmaSnp2","fProfPtSigmaSnp2",fgkN1PtBins,bins1Pt);
 670   fHistList->Add(fProfPtSigmaSnp2);
 671
 672   fProfPtSigmaTgl2 = new TProfile("fProfPtSigmaTgl2","fProfPtSigmaTgl2",fgkN1PtBins,bins1Pt);
 673   fHistList->Add(fProfPtSigmaTgl2);
 674
 675   fProfPtSigma1Pt2 = new TProfile("fProfPtSigma1Pt2","fProfPtSigma1Pt2",fgkN1PtBins,bins1Pt);
 676   fHistList->Add(fProfPtSigma1Pt2);
 677
 678   fProfPtSigma1Pt = new TProfile("fProfPtSigma1Pt","fProfPtSigma1Pt;p_{T};#sigma(1/p_{T})",fgkNPtBins,binsPt);
 679   fHistList->Add(fProfPtSigma1Pt);
 680
 681   fProfPtPtSigma1Pt = new TProfile("fProfPtPtSigma1Pt","fProfPtPtSigma1Pt;p_{T};p_{T}#sigma(1/p_{T})",fgkNPtBins,binsPt);
 682   fHistList->Add(fProfPtPtSigma1Pt);
 683
 684   TH1::AddDirectory(oldStatus);
 685
 686   PostData(1, fHistList);
 687
 688   if(binsPhi)               delete [] binsPhi;
 689   if(binsPt)                delete [] binsPt;
 690   if(binsNClustersTPC)      delete [] binsNClustersTPC;
 691   if(binsDCA2D)             delete [] binsDCA2D;
 692   if(binsDCAZ)              delete [] binsDCAZ;
 693   if(binsNPointITS)         delete [] binsNPointITS;
 694   if(binsNSigmaToVertex)    delete [] binsNSigmaToVertex;
 695   if(binsChi2C)             delete [] binsChi2C;
 696   if(binsEta)               delete [] binsEta;
 697   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
 698   if(binsUncertainty1Pt)    delete [] binsUncertainty1Pt;
 699   if(binsChi2PerClus)       delete [] binsChi2PerClus;
 700   if(binsChi2PerClus)       delete [] binsNCrossedRowsNClusF;
 701   if(bins1Pt)               delete [] bins1Pt;
 702   if(binsSigmaY2)           delete [] binsSigmaY2;
 703   if(binsSigmaZ2)           delete [] binsSigmaZ2;
 704   if(binsSigmaSnp2)         delete [] binsSigmaSnp2;
 705   if(binsSigmaTgl2)         delete [] binsSigmaTgl2;
 706   if(binsSigma1Pt2)         delete [] binsSigma1Pt2;
 707 }
 708
 709 //________________________________________________________________________
 710 Bool_t AliPWG4HighPtTrackQA::SelectEvent() {
 711   //
 712   // Decide if event should be selected for analysis
 713   //
 714
 715   // Checks following requirements:
 716   // - fEvent available
 717   // - trigger info from AliPhysicsSelection
 718   // - MCevent available
 719   // - number of reconstructed tracks > 1
 720   // - primary vertex reconstructed
 721   // - z-vertex < 10 cm
 722   // - centrality in case of PbPb
 723
 724   Bool_t selectEvent = kTRUE;
 725
 726   //fEvent object available?
 727   if (!fEvent) {
 728     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
 729     fNEventReject->Fill("noAliVEvent",1);
 730     selectEvent = kFALSE;
 731     return selectEvent;
 732   }
 733
 734   //Check if number of reconstructed tracks is larger than 1
 735   if(!fEvent->GetNumberOfTracks() || fEvent->GetNumberOfTracks()<2)  {
 736     fNEventReject->Fill("NTracks<2",1);
 737     selectEvent = kFALSE;
 738     return selectEvent;
 739   }
 740
 741   //Check if vertex is reconstructed
 742   if(fDataType==kESD&&dynamic_cast<AliESDEvent*>(fEvent)) {
 743     fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexTracks();
 744
 745     if (!fVtx || !fVtx->GetStatus())
 746       fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexSPD();
 747
 748     if(!fVtx) {
 749       fNEventReject->Fill("noVTX",1);
 750       selectEvent = kFALSE;
 751       return selectEvent;
 752     }
 753
 754     if(!fVtx->GetStatus()) {
 755       fNEventReject->Fill("VtxStatus",1);
 756       selectEvent = kFALSE;
 757       return selectEvent;
 758     }
 759
 760     // Need vertex cut
 761     if(fVtx->GetNContributors()<2) {
 762       fNEventReject->Fill("NCont<2",1);
 763       selectEvent = kFALSE;
 764       return selectEvent;
 765     }
 766
 767     //Check if z-vertex < 10 cm
 768     double primVtx[3];
 769     fVtx->GetXYZ(primVtx);
 770     if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 771       fNEventReject->Fill("ZVTX>10",1);
 772       selectEvent = kFALSE;
 773       return selectEvent;
 774     }
 775   }
 776   else if(fDataType==kAOD&&dynamic_cast<AliAODEvent*>(fEvent)) {
 777     const AliAODVertex *vtx = ((AliAODEvent*)fEvent)->GetPrimaryVertexSPD();
 778     if(!vtx) {
 779       fNEventReject->Fill("noVTX",1);
 780       selectEvent = kFALSE;
 781       return selectEvent;
 782     }
 783
 784     // Need vertex cut
 785     if(vtx->GetNContributors()<2) {
 786       fNEventReject->Fill("NCont<2",1);
 787       selectEvent = kFALSE;
 788       return selectEvent;
 789     }
 790
 791     //Check if z-vertex < 10 cm
 792     double primVtx[3];
 793     vtx->GetXYZ(primVtx);
 794     if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 795       fNEventReject->Fill("ZVTX>10",1);
 796       selectEvent = kFALSE;
 797       return selectEvent;
 798     }
 799
 800   }
 801
 802   //Centrality selection should only be done in case of PbPb
 803   if(IsPbPb()) {
 804     Float_t cent = 0.;
 805     if(fCentClass!=CalculateCentrality(fEvent) && fCentClass!=10) {
 806       fNEventReject->Fill("cent",1);
 807       selectEvent = kFALSE;
 808       return selectEvent;
 809     }
 810     else {
 811       if(fDataType==kESD) {
 812         if(dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()) {
 813           cent = dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()->GetCentralityPercentile("V0M");
 814         }
 815       }
 816       else if(fDataType==kAOD) {
 817         if(dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality())
 818           cent = dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality();
 819        }
 820       if(cent>90.) {
 821         fNEventReject->Fill("cent>90",1);
 822         selectEvent = kFALSE;
 823         return selectEvent;
 824       }
 825       fh1Centrality->Fill(cent);
 826     }
 827   }
 828
 829   return selectEvent;
 830
 831 }
 832
 833 //________________________________________________________________________
 834 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliVEvent *ev){
 835   //
 836   // Get centrality from ESD or AOD
 837   //
 838
 839   if(fDataType==kESD)
 840     return CalculateCentrality(dynamic_cast<AliESDEvent*>(ev));
 841   else if(fDataType==kAOD)
 842     return CalculateCentrality(dynamic_cast<AliAODEvent*>(ev));
 843   else
 844     return 5;
 845 }
 846
 847 //________________________________________________________________________
 848 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliESDEvent *esd){
 849   //
 850   // Get centrality from ESD
 851   //
 852
 853   Float_t cent = -1;
 854
 855   if(esd){
 856     if(esd->GetCentrality()){
 857       cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
 858       if(fDebug>3) printf("centrality: %f\n",cent);
 859     }
 860   }
 861
 862   return GetCentralityClass(cent);
 863
 864 }
 865
 866 //________________________________________________________________________
 867 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(const AliAODEvent *aod){
 868   //
 869   // Get centrality from AOD
 870   //
 871
 872   if(!aod) return 5;
 873   Float_t cent = aod->GetHeader()->GetCentrality();
 874   if(fDebug>3) printf("centrality: %f\n",cent);
 875
 876   return GetCentralityClass(cent);
 877
 878 }
 879
 880 //________________________________________________________________________
 881 Int_t AliPWG4HighPtTrackQA::GetCentralityClass(Float_t cent) const {
 882   //
 883   // Get centrality class
 884   //
 885
 886   if(cent<0)  return 5; // OB - cent sometimes negative
 887   if(cent>80) return 4;
 888   if(cent>50) return 3;
 889   if(cent>30) return 2;
 890   if(cent>10) return 1;
 891   return 0;
 892
 893 }
 894
 895 //________________________________________________________________________
 896 void AliPWG4HighPtTrackQA::UserExec(Option_t *) {
 897   // Main loop
 898   // Called for each event
 899   AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserExec \n"));
 900
 901   fEvent = InputEvent();
 902   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 903
 904   // All events without selection
 905   fNEventAll->Fill(0.);
 906
 907   if(!SelectEvent()) {
 908     // Post output data
 909     PostData(1, fHistList);
 910     return;
 911   }
 912
 913
 914   //Need to keep track of selected events
 915   fNEventSel->Fill(0.);
 916
 917   fVariables = new TArrayF(fNVariables);
 918
 919   if(fDataType==kESD) DoAnalysisESD();
 920   if(fDataType==kAOD) DoAnalysisAOD();
 921
 922   //Delete old fVariables
 923   if(fVariables) delete fVariables;
 924
 925   // Post output data
 926   PostData(1, fHistList);
 927
 928 }
 929
 930 //________________________________________________________________________
 931 void AliPWG4HighPtTrackQA::DoAnalysisESD() {
 932   //
 933   // Run analysis on ESD
 934   //
 935
 936   if(!fESD) {
 937     PostData(1, fHistList);
 938     return;
 939   }
 940
 941   // ---- Get MC Header information (for MC productions in pThard bins) ----
 942   Double_t ptHard = 0.;
 943   Double_t nTrials = 1; // trials for MC trigger weight for real data
 944
 945   AliMCEventHandler *eventHandlerMC = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 946   if (eventHandlerMC) {
 947
 948     if(eventHandlerMC->MCEvent()){
 949       AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(eventHandlerMC->MCEvent());
 950       if(pythiaGenHeader){
 951         nTrials = pythiaGenHeader->Trials();
 952         ptHard  = pythiaGenHeader->GetPtHard();
 953
 954         fh1PtHard->Fill(ptHard);
 955         fh1PtHardTrials->Fill(ptHard,nTrials);
 956
 957         fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 958       }
 959     }
 960   }
 961
 962   Int_t nTracks = fESD->GetNumberOfTracks();
 963   AliDebug(2,Form("nTracks ESD%d", nTracks));
 964
 965   /*
 966     Variables to be put in fVariables
 967     0: pt
 968     1: phi
 969     2: eta
 970     3: dca2D
 971     4: dcaZ
 972     5: nClustersTPC
 973     6: nPointITS
 974     7: chi2C
 975     8: nSigmaToVertex
 976     9: trackLengthTPC
 977     10: chi2PerClusterTPC
 978     11: #crossed rows
 979     12: (#crossed rows)/(#findable clusters)
 980     13: SigmaY2
 981     14: SigmaZ2
 982     15: SigmaSnp2
 983     16: SigmaTgl2
 984     17: Sigma1Pt2
 985     18: NClustersTPCIter1
 986     19: Chi2TPCIter1
 987     20: nClustersTPCShared
 988     21: Golden Chi2 - global vs TPC constrained
 989     22: Chi2 between global and global constrained
 990   */
 991
 992   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 993     fh1NTracksAll->Fill(0.);
 994
 995     //Get track for analysis
 996     AliESDtrack *track = 0x0;
 997     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
 998     if(!esdtrack) {
 999       fh1NTracksReject->Fill("noESDtrack",1);
1000       continue;
1001     }
1002     AliESDtrack *origtrack = new AliESDtrack(*esdtrack);
1003     if(!origtrack)
1004       continue;
1005
1006     if(fTrackType==4) {
1007       if (!(fTrackCuts->AcceptTrack(esdtrack))) {
1008         fh1NTracksReject->Fill("trackCuts",1);
1009         if(origtrack) delete origtrack;
1010         continue;
1011       }
1012     }
1013
1014     if(fTrackType==1)
1015       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1016     else if(fTrackType==2 || fTrackType==4) {
1017       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1018       if(!track) {
1019         fh1NTracksReject->Fill("noTPConly",1);
1020         if(origtrack) delete origtrack;
1021         continue;
1022       }
1023       AliExternalTrackParam exParam;
1024       Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1025       if( !relate ) {
1026         fh1NTracksReject->Fill("relate",1);
1027         if(track) delete track;
1028         if(origtrack) delete origtrack;
1029         continue;
1030       }
1031       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1032     }
1033     else if(fTrackType==5 || fTrackType==6) {
1034       if(fTrackCuts->AcceptTrack(esdtrack)) {
1035         if(origtrack) delete origtrack;
1036         continue;
1037       }
1038       else {
1039         if( !(fTrackCutsITSLoose->AcceptTrack(esdtrack)) && fTrackCutsTPConly->AcceptTrack(esdtrack) ) {
1040
1041           if(fTrackType==5) {
1042             //use TPConly constrained track
1043             track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1044             if(!track) {
1045               fh1NTracksReject->Fill("noTPConly",1);
1046               if(origtrack) delete origtrack;
1047               continue;
1048             }
1049             AliExternalTrackParam exParam;
1050             Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1051             if( !relate ) {
1052               fh1NTracksReject->Fill("relate",1);
1053               if(track) delete track;
1054               if(origtrack) delete origtrack;
1055               continue;
1056             }
1057             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1058           }
1059           else if(fTrackType==6) {
1060             //use global constrained track
1061             track = new AliESDtrack(*esdtrack);
1062             track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
1063
1064           }
1065         }
1066       }
1067     }
1068     else if(fTrackType==7) {
1069       //use global constrained track
1070       track = new AliESDtrack(*esdtrack);
1071     }
1072     else
1073       track = esdtrack;
1074
1075     if(!track) {
1076       if(origtrack) delete origtrack;
1077       continue;
1078     }
1079
1080     if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
1081       //Cut on chi2 of constrained fit
1082       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
1083         fh1NTracksReject->Fill("chi2",1);
1084         if(track) delete track;
1085         if(origtrack) delete origtrack;
1086         continue;
1087       }
1088     }
1089
1090     fPtAll->Fill(track->Pt());
1091
1092     if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
1093       fh1NTracksReject->Fill("trackCuts",1);
1094       if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
1095         if(track) delete track;
1096       }
1097       if(origtrack) delete origtrack;
1098       continue;
1099     }
1100
1101     if(fTrackType==7) {
1102       if(fTrackCutsITSLoose ) {
1103         if(fTrackCutsITSLoose->AcceptTrack(track) ) {
1104           if(track) delete track;
1105           if(origtrack) delete origtrack;
1106           continue;
1107         }
1108       }
1109
1110       if(esdtrack->GetConstrainedParam())
1111         track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
1112     }
1113
1114     if(!track) {
1115       if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
1116         if(track) delete track;
1117       }
1118       if(origtrack) delete origtrack;
1119       continue;
1120     }
1121
1122     fh1NTracksSel->Fill(0.);
1123
1124     fVariables->Reset(0.);
1125
1126     fVariables->SetAt(track->Pt(),0);
1127     fVariables->SetAt(track->Phi(),1);
1128     fVariables->SetAt(track->Eta(),2);
1129
1130     Float_t dca2D = 0.;
1131     Float_t dcaz  = 0.;
1132
1133     if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
1134       track->GetImpactParametersTPC(dca2D,dcaz); //TPConly
1135     }
1136     else
1137       track->GetImpactParameters(dca2D,dcaz);    //Global
1138
1139     fVariables->SetAt(dca2D,3);
1140     fVariables->SetAt(dcaz,4);
1141
1142     fVariables->SetAt((float)track->GetTPCNcls(),5);
1143
1144     Int_t nPointITS = 0;
1145     fITSClusterMap = track->GetITSClusterMap();
1146     UChar_t itsMap = track->GetITSClusterMap();
1147     for (Int_t i=0; i < 6; i++) {
1148       if (itsMap & (1 << i))
1149         nPointITS ++;
1150     }
1151     fVariables->SetAt((float)nPointITS,6);
1152     Float_t chi2C = (float)track->GetConstrainedChi2();
1153     if(fTrackType==1 || fTrackType==2 || fTrackType==4)
1154       chi2C = (float)track->GetConstrainedChi2TPC();
1155     fVariables->SetAt(chi2C,7);
1156     fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
1157
1158     fVariables->SetAt(GetTrackLengthTPC(track),9);
1159
1160     if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
1161
1162     //fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
1163     fVariables->SetAt(track->GetTPCCrossedRows(),11); //#crossed rows
1164
1165     Float_t crossedRowsTPCNClsF = 1.;//track->GetTPCClusterInfo(2,0);
1166     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
1167     fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
1168     fVariables->SetAt(track->GetSigmaY2(),13);
1169     fVariables->SetAt(track->GetSigmaZ2(),14);
1170     fVariables->SetAt(track->GetSigmaSnp2(),15);
1171     fVariables->SetAt(track->GetSigmaTgl2(),16);
1172     fVariables->SetAt(track->GetSigma1Pt2(),17);
1173
1174     fVariables->SetAt(track->GetTPCNclsIter1(),18);
1175     fVariables->SetAt(track->GetTPCchi2Iter1(),19);
1176
1177     fVariables->SetAt(track->GetTPCnclsS(),20);
1178
1179     Float_t chi2Gold = origtrack->GetChi2TPCConstrainedVsGlobal(fVtx);//GetGoldenChi2(origtrack);
1180     Float_t chi2GGC  = GetGGCChi2(origtrack);
1181
1182     fVariables->SetAt(chi2Gold,21);
1183     fVariables->SetAt(chi2GGC,22);
1184
1185     fVariables->SetAt(GetTPCClusterInfoFitMap(track,2,1),23);
1186     Float_t crossedRowsTPCNClsFFit = 1.;
1187     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsFFit = fVariables->At(23)/track->GetTPCNclsF();
1188     fVariables->SetAt(crossedRowsTPCNClsFFit,24);
1189
1190     FillHistograms();
1191
1192     //      int mult = fTrackCuts->CountAcceptedTracks(fESD);
1193
1194     if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
1195       if(track) delete track;
1196     }
1197     if(origtrack) delete origtrack;
1198
1199   }//track loop
1200
1201 }
1202
1203 //________________________________________________________________________
1204 void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
1205   //
1206   // Do QA on AOD input
1207   //
1208   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1209   if(!aod)return;
1210   AliExternalTrackParam *exParam = new  AliExternalTrackParam();
1211   for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
1212
1213     AliAODTrack *aodtrack = aod->GetTrack(iTrack);
1214     if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
1215
1216     fVariables->Reset(0.);
1217
1218     fVariables->SetAt(aodtrack->Pt(),0);
1219     fVariables->SetAt(aodtrack->Phi(),1);
1220     fVariables->SetAt(aodtrack->Eta(),2);
1221
1222     Double_t dca[2] = {1e6,1e6};
1223     Double_t covar[3] = {1e6,1e6,1e6};
1224     if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
1225       fVariables->SetAt(dca[0],3);
1226       fVariables->SetAt(dca[1],4);
1227     }
1228
1229     fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
1230     fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
1231     fVariables->SetAt(aodtrack->Chi2perNDF(),7);
1232     fVariables->SetAt(0.,8);
1233     fVariables->SetAt(GetTrackLengthTPC(aodtrack),9);
1234     fVariables->SetAt(aodtrack->Chi2perNDF(),10);
1235     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2,1),11);
1236     Float_t crossedRowsTPCNClsF = 0.;
1237     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
1238     fVariables->SetAt(crossedRowsTPCNClsF,12);
1239
1240     //get covariance matrix
1241     Double_t cov[21] = {0,};
1242     aodtrack->GetCovMatrix(cov);
1243     Double_t pxpypz[3] = {0,};
1244     aodtrack->PxPyPz(pxpypz);
1245     Double_t xyz[3] = {0,};
1246     aodtrack->GetXYZ(xyz);
1247     Short_t sign = aodtrack->Charge();
1248     exParam->Set(xyz,pxpypz,cov,sign);
1249
1250     fVariables->SetAt(exParam->GetSigmaY2(),13);
1251     fVariables->SetAt(exParam->GetSigmaZ2(),14);
1252     fVariables->SetAt(exParam->GetSigmaSnp2(),15);
1253     fVariables->SetAt(exParam->GetSigmaTgl2(),16);
1254     fVariables->SetAt(exParam->GetSigma1Pt2(),17);
1255
1256     fVariables->SetAt(0.,18);
1257     fVariables->SetAt(0.,19);
1258
1259     TBits sharedClusterMap = aodtrack->GetTPCSharedMap();
1260     fVariables->SetAt(sharedClusterMap.CountBits(),20);
1261
1262     fVariables->SetAt(0.,21); //not available in AOD
1263     fVariables->SetAt(0.,22); //not available in AOD
1264
1265     fVariables->SetAt(0.,23); //not available in AOD
1266     fVariables->SetAt(0.,24); //not available in AOD
1267
1268     fPtAll->Fill(fVariables->At(0));
1269
1270     FillHistograms();
1271
1272   }
1273
1274 }
1275
1276 //________________________________________________________________________
1277 void AliPWG4HighPtTrackQA::FillHistograms() {
1278   //
1279   // Fill all QA histograms
1280   //
1281
1282   fPtSel->Fill(fVariables->At(0));
1283   fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1284   fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1285   fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1286   fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1287   fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1288   fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1289
1290
1291   fPtNClustersTPCIter1->Fill(fVariables->At(0),fVariables->At(18));
1292   fPtNClustersTPCIter1Phi->Fill(fVariables->At(0),fVariables->At(18),fVariables->At(1));
1293   fPtNClustersTPCShared->Fill(fVariables->At(0),fVariables->At(20));
1294   if(fVariables->At(5)>0.)
1295     fPtNClustersTPCSharedFrac->Fill(fVariables->At(0),fVariables->At(20)/fVariables->At(5));
1296
1297   if(fVariables->At(18)>0.)
1298     fPtChi2PerClusterTPCIter1->Fill(fVariables->At(0),fVariables->At(19)/fVariables->At(18));
1299
1300   fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1301   fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1302   fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1303   fPtRelUncertainty1PtNClus->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(5));
1304   fPtRelUncertainty1PtNClusIter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(18));
1305   fPtRelUncertainty1PtNPointITS->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(6));
1306
1307   fPtRelUncertainty1PtITSClusterMap->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),(int)fITSClusterMap);
1308
1309   fPtRelUncertainty1PtChi2->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(10));
1310   if(fVariables->At(18)>0.)
1311     fPtRelUncertainty1PtChi2Iter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(19)/fVariables->At(18));
1312   fPtRelUncertainty1PtPhi->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(1));
1313
1314   fPtUncertainty1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1315   fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1316   fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1317   fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1318   fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1319   fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1320
1321   fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1322   fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1323   fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1324   fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1325   fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1326   fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1327   fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1328
1329   fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1330   fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1331   fPtNCrossedRowsPhi->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(1));
1332   fPtNCrossedRowsNClusFPhi->Fill(fVariables->At(0),fVariables->At(12),fVariables->At(1));
1333   fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1334
1335   fPtChi2Gold->Fill(fVariables->At(0),fVariables->At(21));
1336   fPtChi2GGC->Fill(fVariables->At(0),fVariables->At(22));
1337
1338   fPtChi2GoldPhi->Fill(fVariables->At(0),fVariables->At(21),fVariables->At(1));
1339   fPtChi2GGCPhi->Fill(fVariables->At(0),fVariables->At(22),fVariables->At(1));
1340
1341   fChi2GoldChi2GGC->Fill(fVariables->At(21),fVariables->At(22));
1342
1343   fPtNCrossedRowsFit->Fill(fVariables->At(0),fVariables->At(23));
1344   fPtNCrossedRowsFitPhi->Fill(fVariables->At(0),fVariables->At(23),fVariables->At(1));
1345   fPtNCrossedRowsNClusFFitPhi->Fill(fVariables->At(0),fVariables->At(24),fVariables->At(1));
1346   fNCrossedRowsNCrossedRowsFit->Fill(fVariables->At(11),fVariables->At(23));
1347
1348   fNClustersNCrossedRows->Fill(fVariables->At(5),fVariables->At(11));
1349   fNClustersNCrossedRowsFit->Fill(fVariables->At(5),fVariables->At(23));
1350
1351   fPtRelUncertainty1PtNCrossedRows->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(11));
1352   fPtRelUncertainty1PtNCrossedRowsFit->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(23));
1353
1354 }
1355
1356 //________________________________________________________________________
1357 Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1358   //
1359   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1360   // This is to called in Notify and should provide the path to the AOD/ESD file
1361   // Copied from AliAnalysisTaskJetSpectrum2
1362   //
1363
1364   TString file(currFile);
1365   fXsec = 0;
1366   fTrials = 1;
1367
1368   if(file.Contains("root_archive.zip#")){
1369     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1370     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1371     file.Replace(pos+1,20,"");
1372   }
1373   else {
1374     // not an archive take the basename....
1375     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1376   }
1377   //  Printf("%s",file.Data());
1378
1379
1380   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
1381   if(!fxsec){
1382     // next trial fetch the histgram file
1383     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1384     if(!fxsec){
1385         // not a severe condition but inciate that we have no information
1386       return kFALSE;
1387     }
1388     else{
1389       // find the tlist we want to be independtent of the name so use the Tkey
1390       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1391       if(!key){
1392         fxsec->Close();
1393         return kFALSE;
1394       }
1395       TList *list = dynamic_cast<TList*>(key->ReadObj());
1396       if(!list){
1397         fxsec->Close();
1398         return kFALSE;
1399       }
1400       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1401       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1402       fxsec->Close();
1403     }
1404   } // no tree pyxsec.root
1405   else {
1406     TTree *xtree = (TTree*)fxsec->Get("Xsection");
1407     if(!xtree){
1408       fxsec->Close();
1409       return kFALSE;
1410     }
1411     UInt_t   ntrials  = 0;
1412     Double_t  xsection  = 0;
1413     xtree->SetBranchAddress("xsection",&xsection);
1414     xtree->SetBranchAddress("ntrials",&ntrials);
1415     xtree->GetEntry(0);
1416     fTrials = ntrials;
1417     fXsec = xsection;
1418     fxsec->Close();
1419   }
1420   return kTRUE;
1421 }
1422 //________________________________________________________________________
1423 Bool_t AliPWG4HighPtTrackQA::Notify()
1424 {
1425   //
1426   // Implemented Notify() to read the cross sections
1427   // and number of trials from pyxsec.root
1428   // Copied from AliAnalysisTaskJetSpectrum2
1429   //
1430
1431   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1432   Float_t xsection = 0;
1433   Float_t ftrials  = 1;
1434
1435   fAvgTrials = 1;
1436   if(tree){
1437     TFile *curfile = tree->GetCurrentFile();
1438     if (!curfile) {
1439       Error("Notify","No current file");
1440       return kFALSE;
1441     }
1442     if(!fh1Xsec||!fh1Trials){
1443       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1444       return kFALSE;
1445     }
1446     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1447     fh1Xsec->Fill("<#sigma>",xsection);
1448     // construct a poor man average trials
1449     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1450     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1451   }
1452   return kTRUE;
1453 }
1454
1455 //________________________________________________________________________
1456 AliGenPythiaEventHeader*  AliPWG4HighPtTrackQA::GetPythiaEventHeader(const AliMCEvent *mcEvent){
1457
1458   if(!mcEvent)return 0;
1459   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1460   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1461   if(!pythiaGenHeader){
1462     // cocktail ??
1463     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1464
1465     if (!genCocktailHeader) {
1466       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1467       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1468       return 0;
1469     }
1470     TList* headerList = genCocktailHeader->GetHeaders();
1471     for (Int_t i=0; i<headerList->GetEntries(); i++) {
1472       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1473       if (pythiaGenHeader)
1474         break;
1475     }
1476     if(!pythiaGenHeader){
1477       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1478       return 0;
1479     }
1480   }
1481   return pythiaGenHeader;
1482
1483 }
1484
1485 //_______________________________________________________________________
1486 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(const AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1487 {
1488   //MV: copied from AliESDtrack since method is not available in AliAODTrack
1489
1490   //
1491   // TPC cluster information
1492   // type 0: get fraction of found/findable clusters with neighbourhood definition
1493   //      1: findable clusters with neighbourhood definition
1494   //      2: found clusters
1495   //
1496   // definition of findable clusters:
1497   //            a cluster is defined as findable if there is another cluster
1498   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1499   //           effects with a very simple algorithm.
1500   //
1501
1502   TBits fTPCClusterMap = tr->GetTPCClusterMap();
1503   if (type==2) return fTPCClusterMap.CountBits();
1504
1505   Int_t found=0;
1506   Int_t findable=0;
1507   Int_t last=-nNeighbours;
1508
1509   for (Int_t i=row0; i<row1; ++i){
1510     //look to current row
1511     if (fTPCClusterMap[i]) {
1512       last=i;
1513       ++found;
1514       ++findable;
1515       continue;
1516     }
1517     //look to nNeighbours before
1518     if ((i-last)<=nNeighbours) {
1519       ++findable;
1520       continue;
1521     }
1522     //look to nNeighbours after
1523     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1524       if (fTPCClusterMap[j]){
1525         ++findable;
1526         break;
1527       }
1528     }
1529   }
1530   if (type==1) return findable;
1531
1532   if (type==0){
1533     Float_t fraction=0;
1534     if (findable>0)
1535       fraction=(Float_t)found/(Float_t)findable;
1536     else
1537       fraction=0;
1538     return fraction;
1539   }
1540   return 0;  // undefined type - default value
1541 }
1542
1543 //_______________________________________________________________________
1544 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfoFitMap(const AliESDtrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1545 {
1546   //
1547   // TPC cluster information from fit map
1548   // type 0: get fraction of found/findable clusters with neighbourhood definition
1549   //      1: findable clusters with neighbourhood definition
1550   //      2: found clusters
1551   //
1552   // definition of findable clusters:
1553   //            a cluster is defined as findable if there is another cluster
1554   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1555   //           effects with a very simple algorithm.
1556   //
1557
1558   TBits fTPCFitMap = tr->GetTPCFitMap();
1559   if (type==2) return fTPCFitMap.CountBits();
1560
1561   Int_t found=0;
1562   Int_t findable=0;
1563   Int_t last=-nNeighbours;
1564
1565   for (Int_t i=row0; i<row1; ++i){
1566     //look to current row
1567     if (fTPCFitMap[i]) {
1568       last=i;
1569       ++found;
1570       ++findable;
1571       continue;
1572     }
1573     //look to nNeighbours before
1574     if ((i-last)<=nNeighbours) {
1575       ++findable;
1576       continue;
1577     }
1578     //look to nNeighbours after
1579     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1580       if (fTPCFitMap[j]){
1581         ++findable;
1582         break;
1583       }
1584     }
1585   }
1586   if (type==1) return findable;
1587
1588   if (type==0){
1589     Float_t fraction=0;
1590     if (findable>0)
1591       fraction=(Float_t)found/(Float_t)findable;
1592     else
1593       fraction=0;
1594     return fraction;
1595   }
1596   return 0;  // undefined type - default value
1597 }
1598
1599 //_______________________________________________________________________
1600 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(const AliESDtrack *track) const {
1601   //
1602   // returns distance between 1st and last hit in TPC
1603   // distance given in number of padrows
1604   //
1605
1606   TBits fTPCClusterMap = track->GetTPCClusterMap();
1607   int firstHit = 0;
1608   int lastHit = 0;
1609
1610   for(int i=0; i<=159; i++) {
1611     if(fTPCClusterMap[i]>0) firstHit = i;
1612   }
1613   for(int i=159; i>=0; i--) {
1614     if(fTPCClusterMap[i]>0) lastHit = i;
1615   }
1616
1617   Int_t trackLength = lastHit - firstHit;
1618
1619   return trackLength;
1620 }
1621
1622 //_______________________________________________________________________
1623 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(const AliAODTrack *track) const {
1624   //
1625   // returns distance between 1st and last hit in TPC
1626   // distance given in number of padrows
1627   //
1628
1629   TBits fTPCClusterMap = track->GetTPCClusterMap();
1630   int firstHit = 0;
1631   int lastHit = 0;
1632
1633   for(int i=0; i<=159; i++) {
1634     if(fTPCClusterMap[i]>0) firstHit = i;
1635   }
1636   for(int i=159; i>=0; i--) {
1637     if(fTPCClusterMap[i]>0) lastHit = i;
1638   }
1639
1640   Int_t trackLength = lastHit - firstHit;
1641
1642   return trackLength;
1643 }
1644
1645 //_______________________________________________________________________
1646 Float_t AliPWG4HighPtTrackQA::GetGoldenChi2(AliESDtrack *origtrack) {
1647   //
1648   // Return chi2 between global and TPC constrained track
1649   // track should be the global unconstrained track
1650   //
1651
1652   Float_t chi2Gold = 0.;
1653
1654   AliESDtrack *tpcTrack = 0x0;
1655   tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,origtrack->GetID());
1656   if(tpcTrack) {
1657     AliExternalTrackParam exParam;
1658     Bool_t relate = tpcTrack->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1659     if( relate ) {
1660       tpcTrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1661       //          Double_t pTPC[2],covTPC[3];     tpcTrack->PropagateToDCA(fVtx, fESD->GetMagneticField(), 10000,  pTPC, covTPC);
1662     }
1663
1664     tpcTrack->Propagate(origtrack->GetAlpha(), origtrack->GetX(), fESD->GetMagneticField());
1665     chi2Gold = (Float_t)origtrack->GetPredictedChi2(tpcTrack);
1666   }
1667
1668   if(tpcTrack) delete tpcTrack;
1669
1670   return chi2Gold;
1671
1672 }
1673
1674 //_______________________________________________________________________
1675 Float_t AliPWG4HighPtTrackQA::GetGGCChi2(AliESDtrack *origtrack) {
1676   //
1677   // Return chi2 between global and global constrained track
1678   // track should be the global unconstrained track
1679   //
1680
1681   Float_t chi2GGC = 0.;
1682
1683   AliESDtrack *esdtrackC = new AliESDtrack(*origtrack);
1684   if(esdtrackC) {
1685     if(origtrack->GetConstrainedParam()) {
1686       esdtrackC->Set(origtrack->GetConstrainedParam()->GetX(),origtrack->GetConstrainedParam()->GetAlpha(),origtrack->GetConstrainedParam()->GetParameter(),origtrack->GetConstrainedParam()->GetCovariance());
1687       chi2GGC = (Float_t)origtrack->GetPredictedChi2(esdtrackC);
1688     }
1689     delete esdtrackC;
1690   }
1691
1692   return chi2GGC;
1693
1694 }
1695
1696 //________________________________________________________________________
1697 void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1698 {
1699   // The Terminate() function is the last function to be called during
1700   // a query. It always runs on the client, it can be used to present
1701   // the results graphically or save the results to file.
1702
1703 }
1704
1705 #endif