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