PWGJE/AliPWG4HighPtTrackQA.cxx

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