PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.cxx

   1 // **************************************
   2 // Full pp jet task - ESD input only
   3 // Extract the jet spectrum and all the
   4 // systematic uncertainties
   5 // -R. Ma, Mar 2011
   6 // **************************************
   7
   8 #include <TCanvas.h>
   9 #include <TChain.h>
  10 #include <TFormula.h>
  11 #include <TH1.h>
  12 #include <TH2.h>
  13 #include <TH3.h>
  14 #include <TProfile2D.h>
  15 #include <THnSparse.h>
  16 #include <TROOT.h>
  17 #include <TTree.h>
  18 #include <TArrayI.h>
  19 #include <TClonesArray.h>
  20 #include <TRandom3.h>
  21 #include <TFile.h>
  22 #include <TF1.h>
  23
  24 #include "AliAODEvent.h"
  25 #include "AliAODInputHandler.h"
  26 #include "AliAnalysisManager.h"
  27 #include "AliAnalysisTask.h"
  28 #include "AliCentrality.h"
  29 #include "AliAnalysisTaskFullppJet.h"
  30 #include "AliESDEvent.h"
  31 #include "AliESDInputHandler.h"
  32 #include "AliESDtrackCuts.h"
  33 #include "AliVParticle.h"
  34 #include "AliInputEventHandler.h"
  35 #include "AliEMCALGeometry.h"
  36 #include "AliEMCALRecoUtils.h"
  37 #include "TGeoManager.h"
  38 #include "AliMCEvent.h"
  39 #include "AliMCParticle.h"
  40 #include "AliStack.h"
  41 #include "AliGenEventHeader.h"
  42 #include "AliGenPythiaEventHeader.h"
  43 #include "TGeoGlobalMagField.h"
  44
  45 #include "AliAODJet.h"
  46 #include "AliFJWrapper.h"
  47
  48 #include <iostream>
  49 using std::cout;
  50 using std::cerr;
  51 using std::endl;
  52 using std::vector;
  53
  54 ClassImp(AliAnalysisTaskFullppJet)
  55
  56 const Float_t kRadius[3] = {0.4,0.2,0.3};
  57 const Double_t kPI = TMath::Pi();
  58 const Double_t kdRCut[3] = {0.25,0.1,0.15};
  59
  60 //________________________________________________________________________
  61 AliAnalysisTaskFullppJet::AliAnalysisTaskFullppJet() :
  62   AliAnalysisTaskSE("default"),
  63   fVerbosity(0), fEDSFileCounter(0), fNTracksPerChunk(0), fRejectPileup(kFALSE), fRejectExoticTrigger(kTRUE),
  64   fAnaType(0), fPeriod("lhc11a"), fESD(0), fAOD(0), fMC(0), fStack(0), fTrackArray(0x0), fClusterArray(0x0), fMcPartArray(0x0),
  65   fIsMC(kFALSE), fPhySelForMC(kFALSE), fChargedMC(kFALSE), fXsecScale(0.),
  66   fCentrality(99), fZVtxMax(10),
  67   fTriggerType(-1), fCheckTriggerMask(kTRUE), fIsTPCOnlyVtx(kFALSE),
  68   fIsExoticEvent3GeV(kFALSE), fIsExoticEvent5GeV(kFALSE), fIsEventTriggerBit(kFALSE), fOfflineTrigger(kFALSE), fTriggerMask(0x0),
  69   fGeom(0x0), fRecoUtil(0x0),
  70   fEsdTrackCuts(0x0), fHybridTrackCuts1(0x0), fHybridTrackCuts2(0x0),fTrackCutsType(0), fKinCutType(0), fTrkEtaMax(0.9),
  71   fdEdxMin(73), fdEdxMax(90), fEoverPMin(0.8), fEoverPMax(1.2),
  72   fMatchType(0), fRejectExoticCluster(kTRUE), fRemoveBadChannel(kFALSE), fUseGoodSM(kFALSE),
  73   fStudySubEInHC(kFALSE), fStudyMcOverSubE(kFALSE),
  74   fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0),
  75   fJetNEFMin(0.01), fJetNEFMax(0.98),
  76   fSpotGoodJet(kFALSE), fFindChargedOnlyJet(kFALSE), fFindNeutralOnlyJet(kFALSE),
  77   fCheckTrkEffCorr(kFALSE), fTrkEffCorrCutZ(0.3), fRandomGen(0x0), fRunUE(0), fCheckTPCOnlyVtx(0), fRunSecondaries(0),
  78   fSysJetTrigEff(kFALSE), fVaryJetTrigEff(0),
  79   fSysTrkPtRes(kFALSE), fVaryTrkPtRes(0), fSysTrkEff(kFALSE), fVaryTrkEff(0.),
  80   fSysTrkClsMth(kFALSE), fCutdEta(0.05), fCutdPhi(0.05),
  81   fSysNonLinearity(kFALSE), fNonLinear(0x0), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0),
  82   fSysClusterizer(0),
  83   fNonStdBranch(""), fNonStdFile(""), fAlgorithm("aKt"), fRadius("0.4"), fRecombinationScheme(5),
  84   fConstrainChInEMCal(kFALSE), fRejectNK(kFALSE), fRejectWD(kFALSE), fSmearMC(kFALSE), fTrkPtResData(0x0),
  85   fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventCount(0x0), fhJetEventStat(0x0), fhEventStatTPCVtx(0x0), fhChunkQA(0x0)
  86 {
  87   // Constructor
  88
  89   // Define input and output slots here
  90   // Input slot #0 works with a TChain
  91   DefineInput(0, TChain::Class());
  92   DefineOutput(1, TList::Class());
  93   // Output slot #0 id reserved by the base class for AOD
  94   for(Int_t i=0; i<2; i++)
  95     {
  96       fTrkPtMin[i] = 0.15;
  97       fTrkPtMax[i] = 200;
  98       fClsEtMin[i] = 0.15;
  99       fClsEtMax[i] = 200;
 100
 101       fVertexGenZ[i]       = 0x0;
 102       fEventZ[i]           = 0x0;
 103       fhNTrials[i]         = 0x0;
 104       fhNMatchedTrack[i]   = 0x0;
 105       fhClsE[i]            = 0x0;
 106       for(Int_t k=0; k<2; k++)
 107         {
 108           fhSysClusterE[i][k]     = 0x0;
 109           fhSysNCellVsClsE[i][k]  = 0x0;
 110         }
 111     }
 112
 113   for(Int_t r=0; r<kNRadii; r++)
 114     {
 115       for(Int_t j=0; j<5; j++)
 116         {
 117           fHCOverSubE[r][j]       = 0x0;
 118           fHCOverSubEFrac[r][j]   = 0x0;
 119           for(Int_t k=0; k<2; k++)
 120             {
 121               fhSubClsEVsJetPt[k][r][j] = 0x0;
 122               fhHCTrkPtClean[k][r][j]   = 0x0;
 123               fhHCTrkPtAmbig[k][r][j]   = 0x0;
 124             }
 125
 126           if(j<4) fhSubEVsTrkPt[r][j] = 0x0;
 127
 128           if(j<3)
 129             {
 130               fJetCount[j][r]          = 0x0;
 131               fhNeutralPtInJet[j][r]   = 0x0;
 132               fhTrigNeuPtInJet[j][r]   = 0x0;
 133               fhChargedPtInJet[j][r]   = 0x0;
 134               fhLeadNePtInJet[j][r]    = 0x0;
 135               fhLeadChPtInJet[j][r]    = 0x0;
 136               fJetEnergyFraction[j][r] = 0x0;
 137               fJetNPartFraction[j][r]  = 0x0;
 138             }
 139           if(j<2)
 140             {
 141               fRelTrkCon[j][r]         = 0x0;
 142               fhFcrossVsZleading[j][r] = 0x0;
 143               fhChLeadZVsJetPt[j][r]   = 0x0;
 144               fhJetPtWithTrkThres[j][r]= 0x0;
 145               fhJetPtWithClsThres[j][r]= 0x0;
 146               for(Int_t k=0; k<2; k++)
 147                 {
 148                   fhJetPtVsLowPtCons[j][r][k] = 0x0;
 149                 }
 150               fhJetPtInExoticEvent[j][r] = 0x0;
 151               fhJetInTPCOnlyVtx[j][r]  = 0x0;
 152               fhCorrTrkEffPtBin[j][r] = 0x0;
 153               for(Int_t i=0; i<kNBins; i++)
 154                 {
 155                   fhCorrTrkEffSample[j][r][i] = 0x0;
 156                 }
 157             }
 158
 159           if(r==0 && j<3)
 160             {
 161               for(Int_t k=0; k<2; k++)
 162                 {
 163                   for(Int_t l=0; l<2; l++)
 164                     {
 165                       fhUEJetPtNorm[j][k][l] = 0x0;
 166                       fhUEJetPtVsSumPt[j][k][l] = 0x0;
 167                       fhUEJetPtVsConsPt[j][k][l] = 0x0;
 168                     }
 169                 }
 170             }
 171         }
 172       for(Int_t i=0; i<2; i++)
 173         {
 174           fhNKFracVsJetPt[i][r]     = 0x0;
 175           fhWeakFracVsJetPt[i][r]   = 0x0;
 176           fhJetResponseNK[i][r]     = 0x0;
 177           fhJetResponseWP[i][r]     = 0x0;
 178           fhJetResolutionNK[i][r]   = 0x0;
 179           fhJetResolutionWP[i][r]   = 0x0;
 180           fhJetResponseNKSM[i][r]   = 0x0;
 181           fhJetResponseWPSM[i][r]   = 0x0;
 182           fhJetResolutionNKSM[i][r] = 0x0;
 183           fhJetResolutionWPSM[i][r] = 0x0;
 184         }
 185     }
 186
 187   for(Int_t s=0; s<3; s++)
 188     {
 189       for(Int_t a=0; a<2; a++)
 190         {
 191           for(Int_t r=0; r<kNRadii; r++)
 192             {
 193               fDetJetFinder[s][a][r] = 0x0;
 194               fJetTCA[s][a][r] = 0x0;
 195               if(s==0 && a==0)
 196                 {
 197                   fTrueJetFinder[r] = 0x0;
 198                   fMcTruthAntikt[r] = 0x0;
 199                 }
 200             }
 201         }
 202     }
 203
 204   for(Int_t j=0; j<3; j++)
 205     {
 206       fTrkEffFunc[j] = 0x0;
 207       fhSecondaryResponse[j] = 0x0;
 208     }
 209   for(Int_t i=0; i<10; i++)
 210     {
 211       fTriggerCurve[i] = 0x0;
 212       fTriggerEfficiency[i] = 0x0;
 213     }
 214
 215 }
 216
 217 //________________________________________________________________________
 218 AliAnalysisTaskFullppJet::AliAnalysisTaskFullppJet(const char *name) :
 219   AliAnalysisTaskSE(name),
 220   fVerbosity(0), fEDSFileCounter(0), fNTracksPerChunk(0), fRejectPileup(kFALSE), fRejectExoticTrigger(kTRUE),
 221   fAnaType(0), fPeriod("lhc11a"), fESD(0), fAOD(0), fMC(0), fStack(0), fTrackArray(0x0), fClusterArray(0x0), fMcPartArray(0x0),
 222   fIsMC(kFALSE), fPhySelForMC(kFALSE), fChargedMC(kFALSE), fXsecScale(0.),
 223   fCentrality(99), fZVtxMax(10),
 224   fTriggerType(-1), fCheckTriggerMask(kTRUE), fIsTPCOnlyVtx(kFALSE),
 225   fIsExoticEvent3GeV(kFALSE), fIsExoticEvent5GeV(kFALSE), fIsEventTriggerBit(kFALSE), fOfflineTrigger(kFALSE), fTriggerMask(0x0),
 226   fGeom(0x0), fRecoUtil(0x0),
 227   fEsdTrackCuts(0x0), fHybridTrackCuts1(0x0), fHybridTrackCuts2(0x0), fTrackCutsType(0), fKinCutType(0), fTrkEtaMax(0.9),
 228   fdEdxMin(73), fdEdxMax(90), fEoverPMin(0.8), fEoverPMax(1.2),
 229   fMatchType(0), fRejectExoticCluster(kTRUE), fRemoveBadChannel(kFALSE), fUseGoodSM(kFALSE),
 230   fStudySubEInHC(kFALSE), fStudyMcOverSubE(kFALSE),
 231   fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0),
 232   fJetNEFMin(0.01), fJetNEFMax(0.98),
 233   fSpotGoodJet(kFALSE), fFindChargedOnlyJet(kFALSE), fFindNeutralOnlyJet(kFALSE),
 234   fCheckTrkEffCorr(kFALSE), fTrkEffCorrCutZ(0.3), fRandomGen(0x0), fRunUE(0), fCheckTPCOnlyVtx(0), fRunSecondaries(0),
 235   fSysJetTrigEff(kFALSE), fVaryJetTrigEff(0),
 236   fSysTrkPtRes(kFALSE), fVaryTrkPtRes(0), fSysTrkEff(kFALSE), fVaryTrkEff(0.),
 237   fSysTrkClsMth(kFALSE), fCutdEta(0.05), fCutdPhi(0.05),
 238   fSysNonLinearity(kFALSE), fNonLinear(0x0), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0),
 239   fSysClusterizer(0),
 240   fNonStdBranch(""), fNonStdFile(""), fAlgorithm("aKt"), fRadius("0.4"), fRecombinationScheme(5),
 241   fConstrainChInEMCal(kFALSE), fRejectNK(kFALSE), fRejectWD(kFALSE), fSmearMC(kFALSE), fTrkPtResData(0x0),
 242   fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventCount(0x0), fhJetEventStat(0x0), fhEventStatTPCVtx(0x0), fhChunkQA(0x0)
 243 {
 244   // Constructor
 245
 246   // Define input and output slots here
 247   // Input slot #0 works with a TChain
 248   DefineInput(0, TChain::Class());
 249   DefineOutput(1, TList::Class());
 250   // Output slot #0 id reserved by the base class for AOD
 251   for(Int_t i=0; i<2; i++)
 252     {
 253       fTrkPtMin[i] = 0.15;
 254       fTrkPtMax[i] = 200;
 255       fClsEtMin[i] = 0.15;
 256       fClsEtMax[i] = 200;
 257
 258       fVertexGenZ[i]       = 0x0;
 259       fEventZ[i]           = 0x0;
 260       fhNTrials[i]         = 0x0;
 261       fhNMatchedTrack[i]   = 0x0;
 262       fhClsE[i]            = 0x0;
 263       for(Int_t k=0; k<2; k++)
 264         {
 265           fhSysClusterE[i][k]     = 0x0;
 266           fhSysNCellVsClsE[i][k]  = 0x0;
 267         }
 268     }
 269
 270   for(Int_t r=0; r<kNRadii; r++)
 271     {
 272       for(Int_t j=0; j<5; j++)
 273         {
 274           fHCOverSubE[r][j]       = 0x0;
 275           fHCOverSubEFrac[r][j]   = 0x0;
 276           for(Int_t k=0; k<2; k++)
 277             {
 278               fhSubClsEVsJetPt[k][r][j] = 0x0;
 279               fhHCTrkPtClean[k][r][j]   = 0x0;
 280               fhHCTrkPtAmbig[k][r][j]   = 0x0;
 281             }
 282
 283           if(j<4) fhSubEVsTrkPt[r][j] = 0x0;
 284
 285           if(j<3)
 286             {
 287               fJetCount[j][r]          = 0x0;
 288               fhNeutralPtInJet[j][r]   = 0x0;
 289               fhTrigNeuPtInJet[j][r]   = 0x0;
 290               fhChargedPtInJet[j][r]   = 0x0;
 291               fhLeadNePtInJet[j][r]    = 0x0;
 292               fhLeadChPtInJet[j][r]    = 0x0;
 293               fJetEnergyFraction[j][r] = 0x0;
 294               fJetNPartFraction[j][r]  = 0x0;
 295             }
 296           if(j<2)
 297             {
 298               fRelTrkCon[j][r]         = 0x0;
 299               fhFcrossVsZleading[j][r] = 0x0;
 300               fhChLeadZVsJetPt[j][r]   = 0x0;
 301               fhJetPtWithTrkThres[j][r]= 0x0;
 302               fhJetPtWithClsThres[j][r]= 0x0;
 303               for(Int_t k=0; k<2; k++)
 304                 {
 305                   fhJetPtVsLowPtCons[j][r][k] = 0x0;
 306                 }
 307               fhJetPtInExoticEvent[j][r] = 0x0;
 308               fhJetInTPCOnlyVtx[j][r]  = 0x0;
 309               fhCorrTrkEffPtBin[j][r] = 0x0;
 310               for(Int_t i=0; i<kNBins; i++)
 311                 {
 312                   fhCorrTrkEffSample[j][r][i] = 0x0;
 313                 }
 314             }
 315
 316           if(r==0 && j<3)
 317             {
 318               for(Int_t k=0; k<2; k++)
 319                 {
 320                   for(Int_t l=0; l<2; l++)
 321                     {
 322                       fhUEJetPtNorm[j][k][l] = 0x0;
 323                       fhUEJetPtVsSumPt[j][k][l] = 0x0;
 324                       fhUEJetPtVsConsPt[j][k][l] = 0x0;
 325                     }
 326                 }
 327             }
 328         }
 329       for(Int_t i=0; i<2; i++)
 330         {
 331           fhNKFracVsJetPt[i][r]     = 0x0;
 332           fhWeakFracVsJetPt[i][r]   = 0x0;
 333           fhJetResponseNK[i][r]     = 0x0;
 334           fhJetResponseWP[i][r]     = 0x0;
 335           fhJetResolutionNK[i][r]   = 0x0;
 336           fhJetResolutionWP[i][r]   = 0x0;
 337           fhJetResponseNKSM[i][r]   = 0x0;
 338           fhJetResponseWPSM[i][r]   = 0x0;
 339           fhJetResolutionNKSM[i][r] = 0x0;
 340           fhJetResolutionWPSM[i][r] = 0x0;
 341         }
 342     }
 343
 344   for(Int_t s=0; s<3; s++)
 345     {
 346       for(Int_t a=0; a<2; a++)
 347         {
 348           for(Int_t r=0; r<kNRadii; r++)
 349             {
 350               fDetJetFinder[s][a][r] = 0x0;
 351               fJetTCA[s][a][r] = 0x0;
 352               if(s==0 && a==0)
 353                 {
 354                   fTrueJetFinder[r] = 0x0;
 355                   fMcTruthAntikt[r] = 0x0;
 356                 }
 357             }
 358         }
 359     }
 360
 361   for(Int_t j=0; j<3; j++)
 362     {
 363       fTrkEffFunc[j] = 0x0;
 364       fhSecondaryResponse[j] = 0x0;
 365     }
 366   for(Int_t i=0; i<10; i++)
 367     {
 368       fTriggerCurve[i] = 0x0;
 369       fTriggerEfficiency[i] = 0x0;
 370     }
 371 }
 372
 373 //________________________________________________________________________
 374 AliAnalysisTaskFullppJet::~AliAnalysisTaskFullppJet()
 375 {
 376   //Destructor
 377
 378   if(fEsdTrackCuts) delete fEsdTrackCuts;
 379   if(fHybridTrackCuts1) delete fHybridTrackCuts1;
 380   if(fHybridTrackCuts2) delete fHybridTrackCuts2;
 381   if(fOutputList)
 382     { fOutputList->Delete(); delete fOutputList;}
 383   for(Int_t s=0; s<3; s++)
 384     {
 385       for(Int_t a=0; a<2; a++)
 386         {
 387           for(Int_t r=0; r<kNRadii; r++)
 388             {
 389               if(fDetJetFinder[s][a][r]) delete fDetJetFinder[s][a][r];
 390               if(fJetTCA[s][a][r]) { fJetTCA[s][a][r]->Delete(); delete fJetTCA[s][a][r]; }
 391               if(s==0 && a==0)
 392                 {
 393                   if(fTrueJetFinder[r]) delete fTrueJetFinder[r];
 394                   if(fMcTruthAntikt[r]) { fMcTruthAntikt[r]->Delete(); delete fMcTruthAntikt[r]; }
 395                 }
 396             }
 397         }
 398     }
 399   if(fRandomGen) delete fRandomGen;
 400   for(Int_t i=0; i<3; i++)
 401     {
 402       if(fTrkEffFunc[i]) delete fTrkEffFunc[i];
 403       if(fhSecondaryResponse[i]) delete fhSecondaryResponse[i];
 404     }
 405   for(Int_t i=0; i<10; i++)
 406     {
 407       if(fTriggerEfficiency[i]) delete fTriggerEfficiency[i];
 408       if(fTriggerCurve[i]) delete fTriggerCurve[i];
 409     }
 410   for(Int_t r=0; r<kNRadii; r++)
 411     {
 412       for(Int_t j=0; j<2; j++)
 413         {
 414           if(fhCorrTrkEffPtBin[j][r]) delete fhCorrTrkEffPtBin[j][r];
 415           for(Int_t i=0; i<kNBins; i++)
 416             {
 417               if(fhCorrTrkEffSample[j][r][i]) delete fhCorrTrkEffSample[j][r][i];
 418             }
 419         }
 420     }
 421   if(fTrackArray)   { fTrackArray->Delete(); delete fTrackArray; }
 422   if(fClusterArray) { fClusterArray->Delete(); delete fClusterArray; }
 423   if(fMcPartArray)  { fMcPartArray->Reset(); delete fMcPartArray; }
 424
 425   if(fRecoUtil) delete fRecoUtil;
 426   if(fClusterEResolution) delete fClusterEResolution;
 427   if(fNonLinear) delete fNonLinear;
 428   if(fTrkPtResData) delete fTrkPtResData;
 429 }
 430
 431 //________________________________________________________________________
 432 Bool_t AliAnalysisTaskFullppJet::Notify()
 433 {
 434   //
 435   // Fill the number of tracks per chunk
 436   //
 437
 438   fhChunkQA->SetBinContent(fEDSFileCounter,fNTracksPerChunk);
 439   fNTracksPerChunk = 0;
 440   fEDSFileCounter++;
 441   return kTRUE;
 442 }
 443
 444 //________________________________________________________________________
 445 void AliAnalysisTaskFullppJet::UserCreateOutputObjects()
 446 {
 447   // Create histograms
 448   // Called once
 449   //
 450
 451   if(fRunUE) fFindChargedOnlyJet = kTRUE;
 452
 453   const Int_t nTrkPtBins = 100;
 454   const Float_t lowTrkPtBin=0, upTrkPtBin=100.;
 455
 456   const Int_t nbins = 220;
 457   Double_t xbins[221];
 458   for(Int_t i=0; i<nbins+1; i++)
 459     xbins[i] = i;
 460
 461   OpenFile(1);
 462   fOutputList = new TList();
 463   fOutputList->SetOwner(1);
 464
 465   fhJetEventStat = new TH1F("fhJetEventStat","Event statistics for jet analysis",12,0,12);
 466   fhJetEventStat->GetXaxis()->SetBinLabel(1,"ALL");
 467   fhJetEventStat->GetXaxis()->SetBinLabel(2,"MB");
 468   fhJetEventStat->GetXaxis()->SetBinLabel(3,"MB+vtx+10cm");
 469   fhJetEventStat->GetXaxis()->SetBinLabel(4,"EMC");
 470   fhJetEventStat->GetXaxis()->SetBinLabel(5,"EMC+vtx+10cm");
 471   fhJetEventStat->GetXaxis()->SetBinLabel(6,"MB+vtx");
 472   fhJetEventStat->GetXaxis()->SetBinLabel(7,"EMC+vtx");
 473   fhJetEventStat->GetXaxis()->SetBinLabel(8,"TriggerBit");
 474   fhJetEventStat->GetXaxis()->SetBinLabel(9,"LED");
 475   fhJetEventStat->GetXaxis()->SetBinLabel(10,"MB+TVtx+10cm");
 476   fhJetEventStat->GetXaxis()->SetBinLabel(11,"EMC+TVtx+10cm");
 477   fhJetEventStat->GetXaxis()->SetBinLabel(12,"ALL-Pileup");
 478   fOutputList->Add(fhJetEventStat);
 479
 480   fhJetEventCount = new TH1F("fhJetEventCount","Event statistics for jet analysis",12,0,12);
 481   fhJetEventCount->GetXaxis()->SetBinLabel(1,"ALL");
 482   fhJetEventCount->GetXaxis()->SetBinLabel(2,"MB");
 483   fhJetEventCount->GetXaxis()->SetBinLabel(3,"MB-pileup");
 484   fhJetEventCount->GetXaxis()->SetBinLabel(4,"MB+Vtx");
 485   fhJetEventCount->GetXaxis()->SetBinLabel(5,"MB+Vtx+10cm");
 486   fhJetEventCount->GetXaxis()->SetBinLabel(6,"EMC");
 487   fhJetEventCount->GetXaxis()->SetBinLabel(7,"EMC-pileup");
 488   fhJetEventCount->GetXaxis()->SetBinLabel(8,"EMC+Vtx");
 489   fhJetEventCount->GetXaxis()->SetBinLabel(9,"EMC+Vtx+10cm");
 490   fhJetEventCount->GetXaxis()->SetBinLabel(10,"Good EMC");
 491   fOutputList->Add(fhJetEventCount);
 492
 493   if(fCheckTPCOnlyVtx)
 494     {
 495       fhEventStatTPCVtx = new TH1F("fhEventStatTPCVtx","Event statistics for TPC only vertex",9,0,9);
 496       fhEventStatTPCVtx->GetXaxis()->SetBinLabel(1,"FastOnly");
 497       fhEventStatTPCVtx->GetXaxis()->SetBinLabel(2,"FastOnly+PVtx");
 498       fhEventStatTPCVtx->GetXaxis()->SetBinLabel(3,"FastOnly+TVtx");
 499       fhEventStatTPCVtx->GetXaxis()->SetBinLabel(4,"MB");
 500       fhEventStatTPCVtx->GetXaxis()->SetBinLabel(5,"MB+PVtx");
 501       fhEventStatTPCVtx->GetXaxis()->SetBinLabel(6,"MB+TVtx");
 502       fhEventStatTPCVtx->GetXaxis()->SetBinLabel(7,"EMC");
 503       fhEventStatTPCVtx->GetXaxis()->SetBinLabel(8,"EMC+PVtx");
 504       fhEventStatTPCVtx->GetXaxis()->SetBinLabel(9,"EMC+TVtx");
 505       fOutputList->Add(fhEventStatTPCVtx);
 506     }
 507
 508   fhChunkQA = new TH1F("fhChunkQA","# of hybrid tracks per chunk",200,0,200);
 509   fOutputList->Add(fhChunkQA);
 510
 511   const Int_t dim1 = 3;
 512   Int_t nBins1[dim1]     = {200,50,110};
 513   Double_t lowBin1[dim1] = {0,0,0,};
 514   Double_t upBin1[dim1]  = {200,0.5,1.1};
 515
 516   const Int_t dim2 = 3;
 517   Int_t nBins2[dim2]     = {200,50,50};
 518   Double_t lowBin2[dim2] = {0,0,0,};
 519   Double_t upBin2[dim2]  = {200,0.5,50};
 520
 521   const char* triggerName[3] = {"MB","EMC","MC"};
 522   const char* triggerTitle[3] = {"MB","EMCal-trigger","MC true"};
 523   const char* fraction[5] = {"MIP","30","50","70","100"};
 524   const char* exotic[2] = {"3GeV","5GeV"};
 525   const char* vertexType[2] = {"All MB","MB with vertex"};
 526   const char *vertexName[2] = {"All","Vertex"};
 527   const char *clusterizerName[2] = {"before","after"};
 528   const char *UEName[2] = {"charged","charged_neutral"};
 529   const char *UETitle[2] = {"charged","charged+neutral"};
 530   const char *UEEventName[2] = {"LeadingJet","Back-To-Back"};
 531
 532   if(fIsMC)
 533     {
 534       for(Int_t i=0; i<2; i++)
 535         {
 536           fhNTrials[i] = new TH1F(Form("MC_%s_fhNTrials",triggerName[i]),Form("MC-%s: # of trials",triggerName[i]),1,0,1);
 537           fOutputList->Add(fhNTrials[i]);
 538
 539           fVertexGenZ[i] = new TH1F(Form("%s_fVertexGenZ",vertexName[i]),Form("Distribution of vertex z (%s); z (cm)",vertexType[i]),60,-30,30);
 540           fOutputList->Add(fVertexGenZ[i]);
 541         }
 542     }
 543
 544   for(Int_t i=0; i<3; i++)
 545     {
 546       if(!fIsMC && i==2) continue;
 547
 548       if(fSaveQAHistos)
 549         {
 550           if(i<2)
 551             {
 552               fEventZ[i] = new TH1F(Form("%s_fEventZ",triggerName[i]),Form("%s: Distribution of vertex z; z (cm)",triggerTitle[i]),60,-30,30);
 553               fOutputList->Add(fEventZ[i]);
 554             }
 555
 556           for(Int_t r=0; r<kNRadii; r++)
 557             {
 558               if(!fRadius.Contains("0.4") && r==0) continue;
 559               if(!fRadius.Contains("0.2") && r==1) continue;
 560               if(!fRadius.Contains("0.3") && r==2) continue;
 561
 562               fJetCount[i][r] = new TH1F(Form("%s_fJetCount_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins);
 563               fOutputList->Add(fJetCount[i][r]);
 564
 565               fhNeutralPtInJet[i][r] = new TH2F(Form("%s_fhNeutralPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of neutral constituents vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,nTrkPtBins*10,lowTrkPtBin,upTrkPtBin);
 566               fOutputList->Add(fhNeutralPtInJet[i][r]);
 567
 568               fhTrigNeuPtInJet[i][r] = new TH2F(Form("%s_fhTrigNeuPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of triggered neutral constituents vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,nTrkPtBins*10,lowTrkPtBin,upTrkPtBin);
 569               fOutputList->Add(fhTrigNeuPtInJet[i][r]);
 570
 571               fhChargedPtInJet[i][r] = new TH2F(Form("%s_fhChargedPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of charged constituents vs jet p_{T} in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,nTrkPtBins*10,lowTrkPtBin,upTrkPtBin);
 572               fOutputList->Add(fhChargedPtInJet[i][r]);
 573
 574               fhLeadNePtInJet[i][r] = new TH2F(Form("%s_fhLeadNePtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of leading neutral constituent vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T}^{ne} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,100);
 575               fOutputList->Add(fhLeadNePtInJet[i][r]);
 576
 577               fhLeadChPtInJet[i][r] = new TH2F(Form("%s_fhLeadChPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of leading charged constituent vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T}^{ch} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,100);
 578               fOutputList->Add(fhLeadChPtInJet[i][r]);
 579
 580               fhJetPtVsZ[i][r] = new TH3F(Form("%s_fhJetPtVsZ_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} vs Z_{h} vs constituent type in EMCal (R=%1.1f);p_{T}^{jet} (GeV/c); Z_{h};constituent type",triggerTitle[i],kRadius[r]),200,0,200,110,-0.05,1.05,4,0,4);
 581               fOutputList->Add(fhJetPtVsZ[i][r]);
 582
 583               fJetEnergyFraction[i][r] = new THnSparseF(Form("%s_fJetEnergyFraction_%1.1f",triggerName[i],kRadius[r]),Form("%s: Jet p_{T} vs radius vs energy fraction in EMCal (R=%1.1f,Z<0.98); p_{T};R;Fraction",triggerName[i],kRadius[r]),dim1,nBins1,lowBin1,upBin1);
 584               fOutputList->Add(fJetEnergyFraction[i][r]);
 585
 586               fJetNPartFraction[i][r] = new THnSparseF(Form("%s_fJetNPartFraction_%1.1f",triggerName[i],kRadius[r]),Form("%s: Jet p_{T} vs radius vs NPart in EMCal (R=%1.1f,Z<0.98); p_{T};R;NPart",triggerName[i],kRadius[r]),dim2,nBins2,lowBin2,upBin2);
 587               fOutputList->Add(fJetNPartFraction[i][r]);
 588
 589               if(i<2)
 590                 {
 591                   fhJetPtInExoticEvent[i][r] = new TH1F(Form("EMC_fhJetPtInExoticEvent_%1.1f_%s",kRadius[r],exotic[i]),Form("EMC: jet p_{T} in events with exotic cluster > %s (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c)",exotic[i],kRadius[r]),nbins,xbins);
 592                   fOutputList->Add(fhJetPtInExoticEvent[i][r]);
 593
 594                   fRelTrkCon[i][r] = new TH3F(Form("%s_fRelTrkCon_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} vs (sum p_{T}^{ch})/p_{T}^{jet} vs track class in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); (sum p_{T}^{ch})/p_{T}^{jet}; track class",triggerTitle[i],kRadius[r]),200,0,200,110,-0.05,1.05,3,0,3);
 595                   fOutputList->Add(fRelTrkCon[i][r]);
 596
 597                   fhFcrossVsZleading[i][r] = new TH3F(Form("%s_fhFcrossVsZleading_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} vs F_{cross} vs Z_{leading}^{ne} in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c);F_{cross};Z_{leading}^{ne}",triggerTitle[i],kRadius[r]),200,0,200,55,-0.05,1.05,55,-0.05,1.05);
 598                   fOutputList->Add(fhFcrossVsZleading[i][r]);
 599
 600                   fhChLeadZVsJetPt[i][r]   = new TH2F(Form("%s_fhChLeadZVsJetPt_%1.1f",triggerName[i],kRadius[r]),Form("%s: Z of leading charged constituent vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); Z_{leading}^{ch}",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,1);
 601                   fOutputList->Add(fhChLeadZVsJetPt[i][r]);
 602
 603                   fhJetPtWithTrkThres[i][r]   = new TH2F(Form("%s_fhJetPtWithTrkThres_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of jets containing tracks above certain threshold (15,25,40 GeV/c) (EMCal, R=%1.1f,Z<0.98);Threshold type;p_{T}^{jet} (GeV/c)",triggerTitle[i],kRadius[r]),3,0,3,nbins,xbins);
 604                   fOutputList->Add(fhJetPtWithTrkThres[i][r]);
 605
 606                   fhJetPtWithClsThres[i][r]   = new TH2F(Form("%s_fhJetPtWithClsThres_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of jets containing clusters above certain threshold (15,25,40 GeV) (EMCal, R=%1.1f,Z<0.98);Threshold type;p_{T}^{jet} (GeV/c)",triggerTitle[i],kRadius[r]),3,0,3,nbins,xbins);
 607                   fOutputList->Add(fhJetPtWithClsThres[i][r]);
 608
 609                   fhJetPtVsLowPtCons[i][r][0]   = new TH2F(Form("%s_fhJetPtVsLowPtCons_150-300MeV_%1.1f",triggerName[i],kRadius[r]),Form("%s: energy carried by constituents in 150-300MeV (EMCal, R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c);p_{T,em}^{low} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,1);
 610                   fOutputList->Add(fhJetPtVsLowPtCons[i][r][0]);
 611
 612                   fhJetPtVsLowPtCons[i][r][1]   = new TH2F(Form("%s_fhJetPtVsLowPtCons_300-500MeV_%1.1f",triggerName[i],kRadius[r]),Form("%s: energy carried by constituents in 300-500MeV (EMCal, R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c);p_{T,em}^{low} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,1);
 613                   fOutputList->Add(fhJetPtVsLowPtCons[i][r][1]);
 614
 615                   if(fCheckTPCOnlyVtx)
 616                     {
 617                       fhJetInTPCOnlyVtx[i][r] = new TH3F(Form("%s_fhJetInTPCOnlyVtx_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet pt in events with only TPC vertex (Full, R=%1.1f);p_{T}^{jet} (GeV/c);#phi;#eta",triggerTitle[i],kRadius[r]),20,0,100,36,0,360,20,-1,1);
 618                       fOutputList->Add(fhJetInTPCOnlyVtx[i][r]);
 619                     }
 620
 621                   if(fStudySubEInHC)
 622                     {
 623                       for(Int_t k=0; k<5; k++)
 624                         {
 625                           fhSubClsEVsJetPt[i][r][k] = new TH2F(Form("%s_fhSubClsEVsJetPt_%s_%1.1f",triggerName[i],fraction[k],kRadius[r]),Form("%s: relative %s%% subtracted cluster Et vs jet pt (R=%1.1f);jet p_{T} (GeV/c);E_{t}^{sub}/p_{T,jet}",triggerTitle[i],fraction[k], kRadius[r]),nbins,xbins,50,0,0.5);
 626                           fOutputList->Add(fhSubClsEVsJetPt[i][r][k]);
 627
 628                           fhHCTrkPtClean[i][r][k] = new TH2F(Form("%s_fhHCTrkPtClean_%s_%1.1f",triggerName[i],fraction[k],kRadius[r]),Form("%s: sum of track p_{T} that are cleanly subtracted  vs jet pt (%s%%, R=%1.1f);jet p_{T} (GeV/c);#sum(p_{T,trk}^{clean})/p_{T,jet}",triggerTitle[i],fraction[k], kRadius[r]),nbins,xbins,100,-0.005,0.995);
 629                           fOutputList->Add(fhHCTrkPtClean[i][r][k]);
 630
 631                           fhHCTrkPtAmbig[i][r][k] = new TH2F(Form("%s_fhHCTrkPtAmbig_%s_%1.1f",triggerName[i],fraction[k],kRadius[r]),Form("%s: sum of track p_{T} that are ambiguously subtracted vs jet pt (%s%%, R=%1.1f);jet p_{T} (GeV/c);#sum(p_{T,trk}^{ambig})/p_{T,jet}",triggerTitle[i],fraction[k], kRadius[r]),nbins,xbins,100,-0.005,0.995);
 632                           fOutputList->Add(fhHCTrkPtAmbig[i][r][k]);
 633                         }
 634                     }
 635                 }
 636             }
 637           if(i<2)
 638             {
 639               fhNMatchedTrack[i] = new TH1F(Form("%s_fhNMatchedTrack",triggerName[i]),Form("%s: # of matched tracks per cluster; N_{mth}",triggerTitle[i]),5,0,5);
 640               fOutputList->Add(fhNMatchedTrack[i]);
 641
 642               for(Int_t j=0; j<4; j++)
 643                 {
 644                   fhSubEVsTrkPt[i][j] = new TH2F(Form("%s_fhSubEVsTrkPt_%s",triggerName[i],fraction[j+1]),Form("%s: fractional subtracted energy (%s%% HC);#sum(p_{ch,T}^{mth}) (GeV/c);E_{sub}/#sum(P_{ch}^{mth})",triggerTitle[i],fraction[j+1]),50,0,50,110,0,1.1);
 645                   fOutputList->Add(fhSubEVsTrkPt[i][j]);
 646                 }
 647             }
 648         }
 649
 650       if(fRunUE)
 651         {
 652           for(Int_t k=0; k<2; k++)
 653             {
 654               for(Int_t l=0; l<2; l++)
 655                 {
 656                   fhUEJetPtNorm[i][k][l] = new TH1F(Form("%s_fhUEJetPtNorm_%s_%s",triggerName[i],UEName[k],UEEventName[l]),Form("%s: leading jet p_{T} in TPC (%s in %s event);p_{T,jet}^{ch} (GeV/c)",triggerTitle[i],UETitle[k], UEEventName[l]),nbins,xbins);
 657                   fOutputList->Add(fhUEJetPtNorm[i][k][l]);
 658
 659                   fhUEJetPtVsSumPt[i][k][l] = new TH2F(Form("%s_fhUEJetPtVsSumPt_%s_%s",triggerName[i],UEName[k],UEEventName[l]),Form("%s: leading jet p_{T} vs underlying event contribution (R=0.4,%s in %s event);p_{T,jet}^{ch} (GeV/c);p_{T,UE}^{ch} (GeV/c)",triggerTitle[i],UETitle[k], UEEventName[l]),nbins,xbins,500,0,50);
 660                   fOutputList->Add(fhUEJetPtVsSumPt[i][k][l]);
 661
 662                   fhUEJetPtVsConsPt[i][k][l] = new TH2F(Form("%s_fhUEJetPtVsConsPt_%s_%s",triggerName[i],UEName[k],UEEventName[l]),Form("%s: leading jet p_{T} vs constituent pt in UE (R=0.4,%s in %s event);p_{T,jet}^{ch} (GeV/c);p_{T,cons}^{ch} (GeV/c)",triggerTitle[i],UETitle[k], UEEventName[l]),nbins,xbins,500,0,50);
 663                   fOutputList->Add(fhUEJetPtVsConsPt[i][k][l]);
 664                 }
 665             }
 666         }
 667
 668       if(fSysJetTrigEff)
 669         {
 670           fhClsE[i] = new TH1F(Form("%s_fhClsE",triggerName[i]),Form("%s: cluster E;E (GeV)",triggerTitle[i]),1000,0,100);
 671           fOutputList->Add(fhClsE[i]);
 672         }
 673
 674       if(fSysClusterizer)
 675         {
 676           for(Int_t k=0; k<2; k++)
 677             {
 678               fhSysClusterE[i][k] = new TH1F(Form("%s_fhSysClusterE_%sHC",triggerName[i],clusterizerName[k]),Form("%s: cluster E %s hadronic correction;E (GeV)",triggerTitle[i],clusterizerName[k]),100,0,100);
 679               fOutputList->Add(fhSysClusterE[i][k]);
 680
 681               fhSysNCellVsClsE[i][k] = new TH2F(Form("%s_fhSysNCellVsClsE_%sHC",triggerName[i],clusterizerName[k]),Form("%s: NCell vs cluster E %s hadronic correction;E (GeV);NCell",triggerTitle[i],clusterizerName[k]),100,0,100,50,0,50);
 682               fOutputList->Add(fhSysNCellVsClsE[i][k]);
 683             }
 684         }
 685     }
 686
 687
 688   if(fIsMC)
 689     {
 690       if(fStudyMcOverSubE)
 691         {
 692           for(Int_t r=0; r<kNRadii; r++)
 693             {
 694               if(!fRadius.Contains("0.4") && r==0) continue;
 695               if(!fRadius.Contains("0.2") && r==1) continue;
 696               if(!fRadius.Contains("0.3") && r==2) continue;
 697               for(Int_t i=0; i<5; i++)
 698                 {
 699                   fHCOverSubE[r][i] = new TH2F(Form("%s_HC_over_sub_e_%s_%1.1f",triggerName[2],fraction[i],kRadius[r]),Form("%s: oversubtracted neutral Et by %s%% HC (R=%1.1f);particle jet p_{T} (GeV/c);#DeltaE_{t} (GeV)",triggerName[2],fraction[i],kRadius[r]),nbins,xbins,200,-49.75,50.25);
 700                   fOutputList->Add(fHCOverSubE[r][i]);
 701                   fHCOverSubEFrac[r][i] = new TH2F(Form("%s_HC_over_sub_e_frac_%s_%1.1f",triggerName[2],fraction[i],kRadius[r]),Form("%s: relative oversubtracted neutral Et fraction by %s%% HC (R=%1.1f);jet p_{T} (GeV/c);#DeltaE_{t}/p_{T}^{jet}",triggerName[2],fraction[i],kRadius[r]),nbins,xbins,200,-0.995,1.005);
 702                   fOutputList->Add(fHCOverSubEFrac[r][i]);
 703                 }
 704             }
 705         }
 706       if(fRunSecondaries && fAnaType==0)
 707         {
 708           for(Int_t i=0; i<2; i++)
 709             {
 710               for(Int_t r=0; r<3; r++)
 711                 {
 712                   fhNKFracVsJetPt[i][r]  = new TH2F(Form("%s_fhNKFracVsJetPt_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: energy fraction carried by n,k^{0}_{L} vs jet p_{T} (R=%1.1f,|#eta|<%1.1f);p_{T,jet} (GeV/c);fraction",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,0,1);
 713                   fOutputList->Add(fhNKFracVsJetPt[i][r]);
 714
 715                   fhWeakFracVsJetPt[i][r]  = new TH2F(Form("%s_fhWeakFracVsJetPt_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: energy fraction carried by k^{0}_{S},hyperon vs jet p_{T} (R=%1.1f,|#eta|<%1.1f);p_{T,jet} (GeV/c);fraction",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,0,1);
 716                   fOutputList->Add(fhWeakFracVsJetPt[i][r]);
 717
 718                   fhJetResponseNK[i][r]  = new TH2F(Form("%s_fhJetResponseNK_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to missing n and k^{0}_{L} (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins);
 719                   fOutputList->Add(fhJetResponseNK[i][r]);
 720
 721                   fhJetResponseWP[i][r]  = new TH2F(Form("%s_fhJetResponseWP_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to k^{0}_{S}, #Lambda and hyperon (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins);
 722                   fOutputList->Add(fhJetResponseWP[i][r]);
 723
 724                   fhJetResolutionNK[i][r] = new TH2F(Form("%s_fhJetResolutionNK_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to missing n and k^{0}_{L}: (p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{rec} vs p_{T,jet}^{rec} (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T}/p_{T}",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,-0.995,1.005);
 725                   fOutputList->Add(fhJetResolutionNK[i][r]);
 726
 727                   fhJetResolutionWP[i][r] = new TH2F(Form("%s_fhJetResolutionWP_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to k^{0}_{S}, #Lambda and hyperon: (p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{rec} vs p_{T,jet}^{rec} (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T}/p_{T}",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,-0.995,1.005);
 728                   fOutputList->Add(fhJetResolutionWP[i][r]);
 729
 730                   fhJetResponseNKSM[i][r]  = new TH2F(Form("%s_fhJetResponseNKSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to missing n and k^{0}_{L} via matching (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins);
 731                   fOutputList->Add(fhJetResponseNKSM[i][r]);
 732
 733                   fhJetResponseWPSM[i][r]  = new TH2F(Form("%s_fhJetResponseWPSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to k^{0}_{S}, #Lambda and hyperon via matching(R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins);
 734                   fOutputList->Add(fhJetResponseWPSM[i][r]);
 735
 736                   fhJetResolutionNKSM[i][r] = new TH3F(Form("%s_fhJetResolutionNKSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet resolution due to missing n and k^{0}_{L} via matching (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T,jet}/p_{T,jet}^{rec};dR",triggerName[2],kRadius[r],i*0.5+0.5),220,0,220,200,-0.995,1.005,100,0,1);
 737                   fOutputList->Add(fhJetResolutionNKSM[i][r]);
 738
 739                   fhJetResolutionWPSM[i][r] = new TH3F(Form("%s_fhJetResolutionWPSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet resolution due tok^{0}_{S}, #Lambda and hyperon via matching (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T,jet}/p_{T,jet}^{rec};dR",triggerName[2],kRadius[r],i*0.5+0.5),220,0,220,200,-0.995,1.005,100,0,1);
 740                   fOutputList->Add(fhJetResolutionWPSM[i][r]);
 741                 }
 742             }
 743         }
 744     }
 745
 746   printf("\n=======================================\n");
 747   printf("===== Jet task configuration ==========\n");
 748
 749   if(fNonStdBranch.Length()!=0)
 750     {
 751       const char* species[3] = {"in","ch","ne"};
 752       const char* algorithm[2] = {"akt","kt"};
 753       const char* radii[kNRadii] = {"04","02","03"};
 754       for(Int_t s=0; s<3; s++)
 755         {
 756           if(!fFindChargedOnlyJet && s==1) continue;
 757           if(!fFindNeutralOnlyJet && s==2) continue;
 758           for(Int_t a=0; a<2; a++)
 759             {
 760               if(!fAlgorithm.Contains("aKt") && a==0) continue;
 761               if(!fAlgorithm.Contains("kt")  && a==1) continue;
 762               for(Int_t r=0; r<kNRadii; r++)
 763                 {
 764                   if(!fRadius.Contains("0.4") && r==0) continue;
 765                   if(!fRadius.Contains("0.2") && r==1) continue;
 766                   if(!fRadius.Contains("0.3") && r==2) continue;
 767                   if(fAnaType==0)
 768                     {
 769                       fJetTCA[s][a][r] = new TClonesArray("AliAODJet",0);
 770                       fJetTCA[s][a][r]->SetName(Form("Jet_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()));
 771                       AddAODBranch("TClonesArray",&fJetTCA[s][a][r],fNonStdFile.Data());
 772                       printf("Add branch: Jet_%s_%s_%s_%s\n",species[s],algorithm[a],radii[r],fNonStdBranch.Data());
 773                     }
 774
 775                   fDetJetFinder[s][a][r] = new AliFJWrapper(Form("DetJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()),Form("DetJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()));
 776                   if(a==0) fDetJetFinder[s][a][r]->SetAlgorithm(fastjet::antikt_algorithm);
 777                   if(a==1) fDetJetFinder[s][a][r]->SetAlgorithm(fastjet::kt_algorithm);
 778                   if(fRecombinationScheme==0)     fDetJetFinder[s][a][r]->SetRecombScheme(fastjet::E_scheme);
 779                   fDetJetFinder[s][a][r]->SetR(kRadius[r]);
 780                   fDetJetFinder[s][a][r]->SetMaxRap(0.9);
 781                   fDetJetFinder[s][a][r]->Clear();
 782
 783                   if(s==0 && a==0)
 784                     {
 785                       if(fIsMC && fNonStdBranch.Contains("Baseline",TString::kIgnoreCase))
 786                         {
 787                           if(fAnaType==0)
 788                             {
 789                               fMcTruthAntikt[r] = new TClonesArray("AliAODJet",0);
 790                               fMcTruthAntikt[r]->SetName(Form("Jet_mc_truth_in_akt_%s_%s",radii[r],fNonStdBranch.Data()));
 791                               AddAODBranch("TClonesArray",&fMcTruthAntikt[r],fNonStdFile.Data());
 792                               printf("Add branch: Jet_mc_truth_in_akt_%s_%s\n",radii[r],fNonStdBranch.Data());
 793                             }
 794
 795                           fTrueJetFinder[r] = new AliFJWrapper(Form("TrueJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()),Form("TrueJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()));
 796                           fTrueJetFinder[r]->SetAlgorithm(fastjet::antikt_algorithm);
 797                           fTrueJetFinder[r]->SetR(kRadius[r]);
 798                           fTrueJetFinder[r]->SetMaxRap(0.9);
 799                           if(fRecombinationScheme==0) fTrueJetFinder[r]->SetRecombScheme(fastjet::E_scheme);
 800                           fTrueJetFinder[r]->Clear();
 801                         }
 802                     }
 803                 }
 804             }
 805         }
 806     }
 807
 808   fRandomGen = new TRandom3(0);
 809   if(fCheckTrkEffCorr && fAnaType==0)
 810     {
 811       TFile f ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TrkEffFit.root","read");
 812       for(Int_t i=0; i<3; i++)
 813         {
 814           fTrkEffFunc[i] = new TF1(*((TF1*)f.Get(Form("Trk_eff_fit_%d",i+1))));
 815         }
 816       f.Close();
 817
 818       if(fPeriod.CompareTo("lhc11a",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11aJet",TString::kIgnoreCase)==0)
 819         {
 820           TFile f1 ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TrkEffSampling.root","read");
 821           for(Int_t j=0; j<2; j++)
 822             {
 823               Int_t tmp = j;
 824               if(fPeriod.CompareTo("lhc11aJet",TString::kIgnoreCase)==0) tmp = 2;
 825               for(Int_t r=0; r<kNRadii; r++)
 826                 {
 827                   fhCorrTrkEffPtBin[j][r] = new TH1F(*((TH1F*)f1.Get(Form("%s_%s_NTrackPerPtBin_%1.1f",fPeriod.Data(),triggerName[tmp],kRadius[r]))));
 828                   for(Int_t i=0; i<kNBins; i++)
 829                     {
 830                       fhCorrTrkEffSample[j][r][i] = new TH1F(*((TH1F*)f1.Get(Form("%s_%s_ChTrkPt_Bin%d_%1.1f",fPeriod.Data(),triggerName[tmp],i+1,kRadius[r]))));
 831                     }
 832                 }
 833             }
 834           f1.Close();
 835         }
 836     }
 837
 838   if(fRunSecondaries && fAnaType==0)
 839     {
 840       const char *secondaryName[3] = {"k0S","lamda","hyperon"};
 841       TFile f2 ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/SecondaryResponse.root","read");
 842       for(Int_t j=0; j<3; j++)
 843         {
 844           fhSecondaryResponse[j] = new TH2F(*((TH2F*)f2.Get(Form("DetectorReponse_%s",secondaryName[j]))));
 845         }
 846     }
 847
 848   if(fCheckTriggerMask)
 849     {
 850       TString name = "TriggerCurve.root";
 851       if(fAnaType==0) name = "/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TriggerCurve.root";
 852       TFile f (name.Data(),"read");
 853       fTriggerMask = new TH2F(*((TH2F*)f.Get("lhc11a_TriggerMask")));
 854       if(fOfflineTrigger)
 855         {
 856           for(Int_t i=0; i<10; i++)
 857             {
 858               fTriggerEfficiency[i] = new TF1(*((TF1*)f.Get(Form("lhc11a_TriggerEfficiency_SM%d_fit",i))));
 859               fTriggerCurve[i] = new TH1D(*((TH1D*)f.Get(Form("lhc11a_TriggerCurve_SM%d",i))));
 860             }
 861         }
 862       f.Close();
 863     }
 864
 865   fClusterEResolution = new TF1("fClusterEResolution","sqrt([0]^2+[1]^2*x+([2]*x)^2)*0.01");
 866   fClusterEResolution->SetParameters(4.35,9.07,1.63);
 867
 868   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
 869   if(!fGeom)
 870     {
 871       AliError("No EMCal geometry is available!");
 872       return;
 873     }
 874   fRecoUtil = new AliEMCALRecoUtils();
 875   if(fRejectExoticCluster || fRejectExoticTrigger)
 876     fRecoUtil->SwitchOnRejectExoticCluster();
 877   else
 878     {
 879       fRecoUtil->SwitchOffRejectExoticCluster();
 880       fRecoUtil->SwitchOffRejectExoticCell();
 881     }
 882   if(fRemoveBadChannel)
 883     {
 884       fRecoUtil->SwitchOnBadChannelsRemoval();
 885       // Remove the problematic SM4 region due to wrong event sequence
 886       for(Int_t ieta=0; ieta<36; ieta++)
 887         for(Int_t iphi=0; iphi<8; iphi++)
 888           fRecoUtil->SetEMCALChannelStatus(4,ieta,iphi);
 889     }
 890   else
 891     fRecoUtil->SwitchOffBadChannelsRemoval();
 892
 893   fRecoUtil->SetNonLinearityFunction(AliEMCALRecoUtils::kBeamTestCorrected);
 894   if(fSysNonLinearity)
 895     {
 896       fNonLinear = new TF1("TB_oldBest","([0])*(1./(1.+[1]*exp(-x/[2]))*1./(1.+[3]*exp((x-[4])/[5])))*1/([6])",0.1,110.);
 897       fNonLinear->SetParameters(0.99078, 0.161499, 0.655166, 0.134101, 163.282, 23.6904, 0.978);
 898     }
 899   if(fSysTrkClsMth)
 900     {
 901       fRecoUtil->SwitchOnCutEtaPhiSeparate();
 902       fRecoUtil->SetCutEta(fCutdEta);
 903       fRecoUtil->SetCutPhi(fCutdPhi);
 904     }
 905
 906   fTrackArray   = new TObjArray();
 907   fTrackArray->SetOwner(1);
 908   fClusterArray = new TObjArray();
 909   fClusterArray->SetOwner(1);
 910   fMcPartArray = new TArrayI();
 911
 912
 913   //error calculation in THnSparse
 914   Int_t nObj = fOutputList->GetEntries();
 915   for(Int_t i=0; i<nObj; i++)
 916     {
 917       TObject *obj = (TObject*) fOutputList->At(i);
 918       if (obj->IsA()->InheritsFrom( "THnSparse" ))
 919         {
 920           THnSparseF *hn = (THnSparseF*)obj;
 921           hn->Sumw2();
 922         }
 923     }
 924
 925
 926   PrintConfig();
 927   BookHistos();
 928   PostData(1, fOutputList);
 929 }
 930
 931
 932
 933 //________________________________________________________________________
 934 void AliAnalysisTaskFullppJet::BookHistos()
 935 {
 936   // Book histograms.
 937
 938   if(fVerbosity>10) printf("[i] Booking histograms \n");
 939   return;
 940 }
 941
 942 //________________________________________________________________________
 943 void AliAnalysisTaskFullppJet::UserExec(Option_t *)
 944 {
 945   //
 946   // Main loop, called for each event.
 947   //
 948
 949   // Get event pointers, check for signs of life
 950   Double_t vtxTrueZ = -100;
 951   if(fIsMC)
 952     {
 953       fMC = MCEvent();
 954       if (!fMC)
 955         {
 956           printf("ERROR: Could not retrieve MC event");
 957           return;
 958         }
 959       fStack = fMC->Stack();
 960       TParticle *particle = fStack->Particle(0);
 961       if(particle)
 962         {
 963           vtxTrueZ = particle->Vz(); // True vertex z in MC events
 964           if(fVerbosity>10) printf("Generated vertex coordinate: (x,y,z) = (%4.2f, %4.2f, %4.2f)\n", particle->Vx(), particle->Vy(), particle->Vz());
 965         }
 966     }
 967
 968   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 969   if (!fESD)
 970     {
 971       AliError("fESD is not available");
 972       return;
 973     }
 974
 975   fhJetEventStat->Fill(0.5);
 976   fhJetEventCount->Fill(0.5);
 977   if(fIsMC)
 978     {
 979       GetMCInfo();
 980       if(fVertexGenZ[0]) fVertexGenZ[0]->Fill(vtxTrueZ);
 981     }
 982
 983   // Centrality, vertex, other event variables...
 984   UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 985   if (trigger==0)  return;
 986   if(fCheckTPCOnlyVtx) CheckTPCOnlyVtx(trigger);
 987   if(!fIsMC)
 988     {
 989       if (trigger & AliVEvent::kFastOnly) return;  // Reject fast trigger cluster
 990       if(fPeriod.CompareTo("lhc11a",TString::kIgnoreCase)==0)
 991         {
 992           if (trigger & AliVEvent::kMB)       fTriggerType = 0;
 993           else if(trigger & AliVEvent::kEMC1) fTriggerType = 1;
 994           else fTriggerType = -1;
 995         }
 996       else if (fPeriod.CompareTo("lhc11c",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11d",TString::kIgnoreCase)==0)
 997         {
 998           if (trigger & AliVEvent::kINT7)     fTriggerType = 0;
 999           else if(trigger & AliVEvent::kEMC7) fTriggerType = 1;
1000           else fTriggerType = -1;
1001         }
1002       else
1003         {
1004           return;
1005         }
1006     }
1007   else
1008     {
1009       if(!fPhySelForMC) fTriggerType = 0;
1010       else if (trigger & AliVEvent::kAnyINT) fTriggerType = 0;
1011       else fTriggerType = -1;
1012
1013       if(fOfflineTrigger)
1014         {
1015           RunOfflineTrigger();
1016           if(fIsEventTriggerBit) fTriggerType = 1;
1017         }
1018     }
1019
1020   if(fTriggerType==-1)
1021     {
1022       if(fVerbosity>10) printf("Error: worng trigger type %s\n",(fESD->GetFiredTriggerClasses()).Data());
1023       return;
1024     }
1025
1026   fhJetEventCount->Fill(1.5+fTriggerType*4);
1027   if(fRejectPileup && fESD->IsPileupFromSPD() ) return; // reject pileup
1028   fhJetEventStat->Fill(11.5);
1029   fhJetEventCount->Fill(2.5+fTriggerType*4);
1030
1031
1032   fIsTPCOnlyVtx = 0;
1033   // Reject LED events
1034   if (IsLEDEvent())
1035     {
1036       fhJetEventStat->Fill(8.5);
1037       return;
1038     }
1039
1040   fhJetEventStat->Fill(1.5+fTriggerType*2);
1041
1042   // Check if primary vertex exists
1043   if (!HasPrimaryVertex()) return;
1044   fhJetEventStat->Fill(5.5+fTriggerType);
1045   fhJetEventCount->Fill(3.5+fTriggerType*4);
1046
1047   const AliESDVertex* vtx = fESD->GetPrimaryVertex();
1048   Double_t zVertex    = vtx->GetZ();
1049   if(fEventZ[fTriggerType]) fEventZ[fTriggerType]->Fill(zVertex);
1050   if(fVertexGenZ[1]) fVertexGenZ[1]->Fill(vtxTrueZ);
1051
1052   // Check if |Z_vtx|<10cm
1053   if( TMath::Abs(zVertex) > fZVtxMax ) return;
1054   fhJetEventCount->Fill(4.5+fTriggerType*4);
1055
1056   // Check if only TPC vertex exists primitive
1057   fIsTPCOnlyVtx = IsTPCOnlyVtx();
1058
1059   if(!fIsMC && fTriggerType==1)
1060     {
1061       // Check if event has valid trigger bit
1062       CheckEventTriggerBit();
1063       if(fIsEventTriggerBit)
1064         {
1065           fhJetEventStat->Fill(7.5);
1066           fhJetEventCount->Fill(9.5);
1067         }
1068       else return;
1069     }
1070   fhJetEventStat->Fill(2.5+fTriggerType*2);
1071   if(fIsTPCOnlyVtx) fhJetEventStat->Fill(9.5+fTriggerType);
1072
1073   // Check if event contains exotic clusters
1074   CheckExoticEvent();
1075
1076   // Write jet tree into AOD output in local analysis mode
1077   if(fAnaType==0) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
1078
1079   // Clean up arrays
1080   for(Int_t s=0; s<3; s++)
1081     {
1082       for(Int_t a=0; a<2; a++)
1083         {
1084           for(Int_t r=0; r<kNRadii; r++)
1085             {
1086               if(fJetTCA[s][a][r])       fJetTCA[s][a][r]->Delete();
1087               if(fDetJetFinder[s][a][r]) fDetJetFinder[s][a][r]->Clear();
1088
1089               if(s==0 && a==0)
1090                 {
1091                   if(fMcTruthAntikt[r]) fMcTruthAntikt[r]->Delete();
1092                   if(fTrueJetFinder[r]) fTrueJetFinder[r]->Clear();
1093                 }
1094             }
1095         }
1096     }
1097
1098   if(fVerbosity>5) printf("# of jets after clear: %d\n",fJetTCA[0][0][0]->GetEntries());
1099
1100   if(fTrackArray) fTrackArray->Delete();
1101   if(fClusterArray) fClusterArray->Delete();
1102   fMcPartArray->Reset();
1103
1104
1105   // get the tracks and fill the input vector for the jet finders
1106   GetESDTrax();
1107
1108   // get EMCal clusters and fill the input vector for the jet finders
1109   GetESDEMCalClusters();
1110
1111   for(Int_t s=0; s<3; s++)
1112     {
1113       for(Int_t a=0; a<2; a++)
1114         {
1115           for(Int_t r=0; r<kNRadii; r++)
1116             {
1117               //Detector jets
1118               if(fDetJetFinder[s][a][r])
1119                 {
1120                   FindDetJets(s,a,r);
1121                   if(fJetTCA[s][a][r]) FillAODJets(fJetTCA[s][a][r], fDetJetFinder[s][a][r], 0);
1122                   if(s==0 && a==0 && fNonStdBranch.Contains("Baseline",TString::kIgnoreCase))
1123                     {
1124                       if(fSaveQAHistos)    AnalyzeJets(fDetJetFinder[s][a][r],fTriggerType, r); // run analysis on jets
1125                     }
1126                 }
1127             }
1128         }
1129     }
1130   if(fRunUE && fDetJetFinder[1][0][0]) RunAnalyzeUE(fDetJetFinder[1][0][0],fTriggerType,0); // Run analysis on underlying event
1131
1132   if(fIsMC)
1133     {
1134       for(Int_t r=0; r<kNRadii; r++)
1135         if(fTrueJetFinder[r]) ProcessMC(r); // find particle level jets
1136     }
1137
1138   // Fast Jet calls END --------------------------------------------------------
1139
1140   PostData(1, fOutputList);
1141   return;
1142 }
1143
1144
1145 //________________________________________________________________________
1146 void AliAnalysisTaskFullppJet::FindDetJets(const Int_t s, const Int_t a, const Int_t r)
1147 {
1148   //
1149   // Jet finding is happening here
1150   //
1151
1152   Bool_t isCh = kTRUE;
1153   Bool_t isNe = kTRUE;
1154   if(s==1) isNe = kFALSE;
1155   if(s==2) isCh = kFALSE;
1156
1157   if(isCh)
1158     {
1159       // Feed in charged tracks
1160       Int_t countTracks = 0;
1161       Int_t ntracks = fTrackArray->GetEntriesFast();
1162       for(Int_t it=0; it<ntracks; it++)
1163         {
1164           AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it);
1165           if(!t) continue;
1166           countTracks++;
1167           Double_t e = t->P();
1168           Double_t pt = t->Pt();
1169           if(s==1 && fRunUE && pt<1) continue;
1170           if(fRecombinationScheme==0) e = TMath::Sqrt(t->P()*t->P()+0.139*0.139);
1171           if(fSysTrkPtRes) pt += GetSmearedTrackPt(t);
1172           Double_t phi = t->Phi();
1173           Double_t eta = t->Eta();
1174           Double_t px = pt*TMath::Cos(phi);
1175           Double_t py = pt*TMath::Sin(phi);
1176           if(fConstrainChInEMCal && ( TMath::Abs(eta)>0.7 || phi>kPI || phi<TMath::DegToRad()*80) ) continue;
1177           fDetJetFinder[s][a][r]->AddInputVector(px,py,t->Pz(), e, it+1);
1178           if(fVerbosity>10) printf("%d: m_{ch}=%f\n",it+1,t->P()*t->P()-t->Px()*t->Px()-t->Py()*t->Py()-t->Pz()*t->Pz());
1179         }
1180       if(fVerbosity>5) printf("[i] # of tracks filled: %d\n",countTracks);
1181     }
1182   if(isNe)
1183     {
1184       // Feed in EMCal clusters
1185       Double_t vertex[3] = {0, 0, 0};
1186       fESD->GetVertex()->GetXYZ(vertex) ;
1187       TLorentzVector gamma;
1188
1189       Int_t countClusters = 0;
1190       Int_t nclusters = fClusterArray->GetEntriesFast();
1191       for (Int_t ic = 0; ic < nclusters; ic++)
1192         {
1193           AliESDCaloCluster * cl = (AliESDCaloCluster *)fClusterArray->At(ic);
1194           if(!cl) continue;
1195           cl->GetMomentum(gamma, vertex);
1196           countClusters++;
1197           Double_t e = gamma.P();
1198           if(fRecombinationScheme==0) e = TMath::Sqrt(gamma.P()*gamma.P()+0.139*0.139);
1199           fDetJetFinder[s][a][r]->AddInputVector(gamma.Px(), gamma.Py(), gamma.Pz(), e, (ic+1)*(-1));
1200         }
1201       if(fVerbosity>5) printf("[i] # of clusters filled: %d\n",countClusters);
1202     }
1203   // Run jet finding
1204   fDetJetFinder[s][a][r]->Run();
1205 }
1206
1207
1208 //________________________________________________________________________
1209 void AliAnalysisTaskFullppJet::FillAODJets(TClonesArray *fJetArray, AliFJWrapper *jetFinder, const Bool_t isTruth)
1210 {
1211   //
1212   // Fill the found jets into AOD output
1213   // Only consider jets pointing to EMCal and with pt above 1GeV/c
1214   //
1215
1216   Int_t radiusIndex = 0;
1217   TString arrayName = fJetArray->GetName();
1218   if(arrayName.Contains("_02_"))  radiusIndex = 1;
1219   if(arrayName.Contains("_03_"))  radiusIndex = 2;
1220   std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1221   if(fVerbosity>5 && radiusIndex==0) printf("[i] # of jets in %s : %d\n",fJetArray->GetName(),(Int_t)jetsIncl.size());
1222   AliAODJet *jet = 0;
1223   Int_t jetCount = 0;
1224   for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1225     {
1226       if(fVerbosity>10) printf("fastjet: eta=%f, phi=%f\n",jetsIncl[ij].eta(),jetsIncl[ij].phi()*TMath::RadToDeg());
1227       if(!IsGoodJet(jetsIncl[ij],0)) continue;
1228
1229       AliAODJet tmpJet (jetsIncl[ij].px(), jetsIncl[ij].py(), jetsIncl[ij].pz(), jetsIncl[ij].E());
1230       jet = new ((*fJetArray)[jetCount]) AliAODJet(tmpJet);
1231       jetCount++;
1232       if(fVerbosity>10 && radiusIndex==0) printf("AOD jet: ij=%d, pt=%f, eta=%f, phi=%f\n",ij, jet->Pt(), jet->Eta(),jet->Phi()*TMath::RadToDeg());
1233
1234       jet->GetRefTracks()->Clear();
1235       std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
1236       Double_t totE=0, totPt=0, totChPt=0, leadChPt=0, neE=0, totNePt=0, leadNePt=0, leadPt=0;
1237       Double_t leadTrkType=0, nChPart = 0, nNePart = 0;
1238       Bool_t isHighPtTrigger = kFALSE, isTriggering = kFALSE, isHyperTrack = kFALSE;
1239       Int_t leadIndex = -1;
1240       for(UInt_t ic=0; ic<constituents.size(); ic++)
1241         {
1242           if(fVerbosity>10 && radiusIndex==0) printf("ic=%d: user_index=%d, E=%f\n",ic,constituents[ic].user_index(),constituents[ic].E());
1243           if(constituents[ic].user_index()<0)
1244             {
1245               nNePart ++;
1246               totNePt += constituents[ic].perp();
1247               if(constituents[ic].perp()>leadNePt)
1248                 leadNePt = constituents[ic].perp();
1249
1250               neE += constituents[ic].E();
1251               if(constituents[ic].perp()>fClsEtMax[fTriggerType])
1252                 isHighPtTrigger = kTRUE;
1253               if(!isTruth)
1254                 {
1255                   AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
1256                   if(cluster->Chi2()>0.5) isTriggering = kTRUE;
1257                 }
1258             }
1259           else
1260             {
1261               totChPt += constituents[ic].perp();
1262               nChPart ++;
1263               if(constituents[ic].perp()>leadChPt)
1264                 {
1265                   leadChPt = constituents[ic].perp();
1266                   if(!isTruth)
1267                     {
1268                       AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1269                       leadTrkType = track->GetTRDQuality();
1270                       if(track->GetIntegratedLength()>500) isHyperTrack = kTRUE;
1271                     }
1272                 }
1273               if(constituents[ic].perp()>fTrkPtMax[fTriggerType])
1274                 isHighPtTrigger = kTRUE;
1275             }
1276           TParticle part;
1277           part.SetMomentum(constituents[ic].px(),constituents[ic].py(),constituents[ic].pz(),constituents[ic].E());
1278           jet->AddTrack(&part); //The references are not usable, this line is aimed to count the number of contituents
1279           totE  += constituents[ic].E();
1280           totPt += constituents[ic].perp();
1281           if(constituents[ic].perp()>leadPt)
1282             {
1283               leadPt = constituents[ic].perp();
1284               leadIndex = ic;
1285             }
1286         }
1287       if(fIsEventTriggerBit)   jet->SetTrigger(AliAODJet::kTRDTriggered);
1288       if(isHighPtTrigger)  jet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1289       if(fTriggerType==1)  jet->SetTrigger(AliAODJet::kEMCALTriggered);
1290       if(GetLeadingZ(ij,jetFinder) > 0.98 )  jet->SetTrigger(AliAnalysisTaskFullppJet::kHighZ);
1291       if(isTriggering) jet->SetTrigger(AliAnalysisTaskFullppJet::kTrigger);
1292       if(isHyperTrack) jet->SetTrigger(AliAnalysisTaskFullppJet::kSuspicious);
1293       if(fIsTPCOnlyVtx) jet->SetTrigger(AliAnalysisTaskFullppJet::kTPCOnlyVtx);
1294       if(constituents[leadIndex].user_index()>0) jet->SetTrigger(AliAnalysisTaskFullppJet::kLeadCh);
1295
1296       if(jetsIncl[ij].E()>0)  jet->SetNEF(neE/jetsIncl[ij].E());
1297       jet->SetPtLeading(leadPt);
1298       jet->SetPtSubtracted(leadTrkType,0);
1299
1300       Double_t  effAErrCh = 0, effAErrNe = leadTrkType;
1301       Double_t chBgEnergy = 10, neBgEnergy = nNePart;
1302       if(!isTruth)
1303         {
1304           effAErrCh = GetMeasuredJetPtResolution(ij,jetFinder);
1305           if(fIsMC && constituents[leadIndex].user_index()>0)
1306             {
1307               AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[leadIndex].user_index()-1);
1308               Double_t pt = track->Pt();
1309               Int_t ipart = track->GetLabel();
1310               if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
1311                 {
1312                   AliVParticle* vParticle = fMC->GetTrack(ipart);
1313                   if(vParticle)
1314                     {
1315                       Double_t truePt = vParticle->Pt();
1316                       chBgEnergy = (truePt-pt)/truePt;
1317                     }
1318                 }
1319             }
1320         }
1321
1322       jet->SetEffArea(leadPt, jetFinder->GetJetArea(ij),effAErrCh,effAErrNe);
1323       jet->SetBgEnergy(chBgEnergy,neBgEnergy);
1324       if(fVerbosity>10 && isTruth) cout<<jet->ErrorEffectiveAreaCharged()<<"  "<<jet->ErrorEffectiveAreaNeutral()<<endl;
1325       if(fVerbosity>10) cout<<"jet pt="<<jetsIncl[ij].perp()<<", nef="<<jet->GetNEF()<<", trk eff corr="<<chBgEnergy<<endl;
1326
1327       if(fVerbosity>5)
1328         printf("# of ref tracks: %d = %d, and nef=%f\n",jet->GetRefTracks()->GetEntries(), (Int_t)constituents.size(), jet->GetNEF());
1329
1330       // For catch good high-E jets
1331       if(fSpotGoodJet && !isTruth)
1332         {
1333           if(jetsIncl[ij].perp()>100)
1334             {
1335               printf("\n\n--- HAHAHA: High pt jets ---\n");
1336               printf("File: %s, event = %d, pt=%f\n",CurrentFileName(),(Int_t)Entry(),jetsIncl[ij].perp());
1337               printf("%s , pt < %f\n", fJetArray->GetName(), fTrkPtMax[1]);
1338               printf("Jet: E=%f, eta=%f, phi=%f, # of constituents=%d, nef=%f\n",jetsIncl[ij].E(),jetsIncl[ij].eta(), jetsIncl[ij].phi()*TMath::RadToDeg(), (Int_t)constituents.size(),jet->GetNEF());
1339               for(UInt_t ic=0; ic<constituents.size(); ic++)
1340                 {
1341                   if(constituents[ic].user_index()<0)
1342                     {
1343                       AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
1344                       printf("id = %d, cluster with pt = %f, ncell = %d\n",ic,constituents[ic].perp(), cluster->GetNCells());
1345                     }
1346                   else
1347                     {
1348                       AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1349                       printf("id = %d, track with pt = %f, track class = %d, originalPt = %f\n",ic,constituents[ic].perp(),(Int_t)track->GetTRDQuality(), track->GetIntegratedLength());
1350
1351                     }
1352                 }
1353               printf("==============================\n\n");
1354             }
1355         }
1356       // End of catching good high-E jets
1357     }
1358   if(fVerbosity>5) printf("%s has %d jets\n",fJetArray->GetName(), fJetArray->GetEntries());
1359 }
1360
1361 //________________________________________________________________________
1362 Int_t AliAnalysisTaskFullppJet::GetParentParton(const Int_t ipart)
1363 {
1364   //
1365   // Find the parent parton for a given MC particle by tracing back the parent
1366   // DUMMY; NEED improvement
1367
1368   return ipart;
1369 }
1370
1371
1372 //________________________________________________________________________
1373 Bool_t AliAnalysisTaskFullppJet::IsGoodJet(fastjet::PseudoJet jet, Double_t rad)
1374 {
1375   //
1376   // Check if the jet pt and direction fulfill the requirement
1377   //
1378
1379   if(jet.perp()<1) return kFALSE;
1380   if(TMath::Abs(jet.eta())>(0.7-rad)) return kFALSE;
1381   if(jet.phi() < (80*TMath::DegToRad()+rad) || jet.phi() > (180*TMath::DegToRad()-rad) ) return kFALSE;
1382   return kTRUE;
1383 }
1384
1385
1386 //________________________________________________________________________
1387 void AliAnalysisTaskFullppJet::ProcessMC(const Int_t r)
1388 {
1389   //
1390   // Main function for jet finding in MC
1391   //
1392
1393   Int_t npart = fMC->GetNumberOfTracks();
1394   fMcPartArray->Set(npart);
1395   Int_t countPart = 0;
1396   for(Int_t ipart=0; ipart<npart; ipart++)
1397     {
1398       AliVParticle* vParticle = fMC->GetTrack(ipart);
1399       if(!IsGoodMcPartilce(vParticle,ipart)) continue;
1400
1401       Int_t pdgCode = vParticle->PdgCode();
1402       if( fRejectNK && (pdgCode==130 || pdgCode==2112) ) continue;
1403
1404       if( fRejectWD && (pdgCode==310 || pdgCode==3112 || pdgCode==3122 || pdgCode==3222 || pdgCode==3312 || pdgCode==3322 || pdgCode==3334)) continue;
1405
1406       if( fChargedMC && vParticle->Charge()==0 ) continue;
1407
1408       Int_t index = 1;
1409       if(vParticle->Charge()==0) { index=-1; }
1410       fMcPartArray->AddAt(ipart, countPart);
1411       countPart++;
1412
1413       fTrueJetFinder[r]->AddInputVector(vParticle->Px(), vParticle->Py(), vParticle->Pz(), vParticle->P(), (ipart+1)*index);
1414       if(fVerbosity>10)
1415         printf("Input particle: ipart=%d, pdg=%d, species=%s, charge=%d, E=%4.3f\n",(ipart+1)*index,pdgCode,vParticle->GetName(), vParticle->Charge(),vParticle->P());
1416     }
1417   fMcPartArray->Set(countPart);
1418   fTrueJetFinder[r]->Run();
1419
1420   if(fMcTruthAntikt[r]) FillAODJets(fMcTruthAntikt[r], fTrueJetFinder[r], 1);
1421   if(fSaveQAHistos)     AnalyzeJets(fTrueJetFinder[r], 2, r);
1422   if(fRunUE && r==0)    RunAnalyzeUE(fTrueJetFinder[r], 2, 1);
1423
1424   // Run analysis on secondary particles
1425   if(fRunSecondaries && fAnaType==0)
1426     {
1427       for(Int_t i=0; i<2; i++)
1428         {
1429           AnalyzeSecondaryContribution(fTrueJetFinder[r],r,i);
1430           AnalyzeSecondaryContributionViaMatching(fTrueJetFinder[r],r,0,i);
1431           AnalyzeSecondaryContributionViaMatching(fTrueJetFinder[r],r,1,i);
1432         }
1433     }
1434 }
1435
1436
1437 //________________________________________________________________________
1438 void AliAnalysisTaskFullppJet::AnalyzeJets(AliFJWrapper *jetFinder, const Int_t type,  const Int_t r)
1439 {
1440   //
1441   // Fill all the QA plots for jets
1442   // Especailly all the constituents plot must be filled here since the information
1443   // is not available in the output AOD
1444
1445   Double_t vertex[3] = {0, 0, 0};
1446   fESD->GetVertex()->GetXYZ(vertex) ;
1447   TLorentzVector gamma;
1448
1449   const Int_t nBins = fJetEnergyFraction[type][r]->GetAxis(1)->GetNbins();
1450   Float_t radiusCut[nBins];
1451   Float_t eFraction[nBins];
1452   Int_t   nPart[nBins];
1453   for(Int_t i=0; i<nBins; i++)
1454     {
1455       radiusCut[i] = (fJetEnergyFraction[type][r]->GetAxis(1)->GetBinWidth(1)) * (i+1);
1456       eFraction[i] = 0.;
1457       nPart[nBins] = 0;
1458     }
1459   std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1460
1461   for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1462     {
1463       if(type<2 && fCheckTPCOnlyVtx && fIsTPCOnlyVtx)
1464         {
1465           fhJetInTPCOnlyVtx[type][r]->Fill(jetsIncl[ij].perp(),jetsIncl[ij].phi()*TMath::RadToDeg(),jetsIncl[ij].eta());
1466         }
1467       if(!IsGoodJet(jetsIncl[ij],kRadius[r])) continue; // Fidiual cut
1468       Float_t jetEta = jetsIncl[ij].eta();
1469       Float_t jetPhi = jetsIncl[ij].phi();
1470       Float_t jetE   = jetsIncl[ij].E();
1471       Float_t jetPt  = jetsIncl[ij].perp();
1472
1473       std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
1474       Double_t neE=0, leadingZ = 0, maxPt = 0;
1475       Int_t constituentType = -1, leadingIndex = 0;
1476       for(UInt_t ic=0; ic<constituents.size(); ic++)
1477         {
1478           if(constituents[ic].perp()>maxPt)
1479             {
1480               maxPt = constituents[ic].perp();
1481               leadingIndex = constituents[ic].user_index();
1482             }
1483
1484           if(constituents[ic].user_index()<0)
1485             {
1486               neE += constituents[ic].E();
1487               constituentType = 3;
1488             }
1489           else
1490             {
1491               if(type==2) constituentType = 0;
1492               else
1493                 {
1494                   AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1495                   constituentType = (Int_t)track->GetTRDQuality();
1496                 }
1497             }
1498           Double_t cz = GetZ(constituents[ic].px(),constituents[ic].py(),constituents[ic].pz(),jetsIncl[ij].px(),jetsIncl[ij].py(),jetsIncl[ij].pz());
1499           fhJetPtVsZ[type][r]->Fill(jetPt,cz, constituentType);
1500           if(cz>leadingZ) leadingZ = cz;
1501         }
1502
1503       if(type<2 && leadingIndex<0)
1504         {
1505           AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(leadingIndex*(-1)-1);
1506           fhFcrossVsZleading[type][r]->Fill(jetPt,GetExoticEnergyFraction(cluster),leadingZ);
1507         }
1508
1509       if(leadingZ > 0.98) continue;  // Z cut
1510
1511       fJetCount[type][r]->Fill(jetPt);
1512       if(type==1 && fIsExoticEvent3GeV) fhJetPtInExoticEvent[0][r]->Fill(jetPt);
1513       if(type==1 && fIsExoticEvent5GeV) fhJetPtInExoticEvent[1][r]->Fill(jetPt);
1514
1515       Double_t totTrkPt[3] = {0.,0.,0.};
1516       Double_t chPt = 0;
1517       for(Int_t i=0; i<nBins; i++) { eFraction[i] = 0.; nPart[i] = 0;}
1518       Double_t mcSubE=0.;
1519       Double_t subClsE[5] = {0,0,0,0,0};
1520       Double_t subTrkPtClean[5] = {0,0,0,0,0};
1521       Double_t subTrkPtAmbig[5] = {0,0,0,0,0};
1522       Double_t fraction[5] = {270,0.3,0.5,0.7,1};
1523       Double_t leadChPt=0., leadNePt=0.;
1524       Int_t leadChIndex = -1;
1525       for(UInt_t ic=0; ic<constituents.size(); ic++)
1526         {
1527           Float_t partEta = constituents[ic].eta();
1528           Float_t partPhi = constituents[ic].phi();
1529           Float_t partE   = constituents[ic].E();
1530           Float_t partPt  = constituents[ic].perp();
1531
1532           if(constituents[ic].user_index()<0)
1533             {
1534               fhNeutralPtInJet[type][r]->Fill(jetPt, partPt);
1535               if(partPt>leadNePt)
1536                 leadNePt = partPt;
1537               if(type<2)
1538                 {
1539                   AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
1540                   Double_t clsE = cluster->E();
1541                   if(cluster->Chi2()>0.5) fhTrigNeuPtInJet[type][r]->Fill(jetPt, partPt);
1542                   if(fStudySubEInHC || fStudyMcOverSubE)
1543                     {
1544                       cluster->GetMomentum(gamma, vertex);
1545                       Double_t sinTheta = TMath::Sqrt(1-TMath::Power(gamma.CosTheta(),2));
1546
1547                       Double_t subEtmp = cluster->GetDispersion();
1548                       Double_t mipETmp = cluster->GetEmcCpvDistance();
1549                       mcSubE += cluster->GetDistanceToBadChannel()*sinTheta;
1550                       subClsE[0]   += ((mipETmp>clsE)?clsE:mipETmp)*sinTheta;
1551                       for(Int_t j=1; j<5; j++)
1552                         {
1553                           subClsE[j]   += ((fraction[j]*subEtmp>clsE)?clsE:fraction[j]*subEtmp)*sinTheta;
1554                         }
1555                     }
1556                 }
1557             }
1558           else
1559             {
1560               if(partPt>leadChPt) {leadChPt = partPt;leadChIndex=ic;}
1561               fhChargedPtInJet[type][r]->Fill(jetPt, partPt);
1562               chPt += constituents[ic].perp();
1563               if(type<2)
1564                 {
1565                   AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1566                   Int_t trkType = (Int_t)track->GetTRDQuality();
1567                   totTrkPt[trkType] += partPt;
1568                   Int_t clusterIndex = track->GetEMCALcluster();
1569                   Int_t clusterPos = -1;
1570                   Bool_t isExist = kFALSE;
1571                   for(Int_t j=0; j<fClusterArray->GetEntriesFast(); j++)
1572                     {
1573                       AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(j);
1574                       if( clusterIndex == cluster->GetID() )
1575                         {
1576                           isExist = kTRUE;
1577                           clusterPos = j;
1578                           break;
1579                         }
1580                     }
1581                   if(isExist)
1582                     {
1583                       AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(clusterPos);
1584                       Double_t subEtmp = cluster->GetDispersion();
1585                       Double_t mipETmp = cluster->GetEmcCpvDistance();
1586                       for(Int_t k=0; k<5; k++)
1587                         {
1588                           if(k==0)
1589                             {
1590                               if(mipETmp>cluster->E()) subTrkPtClean[k] += partPt;
1591                               else             subTrkPtAmbig[k] += partPt;
1592                             }
1593                           else
1594                             {
1595                               if(fraction[k]*subEtmp>cluster->E()) subTrkPtClean[k] += partPt;
1596                               else                         subTrkPtAmbig[k] += partPt;
1597                             }
1598                         }
1599                     }
1600                 }
1601             }
1602
1603           Float_t dR = TMath::Sqrt( (partEta-jetEta)*(partEta-jetEta) + (partPhi-jetPhi)*(partPhi-jetPhi) );
1604           for(Int_t i=0; i<nBins; i++)
1605             {
1606               if( dR < radiusCut[i] )
1607                 {
1608                   eFraction[i] += partE;
1609                   nPart[i]++;
1610                 }
1611             }
1612         }
1613
1614       fhLeadNePtInJet[type][r]->Fill(jetPt, leadNePt);
1615       fhLeadChPtInJet[type][r]->Fill(jetPt, leadChPt);
1616       if(type<2 && leadChIndex>-1)
1617         {
1618           Double_t cz = GetZ(constituents[leadChIndex].px(),constituents[leadChIndex].py(),constituents[leadChIndex].pz(),jetsIncl[ij].px(),jetsIncl[ij].py(),jetsIncl[ij].pz());
1619           fhChLeadZVsJetPt[type][r]->Fill(jetPt, cz);
1620         }
1621
1622       if((fStudySubEInHC||fStudyMcOverSubE) && type<2)
1623         {
1624           for(Int_t i=0; i<5; i++)
1625             {
1626               if(fStudySubEInHC)
1627                 {
1628                   fhSubClsEVsJetPt[type][r][i]->Fill(jetPt-subClsE[4],subClsE[i]/(jetPt-subClsE[4]));
1629                   fhHCTrkPtClean[type][r][i]->Fill(jetPt-subClsE[4],subTrkPtClean[i]/(jetPt-subClsE[4]));
1630                   fhHCTrkPtAmbig[type][r][i]->Fill(jetPt-subClsE[4],subTrkPtAmbig[i]/(jetPt-subClsE[4]));
1631                 }
1632               if(type==0 && fIsMC && fStudyMcOverSubE)
1633                 {
1634                   fHCOverSubE[r][i]->Fill(jetPt-subClsE[4],subClsE[i]-mcSubE);
1635                   fHCOverSubEFrac[r][i]->Fill(jetPt-subClsE[4],(subClsE[i]-mcSubE)/(jetPt-subClsE[4]));
1636                 }
1637             }
1638         }
1639       if(type<2 && chPt>0)
1640         {
1641           for(Int_t i=0; i<3; i++)
1642             {
1643               fRelTrkCon[type][r]->Fill(jetPt,totTrkPt[i]/chPt,i);
1644             }
1645         }
1646
1647
1648       for(Int_t ibin=0; ibin<nBins; ibin++)
1649         {
1650           Double_t fill1[3] = {jetPt,radiusCut[ibin]-0.005,eFraction[ibin]/jetE};
1651           fJetEnergyFraction[type][r]->Fill(fill1);
1652           Double_t fill2[3] = {jetPt,radiusCut[ibin]-0.005,nPart[ibin]};
1653           fJetNPartFraction[type][r]->Fill(fill2);
1654         }
1655
1656       // Get the jet pt containing tracks or clusters above some threshold
1657       if(type<2)
1658         {
1659           Double_t thres[3] = {15,25,40};
1660           Int_t okCh[3] = {0,0,0};
1661           Int_t okNe[3] = {0,0,0};
1662           Double_t lowPt[2] = {0.,0.};
1663           for(UInt_t ic=0; ic<constituents.size(); ic++)
1664             {
1665               Float_t partPt  = constituents[ic].perp();
1666
1667               if(partPt<0.3) lowPt[0] += partPt;
1668               else if(partPt<0.5) lowPt[1] += partPt;
1669
1670               if(constituents[ic].user_index()>0)
1671                 {
1672                   for(Int_t it=0; it<3; it++)
1673                     {
1674                       if(partPt>thres[it]) okCh[it] = 1;
1675                     }
1676                 }
1677               else
1678                 {
1679                   for(Int_t icl=0; icl<3; icl++)
1680                     {
1681                       if(partPt>thres[icl]) okNe[icl] = 1;
1682                     }
1683                 }
1684             }
1685           for(Int_t i=0; i<3; i++)
1686             {
1687               if(okCh[i]==1)
1688                 fhJetPtWithTrkThres[type][r]->Fill(i,jetPt);
1689               if(okNe[i]==1)
1690                 fhJetPtWithClsThres[type][r]->Fill(i,jetPt);
1691             }
1692           for(Int_t k=0; k<2; k++) fhJetPtVsLowPtCons[type][r][k]->Fill(jetPt,lowPt[k]/jetPt);
1693         }
1694     }
1695 }
1696
1697 //________________________________________________________________________
1698 void AliAnalysisTaskFullppJet::AnalyzeSecondaryContribution(AliFJWrapper *jetFinder, const Int_t r, const Int_t etaCut)
1699 {
1700   //
1701   // Analyze secondaries
1702   //
1703
1704   std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1705   for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1706     {
1707       Float_t jetEta = jetsIncl[ij].eta();
1708       Float_t jetPt  = jetsIncl[ij].perp();
1709       if(TMath::Abs(jetEta)>0.5*etaCut+0.5 || jetPt<1) continue;
1710       std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
1711       Double_t dNKPt = 0, dWPPt = 0, weakPt=0, nkPt=0;
1712       Int_t type = -1;
1713       for(UInt_t ic=0; ic<constituents.size(); ic++)
1714         {
1715           Int_t ipart = TMath::Abs(constituents[ic].user_index())-1;
1716           AliMCParticle* mcParticle = (AliMCParticle*)fMC->GetTrack(ipart);
1717           if(!mcParticle) continue;
1718           Int_t pdg = mcParticle->PdgCode();
1719           if(pdg==310 || (pdg<=3400 && pdg>=3100)) weakPt += mcParticle->Pt();
1720           if(pdg==130 || pdg==2112)  nkPt += mcParticle->Pt();
1721
1722           // Weak decaying particles
1723           if(pdg==310) type=0;
1724           else if (pdg==3122) type=1;
1725           else if (pdg<=3400 && pdg>=3100) type=2;
1726           else type=-1;
1727           if(type>-1)
1728             {
1729               Int_t binx = fhSecondaryResponse[type]->GetXaxis()->FindFixBin(mcParticle->Pt());
1730               TH1F *htmp = (TH1F*)fhSecondaryResponse[type]->ProjectionY(Form("hpro_%d_%d_%d",(Int_t)Entry(),ij,ic),binx,binx);
1731               Double_t pro = htmp->GetRandom();
1732               dWPPt += pro*mcParticle->Pt() - mcParticle->Pt();
1733               delete htmp;
1734             }
1735
1736           // Missing neutron and K0L
1737           if(pdg==130 || pdg==2112) dNKPt -= mcParticle->Pt();
1738         }
1739
1740       fhNKFracVsJetPt[etaCut][r]->Fill(jetPt,nkPt/jetPt);
1741       fhWeakFracVsJetPt[etaCut][r]->Fill(jetPt,weakPt/jetPt);
1742
1743       Double_t jetNewWPPt = jetPt + dWPPt;
1744       fhJetResponseWP[etaCut][r]->Fill(jetPt,jetNewWPPt);
1745       fhJetResolutionWP[etaCut][r]->Fill(jetPt,(jetPt-jetNewWPPt)/jetPt);
1746
1747       Double_t jetNewNKPt = jetPt + dNKPt;
1748       fhJetResponseNK[etaCut][r]->Fill(jetPt,jetNewNKPt);
1749       fhJetResolutionNK[etaCut][r]->Fill(jetPt,(jetPt-jetNewNKPt)/jetPt);
1750     }
1751 }
1752
1753
1754 //________________________________________________________________________
1755 void AliAnalysisTaskFullppJet::AnalyzeSecondaryContributionViaMatching(AliFJWrapper *jetFinder, const Int_t r, const Int_t type, const Int_t etaCut)
1756 {
1757   //
1758   // Estimate the contribution of missing energies via matching
1759   // type = 0 -> NK
1760   // type = 1 -> WP
1761   //
1762
1763   if(type!=0 && type!=1) return;
1764
1765   AliFJWrapper fj("fj","fj");
1766   fj.CopySettingsFrom(*jetFinder);
1767
1768   Int_t npart = fMC->GetNumberOfTracks();
1769   Int_t pType = -1;
1770   for(Int_t ipart=0; ipart<npart; ipart++)
1771     {
1772       AliVParticle* vParticle = fMC->GetTrack(ipart);
1773       if(!IsGoodMcPartilce(vParticle,ipart)) continue;
1774       Int_t pdgCode = vParticle->PdgCode();
1775       Double_t pt = vParticle->Pt();
1776       if( type==0 && (pdgCode==130 || pdgCode==2112) ) continue; // Reject NK
1777       if(type==1)
1778         {
1779           if(pdgCode==310) pType=0;
1780           else if (pdgCode==3122) pType=1;
1781           else if (pdgCode<=3400 && pdgCode>=3100) pType=2;
1782           else pType=-1;
1783           if(pType>-1)
1784             {
1785               Int_t binx = fhSecondaryResponse[pType]->GetXaxis()->FindFixBin(vParticle->Pt());
1786               TH1F *htmp = (TH1F*)fhSecondaryResponse[pType]->ProjectionY(Form("hpro_%d_%d",(Int_t)Entry(),ipart),binx,binx);
1787               Double_t pro = htmp->GetRandom();
1788               pt = pro * vParticle->Pt();
1789               delete htmp;
1790             }
1791         }
1792       Int_t index = 1;
1793       if(vParticle->Charge()==0) { index=-1; }
1794       fj.AddInputVector(vParticle->Px()*pt/vParticle->Pt(), vParticle->Py()*pt/vParticle->Pt(), vParticle->Pz(), vParticle->P(), (ipart+1)*index);
1795     }
1796
1797   fj.Run();
1798   std::vector<fastjet::PseudoJet> jets_incl_mc = fj.GetInclusiveJets();
1799   std::vector<fastjet::PseudoJet> jets_incl_mc_std = jetFinder->GetInclusiveJets();
1800
1801   for(UInt_t ij=0; ij<jets_incl_mc.size(); ij++)
1802     {
1803       Float_t jetEta = jets_incl_mc[ij].eta();
1804       Float_t jetPt  = jets_incl_mc[ij].perp();
1805       if(TMath::Abs(jetEta)>0.5*etaCut+0.5 || jetPt<5) continue;
1806       Double_t dEta, dPhi;
1807       Int_t index = FindSpatialMatchedJet(jets_incl_mc[ij], jetFinder, dEta, dPhi, 1);
1808       if(index>-1)
1809         {
1810           Double_t jetPtStd = jets_incl_mc_std[index].perp();
1811           Double_t dR = TMath::Sqrt(dEta*dEta+dPhi*dPhi);
1812           if(type==0) fhJetResolutionNKSM[etaCut][r]->Fill(jetPtStd,(jetPtStd-jetPt)/jetPtStd,dR);
1813           if(type==1) fhJetResolutionWPSM[etaCut][r]->Fill(jetPtStd,(jetPtStd-jetPt)/jetPtStd,dR);
1814
1815           if(dR<kdRCut[r])
1816             {
1817               if(type==0) fhJetResponseNKSM[etaCut][r]->Fill(jetPtStd,jetPt);
1818               if(type==1) fhJetResponseWPSM[etaCut][r]->Fill(jetPtStd,jetPt);
1819             }
1820         }
1821     }
1822 }
1823
1824 //________________________________________________________________________
1825 void AliAnalysisTaskFullppJet::RunAnalyzeUE(AliFJWrapper *jetFinder, const Int_t type, const Bool_t isMCTruth)
1826 {
1827   //
1828   // Run analysis to estimate the underlying event
1829   // Adapted from Oliver
1830
1831   std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1832   Double_t leadpt=0, leadPhi = -999, leadEta = -999;
1833   Int_t leadIndex = -1;
1834   for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1835     {
1836       Double_t jetEta = jetsIncl[ij].eta();
1837       Double_t jetPt = jetsIncl[ij].perp();
1838       Double_t jetPhi = jetsIncl[ij].phi();
1839       if(leadpt<jetPt)
1840         {
1841           leadpt = jetPt;
1842           leadPhi = jetPhi;
1843           leadEta = jetEta;
1844           leadIndex = ij;
1845         }
1846     }
1847   if(leadpt<1 || TMath::Abs(leadEta)>0.5 ) return;
1848
1849   Int_t backIndex = -1;
1850   Bool_t isBackToBack = kFALSE;
1851   for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1852     {
1853       Double_t dPhi = TMath::Abs(jetsIncl[ij].phi()-leadPhi);
1854       if(dPhi > kPI) dPhi = 2*kPI - dPhi;
1855       if(dPhi>150*TMath::DegToRad() && jetsIncl[ij].perp()/leadpt > 0.8)
1856         {
1857           backIndex = ij;
1858           isBackToBack = kTRUE;
1859         }
1860     }
1861
1862   for(Int_t ij=0; ij<(Int_t)jetsIncl.size(); ij++)
1863     {
1864       if(ij==leadIndex || ij==backIndex) continue;
1865       if(TMath::Abs(jetsIncl[ij].eta())>0.5) continue;
1866       if(jetsIncl[ij].perp()>15) isBackToBack = kFALSE;
1867     }
1868
1869
1870   Double_t axis[2] = {leadPhi+0.5 * kPI, leadPhi-0.5 * kPI};
1871   if(axis[0]>2*kPI) axis[0] -= 2*kPI;
1872   if(axis[1]<0) axis[1] += 2*kPI;
1873
1874   fhUEJetPtNorm[type][0][0]->Fill(leadpt);
1875   if(isBackToBack)  fhUEJetPtNorm[type][0][1]->Fill(leadpt);
1876
1877   if(!isMCTruth)
1878     {
1879       Int_t ntracks = fTrackArray->GetEntriesFast();
1880       Double_t vertex[3] = {0, 0, 0};
1881       fESD->GetVertex()->GetXYZ(vertex) ;
1882       TLorentzVector gamma;
1883       Int_t nclusters = fClusterArray->GetEntriesFast();
1884
1885       for(Int_t j=0; j<2; j++)
1886         {
1887           Double_t ueCh = 0, ueChNe = 0.;
1888           for(Int_t it=0; it<ntracks; it++)
1889             {
1890               AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it);
1891               if(!t) continue;
1892               Double_t dPhi = TMath::Abs(axis[j]-t->Phi());
1893               Double_t dEta = TMath::Abs(leadEta-t->Eta());
1894               if(dPhi > kPI) dPhi = 2*kPI - dPhi;
1895               if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4)
1896                 {
1897                   fhUEJetPtVsConsPt[type][0][0]->Fill(leadpt,t->Pt());
1898                   if(isBackToBack) fhUEJetPtVsConsPt[type][0][1]->Fill(leadpt,t->Pt());
1899                   ueCh += t->Pt();
1900                   if(TMath::Abs(leadEta-0.4)<0.7 && axis[j]>80*TMath::DegToRad()+0.4 && axis[j]<180*TMath::DegToRad()-0.4)
1901                     {
1902                       fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,t->Pt());
1903                       if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,t->Pt());
1904                       ueChNe  += t->Pt();
1905                     }
1906                 }
1907             }
1908           fhUEJetPtVsSumPt[type][0][0]->Fill(leadpt,ueCh);
1909           if(isBackToBack) fhUEJetPtVsSumPt[type][0][1]->Fill(leadpt,ueCh);
1910
1911           if(!(TMath::Abs(leadEta-0.4)<0.7 && axis[j]>80*TMath::DegToRad()+0.4 && axis[j]<180*TMath::DegToRad()-0.4)) continue;
1912           fhUEJetPtNorm[type][1][0]->Fill(leadpt);
1913           if(isBackToBack) fhUEJetPtNorm[type][1][1]->Fill(leadpt);
1914           for(Int_t ic=0; ic<nclusters; ic++)
1915             {
1916               AliESDCaloCluster * cl = (AliESDCaloCluster *)fClusterArray->At(ic);
1917               if(!cl) continue;
1918               cl->GetMomentum(gamma, vertex);
1919               Double_t clsPhi = gamma.Phi();
1920               if(clsPhi<0) clsPhi += 2*kPI;
1921               Double_t dPhi = TMath::Abs(axis[j]-clsPhi);
1922               Double_t dEta = TMath::Abs(leadEta-gamma.Eta());
1923               if(dPhi > kPI) dPhi = 2*kPI - dPhi;
1924               if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4)
1925                 {
1926                   ueChNe  += gamma.Pt();
1927                   fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,gamma.Pt());
1928                   if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,gamma.Pt());
1929                 }
1930             }
1931           fhUEJetPtVsSumPt[type][1][0]->Fill(leadpt,ueChNe);
1932           if(isBackToBack)  fhUEJetPtVsSumPt[type][1][1]->Fill(leadpt,ueChNe);
1933         }
1934     }
1935   else
1936     {
1937       fhUEJetPtNorm[type][1][0]->Fill(leadpt);
1938       if(isBackToBack) fhUEJetPtNorm[type][1][1]->Fill(leadpt);
1939       Int_t npart = fMcPartArray->GetSize();
1940       for(Int_t j=0; j<2; j++)
1941         {
1942           Double_t ueCh = 0, ueChNe = 0., ueChNe2 = 0.;
1943           for(Int_t ipos=0; ipos<npart; ipos++)
1944             {
1945               AliVParticle* vParticle = fMC->GetTrack(fMcPartArray->At(ipos));
1946               if(!vParticle) continue;
1947               Double_t dPhi = TMath::Abs(axis[j]-vParticle->Phi());
1948               Double_t dEta = TMath::Abs(leadEta-vParticle->Eta());
1949               if(dPhi > kPI) dPhi = 2*kPI - dPhi;
1950               if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4)
1951                 {
1952                   Double_t pt = vParticle->Pt();
1953                   ueChNe += pt;
1954                   fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,pt);
1955                   if(isBackToBack)  fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,pt);
1956                   if( vParticle->Charge()!=0 )
1957                     {
1958                       fhUEJetPtVsConsPt[type][0][0]->Fill(leadpt,pt);
1959                       if(isBackToBack) fhUEJetPtVsConsPt[type][0][1]->Fill(leadpt,pt);
1960                       ueCh += pt;
1961                       ueChNe2 += pt;
1962                     }
1963                   else
1964                     {
1965                       if(pt>0.5) ueChNe2 += pt;
1966                     }
1967                 }
1968             }
1969           fhUEJetPtVsSumPt[type][0][0]->Fill(leadpt,ueCh);
1970           fhUEJetPtVsSumPt[type][1][0]->Fill(leadpt,ueChNe);
1971           if(isBackToBack) fhUEJetPtVsSumPt[type][0][1]->Fill(leadpt,ueCh);
1972           if(isBackToBack) fhUEJetPtVsSumPt[type][1][1]->Fill(leadpt,ueChNe);
1973         }
1974     }
1975 }
1976
1977 //________________________________________________________________________
1978 Bool_t AliAnalysisTaskFullppJet::IsGoodJet(AliAODJet *jet, Double_t rad)
1979 {
1980   //
1981   // Check if it is a good jet
1982   //
1983
1984   if(jet->Pt()<1) return kFALSE;
1985   if(TMath::Abs(jet->Eta())>(0.7-rad)) return kFALSE;
1986   if(jet->Phi() < (80*TMath::DegToRad()+rad) || jet->Phi() > (180*TMath::DegToRad()-rad) ) return kFALSE;
1987   return kTRUE;
1988 }
1989
1990
1991 //________________________________________________________________________
1992 Double_t AliAnalysisTaskFullppJet::GetLeadingZ(const Int_t jetIndex, AliFJWrapper *jetFinder)
1993 {
1994   //
1995   // Get the leading z of the input jet
1996   //
1997
1998   Double_t z = 0;
1999   std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
2000   std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
2001   Int_t index = -1;
2002   Double_t maxPt = 0;
2003   for(UInt_t ic=0; ic<constituents.size(); ic++)
2004     {
2005       if(constituents[ic].perp()>maxPt)
2006         {
2007           maxPt = constituents[ic].perp();
2008           index = ic;
2009         }
2010     }
2011   if(index>-1)
2012     z = GetZ(constituents[index].px(),constituents[index].py(),constituents[index].pz(),jetsIncl[jetIndex].px(),jetsIncl[jetIndex].py(),jetsIncl[jetIndex].pz());
2013   return z;
2014 }
2015
2016 //________________________________________________________________________
2017 Double_t AliAnalysisTaskFullppJet::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
2018 {
2019   //
2020   // Get the z of a constituent inside of a jet
2021   //
2022
2023   return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/(jetPx*jetPx+jetPy*jetPy+jetPz*jetPz);
2024 }
2025
2026 //________________________________________________________________________
2027 Double_t AliAnalysisTaskFullppJet::GetMeasuredJetPtResolution(const Int_t jetIndex, AliFJWrapper *jetFinder)
2028 {
2029   //
2030   // Get jet energy resoultion due to intrinsic detector effects
2031   //
2032
2033   Double_t jetSigma2 = 0;
2034   std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
2035   for(UInt_t ic=0; ic<constituents.size(); ic++)
2036     {
2037       if(constituents[ic].user_index()>0)
2038         {
2039           AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
2040           jetSigma2 += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
2041         }
2042       else
2043         {
2044           AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
2045           jetSigma2 += TMath::Power(cluster->GetTOF(),2);
2046         }
2047     }
2048   return TMath::Sqrt(jetSigma2);
2049 }
2050
2051
2052
2053 //________________________________________________________________________
2054 Double_t AliAnalysisTaskFullppJet::GetJetMissPtDueToTrackingEfficiency(const Int_t jetIndex, AliFJWrapper *jetFinder, const Int_t radiusIndex)
2055 {
2056   //
2057   // Correct the tracking inefficiency explicitly
2058   //
2059
2060   Double_t misspt = 0;
2061   std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
2062   Double_t jetPt = jetsIncl[jetIndex].perp();
2063   if(!fhCorrTrkEffPtBin[fTriggerType][radiusIndex])
2064     {
2065       printf("Warning: can't get the mean # of tracks per jet with pt=%f in: trigger=%d, radiusIndex=%d\n",jetPt,fTriggerType,radiusIndex);
2066       return 0;
2067     }
2068   Int_t ibin = fhCorrTrkEffPtBin[fTriggerType][radiusIndex]->FindFixBin(jetPt);
2069   if(!fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1] || fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1]->Integral()<0.001)
2070     {
2071       printf("Warning: no sampling distrubtion for jet with pt=%f\n",jetPt);
2072       return 0;
2073     }
2074   Double_t ntrack = fhCorrTrkEffPtBin[fTriggerType][radiusIndex]->GetBinContent(ibin);
2075   Int_t nTrk = (Int_t) ntrack;
2076   Double_t res = ntrack-nTrk;
2077   Double_t pro1 = fRandomGen->Uniform();
2078   if(pro1<res) nTrk++;
2079   for(Int_t itry=0; itry<nTrk; itry++)
2080     {
2081       Double_t trkPt = fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1]->GetRandom();
2082       if(trkPt/jetPt>fTrkEffCorrCutZ) continue;
2083       Double_t eff = GetTrkEff(trkPt);
2084       Double_t pro = fRandomGen->Uniform();
2085       if(pro>eff) misspt += trkPt;
2086     }
2087   return misspt;
2088 }
2089
2090
2091 //________________________________________________________________________
2092 Bool_t AliAnalysisTaskFullppJet::IsGoodMcPartilce(const AliVParticle* vParticle, const Int_t ipart)
2093 {
2094   //
2095   // Select good primary particles to feed into jet finder
2096   //
2097
2098   if(!vParticle) return kFALSE;
2099   if(!fMC->IsPhysicalPrimary(ipart)) return kFALSE;
2100   if (TMath::Abs(vParticle->Eta())>2) return kFALSE;
2101   return kTRUE;
2102 }
2103
2104 //________________________________________________________________________
2105 Int_t AliAnalysisTaskFullppJet::FindSpatialMatchedJet(fastjet::PseudoJet jet, AliFJWrapper *jetFinder, Double_t &dEta, Double_t &dPhi, Double_t maxR)
2106 {
2107   //
2108   // Find spatially matched detector and particle jets
2109   //
2110
2111   Int_t index=-1;
2112   dEta=-999, dPhi=-999;
2113   std::vector<fastjet::PseudoJet> jets_incl = jetFinder->GetInclusiveJets();
2114   for(UInt_t ij=0; ij<jets_incl.size(); ij++)
2115     {
2116       if(jets_incl[ij].perp()<5) continue;
2117       if(TMath::Abs(jets_incl[ij].eta())>1) continue;
2118       Double_t tmpR = TMath::Sqrt( TMath::Power(jet.eta()-jets_incl[ij].eta(), 2) + TMath::Power(jet.phi()-jets_incl[ij].phi(), 2) );
2119       if(tmpR<maxR)
2120         {
2121           maxR=tmpR;
2122           index=ij;
2123           dEta = jet.eta()-jets_incl[ij].eta();
2124           dPhi = jet.phi()-jets_incl[ij].phi();
2125         }
2126     }
2127   return index;
2128 }
2129
2130 //________________________________________________________________________
2131 Int_t AliAnalysisTaskFullppJet::FindEnergyMatchedJet(AliFJWrapper *jetFinder1, const Int_t index1, AliFJWrapper *jetFinder2, const Double_t fraction)
2132 {
2133   //
2134   // Find matched detector and particle jets based on shared constituents
2135   //
2136
2137   Int_t matchedIndex = -1;
2138   std::vector<fastjet::PseudoJet> jetsIncl1 = jetFinder1->GetInclusiveJets();
2139   std::vector<fastjet::PseudoJet> jetsIncl2 = jetFinder2->GetInclusiveJets();
2140   std::vector<fastjet::PseudoJet> constituents1 = jetFinder1->GetJetConstituents(index1);
2141   Double_t jetPt1 = jetsIncl1[index1].perp();
2142   if(jetPt1<0) return matchedIndex;
2143
2144   for(UInt_t ij2=0; ij2<jetsIncl2.size(); ij2++)
2145     {
2146       Double_t jetPt2  = jetsIncl2[ij2].perp();
2147       if(jetPt2<0) return matchedIndex;
2148       std::vector<fastjet::PseudoJet> constituents2 = jetFinder2->GetJetConstituents(ij2);
2149       Double_t sharedPt1 = 0., sharedPt2 = 0.;
2150       for(UInt_t ic2=0; ic2<constituents2.size(); ic2++)
2151         {
2152           Int_t mcLabel = constituents2[ic2].user_index()-1;
2153           Double_t consPt2 = constituents2[ic2].perp();
2154           for(UInt_t ic1=0; ic1<constituents1.size(); ic1++)
2155             {
2156               Double_t consPt1 = constituents1[ic1].perp();
2157               if(constituents1[ic1].user_index()>0)
2158                 {
2159                   AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents1[ic1].user_index()-1);
2160                   if(track->GetLabel()==mcLabel)  {sharedPt2 += consPt2; sharedPt1 += consPt1; cout<<"Found a matched track"<<endl;break;}
2161                 }
2162               else
2163                 {
2164                   AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents1[ic1].user_index()*(-1)-1);
2165                   if(cluster->GetLabel()==mcLabel) {sharedPt2 += consPt2; sharedPt1 += consPt1; cout<<"Found a matched cluster"<<endl;break;}
2166                 }
2167             }
2168         }
2169       cout<<sharedPt1/jetPt1<<"  "<<sharedPt2/jetPt2<<endl;
2170       if(sharedPt1/jetPt1 > fraction && sharedPt2/jetPt2 > fraction)
2171         {
2172           matchedIndex = ij2;
2173           break;
2174         }
2175     }
2176   return matchedIndex;
2177 }
2178
2179 //________________________________________________________________________
2180 void AliAnalysisTaskFullppJet::GetESDTrax()
2181 {
2182   //
2183   // Get esd tracks.
2184   //
2185
2186   Int_t nTrax = fESD->GetNumberOfTracks();
2187   if(fVerbosity>5) printf("[i] # of tracks in event: %d\n",nTrax);
2188
2189   for (Int_t i = 0; i < nTrax; ++i)
2190     {
2191       AliESDtrack* esdtrack = fESD->GetTrack(i);
2192       if (!esdtrack)
2193         {
2194           AliError(Form("Couldn't get ESD track %d\n", i));
2195           continue;
2196         }
2197
2198       AliESDtrack *newtrack = GetAcceptTrack(esdtrack);
2199       if(!newtrack) continue;
2200
2201       Double_t pt  = newtrack->Pt();
2202       Double_t eta = newtrack->Eta();
2203       if (pt < fTrkPtMin[fTriggerType] || pt > fTrkPtMax[fTriggerType] || TMath::Abs(eta) > fTrkEtaMax)
2204         {
2205           delete newtrack;
2206           continue;
2207         }
2208
2209       if(fSysTrkEff)
2210         {
2211           Double_t rand = fRandomGen->Uniform();
2212           if(rand<fVaryTrkEff)
2213             {
2214               delete newtrack;
2215               continue;
2216             }
2217         }
2218       newtrack->SetIntegratedLength(esdtrack->Pt());
2219       fTrackArray->Add(newtrack);
2220     }
2221   if(fVerbosity>5) printf("[i] # of tracks in event: %d\n", fTrackArray->GetEntries());
2222   fNTracksPerChunk += fTrackArray->GetEntries();
2223 }
2224
2225 //
2226 //________________________________________________________________________
2227 //
2228 AliESDtrack *AliAnalysisTaskFullppJet::GetAcceptTrack(AliESDtrack *esdtrack)
2229 {
2230   //
2231   // Get the hybrid tracks
2232   //
2233
2234   AliESDtrack *newTrack = 0x0;
2235
2236   if(fTrackCutsType==0 || fTrackCutsType==3)
2237     {
2238       if(fEsdTrackCuts->AcceptTrack(esdtrack))
2239         {
2240           newTrack = new AliESDtrack(*esdtrack);
2241           newTrack->SetTRDQuality(0);
2242         }
2243       else if(fHybridTrackCuts1->AcceptTrack(esdtrack))
2244         {
2245           if(esdtrack->GetConstrainedParam())
2246             {
2247               newTrack = new AliESDtrack(*esdtrack);
2248               const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
2249               newTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
2250               newTrack->SetTRDQuality(1);
2251             }
2252           else
2253             return 0x0;
2254         }
2255       else if(fHybridTrackCuts2->AcceptTrack(esdtrack))
2256         {
2257           if(esdtrack->GetConstrainedParam())
2258             {
2259               newTrack = new AliESDtrack(*esdtrack);
2260               const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
2261               newTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
2262               newTrack->SetTRDQuality(2);
2263             }
2264           else
2265             return 0x0;
2266         }
2267       else
2268         {
2269           return 0x0;
2270         }
2271     }
2272   else if(fTrackCutsType==1)
2273     {
2274       if(fEsdTrackCuts->AcceptTrack(esdtrack))
2275         {
2276           newTrack = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(fESD),esdtrack->GetID());// use TPC only tracks with non default SPD vertex
2277           if(!newTrack) return 0x0;
2278           AliExternalTrackParam exParam;
2279           Bool_t relate = newTrack->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam); //constrain to SPD vertex
2280           if( !relate )
2281             {
2282               delete newTrack;
2283               return 0x0;
2284             }
2285           newTrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
2286           newTrack->SetTRDQuality(1);
2287         }
2288       else
2289         {
2290           return 0x0;
2291         }
2292     }
2293   else if (fTrackCutsType==2)
2294     {
2295       if(fEsdTrackCuts->AcceptTrack(esdtrack))
2296         {
2297           newTrack = new AliESDtrack(*esdtrack);
2298           newTrack->SetTRDQuality(0);
2299         }
2300       else
2301         return 0x0;
2302     }
2303   else
2304     {
2305       printf("Unknown track cuts type: %d\n",fTrackCutsType);
2306       return 0x0;
2307     }
2308
2309   return newTrack;
2310 }
2311
2312 //________________________________________________________________________
2313 void AliAnalysisTaskFullppJet::GetESDEMCalClusters()
2314 {
2315   //
2316   // Get emcal clusters - selected
2317   //
2318
2319   if(fSysTrkClsMth)
2320     {
2321       if (!TGeoGlobalMagField::Instance()->GetField()) fESD->InitMagneticField();
2322       fRecoUtil->FindMatches(fESD,0x0,fGeom);
2323       fRecoUtil->SetClusterMatchedToTrack(fESD);
2324       fRecoUtil->SetTracksMatchedToCluster(fESD);
2325     }
2326
2327   const Int_t nCaloClusters = fESD->GetNumberOfCaloClusters();
2328   for(Int_t i = 0 ; i < nCaloClusters; i++)
2329     {
2330       AliESDCaloCluster * cl = (AliESDCaloCluster *) fESD->GetCaloCluster(i);
2331       if(!IsGoodCluster(cl)) continue;
2332       AliESDCaloCluster *newCluster = new AliESDCaloCluster(*cl);
2333       if (newCluster->E() < fClsEtMin[fTriggerType] || newCluster->E() > fClsEtMax[fTriggerType])
2334         {delete newCluster; continue;}
2335
2336       // Absolute scale
2337       if(fPeriod.Contains("lhc11a",TString::kIgnoreCase))
2338         {
2339           Double_t newE = newCluster->E() * 1.02;
2340           newCluster->SetE(newE);
2341         }
2342
2343       // Trigger efficiency systematic uncertainty
2344       if(fSysJetTrigEff)
2345         {
2346           Double_t newE = newCluster->E() * (1+fVaryJetTrigEff);
2347           newCluster->SetE(newE);
2348           if(fTriggerType==0) fhClsE[fTriggerType]->Fill(newCluster->E());
2349           if(fTriggerType==1)
2350             {
2351               if(fPeriod.Contains("lhc12a15a",TString::kIgnoreCase)) fhClsE[0]->Fill(newCluster->E());
2352               if(newCluster->Chi2()>0.5) fhClsE[fTriggerType]->Fill(newCluster->E());
2353             }
2354         }
2355
2356       // Cluster energy scale systematic uncertainty
2357       if(fSysClusterEScale)
2358         {
2359           Double_t newE = newCluster->E() * (1+fVaryClusterEScale);
2360           newCluster->SetE(newE);
2361         }
2362
2363       // Cluster energy resolution systematic uncertainty
2364       if(fSysClusterERes)
2365         {
2366           Double_t oldE = newCluster->E();
2367           Double_t resolution = fClusterEResolution->Eval(oldE);
2368           Double_t smear = resolution * TMath::Sqrt((1+fVaryClusterERes)*(1+fVaryClusterERes)-1);
2369           Double_t newE = oldE + fRandomGen->Gaus(0, smear);
2370           newCluster->SetE(newE);
2371         }
2372
2373       // non-linearity systematic uncertainty
2374       if(fSysNonLinearity)
2375         {
2376           Double_t oldE = newCluster->E();
2377           Double_t newE = oldE/fNonLinear->Eval(oldE) * 1.012;
2378           newCluster->SetE(newE);
2379         }
2380
2381       // clusterizer systematic uncertainty
2382       if(fSysClusterizer)
2383         {
2384           fhSysClusterE[fTriggerType][0]->Fill(newCluster->E());
2385           fhSysNCellVsClsE[fTriggerType][0]->Fill(newCluster->E(),newCluster->GetNCells());
2386         }
2387
2388       // hadronic correction
2389       Double_t clsE = newCluster->E();
2390       Double_t subE = 0., eRes = 0., mcSubE = 0, MIPE = 0.;
2391       if(fElectronRejection || fHadronicCorrection)
2392         subE= SubtractClusterEnergy(newCluster,  eRes, MIPE, mcSubE);
2393
2394       if(!fStudySubEInHC && !fStudyMcOverSubE)  clsE -= subE;
2395       if(clsE<0) {delete newCluster; continue;}
2396       newCluster->SetE(clsE);
2397       newCluster->SetTOF(eRes);
2398       newCluster->SetDispersion(subE);
2399       newCluster->SetEmcCpvDistance(MIPE);
2400       newCluster->SetDistanceToBadChannel(mcSubE);
2401       fClusterArray->Add(newCluster);
2402
2403       // clusterizer systematic uncertainty
2404       if(fSysClusterizer)
2405         {
2406           fhSysClusterE[fTriggerType][1]->Fill(newCluster->E());
2407           fhSysNCellVsClsE[fTriggerType][1]->Fill(newCluster->E(),newCluster->GetNCells());
2408         }
2409     }
2410
2411   if(fVerbosity>5) printf("[i] # of EMCal clusters in event: %d\n", fClusterArray->GetEntries());
2412
2413 }
2414
2415
2416 //________________________________________________________________________
2417 Bool_t AliAnalysisTaskFullppJet::IsGoodCluster(AliESDCaloCluster *cluster)
2418 {
2419   //
2420   // Select good clusters
2421   //
2422
2423   if(!cluster) return kFALSE;
2424   if (!cluster->IsEMCAL()) return kFALSE;
2425   if(fRejectExoticCluster && fRecoUtil->IsExoticCluster(cluster, (AliVCaloCells*)fESD->GetEMCALCells())) return kFALSE;
2426   if(fRemoveBadChannel && fRecoUtil->ClusterContainsBadChannel(fGeom, cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
2427   return kTRUE;
2428 }
2429
2430
2431 //________________________________________________________________________
2432 Double_t AliAnalysisTaskFullppJet::SubtractClusterEnergy(AliESDCaloCluster *cluster, Double_t &eRes, Double_t &MIPE, Double_t &mcSubE)
2433 {
2434   //
2435   // Hadronic correction
2436   //
2437
2438   mcSubE = 0;
2439   eRes = 0;
2440   MIPE = 0;
2441
2442   eRes += TMath::Power(fClusterEResolution->Eval(cluster->E()),2);
2443   Double_t subE = 0., sumTrkPt = 0., sumTrkP = 0.;
2444   Int_t nTrack = 0;
2445   TArrayI *matched = cluster->GetTracksMatched();
2446   if(!matched) return 0;
2447   for(Int_t im=0; im<matched->GetSize(); im++)
2448     {
2449       Int_t trkIndex = matched->At(im);
2450       if(trkIndex<0 || trkIndex>=fESD->GetNumberOfTracks()) continue;
2451       Bool_t isSelected = kFALSE;
2452       Int_t index = -1;
2453       for(Int_t j=0; j<fTrackArray->GetEntriesFast(); j++)
2454         {
2455           AliESDtrack *tr = (AliESDtrack*)fTrackArray->At(j);
2456           if( trkIndex == tr->GetID() )
2457             {
2458               isSelected = kTRUE;
2459               index = j;
2460               break;
2461             }
2462         }
2463       if(!isSelected) continue;
2464       nTrack++;
2465       AliESDtrack *track = (AliESDtrack*)fTrackArray->At(index);
2466       Double_t trkP = track->P();
2467       sumTrkPt += track->Pt(); sumTrkP += track->P();
2468       if(fSysTrkPtRes) trkP = TMath::Sqrt(TMath::Power(track->Pt()+GetSmearedTrackPt(track),2)+TMath::Power(track->Pz(),2));
2469       if(IsElectron(track,cluster->E())) //electrons
2470         {
2471           if(fElectronRejection)
2472             {
2473               subE+= trkP;
2474               MIPE+= trkP;
2475               eRes += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
2476             }
2477         }
2478       else //hadrons
2479         {
2480           if(fHadronicCorrection)
2481             {
2482               MIPE += (trkP>0.27)?0.27:trkP; //MIP correction
2483               if(fFractionHC>2)
2484                 {
2485                   if(trkP>0.27)
2486                     {
2487                       subE+=0.27;
2488                     }
2489                   else
2490                     {
2491                       subE+=trkP;
2492                       eRes += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
2493                     }
2494                 }
2495               else
2496                 {
2497                   if(trkP>fHCLowerPtCutMIP) subE += 0.27;
2498                   else subE+=fFractionHC*trkP;
2499                   eRes += TMath::Power(fFractionHC*track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
2500                 }
2501             }
2502         }
2503     }
2504
2505   if(fSaveQAHistos) fhNMatchedTrack[fTriggerType]->Fill(nTrack);
2506   if(nTrack>0)
2507     {
2508       Double_t fraction[4] = {0.3,0.5,0.7,1.0};
2509       for(Int_t j=0; j<4; j++)
2510         {
2511           Double_t subETmp = sumTrkP*fraction[j];
2512           if(subETmp>cluster->E()) subETmp = cluster->E();
2513           if(fSaveQAHistos) fhSubEVsTrkPt[fTriggerType][j]->Fill(sumTrkPt,subETmp/sumTrkP);
2514         }
2515     }
2516
2517   eRes = TMath::Sqrt(eRes);
2518
2519   if(fIsMC && nTrack>0)
2520     {
2521       Double_t neutralE = 0;
2522       TArrayI* labels = cluster->GetLabelsArray();
2523       if(labels)
2524         {
2525           for(Int_t il=0; il<labels->GetSize(); il++)
2526             {
2527               Int_t ipart = labels->At(il);
2528               if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
2529                 {
2530                   AliVParticle* vParticle = fMC->GetTrack(ipart);
2531                   if(vParticle->Charge()==0)
2532                     {
2533                       neutralE += vParticle->E();
2534                     }
2535                 }
2536             }
2537         }
2538       mcSubE = cluster->E() - neutralE;
2539       if(mcSubE<0)
2540         mcSubE=0;
2541     }
2542
2543
2544   return subE;
2545 }
2546
2547
2548 //
2549 //________________________________________________________________________
2550 //
2551 Double_t AliAnalysisTaskFullppJet::GetExoticEnergyFraction(AliESDCaloCluster *cluster)
2552 {
2553   //
2554   // Exotic fraction: f_cross
2555   // Adpated from AliEMCalRecoUtils
2556
2557   if(!cluster) return -1;
2558   AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells();
2559   if(!cells)  return -1;
2560
2561   // Get highest energy tower
2562   Int_t iSupMod = -1, absID = -1, ieta = -1, iphi = -1,iTower = -1, iIphi = -1, iIeta = -1;
2563   Bool_t share = kFALSE;
2564   fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, iSupMod, ieta, iphi, share);
2565   fGeom->GetCellIndex(absID,iSupMod,iTower,iIphi,iIeta);
2566   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
2567
2568   Int_t absID1 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi+1, ieta);
2569   Int_t absID2 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi-1, ieta);
2570   Int_t absID3 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi, ieta+1);
2571   Int_t absID4 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi, ieta-1);
2572
2573   Float_t  ecell  = 0, ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
2574   Double_t tcell  = 0, tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
2575   Bool_t   accept = 0, accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0;
2576   const Int_t bc  = 0;
2577
2578   accept  = fRecoUtil->AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
2579
2580   if(!accept) return -1;
2581
2582   if(ecell < 0.5) return -1;
2583
2584   accept1 = fRecoUtil->AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
2585   accept2 = fRecoUtil->AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
2586   accept3 = fRecoUtil->AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
2587   accept4 = fRecoUtil->AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
2588
2589   const Double_t exoticCellDiffTime = 1e6;
2590   if(TMath::Abs(tcell-tcell1)*1.e9 > exoticCellDiffTime) ecell1 = 0 ;
2591   if(TMath::Abs(tcell-tcell2)*1.e9 > exoticCellDiffTime) ecell2 = 0 ;
2592   if(TMath::Abs(tcell-tcell3)*1.e9 > exoticCellDiffTime) ecell3 = 0 ;
2593   if(TMath::Abs(tcell-tcell4)*1.e9 > exoticCellDiffTime) ecell4 = 0 ;
2594
2595   Float_t eCross = ecell1+ecell2+ecell3+ecell4;
2596
2597   return 1-eCross/ecell;
2598 }
2599
2600
2601 //________________________________________________________________________
2602 void AliAnalysisTaskFullppJet::GetMCInfo()
2603 {
2604   //
2605   // Get # of trials per ESD event
2606   //
2607
2608   AliStack *stack = fMC->Stack();
2609   if(stack)
2610     {
2611       AliGenEventHeader *head = dynamic_cast<AliGenEventHeader*>(fMC->GenEventHeader());
2612       if (head == 0x0)
2613        {
2614          AliError("Could not get the event header");
2615          return;
2616        }
2617
2618      AliGenPythiaEventHeader *headPy = dynamic_cast<AliGenPythiaEventHeader*>(head);
2619      if (headPy != 0x0)
2620        {
2621          if (headPy->Trials() > 0)
2622            {
2623              fhNTrials[0]->Fill(0.5,headPy->Trials());
2624
2625            }
2626        }
2627     }
2628 }
2629
2630
2631
2632 //________________________________________________________________________
2633 void AliAnalysisTaskFullppJet::Terminate(Option_t *)
2634 {
2635   //
2636   // Called once at the end of the query
2637   //
2638   Info("Terminate","Terminate");
2639   AliAnalysisTaskSE::Terminate();
2640
2641 }
2642
2643 //
2644 //________________________________________________________________________
2645 //
2646 Int_t AliAnalysisTaskFullppJet::RunOfflineTrigger()
2647 {
2648   //
2649   // Run trigger offline
2650   //
2651
2652   fIsEventTriggerBit = 0;
2653   Int_t isTrigger = 0;
2654   Int_t ncl = fESD->GetNumberOfCaloClusters();
2655   for(Int_t icl=0; icl<ncl; icl++)
2656     {
2657       // Check every cluster
2658       AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2659       if(!IsGoodCluster(cluster)) continue;
2660       Double_t pro = GetOfflineTriggerProbability(cluster);
2661       Double_t rand = fRandomGen->Uniform();
2662       if(rand<pro)
2663         {
2664           isTrigger = 1;
2665           fIsEventTriggerBit = 1;
2666           cluster->SetChi2(1);
2667         }
2668       else
2669         cluster->SetChi2(0);
2670     }
2671   return isTrigger;
2672 }
2673
2674
2675 //
2676 //________________________________________________________________________
2677 //
2678 Double_t AliAnalysisTaskFullppJet::GetOfflineTriggerProbability(AliESDCaloCluster *cluster)
2679 {
2680   //
2681   // Get the probablity of the given cluster to trigger the event
2682   //
2683
2684   Double_t pro = 0;
2685   // Check the trigger mask
2686   AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells();
2687   Int_t iSupMod = -1, absID = -1, ieta = -1, iphi = -1,iTower = -1, iIphi = -1, iIeta = -1;
2688   Bool_t share = kFALSE;
2689   fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, iSupMod, ieta, iphi, share); // Get the position of the most energetic cell in the cluster
2690   fGeom->GetCellIndex(absID,iSupMod,iTower,iIphi,iIeta);
2691   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
2692   // convert co global phi eta
2693   Int_t gphi = iphi + 24*(iSupMod/2);
2694   Int_t geta = ieta + 48*(iSupMod%2);
2695   // get corresponding FALTRO
2696   Int_t fphi = gphi / 2;
2697   Int_t feta = geta / 2;
2698   if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5 && iSupMod>-1) // check the trigger mask
2699     {
2700       Double_t clsE = cluster->E();
2701       if(fSysClusterEScale) clsE = clsE * (1+fVaryClusterEScale); // Used for systematic uncertainty. Not needed for regular analysis
2702       if(clsE>10) pro = 1; // Probability is 1 at high E
2703       else
2704         {
2705           Int_t bin = fTriggerCurve[iSupMod]->FindFixBin(clsE);
2706           pro = fTriggerCurve[iSupMod]->GetBinContent(bin)/fTriggerEfficiency[iSupMod]->Eval(10); // Read the probability from trigger turn-on curves
2707         }
2708     }
2709   return pro;
2710 }
2711
2712
2713
2714 //
2715 //________________________________________________________________________
2716 //
2717 Int_t AliAnalysisTaskFullppJet::GetClusterSuperModule(AliESDCaloCluster *cluster)
2718 {
2719   //
2720   // Return the given cluster supermodule
2721   //
2722
2723   Float_t pos[3];
2724   cluster->GetPosition(pos);
2725   TVector3 clsVec(pos[0],pos[1],pos[2]);
2726
2727   Int_t sMod=-1;
2728   fGeom->SuperModuleNumberFromEtaPhi(clsVec.Eta(), clsVec.Phi(), sMod);
2729   return sMod;
2730 }
2731
2732
2733 //________________________________________________________________________
2734 Bool_t AliAnalysisTaskFullppJet::IsElectron(AliESDtrack *track, Double_t clsE) const
2735 {
2736   //
2737   // Check if the given track is an electron candidate based on de/dx
2738   //
2739
2740   if(track->GetTPCsignal()<=fdEdxMax && track->GetTPCsignal()>=fdEdxMin && (clsE/track->P())<=fEoverPMax && (clsE/track->P())>=fEoverPMin )
2741     return kTRUE;
2742   else
2743     return kFALSE;
2744 }
2745
2746
2747 //________________________________________________________________________
2748 Double_t AliAnalysisTaskFullppJet::GetTrkEff(Double_t inPt)
2749 {
2750   //
2751   // Get tracking efficiency estimated from simulation
2752   //
2753
2754   Double_t eff = 1;
2755   Double_t ptBound[4] = {0, 0.5, 3.8, 300};
2756
2757   for(Int_t i=0; i<3; i++)
2758     {
2759       if( inPt < ptBound[i+1] && inPt >= ptBound[i])
2760         {
2761           eff = fTrkEffFunc[i]->Eval(inPt);
2762           break;
2763         }
2764     }
2765   return eff;
2766 }
2767
2768 //________________________________________________________________________
2769 Double_t AliAnalysisTaskFullppJet::GetSmearedTrackPt(AliESDtrack *track)
2770 {
2771   //
2772   // Smear track momentum
2773   //
2774
2775   Double_t resolution = track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2());
2776   Double_t smear = resolution*TMath::Sqrt((1+fVaryTrkPtRes)*(1+fVaryTrkPtRes)-1);
2777   return fRandomGen->Gaus(0, smear);
2778
2779 }
2780
2781
2782 //________________________________________________________________________
2783 void AliAnalysisTaskFullppJet::CheckEventTriggerBit()
2784 {
2785   //
2786   // Check if the triggered events have correct trigger bit
2787   //
2788
2789  fIsEventTriggerBit = 0;
2790   // constants
2791   const Int_t nColsModule = 2;
2792   const Int_t nRowsModule = 5;
2793   const Int_t nRowsFeeModule = 24;
2794   const Int_t nColsFeeModule = 48;
2795   const Int_t nColsFaltroModule = 24;
2796   const Int_t nRowsFaltroModule = 12;
2797
2798   // part 1, trigger extraction -------------------------------------
2799   Int_t trigtimes[30], globCol, globRow, ntimes;
2800   Int_t trigger[nColsFaltroModule*nColsModule][nRowsFaltroModule*nRowsModule];
2801
2802   // erase trigger maps
2803   for( Int_t i = 0; i < nColsFaltroModule*nColsModule; i++ )
2804     {
2805       for( Int_t j = 0; j < nRowsFaltroModule*nRowsModule; j++ )
2806         {
2807           trigger[i][j] = 0;
2808         }
2809     }
2810
2811   Int_t fTrigCutLow = 7, fTrigCutHigh = 10, trigInCut = 0;
2812   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
2813   // go through triggers
2814   if( fCaloTrigger->GetEntries() > 0 )
2815     {
2816       // needs reset
2817       fCaloTrigger->Reset();
2818       while( fCaloTrigger->Next() )
2819         {
2820           fCaloTrigger->GetPosition( globCol, globRow );  // get position in global 2x2 tower coordinates
2821           fCaloTrigger->GetNL0Times( ntimes );   // get dimension of time arrays
2822           if( ntimes < 1 ) continue;    // no L0s in this channel. Presence of the channel in the iterator still does not guarantee that L0 was produced!!
2823           fCaloTrigger->GetL0Times( trigtimes );  // get timing array
2824           if(fCheckTriggerMask && fTriggerMask->GetBinContent(globCol+1,globRow+1)<0.5) continue;
2825           trigInCut = 0;
2826           for( Int_t i = 0; i < ntimes; i++ )
2827             {
2828               if( trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )  // check if in cut
2829                 {
2830                   trigInCut = 1;
2831                 }
2832             }
2833           if(trigInCut==1) trigger[globCol][globRow] = 1;
2834         } // calo trigger entries
2835     } // has calo trigger entries
2836
2837
2838
2839   // part 2 go through the clusters here -----------------------------------
2840   Int_t nCell, iCell;
2841   UShort_t *cellAddrs;
2842   Int_t absID = -1, nSupMod=-1, nModule=-1, nIphi=-1, nIeta=-1, iphi=-1, ieta=-1, gphi=-1, geta=-1, feta=-1, fphi=-1;
2843   Bool_t share = kFALSE;
2844   AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells();
2845   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
2846     {
2847       AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2848       if(!cluster || !cluster->IsEMCAL()) continue;
2849       cluster->SetChi2(0);
2850       if(!IsGoodCluster(cluster)) continue;
2851
2852       //Clusters with most energetic cell in dead region can't be triggered
2853       fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, nSupMod, ieta, iphi, share);
2854       fGeom->GetCellIndex(absID,nSupMod,nModule, nIphi, nIeta);
2855       fGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2856       gphi = iphi + nRowsFeeModule*(nSupMod/2);
2857       geta = ieta + nColsFeeModule*(nSupMod%2);
2858       fphi = gphi / 2;
2859       feta = geta / 2;
2860
2861       if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5)
2862         {
2863           nCell = cluster->GetNCells();  // get cluster cells
2864           cellAddrs = cluster->GetCellsAbsId();  // get the cell addresses
2865           for( iCell = 0; iCell < nCell; iCell++ )
2866             {
2867               // get cell position
2868               fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
2869               fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2870
2871               // convert co global phi eta
2872               gphi = iphi + nRowsFeeModule*(nSupMod/2);
2873               geta = ieta + nColsFeeModule*(nSupMod%2);
2874
2875               // get corresponding FALTRO
2876               fphi = gphi / 2;
2877               feta = geta / 2;
2878               // try to match with a triggered
2879               if( trigger[feta][fphi] == 1)
2880                 {
2881                   cluster->SetChi2(1);
2882                   fIsEventTriggerBit = 1;
2883                   //break;
2884                 }
2885             } // cells
2886         }
2887     } // clusters
2888 }
2889
2890
2891 //________________________________________________________________________
2892 void AliAnalysisTaskFullppJet::PrintConfig()
2893 {
2894   //
2895   // Print configuration
2896   //
2897
2898   const char *trackCutName[3] = {"Hybrid tracks","TPCOnly tracks","Golden tracks"};
2899   const char *triggerType[2]  = {"MB","EMC"};
2900   const char *decision[2] = {"no","yes"};
2901   const char *recombination[] = {"E_scheme","pt_scheme","pt2_scheme","Et_scheme","Et2_scheme","BIpt_scheme","BIpt2_scheme"};
2902   const char *type[2] = {"local","grid"};
2903   if(fStudySubEInHC || fStudyMcOverSubE)
2904     {
2905       printf("\n\n=================================\n");
2906       printf("======WARNING: HC is ingored!======\n");
2907       printf("======    NOT for PHYSICS!   ======\n\n");
2908     }
2909   printf("Run period: %s\n",fPeriod.Data());
2910   printf("Reject SPD pileup: %s\n",decision[fRejectPileup]);
2911   printf("Reject exotic triggered events: %s\n",decision[fRejectExoticTrigger]);
2912   printf("Is this MC data: %s\n",decision[fIsMC]);
2913   printf("Only find charged jets in MC: %s\n",decision[fChargedMC]);
2914   printf("Analyze on local or grid? %s\n",type[fAnaType]);
2915   printf("Run offline trigger on MC: %s\n",decision[fOfflineTrigger]);
2916   if(fIsMC)
2917     printf("Is K0 and n included: %s\n",decision[1-fRejectNK]);
2918   printf("Constrain tracks in EMCal acceptance: %s\n",decision[fConstrainChInEMCal]);
2919   printf("Track selection:    %s, |eta| < %2.1f\n",trackCutName[fTrackCutsType], fTrkEtaMax);
2920   for(Int_t i=0; i<2; i++)
2921     {
2922       printf("Track pt cut:       %s -> %2.2f < pT < %2.1f\n",triggerType[i], fTrkPtMin[i], fTrkPtMax[i]);
2923     }
2924   for(Int_t i=0; i<2; i++)
2925     {
2926       printf("Cluster selection:  %s -> %2.2f < Et < %2.1f\n",triggerType[i],fClsEtMin[i], fClsEtMax[i]);
2927     }
2928   printf("Electron selectoin: %2.0f < dE/dx < %2.0f, %1.1f < E/P < %1.1f\n",fdEdxMin, fdEdxMax, fEoverPMin, fEoverPMax);
2929   printf("Reject exotic cluster: %s\n",decision[fRejectExoticCluster]);
2930   printf("Remove problematic region in SM4: %s\n",decision[fRemoveBadChannel]);
2931   printf("Use only good SM (1,2,6,7,8,9) for trigger: %s\n",decision[fUseGoodSM]);
2932   printf("Reject electron: %s\n", decision[fElectronRejection]);
2933   printf("Correct hadron: %s\n",decision[fHadronicCorrection]);
2934   printf("HC fraction: %2.1f up to %2.0f GeV/c\n",fFractionHC,fHCLowerPtCutMIP);
2935   printf("Find charged jets: %s\n", decision[fFindChargedOnlyJet]);
2936   printf("Find netural jets: %s\n", decision[fFindNeutralOnlyJet]);
2937   printf("Find good jets: %s\n",decision[fSpotGoodJet]);
2938   printf("Jet radius: %s\n",fRadius.Data());
2939   printf("Jet recombination scheme: %s\n",recombination[fRecombinationScheme]);
2940   printf("Correct tracking efficiency: %s\n",decision[fCheckTrkEffCorr]);
2941   printf("Save jet QA histos: %s\n",decision[fSaveQAHistos]);
2942   printf("Systematics: jet efficiency: %s with variation %1.0f%%\n",decision[fSysJetTrigEff],fVaryJetTrigEff*100);
2943   printf("Systematics: tracking efficiency: %s with variation %1.0f%%\n",decision[fSysTrkEff],fVaryTrkEff*100);
2944   printf("Systematics: track pt resolution: %s with variation %1.0f%%\n",decision[fSysTrkPtRes],fVaryTrkPtRes*100);
2945   printf("Systematics: track-cluster matching: %s with |dEta|<%2.3f, |dPhi|<%2.3f\n",decision[fSysTrkClsMth],fCutdEta,fCutdPhi);
2946   printf("Systematics: EMCal non-linearity: %s\n",decision[fSysNonLinearity]);
2947   printf("Systematics: EMCal energy scale: %s with uncertainty of %1.0f%%\n",decision[fSysClusterEScale],fVaryClusterEScale*100);
2948   printf("Systematics: EMCal energy resolution: %s with uncertainty of %1.0f%%\n",decision[fSysClusterERes],fVaryClusterERes*100);
2949   printf("Smear lhc12a15a: %s\n",decision[fSmearMC]);
2950   printf("Run UE analysis: %s\n",decision[fRunUE]);
2951   printf("Run secondaries: %s\n",decision[fRunSecondaries]);
2952   printf("=======================================\n\n");
2953 }
2954
2955 //________________________________________________________________________
2956 void AliAnalysisTaskFullppJet::CheckExoticEvent()
2957 {
2958   //
2959   // Check if the event containts exotic clusters
2960   //
2961
2962   Double_t leadingE = 0;
2963   Int_t nCls = fESD->GetNumberOfCaloClusters();
2964   for(Int_t icl=0; icl<nCls; icl++)
2965     {
2966       AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2967       if(!cluster) continue;
2968       if (!cluster->IsEMCAL()) continue;
2969       if(fRecoUtil->IsExoticCluster(cluster, (AliVCaloCells*)fESD->GetEMCALCells()) && cluster->E()>leadingE)
2970         {
2971           leadingE = cluster->E();
2972         }
2973     }
2974   if(leadingE>3) fIsExoticEvent3GeV = kTRUE;
2975   if(leadingE>5) fIsExoticEvent5GeV = kTRUE;
2976 }
2977
2978
2979 //________________________________________________________________________
2980 Bool_t AliAnalysisTaskFullppJet::HasPrimaryVertex() const
2981 {
2982   //
2983   // Check if the primary vertex exists
2984   //
2985   const AliESDVertex* vtx = fESD->GetPrimaryVertex();
2986   if (!vtx || vtx->GetNContributors()<1) return kFALSE;
2987   else return kTRUE;
2988 }
2989
2990
2991 //________________________________________________________________________
2992 Bool_t AliAnalysisTaskFullppJet::IsPrimaryVertexOk() const
2993 {
2994   //
2995   // Check if the event vertex is good
2996   //
2997   return kTRUE;
2998 }
2999
3000 //________________________________________________________________________
3001 Bool_t AliAnalysisTaskFullppJet::IsTPCOnlyVtx() const
3002 {
3003   //
3004   // Check if the event only has valid TPC vertex
3005   //
3006
3007   const Bool_t isPrimaryVtx = HasPrimaryVertex();
3008
3009   Bool_t isGoodVtx = kTRUE;
3010   AliESDVertex* goodvtx = const_cast<AliESDVertex*>(fESD->GetPrimaryVertexTracks());
3011   if(!goodvtx || goodvtx->GetNContributors()<1)
3012     goodvtx = const_cast<AliESDVertex*>(fESD->GetPrimaryVertexSPD()); // SPD vertex
3013   if(!goodvtx || !goodvtx->GetStatus()) isGoodVtx = kFALSE;  // Event rejected
3014
3015   if( isPrimaryVtx && !isGoodVtx )
3016     return kTRUE;
3017   else
3018     return kFALSE;
3019 }
3020
3021 //________________________________________________________________________
3022 void AliAnalysisTaskFullppJet::CheckTPCOnlyVtx(const UInt_t trigger)
3023 {
3024   //
3025   // Check the fraction of accepted events that have only TPC vertex
3026   //
3027   Int_t lTriggerType = -1;
3028   if (trigger & AliVEvent::kMB)       lTriggerType = 1;
3029   if (trigger & AliVEvent::kFastOnly) lTriggerType = 0;
3030   if (trigger & AliVEvent::kEMC1)     lTriggerType = 2;
3031   if(lTriggerType==-1) return;
3032   fhEventStatTPCVtx->Fill(0.5+lTriggerType*3);
3033   if(HasPrimaryVertex()) fhEventStatTPCVtx->Fill(1.5+lTriggerType*3);
3034   if(IsTPCOnlyVtx()) fhEventStatTPCVtx->Fill(2.5+lTriggerType*3);
3035 }
3036
3037 //________________________________________________________________________
3038 Bool_t AliAnalysisTaskFullppJet::IsLEDEvent() const
3039 {
3040   //
3041   // Check if the event is contaminated by LED signal
3042   //
3043
3044   AliESDCaloCells *cells = fESD->GetEMCALCells();
3045   Short_t nCells = cells->GetNumberOfCells();
3046   Int_t nCellCount[2] = {0,0};
3047   for(Int_t iCell=0; iCell<nCells; iCell++)
3048     {
3049       Int_t cellId = cells->GetCellNumber(iCell);
3050       Double_t cellE = cells->GetCellAmplitude(cellId);
3051       Int_t sMod = fGeom->GetSuperModuleNumber(cellId);
3052
3053       if(sMod==3 || sMod==4)
3054         {
3055           if(cellE>0.1)
3056             nCellCount[sMod-3]++;
3057         }
3058     }
3059   Bool_t isLED=kFALSE;
3060
3061   if(fPeriod.CompareTo("lhc11a")==0)
3062     {
3063       if (nCellCount[1] > 100)
3064         isLED = kTRUE;
3065       Int_t runN = fESD->GetRunNumber();
3066       if (runN>=146858 && runN<=146860)
3067         {
3068           if(fTriggerType==0 && nCellCount[0]>=21) isLED=kTRUE;
3069           if(fTriggerType==1 && nCellCount[0]>=35) isLED=kTRUE;
3070         }
3071     }
3072
3073   return isLED;
3074 }