PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.cxx

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