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     23: #crossed rows from fit map
 991     24: (#crossed rows)/(#findable clusters) from fit map
 992   */
 993
 994   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 995     fh1NTracksAll->Fill(0.);
 996
 997     //Get track for analysis
 998     AliESDtrack *track = 0x0;
 999     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
1000     if(!esdtrack) {
1001       fh1NTracksReject->Fill("noESDtrack",1);
1002       continue;
1003     }
1004     AliESDtrack *origtrack = new AliESDtrack(*esdtrack);
1005     if(!origtrack)
1006       continue;
1007
1008     if(fTrackType==4) {
1009       if (!(fTrackCuts->AcceptTrack(esdtrack))) {
1010         fh1NTracksReject->Fill("trackCuts",1);
1011         if(origtrack) delete origtrack;
1012         continue;
1013       }
1014     }
1015
1016     if(fTrackType==1)
1017       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1018     else if(fTrackType==2 || fTrackType==4) {
1019       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1020       if(!track) {
1021         fh1NTracksReject->Fill("noTPConly",1);
1022         if(origtrack) delete origtrack;
1023         continue;
1024       }
1025       AliExternalTrackParam exParam;
1026       Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1027       if( !relate ) {
1028         fh1NTracksReject->Fill("relate",1);
1029         if(track) delete track;
1030         if(origtrack) delete origtrack;
1031         continue;
1032       }
1033       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1034     }
1035     else if(fTrackType==5 || fTrackType==6) {
1036       if(fTrackCuts->AcceptTrack(esdtrack)) {
1037         if(origtrack) delete origtrack;
1038         continue;
1039       }
1040       else {
1041         if( !(fTrackCutsITSLoose->AcceptTrack(esdtrack)) && fTrackCutsTPConly->AcceptTrack(esdtrack) ) {
1042
1043           if(fTrackType==5) {
1044             //use TPConly constrained track
1045             track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
1046             if(!track) {
1047               fh1NTracksReject->Fill("noTPConly",1);
1048               if(origtrack) delete origtrack;
1049               continue;
1050             }
1051             AliExternalTrackParam exParam;
1052             Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1053             if( !relate ) {
1054               fh1NTracksReject->Fill("relate",1);
1055               if(track) delete track;
1056               if(origtrack) delete origtrack;
1057               continue;
1058             }
1059             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1060           }
1061           else if(fTrackType==6) {
1062             //use global constrained track
1063             track = new AliESDtrack(*esdtrack);
1064             track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
1065
1066           }
1067         }
1068       }
1069     }
1070     else if(fTrackType==7) {
1071       //use global constrained track
1072       track = new AliESDtrack(*esdtrack);
1073     }
1074     else
1075       track = esdtrack;
1076
1077     if(!track) {
1078       if(origtrack) delete origtrack;
1079       continue;
1080     }
1081
1082     if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
1083       //Cut on chi2 of constrained fit
1084       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
1085         fh1NTracksReject->Fill("chi2",1);
1086         if(track) delete track;
1087         if(origtrack) delete origtrack;
1088         continue;
1089       }
1090     }
1091
1092     fPtAll->Fill(track->Pt());
1093
1094     if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
1095       fh1NTracksReject->Fill("trackCuts",1);
1096       if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
1097         if(track) delete track;
1098       }
1099       if(origtrack) delete origtrack;
1100       continue;
1101     }
1102
1103     if(fTrackType==7) {
1104       if(fTrackCutsITSLoose ) {
1105         if(fTrackCutsITSLoose->AcceptTrack(track) ) {
1106           if(track) delete track;
1107           if(origtrack) delete origtrack;
1108           continue;
1109         }
1110       }
1111
1112       if(esdtrack->GetConstrainedParam())
1113         track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
1114     }
1115
1116     if(!track) {
1117       if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
1118         if(track) delete track;
1119       }
1120       if(origtrack) delete origtrack;
1121       continue;
1122     }
1123
1124     fh1NTracksSel->Fill(0.);
1125
1126     fVariables->Reset(0.);
1127
1128     fVariables->SetAt(track->Pt(),0);
1129     fVariables->SetAt(track->Phi(),1);
1130     fVariables->SetAt(track->Eta(),2);
1131
1132     Float_t dca2D = 0.;
1133     Float_t dcaz  = 0.;
1134
1135     if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
1136       track->GetImpactParametersTPC(dca2D,dcaz); //TPConly
1137     }
1138     else
1139       track->GetImpactParameters(dca2D,dcaz);    //Global
1140
1141     fVariables->SetAt(dca2D,3);
1142     fVariables->SetAt(dcaz,4);
1143
1144     fVariables->SetAt((float)track->GetTPCNcls(),5);
1145
1146     Int_t nPointITS = 0;
1147     fITSClusterMap = track->GetITSClusterMap();
1148     UChar_t itsMap = track->GetITSClusterMap();
1149     for (Int_t i=0; i < 6; i++) {
1150       if (itsMap & (1 << i))
1151         nPointITS ++;
1152     }
1153     fVariables->SetAt((float)nPointITS,6);
1154     Float_t chi2C = (float)track->GetConstrainedChi2();
1155     if(fTrackType==1 || fTrackType==2 || fTrackType==4)
1156       chi2C = (float)track->GetConstrainedChi2TPC();
1157     fVariables->SetAt(chi2C,7);
1158     fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
1159
1160     fVariables->SetAt(GetTrackLengthTPC(track),9);
1161
1162     if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
1163
1164     //fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
1165     fVariables->SetAt(track->GetTPCCrossedRows(),11); //#crossed rows
1166
1167     Float_t crossedRowsTPCNClsF = 1.;//track->GetTPCClusterInfo(2,0);
1168     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
1169     fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
1170     fVariables->SetAt(track->GetSigmaY2(),13);
1171     fVariables->SetAt(track->GetSigmaZ2(),14);
1172     fVariables->SetAt(track->GetSigmaSnp2(),15);
1173     fVariables->SetAt(track->GetSigmaTgl2(),16);
1174     fVariables->SetAt(track->GetSigma1Pt2(),17);
1175
1176     fVariables->SetAt(track->GetTPCNclsIter1(),18);
1177     fVariables->SetAt(track->GetTPCchi2Iter1(),19);
1178
1179     fVariables->SetAt(track->GetTPCnclsS(),20);
1180
1181     Float_t chi2Gold = origtrack->GetChi2TPCConstrainedVsGlobal(fVtx);//GetGoldenChi2(origtrack);
1182     Float_t chi2GGC  = GetGGCChi2(origtrack);
1183
1184     fVariables->SetAt(chi2Gold,21);
1185     fVariables->SetAt(chi2GGC,22);
1186
1187     fVariables->SetAt(GetTPCClusterInfoFitMap(track,2,1),23);
1188     Float_t crossedRowsTPCNClsFFit = 1.;
1189     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsFFit = fVariables->At(23)/track->GetTPCNclsF();
1190     fVariables->SetAt(crossedRowsTPCNClsFFit,24);
1191
1192     FillHistograms();
1193
1194     //      int mult = fTrackCuts->CountAcceptedTracks(fESD);
1195
1196     if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
1197       if(track) delete track;
1198     }
1199     if(origtrack) delete origtrack;
1200
1201   }//track loop
1202
1203 }
1204
1205 //________________________________________________________________________
1206 void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
1207   //
1208   // Do QA on AOD input
1209   //
1210   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1211   if(!aod)return;
1212   AliExternalTrackParam *exParam = new  AliExternalTrackParam();
1213   for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
1214
1215     AliAODTrack *aodtrack = aod->GetTrack(iTrack);
1216     if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
1217
1218     fVariables->Reset(0.);
1219
1220     fVariables->SetAt(aodtrack->Pt(),0);
1221     fVariables->SetAt(aodtrack->Phi(),1);
1222     fVariables->SetAt(aodtrack->Eta(),2);
1223
1224     Double_t dca[2] = {1e6,1e6};
1225     Double_t covar[3] = {1e6,1e6,1e6};
1226     if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
1227       fVariables->SetAt(dca[0],3);
1228       fVariables->SetAt(dca[1],4);
1229     }
1230
1231     fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
1232     fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
1233     fVariables->SetAt(aodtrack->Chi2perNDF(),7);
1234     fVariables->SetAt(0.,8);
1235     fVariables->SetAt(GetTrackLengthTPC(aodtrack),9);
1236     fVariables->SetAt(aodtrack->Chi2perNDF(),10);
1237     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2,1,0,159,kFALSE),11);
1238     Float_t crossedRowsTPCNClsF = 0.;
1239     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
1240     fVariables->SetAt(crossedRowsTPCNClsF,12);
1241
1242     //get covariance matrix
1243     Double_t cov[21] = {0,};
1244     aodtrack->GetCovMatrix(cov);
1245     Double_t pxpypz[3] = {0,};
1246     aodtrack->PxPyPz(pxpypz);
1247     Double_t xyz[3] = {0,};
1248     aodtrack->GetXYZ(xyz);
1249     Short_t sign = aodtrack->Charge();
1250     exParam->Set(xyz,pxpypz,cov,sign);
1251
1252     fVariables->SetAt(exParam->GetSigmaY2(),13);
1253     fVariables->SetAt(exParam->GetSigmaZ2(),14);
1254     fVariables->SetAt(exParam->GetSigmaSnp2(),15);
1255     fVariables->SetAt(exParam->GetSigmaTgl2(),16);
1256     fVariables->SetAt(exParam->GetSigma1Pt2(),17);
1257
1258     fVariables->SetAt(0.,18); //NClustersTPCIter1
1259     fVariables->SetAt(0.,19); //Chi2TPCIter1
1260
1261     TBits sharedClusterMap = aodtrack->GetTPCSharedMap();
1262     fVariables->SetAt(sharedClusterMap.CountBits(),20);
1263
1264     fVariables->SetAt(0.,21); //not available in AOD golden chi2
1265     fVariables->SetAt(0.,22); //not available in AOD  Chi2 between global and global constrained
1266
1267     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2,1,0,159,kTRUE),23); //not available in AOD #crossed rows from fit map
1268     Float_t crossedRowsTPCNClsFFit = 0.;
1269     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsFFit = fVariables->At(23)/aodtrack->GetTPCNclsF();
1270     fVariables->SetAt(crossedRowsTPCNClsFFit,24); //(#crossed rows)/(#findable clusters) from fit map
1271
1272     fPtAll->Fill(fVariables->At(0));
1273
1274     FillHistograms();
1275
1276   }
1277
1278 }
1279
1280 //________________________________________________________________________
1281 void AliPWG4HighPtTrackQA::FillHistograms() {
1282   //
1283   // Fill all QA histograms
1284   //
1285
1286   fPtSel->Fill(fVariables->At(0));
1287   fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1288   fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1289   fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1290   fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1291   fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1292   fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1293
1294
1295   fPtNClustersTPCIter1->Fill(fVariables->At(0),fVariables->At(18));
1296   fPtNClustersTPCIter1Phi->Fill(fVariables->At(0),fVariables->At(18),fVariables->At(1));
1297   fPtNClustersTPCShared->Fill(fVariables->At(0),fVariables->At(20));
1298   if(fVariables->At(5)>0.)
1299     fPtNClustersTPCSharedFrac->Fill(fVariables->At(0),fVariables->At(20)/fVariables->At(5));
1300
1301   if(fVariables->At(18)>0.)
1302     fPtChi2PerClusterTPCIter1->Fill(fVariables->At(0),fVariables->At(19)/fVariables->At(18));
1303
1304   fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1305   fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1306   fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1307   fPtRelUncertainty1PtNClus->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(5));
1308   fPtRelUncertainty1PtNClusIter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(18));
1309   fPtRelUncertainty1PtNPointITS->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(6));
1310
1311   fPtRelUncertainty1PtITSClusterMap->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),(int)fITSClusterMap);
1312
1313   fPtRelUncertainty1PtChi2->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(10));
1314   if(fVariables->At(18)>0.)
1315     fPtRelUncertainty1PtChi2Iter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(19)/fVariables->At(18));
1316   fPtRelUncertainty1PtPhi->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(1));
1317
1318   fPtUncertainty1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1319   fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1320   fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1321   fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1322   fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1323   fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1324
1325   fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1326   fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1327   fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1328   fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1329   fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1330   fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1331   fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1332
1333   fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1334   fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1335   fPtNCrossedRowsPhi->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(1));
1336   fPtNCrossedRowsNClusFPhi->Fill(fVariables->At(0),fVariables->At(12),fVariables->At(1));
1337   fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1338
1339   fPtChi2Gold->Fill(fVariables->At(0),fVariables->At(21));
1340   fPtChi2GGC->Fill(fVariables->At(0),fVariables->At(22));
1341
1342   fPtChi2GoldPhi->Fill(fVariables->At(0),fVariables->At(21),fVariables->At(1));
1343   fPtChi2GGCPhi->Fill(fVariables->At(0),fVariables->At(22),fVariables->At(1));
1344
1345   fChi2GoldChi2GGC->Fill(fVariables->At(21),fVariables->At(22));
1346
1347   fPtNCrossedRowsFit->Fill(fVariables->At(0),fVariables->At(23));
1348   fPtNCrossedRowsFitPhi->Fill(fVariables->At(0),fVariables->At(23),fVariables->At(1));
1349   fPtNCrossedRowsNClusFFitPhi->Fill(fVariables->At(0),fVariables->At(24),fVariables->At(1));
1350   fNCrossedRowsNCrossedRowsFit->Fill(fVariables->At(11),fVariables->At(23));
1351
1352   fNClustersNCrossedRows->Fill(fVariables->At(5),fVariables->At(11));
1353   fNClustersNCrossedRowsFit->Fill(fVariables->At(5),fVariables->At(23));
1354
1355   fPtRelUncertainty1PtNCrossedRows->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(11));
1356   fPtRelUncertainty1PtNCrossedRowsFit->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(23));
1357
1358 }
1359
1360 //________________________________________________________________________
1361 Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1362   //
1363   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1364   // This is to called in Notify and should provide the path to the AOD/ESD file
1365   // Copied from AliAnalysisTaskJetSpectrum2
1366   //
1367
1368   TString file(currFile);
1369   fXsec = 0;
1370   fTrials = 1;
1371
1372   if(file.Contains("root_archive.zip#")){
1373     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1374     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1375     file.Replace(pos+1,20,"");
1376   }
1377   else {
1378     // not an archive take the basename....
1379     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1380   }
1381   //  Printf("%s",file.Data());
1382
1383
1384   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
1385   if(!fxsec){
1386     // next trial fetch the histgram file
1387     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1388     if(!fxsec){
1389         // not a severe condition but inciate that we have no information
1390       return kFALSE;
1391     }
1392     else{
1393       // find the tlist we want to be independtent of the name so use the Tkey
1394       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1395       if(!key){
1396         fxsec->Close();
1397         return kFALSE;
1398       }
1399       TList *list = dynamic_cast<TList*>(key->ReadObj());
1400       if(!list){
1401         fxsec->Close();
1402         return kFALSE;
1403       }
1404       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1405       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1406       fxsec->Close();
1407     }
1408   } // no tree pyxsec.root
1409   else {
1410     TTree *xtree = (TTree*)fxsec->Get("Xsection");
1411     if(!xtree){
1412       fxsec->Close();
1413       return kFALSE;
1414     }
1415     UInt_t   ntrials  = 0;
1416     Double_t  xsection  = 0;
1417     xtree->SetBranchAddress("xsection",&xsection);
1418     xtree->SetBranchAddress("ntrials",&ntrials);
1419     xtree->GetEntry(0);
1420     fTrials = ntrials;
1421     fXsec = xsection;
1422     fxsec->Close();
1423   }
1424   return kTRUE;
1425 }
1426 //________________________________________________________________________
1427 Bool_t AliPWG4HighPtTrackQA::Notify()
1428 {
1429   //
1430   // Implemented Notify() to read the cross sections
1431   // and number of trials from pyxsec.root
1432   // Copied from AliAnalysisTaskJetSpectrum2
1433   //
1434
1435   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1436   Float_t xsection = 0;
1437   Float_t ftrials  = 1;
1438
1439   fAvgTrials = 1;
1440   if(tree){
1441     TFile *curfile = tree->GetCurrentFile();
1442     if (!curfile) {
1443       Error("Notify","No current file");
1444       return kFALSE;
1445     }
1446     if(!fh1Xsec||!fh1Trials){
1447       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1448       return kFALSE;
1449     }
1450     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1451     fh1Xsec->Fill("<#sigma>",xsection);
1452     // construct a poor man average trials
1453     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1454     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1455   }
1456   return kTRUE;
1457 }
1458
1459 //________________________________________________________________________
1460 AliGenPythiaEventHeader*  AliPWG4HighPtTrackQA::GetPythiaEventHeader(const AliMCEvent *mcEvent){
1461
1462   if(!mcEvent)return 0;
1463   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1464   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1465   if(!pythiaGenHeader){
1466     // cocktail ??
1467     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1468
1469     if (!genCocktailHeader) {
1470       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1471       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1472       return 0;
1473     }
1474     TList* headerList = genCocktailHeader->GetHeaders();
1475     for (Int_t i=0; i<headerList->GetEntries(); i++) {
1476       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1477       if (pythiaGenHeader)
1478         break;
1479     }
1480     if(!pythiaGenHeader){
1481       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1482       return 0;
1483     }
1484   }
1485   return pythiaGenHeader;
1486
1487 }
1488
1489 //_______________________________________________________________________
1490 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
1491 {
1492   //MV: copied from AliESDtrack since method is not available in AliAODTrack
1493
1494   //
1495   // TPC cluster information
1496   // type 0: get fraction of found/findable clusters with neighbourhood definition
1497   //      1: findable clusters with neighbourhood definition
1498   //      2: found clusters
1499   //
1500   // definition of findable clusters:
1501   //            a cluster is defined as findable if there is another cluster
1502   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1503   //           effects with a very simple algorithm.
1504   //
1505
1506   TBits fTPCClusterMap = 0;
1507   if(useFitMap)
1508     fTPCClusterMap = tr->GetTPCFitMap();
1509   else
1510     fTPCClusterMap = tr->GetTPCClusterMap();
1511
1512   if (type==2) return fTPCClusterMap.CountBits();
1513
1514   Int_t found=0;
1515   Int_t findable=0;
1516   Int_t last=-nNeighbours;
1517
1518   for (Int_t i=row0; i<row1; ++i){
1519     //look to current row
1520     if (fTPCClusterMap[i]) {
1521       last=i;
1522       ++found;
1523       ++findable;
1524       continue;
1525     }
1526     //look to nNeighbours before
1527     if ((i-last)<=nNeighbours) {
1528       ++findable;
1529       continue;
1530     }
1531     //look to nNeighbours after
1532     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1533       if (fTPCClusterMap[j]){
1534         ++findable;
1535         break;
1536       }
1537     }
1538   }
1539   if (type==1) return findable;
1540
1541   if (type==0){
1542     Float_t fraction=0;
1543     if (findable>0)
1544       fraction=(Float_t)found/(Float_t)findable;
1545     else
1546       fraction=0;
1547     return fraction;
1548   }
1549   return 0;  // undefined type - default value
1550 }
1551
1552 //_______________________________________________________________________
1553 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfoFitMap(const AliESDtrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1554 {
1555   //
1556   // TPC cluster information from fit map
1557   // type 0: get fraction of found/findable clusters with neighbourhood definition
1558   //      1: findable clusters with neighbourhood definition
1559   //      2: found clusters
1560   //
1561   // definition of findable clusters:
1562   //            a cluster is defined as findable if there is another cluster
1563   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1564   //           effects with a very simple algorithm.
1565   //
1566
1567   TBits fTPCFitMap = tr->GetTPCFitMap();
1568   if (type==2) return fTPCFitMap.CountBits();
1569
1570   Int_t found=0;
1571   Int_t findable=0;
1572   Int_t last=-nNeighbours;
1573
1574   for (Int_t i=row0; i<row1; ++i){
1575     //look to current row
1576     if (fTPCFitMap[i]) {
1577       last=i;
1578       ++found;
1579       ++findable;
1580       continue;
1581     }
1582     //look to nNeighbours before
1583     if ((i-last)<=nNeighbours) {
1584       ++findable;
1585       continue;
1586     }
1587     //look to nNeighbours after
1588     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1589       if (fTPCFitMap[j]){
1590         ++findable;
1591         break;
1592       }
1593     }
1594   }
1595   if (type==1) return findable;
1596
1597   if (type==0){
1598     Float_t fraction=0;
1599     if (findable>0)
1600       fraction=(Float_t)found/(Float_t)findable;
1601     else
1602       fraction=0;
1603     return fraction;
1604   }
1605   return 0;  // undefined type - default value
1606 }
1607
1608 //_______________________________________________________________________
1609 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(const AliESDtrack *track) const {
1610   //
1611   // returns distance between 1st and last hit in TPC
1612   // distance given in number of padrows
1613   //
1614
1615   TBits fTPCClusterMap = track->GetTPCClusterMap();
1616   int firstHit = 0;
1617   int lastHit = 0;
1618
1619   for(int i=0; i<=159; i++) {
1620     if(fTPCClusterMap[i]>0) firstHit = i;
1621   }
1622   for(int i=159; i>=0; i--) {
1623     if(fTPCClusterMap[i]>0) lastHit = i;
1624   }
1625
1626   Int_t trackLength = lastHit - firstHit;
1627
1628   return trackLength;
1629 }
1630
1631 //_______________________________________________________________________
1632 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(const AliAODTrack *track) const {
1633   //
1634   // returns distance between 1st and last hit in TPC
1635   // distance given in number of padrows
1636   //
1637
1638   TBits fTPCClusterMap = track->GetTPCClusterMap();
1639   int firstHit = 0;
1640   int lastHit = 0;
1641
1642   for(int i=0; i<=159; i++) {
1643     if(fTPCClusterMap[i]>0) firstHit = i;
1644   }
1645   for(int i=159; i>=0; i--) {
1646     if(fTPCClusterMap[i]>0) lastHit = i;
1647   }
1648
1649   Int_t trackLength = lastHit - firstHit;
1650
1651   return trackLength;
1652 }
1653
1654 //_______________________________________________________________________
1655 Float_t AliPWG4HighPtTrackQA::GetGoldenChi2(AliESDtrack *origtrack) {
1656   //
1657   // Return chi2 between global and TPC constrained track
1658   // track should be the global unconstrained track
1659   //
1660
1661   Float_t chi2Gold = 0.;
1662
1663   AliESDtrack *tpcTrack = 0x0;
1664   tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,origtrack->GetID());
1665   if(tpcTrack) {
1666     AliExternalTrackParam exParam;
1667     Bool_t relate = tpcTrack->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
1668     if( relate ) {
1669       tpcTrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1670       //          Double_t pTPC[2],covTPC[3];     tpcTrack->PropagateToDCA(fVtx, fESD->GetMagneticField(), 10000,  pTPC, covTPC);
1671     }
1672
1673     tpcTrack->Propagate(origtrack->GetAlpha(), origtrack->GetX(), fESD->GetMagneticField());
1674     chi2Gold = (Float_t)origtrack->GetPredictedChi2(tpcTrack);
1675   }
1676
1677   if(tpcTrack) delete tpcTrack;
1678
1679   return chi2Gold;
1680
1681 }
1682
1683 //_______________________________________________________________________
1684 Float_t AliPWG4HighPtTrackQA::GetGGCChi2(AliESDtrack *origtrack) {
1685   //
1686   // Return chi2 between global and global constrained track
1687   // track should be the global unconstrained track
1688   //
1689
1690   Float_t chi2GGC = 0.;
1691
1692   AliESDtrack *esdtrackC = new AliESDtrack(*origtrack);
1693   if(esdtrackC) {
1694     if(origtrack->GetConstrainedParam()) {
1695       esdtrackC->Set(origtrack->GetConstrainedParam()->GetX(),origtrack->GetConstrainedParam()->GetAlpha(),origtrack->GetConstrainedParam()->GetParameter(),origtrack->GetConstrainedParam()->GetCovariance());
1696       chi2GGC = (Float_t)origtrack->GetPredictedChi2(esdtrackC);
1697     }
1698     delete esdtrackC;
1699   }
1700
1701   return chi2GGC;
1702
1703 }
1704
1705 //________________________________________________________________________
1706 void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1707 {
1708   // The Terminate() function is the last function to be called during
1709   // a query. It always runs on the client, it can be used to present
1710   // the results graphically or save the results to file.
1711
1712 }
1713
1714 #endif