PWGJE/UserTasks/AliAnalysisTaskPPJetSpectra.cxx

   1 #include "AliAnalysisTaskSE.h"
   2 #include "AliAnalysisHelperJetTasks.h"
   3 #include "THnSparse.h"
   4 #include "AliAnalysisManager.h"
   5 #include "AliInputEventHandler.h"
   6 #include "AliAODHandler.h"
   7 #include "AliAODInputHandler.h"
   8 #include "AliAODEvent.h"
   9 #include "AliESDEvent.h"
  10 #include "TAxis.h"
  11 #include "AliAODMCParticle.h"
  12
  13
  14 #include "AliAnalysisTaskPPJetSpectra.h"
  15 #include "AliAnalysisHelperJetTasks.h"
  16
  17 ClassImp(AliAnalysisTaskPPJetSpectra)
  18
  19 //---------------------------------------------------------------------------------------------------
  20 AliAnalysisTaskPPJetSpectra::AliAnalysisTaskPPJetSpectra()
  21  :AliAnalysisTaskSE()
  22   ,fOutputList(0)
  23   ,fESD(0)
  24   ,fAOD(0)
  25   ,fMC(0)
  26   ,fAODJets(0)
  27   ,fAODExtension(0)
  28   ,fNonStdFile("")
  29   ,fDebug(0)
  30   ,fUseMC(kFALSE)
  31   ,fEvtSelectionMask(0)
  32   ,fEventClass(-1)
  33   ,nVtxContCut(2)
  34   ,fVtxZcut(10)
  35   ,fVtxRcut(1)
  36   ,nTrackFilter(272)
  37   ,trackPtMin(0.15)
  38   ,trackPtMax(100.)
  39   ,trackEtaAbsMax(0.9)
  40   ,jetPtMin(0.)
  41   ,jetEtaCut(0.)
  42   ,jetZmax(0.)
  43   ,fhnEvent(0)
  44   ,fhnTracks(0)
  45   ,fhnMC(0)
  46   ,fhnMC2(0)
  47   ,fhnRecJetsNoCut(0)
  48   ,fhnGenJetsNoCut(0)
  49   ,fhnRecJetsCut(0)
  50   ,fhnGenJetsCut(0)
  51   ,fhnRecBckg(0)
  52   ,fhnGenBckg(0)
  53   ,fhnRecJetsTrackUEcor(0)
  54   ,fhnGenJetsTrackUEcor(0)
  55   ,fhnRecJetsBckgUEcor(0)
  56   ,fhnGenJetsBckgUEcor(0)
  57   ,fhnTrackUE(0)
  58   ,fhnParticleUE(0)
  59   ,fhnBckgRecUE(0)
  60   ,fhnBckgGenUE(0)
  61   ,fRecJetBranch("")
  62   ,fGenJetBranch("")
  63   ,fRecBckgBranch("")
  64   ,fGenBckgBranch("")
  65   ,fRecJetR(0.)
  66   ,fGenJetR(0.)
  67   ,fRecBckgR(0.)
  68   ,fGenBckgR(0.)
  69   ,fTrackType(0)
  70   ,fParticleType(0)
  71   ,fhnMatching(0)
  72   ,fhnTrackCorrMatching(0)
  73   ,fhnBckgCorrMatching(0)
  74   ,fhnTrackUEanal(0)
  75   ,kDoUEanalysis(kFALSE)
  76   ,fRejectPileUp(1)
  77 {
  78   if(fDebug) printf("%s: %d  Constructor\n",(char*)__FILE__, __LINE__);
  79   //DefineOutput(1, TList::Class());
  80 }
  81
  82 //---------------------------------------------------------------------------------------------------
  83 AliAnalysisTaskPPJetSpectra::AliAnalysisTaskPPJetSpectra(const char* name):AliAnalysisTaskSE(name)
  84   ,fOutputList(0)
  85   ,fESD(0)
  86   ,fAOD(0)
  87   ,fMC(0)
  88   ,fAODJets(0)
  89   ,fAODExtension(0)
  90   ,fNonStdFile("")
  91   ,fDebug(0)
  92   ,fUseMC(kFALSE)
  93   ,fEvtSelectionMask(0)
  94   ,fEventClass(-1)
  95   ,nVtxContCut(2)
  96   ,fVtxZcut(10)
  97   ,fVtxRcut(1)
  98   ,nTrackFilter(272)
  99   ,trackPtMin(0.15)
 100   ,trackPtMax(100.)
 101   ,trackEtaAbsMax(0.9)
 102   ,jetPtMin(0.)
 103   ,jetEtaCut(0.)
 104   ,jetZmax(0.)
 105   ,fhnEvent(0)
 106   ,fhnTracks(0)
 107   ,fhnMC(0)
 108   ,fhnMC2(0)
 109   ,fhnRecJetsNoCut(0)
 110   ,fhnGenJetsNoCut(0)
 111   ,fhnRecJetsCut(0)
 112   ,fhnGenJetsCut(0)
 113   ,fhnRecBckg(0)
 114   ,fhnGenBckg(0)
 115   ,fhnRecJetsTrackUEcor(0)
 116   ,fhnGenJetsTrackUEcor(0)
 117   ,fhnRecJetsBckgUEcor(0)
 118   ,fhnGenJetsBckgUEcor(0)
 119   ,fhnTrackUE(0)
 120   ,fhnParticleUE(0)
 121   ,fhnBckgRecUE(0)
 122   ,fhnBckgGenUE(0)
 123   ,fRecJetBranch("")
 124   ,fGenJetBranch("")
 125   ,fRecBckgBranch("")
 126   ,fGenBckgBranch("")
 127   ,fRecJetR(0.)
 128   ,fGenJetR(0.)
 129   ,fRecBckgR(0.)
 130   ,fGenBckgR(0.)
 131   ,fTrackType(0)
 132   ,fParticleType(0)
 133   ,fhnMatching(0)
 134   ,fhnTrackCorrMatching(0)
 135   ,fhnBckgCorrMatching(0)
 136   ,fhnTrackUEanal(0)
 137   ,kDoUEanalysis(kFALSE)
 138   ,fRejectPileUp(1)
 139 {
 140   if(fDebug) printf("%s: %d  Constructor\n",(char*)__FILE__, __LINE__);
 141   DefineOutput(1, TList::Class());
 142 }
 143
 144 //---------------------------------------------------------------------------------------------------
 145 void AliAnalysisTaskPPJetSpectra::Init()
 146 {
 147   if(fDebug) printf("%s: %d  Initialize\n", (char*)__FILE__, __LINE__);
 148   return;
 149 }
 150
 151 //---------------------------------------------------------------------------------------------------
 152 Bool_t AliAnalysisTaskPPJetSpectra::Notify()
 153 {
 154
 155   if(fDebug) printf("%s: %d  Notify\n", (char*)__FILE__, __LINE__);
 156   return kTRUE;
 157 }
 158
 159 //---------------------------------------------------------------------------------------------------
 160 void AliAnalysisTaskPPJetSpectra::UserCreateOutputObjects()
 161 {
 162   fOutputList = new TList;
 163   PostData(1,fOutputList);
 164   fOutputList->SetOwner(1);
 165
 166   const Double_t binsPt[] = {0.,2.,4.,6.,8.,10.,12., 14.,16., 18.,20.,24.,28.,32.,38.,44.,50.,58.,66.,76.,86.,100.,114.,128.,150.,250.};
 167   const Int_t nBinsPt = 25;
 168
 169   const Double_t trackPtBins[42] = {0,0.2,0.4,0.6,0.8,1.,1.2,1.4,1.6,1.8,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,8.,10.,12.,14.,16.,18.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90.,95.,100};
 170   const Int_t nTrackPtBins = 41;
 171
 172   if(fDebug) printf("%s: %d  Creating output objects\n",(char*)__FILE__,__LINE__);
 173   // triggered - centrality - z-vtxs - r-vtx - ncontrib - ntracks - IsPileUp?- njetsrec - njetsgen - accepted
 174   const Int_t evtsDim = 7;
 175   Int_t     evtsBins[evtsDim] = {   2,   22,   40,  20,  20,    2,  2};
 176   Double_t  evtsMins[evtsDim] = {-0.5,  -5., -20.,  0.,  0., -0.5, -0.5};
 177   Double_t  evtsMaxs[evtsDim] = { 1.5, 105.,  20.,  2., 40.,  1.5, 1.5};
 178   // pt - eta - phi - is272 - charge
 179   const Int_t trksDim = 7;
 180   Int_t     trksBins[trksDim] = {  nTrackPtBins,  10,             10,    2,   2,   2,    3};
 181   Double_t  trksMins[trksDim] = {            0.,  -1.,            0., -0.5,-0.5,-0.5, -1.5};
 182   Double_t  trksMaxs[trksDim] = {            0.,   1., 2*TMath::Pi(),  1.5, 1.5, 1.5,  1.5};
 183   // pt - eta - phi - AreaCh - AreaN - ptLeadingTrack/ptJet - nTracks
 184   const Int_t jetsDim = 7;
 185   Int_t     jetsBins[jetsDim] = { nBinsPt,  20,            10,3, 5, 20, 30};
 186   Double_t  jetsMins[jetsDim] = {      0., -1.,            0., 0., 0., 0., 0.};
 187   Double_t  jetsMaxs[jetsDim] = {      0.,  1., 2*TMath::Pi(), 3., 5., 1.,30.};
 188   // pt - eta - phi - AreaCH - AreaN - UEdensity - ptLeadingTrack/ptJet - ptOrig
 189   const Int_t jetsUEcorDim = 8;
 190   Int_t     jetsUEcorBins[jetsUEcorDim] = { nBinsPt,  10,            10, 3,  5, 100,  20, nBinsPt};
 191   Double_t  jetsUEcorMins[jetsUEcorDim] = {      0., -1.,            0., 0., 0.,  0.,  0.,  0.};
 192   Double_t  jetsUEcorMaxs[jetsUEcorDim] = {      0.,  1., 2*TMath::Pi(), 3., 5., 10.,  1.,  0.};
 193   // leading jet pt - eta - sumAllPt1 - sumAllPt2 - Rpar
 194   const Int_t ueDim = 5;
 195   Int_t       ueBins[ueDim] = {  nBinsPt, 20, 100, 100, 100};
 196   Double_t    ueMins[ueDim] = {       0., -1,  0.,  0., 0.};
 197   Double_t    ueMaxs[ueDim] = {       0.,  1, 10., 10., 10.};
 198   // ptGen - ptRec - ptFraction - maxDist - bothWays - (1 - rec/gen) - dR - dEta - dPhi - isSamePtBin
 199   const Int_t matchingDim = 10;
 200   Int_t       matchingBins[matchingDim] = {nBinsPt, nBinsPt,  20,    5,    2, 100, 10,  10, 10, 2  };
 201   Double_t    matchingMins[matchingDim] = {     0.,      0.,  0., 0.05, -0.5,   0,  0,   0,  0,-0.5};
 202   Double_t    matchingMaxs[matchingDim] = {     0.,      0.,  1., 0.55,  1.5,   1,0.5, 0.5,0.5, 1.5};
 203
 204   fhnEvent    = new THnSparseF("fhnEvent", "fhnEvent;triggered;centrality;zVtx;rVtx;Ncontrib;IsPileUp;isAccepted", evtsDim, evtsBins, evtsMins, evtsMaxs);
 205   fhnTracks   = new THnSparseF("fhnTracks", "fhnTraks;pT;eta;phi;isFilterMask272;isMask16;isMask256;charge;", trksDim, trksBins, trksMins, trksMaxs);
 206   fhnMC       = new THnSparseF("fhnMC",     "fhnTraks;pT;eta;phi;isFilterMask272;isMask16;isMask256;charge;", trksDim, trksBins, trksMins, trksMaxs);
 207   fhnMC2      = new THnSparseF("fhnMC2",    "fhnTraks;pT;eta;phi;isFilterMask272;isMask16;isMask256;charge;", trksDim, trksBins, trksMins, trksMaxs);
 208
 209   fhnRecJetsNoCut  = new THnSparseF("fhnRecJetsNoCut", "fhnRecJetsNoCut;pT;eta;phi;AreaCH;AreaN;z;nTracks;", jetsDim, jetsBins, jetsMins, jetsMaxs);
 210   fhnGenJetsNoCut  = new THnSparseF("fhnGenJetsNoCut", "fhnGenJetsNoCut;pT;eta;phi;AreaCH;AreaN;z;nTracks;", jetsDim, jetsBins, jetsMins, jetsMaxs);
 211   fhnRecJetsCut    = new THnSparseF("fhnRecJetsCut",     "fhnRecJetsCut;pT;eta;phi;AreaCH;AreaN;z;nTracks;", jetsDim, jetsBins, jetsMins, jetsMaxs);
 212   fhnGenJetsCut    = new THnSparseF("fhnGenJetsCut",     "fhnGenJetsCut;pT;eta;phi;AreaCH;AreaN;z;nTracks;", jetsDim, jetsBins, jetsMins, jetsMaxs);
 213   //UE sparses
 214   fhnTrackUE     = new THnSparseF("fhnTrackUE",      "fhnTrackUE;pT-leading jet;eta;sum1;sum2;avgSum;", ueDim, ueBins, ueMins, ueMaxs);
 215   fhnParticleUE  = new THnSparseF("fhnParticleUE","fhnParticleUE;pT-leading jet;eta;sum1;sum2;avgSum;", ueDim, ueBins, ueMins, ueMaxs);
 216   // pt - eta - phi - AreaCH - AreaN - UEdensity - ptLeadingTrack/ptJet - ptOrig
 217   fhnRecJetsTrackUEcor  = new THnSparseF("fhnRecJetsTrackUEcor", "fhnRecJetsTrackUEcor;pT;eta;phi;AreaCH;AreaN;UE;z;pTorig;", jetsUEcorDim, jetsUEcorBins, jetsUEcorMins, jetsUEcorMaxs);
 218   fhnGenJetsTrackUEcor  = new THnSparseF("fhnGenJetsTrackUEcor", "fhnGenJetsTrackUEcor;pT;eta;phi;AreaCH;AreaN;UE;z;pTorig;", jetsUEcorDim, jetsUEcorBins, jetsUEcorMins, jetsUEcorMaxs);
 219   // ptGen - ptRec - ptFraction - maxDist - bothWays
 220   fhnMatching = new THnSparseF("fhnMatching","fhnMatching;pTgen;pTrec;ptFraction;maxDist;bothWays;1-p_{T,rec}/p_{T,gen}; dR; dEta; dPhi;;isSamePtBin;", matchingDim, matchingBins, matchingMins, matchingMaxs);
 221   fhnTrackCorrMatching = new THnSparseF("fhnTrackCorrMatching","fhnTrackCorrMatching;pTgen;pTrec;ptFraction;maxDist;bothWays;1-p_{T,rec}/p_{T,gen}; dR; dEta; dPhi;isSamePtBin;", matchingDim, matchingBins, matchingMins, matchingMaxs);
 222
 223   fhnTracks->SetBinEdges(0, trackPtBins);
 224   fhnTrackUE->SetBinEdges(0, trackPtBins);
 225   fhnParticleUE->SetBinEdges(0, trackPtBins);
 226
 227   fhnRecJetsNoCut->SetBinEdges(0, binsPt);
 228   fhnGenJetsNoCut->SetBinEdges(0, binsPt);
 229   fhnRecJetsCut->SetBinEdges(0, binsPt);
 230   fhnGenJetsCut->SetBinEdges(0, binsPt);
 231
 232   fhnRecJetsTrackUEcor->SetBinEdges(0,binsPt);
 233   fhnRecJetsTrackUEcor->SetBinEdges(7,binsPt);
 234   fhnGenJetsTrackUEcor->SetBinEdges(0,binsPt);
 235   fhnGenJetsTrackUEcor->SetBinEdges(7,binsPt);
 236
 237   fhnMatching->SetBinEdges(0,binsPt);
 238   fhnMatching->SetBinEdges(1,binsPt);
 239   fhnTrackCorrMatching->SetBinEdges(0,binsPt);
 240   fhnTrackCorrMatching->SetBinEdges(1,binsPt);
 241
 242   fhnEvent->Sumw2();
 243   fhnTracks->Sumw2();
 244   fhnMC->Sumw2();
 245   fhnMC2->Sumw2();
 246   fhnRecJetsNoCut->Sumw2();
 247   fhnGenJetsNoCut->Sumw2();
 248   fhnRecJetsCut->Sumw2();
 249   fhnGenJetsCut->Sumw2();
 250   fhnRecJetsTrackUEcor->Sumw2();
 251   fhnGenJetsTrackUEcor->Sumw2();
 252   fhnTrackUE->Sumw2();
 253   fhnParticleUE->Sumw2();
 254   fhnMatching->Sumw2();
 255   fhnTrackCorrMatching->Sumw2();
 256
 257   fOutputList->Add(fhnEvent);
 258   fOutputList->Add(fhnTracks);
 259   fOutputList->Add(fhnMC);
 260   fOutputList->Add(fhnMC2);
 261   fOutputList->Add(fhnRecJetsNoCut);
 262   fOutputList->Add(fhnGenJetsNoCut);
 263   fOutputList->Add(fhnRecJetsCut);
 264   fOutputList->Add(fhnGenJetsCut);
 265   fOutputList->Add(fhnRecJetsTrackUEcor);
 266   fOutputList->Add(fhnGenJetsTrackUEcor);
 267   fOutputList->Add(fhnTrackUE);
 268   fOutputList->Add(fhnParticleUE);
 269   fOutputList->Add(fhnMatching);
 270   fOutputList->Add(fhnTrackCorrMatching);
 271
 272   const Int_t trackUEdim = 11;
 273   //max track pt - ch_part_dens - ch_pt_dens - avg_pt_track - zero particles -tMAX_part_dens-tMIN_part_dens-tDIF_part_dens-tMAX_pt_dens-tMIN_pt_dens-tDIF_pt_dens
 274   Int_t trackUEbins[trackUEdim] = { 100, 50, 50, 50,   2, 50, 50, 50, 50, 50, 50};
 275   Double_t trackUEmins[trackUEdim] = {  0., 0., 0.,  0.,-0.5, 0., 0., 0., 0., 0., 0.};
 276   Double_t trackUEmaxs[trackUEdim] = {100., 5., 5., 5., 1.5, 5., 5., 5., 5., 5., 5.};
 277
 278   fhnTrackUEanal = new THnSparseF("fhnTrackUEanal","fhnTrackUEanal;PTmax;#rho^{ch.part}_{T};#rho^{pT}_{T};avg-pt/track;zero_particles;#rho^{ch.part}_{transMAX};#rho^{ch.part}_{transMIN};#rho^{ch.part}_{transDIF};#rho^{pT}_{transMAX};#rho^{pT}_{transMIN};#rho^{pT}_{transDIF}", trackUEdim, trackUEbins, trackUEmins, trackUEmaxs);
 279   fhnTrackUEanal->Sumw2();
 280   fOutputList->Add(fhnTrackUEanal);
 281
 282 }
 283
 284 //---------------------------------------------------------------------------------------------------
 285 void AliAnalysisTaskPPJetSpectra::UserExec(Option_t *option)
 286 {
 287   Double_t evtContainer[7];
 288   Bool_t isEventSelected = EventSelection(evtContainer);
 289
 290   TList trackList;
 291   TList particleList;
 292   Int_t nTracks         = GetListOfTracks(fTrackType, &trackList);
 293   Int_t nParticles      = GetListOfTracks(fParticleType, &particleList);
 294
 295  if(fDebug > 10 )
 296     std::cout<< nTracks<< " " << nParticles<<std::endl;
 297
 298   TClonesArray *tcaRecJets = 0;
 299   TClonesArray *tcaGenJets = 0;
 300   TClonesArray *tcaRecBckg = 0;
 301   TClonesArray *tcaGenBckg = 0;
 302
 303   Int_t rJ, gJ, rB, gB;
 304
 305   if(fRecJetBranch.Length() > 0) {
 306     tcaRecJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fRecJetBranch.Data()));
 307     if(!tcaRecJets)                tcaRecJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fRecJetBranch.Data()));
 308     if(!tcaRecJets)                tcaRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fRecJetBranch.Data()));
 309     if(!tcaRecJets)                tcaRecJets = dynamic_cast<TClonesArray*>(fAOD->GetList()->FindObject(fRecJetBranch.Data()));
 310     if(fAODExtension&&!tcaRecJets) tcaRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fRecJetBranch.Data()));
 311     if(tcaRecJets) rJ =  tcaRecJets->GetEntries();
 312     else rJ = -1;
 313   } else rJ = -1;
 314
 315   if(fGenJetBranch.Length() > 0) {
 316     tcaGenJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fGenJetBranch.Data()));
 317     if(!tcaGenJets)                tcaGenJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fGenJetBranch.Data()));
 318     if(!tcaGenJets)                tcaGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fGenJetBranch.Data()));
 319     if(!tcaGenJets)                tcaGenJets = dynamic_cast<TClonesArray*>(fAOD->GetList()->FindObject(fGenJetBranch.Data()));
 320     if(fAODExtension&&!tcaGenJets) tcaGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fGenJetBranch.Data()));
 321     if(tcaGenJets) gJ =  tcaGenJets->GetEntries();
 322     else gJ = -1;
 323   } else gJ = -1;
 324
 325   if(fRecBckgBranch.Length() > 0) {
 326     tcaRecBckg = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fRecBckgBranch.Data()));
 327     if(!tcaRecBckg)                tcaRecBckg = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fRecBckgBranch.Data()));
 328     if(!tcaRecBckg)                tcaRecBckg = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fRecBckgBranch.Data()));
 329     if(!tcaRecBckg)                tcaRecBckg = dynamic_cast<TClonesArray*>(fAOD->GetList()->FindObject(fRecBckgBranch.Data()));
 330     if(fAODExtension&&!tcaRecBckg) tcaRecBckg = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fRecBckgBranch.Data()));
 331     rB =  tcaRecBckg->GetEntries();
 332   } else rB = -1;
 333
 334   if(fGenBckgBranch.Length() > 0) {
 335     tcaGenBckg = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fGenBckgBranch.Data()));
 336     if(!tcaGenBckg)                tcaGenBckg = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fGenBckgBranch.Data()));
 337     if(!tcaGenBckg)                tcaGenBckg = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fGenBckgBranch.Data()));
 338     if(!tcaGenBckg)                tcaGenBckg = dynamic_cast<TClonesArray*>(fAOD->GetList()->FindObject(fGenBckgBranch.Data()));
 339     if(fAODExtension&&!tcaGenBckg) tcaGenBckg = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fGenBckgBranch.Data()));
 340     gB =  tcaGenBckg->GetEntries();
 341   } else gB = -1;
 342
 343   TList jetRecListNoCut;
 344   TList jetGenListNoCut;
 345   Int_t nRecJetsNoCut   = (rJ>0?GetListOfJets(tcaRecJets, &jetRecListNoCut, kFALSE):0);
 346   Int_t nGenJetsNoCut   = (gJ>0?GetListOfJets(tcaGenJets, &jetGenListNoCut, kFALSE):0);
 347
 348   TList jetRecListCut;
 349   TList jetGenListCut;
 350   Int_t nRecJetsCut     = (rJ>0?GetListOfJets(tcaRecJets, &jetRecListCut, kTRUE):0);
 351   Int_t nGenJetsCut     = (gJ>0?GetListOfJets(tcaGenJets, &jetGenListCut, kTRUE):0);
 352
 353   TList bckgRecListCut;
 354   TList bckgGenListCut;
 355   Int_t nRecBckgCut     = (rB>0?GetListOfJets(tcaRecBckg, &bckgRecListCut, kTRUE):0);
 356   Int_t nGenBckgCut     = (gB>0?GetListOfJets(tcaGenBckg, &bckgGenListCut, kTRUE):0);
 357
 358   if(fDebug > 10) std::cout<<nRecJetsNoCut<<" "<<nGenJetsNoCut<<" "<<nRecJetsCut<< " "<< nGenJetsCut<<" "<<nRecBckgCut<<" "<<nGenBckgCut<<std::endl;
 359
 360   evtContainer[6] = (isEventSelected?1:0);
 361   fhnEvent->Fill(evtContainer);
 362
 363   if(!isEventSelected) {
 364     PostData(1,fOutputList);
 365     return;
 366   }
 367
 368   if(kDoUEanalysis) DoUEAnalysis(&trackList, (Double_t)0.15, (Double_t)0.9);
 369
 370   Double_t  recUEdensity = GetUE(&jetRecListCut, &trackList, fRecJetR, fhnTrackUE);
 371   Double_t  genUEdensity = GetUE(&jetGenListCut, &particleList, fGenJetR, fhnParticleUE);
 372
 373   TList  trackUEcorrJetRecListCut;
 374   TList  trackUEcorrJetGenListCut;
 375
 376   CorrectForUE(&jetRecListCut, recUEdensity, &trackUEcorrJetRecListCut, fhnRecJetsTrackUEcor);
 377   CorrectForUE(&jetGenListCut, genUEdensity, &trackUEcorrJetGenListCut, fhnGenJetsTrackUEcor);
 378
 379   FillJetContainer(&jetRecListNoCut, fhnRecJetsNoCut);
 380   FillJetContainer(&jetGenListNoCut, fhnGenJetsNoCut);
 381   FillJetContainer(&jetRecListCut,   fhnRecJetsCut);
 382   FillJetContainer(&jetGenListCut,   fhnGenJetsCut);
 383
 384     if(trackUEcorrJetRecListCut.GetEntries() != 0  && trackUEcorrJetGenListCut.GetEntries() != 0 ) {
 385         MatchJets(kTRUE,  &trackUEcorrJetRecListCut,    &trackUEcorrJetGenListCut, 0.1,  fhnTrackCorrMatching);
 386         MatchJets(kTRUE,  &trackUEcorrJetRecListCut,    &trackUEcorrJetGenListCut, 0.3,  fhnTrackCorrMatching);
 387         MatchJets(kTRUE,  &trackUEcorrJetRecListCut,    &trackUEcorrJetGenListCut, 0.5,  fhnTrackCorrMatching);
 388     }
 389     if(jetRecListCut.GetEntries() != 0  && jetGenListCut.GetEntries() != 0 ) {
 390         MatchJets(kTRUE,  &jetRecListCut,               &jetGenListCut,           0.1,  fhnMatching);
 391         MatchJets(kTRUE,  &jetRecListCut,               &jetGenListCut,           0.3,  fhnMatching);
 392         MatchJets(kTRUE,  &jetRecListCut,               &jetGenListCut,           0.5,  fhnMatching);
 393     }
 394
 395
 396   PostData(1,fOutputList);
 397   return;
 398 }
 399
 400 //---------------------------------------------------------------------------------------------------
 401 void AliAnalysisTaskPPJetSpectra::Terminate(Option_t *)
 402 {
 403   if(fDebug) printf("%s: %d Terminating\n",(char*)__FILE__, __LINE__);
 404   return;
 405 }
 406
 407 //---------------------------------------------------------------------------------------------------
 408 void AliAnalysisTaskPPJetSpectra::SetRecJetBranch(TString s) {
 409   fRecJetBranch = s;
 410   if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fRecJetR = 0.2;
 411   else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fRecJetR = 0.3;
 412   else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fRecJetR = 0.4;
 413   else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fRecJetR = 0.5;
 414   else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fRecJetR = 0.6;
 415   else fRecJetR = 0;
 416   return;
 417 }
 418
 419 //---------------------------------------------------------------------------------------------------
 420 void AliAnalysisTaskPPJetSpectra::SetGenJetBranch(TString s) {
 421   fGenJetBranch = s;
 422   if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fGenJetR = 0.2;
 423   else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fGenJetR = 0.3;
 424   else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fGenJetR = 0.4;
 425   else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fGenJetR = 0.5;
 426   else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fGenJetR = 0.6;
 427   else fGenJetR = 0;
 428   return;
 429 }
 430
 431 //---------------------------------------------------------------------------------------------------
 432 void AliAnalysisTaskPPJetSpectra::SetRecBckgBranch(TString s) {
 433   fRecBckgBranch = s;
 434   if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fRecBckgR = 0.2;
 435   else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fRecBckgR = 0.3;
 436   else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fRecBckgR = 0.4;
 437   else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fRecBckgR = 0.5;
 438   else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fRecBckgR = 0.6;
 439   else fRecBckgR = 0;
 440   return;
 441 }
 442
 443 //---------------------------------------------------------------------------------------------------
 444 void AliAnalysisTaskPPJetSpectra::SetGenBckgBranch(TString s) {
 445   fGenBckgBranch = s;
 446   if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fGenBckgR = 0.2;
 447   else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fGenBckgR = 0.3;
 448   else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fGenBckgR = 0.4;
 449   else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fGenBckgR = 0.5;
 450   else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fGenBckgR = 0.6;
 451   else fGenBckgR = 0;
 452   return;
 453 }
 454
 455 //---------------------------------------------------------------------------------------------------
 456 void AliAnalysisTaskPPJetSpectra::SetVertexCuts(Int_t nCont, Double_t Zcut, Double_t Rcut) {
 457   if(fDebug)
 458     printf("%s: %d Setting Vertex Cuts\n",(char*)__FILE__,__LINE__);
 459   nVtxContCut = nCont;
 460   fVtxZcut = Zcut;
 461   fVtxRcut = Rcut;
 462   return;
 463 }
 464
 465 //---------------------------------------------------------------------------------------------------
 466 void AliAnalysisTaskPPJetSpectra::SetTrackCuts(Double_t ptMin, Double_t ptMax, Double_t etaMax) {
 467   if(fDebug > 1)
 468     printf("%s: %d Track cuts set\n",(char*)__FILE__,__LINE__);
 469   trackPtMin = ptMin;
 470   trackPtMax = ptMax;
 471   trackEtaAbsMax = etaMax;
 472   return;
 473 }
 474
 475 //---------------------------------------------------------------------------------------------------
 476 void AliAnalysisTaskPPJetSpectra::SetJetCuts(Double_t ptMin, Double_t eta, Double_t z) {
 477   jetPtMin = ptMin;
 478   jetEtaCut = eta;
 479   jetZmax = z;
 480   return;
 481 }
 482
 483 //---------------------------------------------------------------------------------------------------
 484 Bool_t AliAnalysisTaskPPJetSpectra::EventSelection(Double_t evtContainer[6]) {
 485   if(fDebug > 1)
 486     printf("%s: %d UserExec analysing event %d.\n", (char*)__FILE__, __LINE__, (Int_t)fEntry);
 487
 488   // Trigger selection
 489   AliInputEventHandler* inputHandler = (AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 490   if(fDebug > 2 )
 491     printf("IsEventSelected: %d eventSelectionMask: %d\n", (Int_t)inputHandler->IsEventSelected(), (Int_t)fEvtSelectionMask);
 492
 493   Bool_t isTriggerSelected = inputHandler->IsEventSelected() & fEvtSelectionMask;
 494   evtContainer[0] = (Int_t)isTriggerSelected;
 495
 496   fESD = (AliESDEvent*)InputEvent();
 497   if(!fESD && fDebug > 2)
 498     printf("%s: %d No ESD event found\n",(char*)__FILE__, __LINE__);
 499
 500   fMC = MCEvent();
 501   if(!fMC && fDebug > 2)
 502     printf("%s: %d No MC event found\n", (char*)__FILE__, __LINE__);
 503
 504
 505   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 506   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
 507     fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
 508     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 509     fAODJets = ((AliAODHandler*)handler)->GetAOD();
 510     if (fDebug > 1)
 511       Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
 512   }
 513   else {
 514     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 515     if( handler && handler->InheritsFrom("AliAODHandler") ) {
 516       fAOD = ((AliAODHandler*)handler)->GetAOD();
 517       fAODJets = fAOD;
 518       if (fDebug > 1)
 519         Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
 520     }
 521   }
 522
 523   if(!fAODJets)
 524   {
 525     if(fDebug > 1)
 526       printf("%s: %d No AOD found\n",(char*)__FILE__,__LINE__);
 527     return kFALSE;
 528   }
 529   handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 530   Float_t centrality = -1;
 531   if(fEventClass > 0)
 532   {
 533     if(handler->InheritsFrom("AliAODHandler")) centrality = fAOD->GetHeader()->GetCentrality();
 534     else if(fESD) centrality = fESD->GetCentrality()->GetCentralityPercentile("V0M");
 535     else centrality = AliAnalysisHelperJetTasks::EventClass();
 536   }
 537   evtContainer[1] = centrality;
 538
 539   // Event selection
 540   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
 541   if(!primVtx) {
 542     evtContainer[2]  =  1e10;
 543     evtContainer[3]  =  1e10;
 544     evtContainer[4]  =  1e10;
 545     return kFALSE;
 546   } else {
 547     primVtx->Print();
 548     evtContainer[2] = primVtx->GetZ();
 549     evtContainer[3] = TMath::Sqrt(primVtx->GetX()*primVtx->GetX() + primVtx->GetY()*primVtx->GetY() );
 550     evtContainer[4] = primVtx->GetNContributors();
 551   }
 552
 553   if(fRejectPileUp && AliAnalysisHelperJetTasks::IsPileUp()){
 554     if (fDebug > 1) Printf("%s:%d SPD pileup: event REJECTED...", (char*)__FILE__,__LINE__);
 555     evtContainer[5] = 0.;
 556     return kFALSE;
 557   }
 558   else
 559     evtContainer[5] = 1.;
 560
 561   if(evtContainer[0] == 0)
 562   {
 563     if(fDebug > 2)
 564       printf("%s: %d Event rejected: Trigger\n",(char*)__FILE__, __LINE__);
 565     return kFALSE;
 566   }
 567   if(evtContainer[1] >= 0 && evtContainer[1] < 10*((Int_t)fEventClass/10) && evtContainer[1] > 10*((Int_t)fEventClass/10 + 1) )
 568   {
 569     if(fDebug > 2)
 570       printf("%s: %d Event rejected: Centrality\n",(char*)__FILE__, __LINE__);
 571     return kFALSE;
 572   }
 573   if(TMath::Abs(evtContainer[2]) > fVtxZcut)
 574   {
 575     if(fDebug > 2)
 576       printf("%s: %d Event rejected: Vertex position (z)\n",(char*)__FILE__, __LINE__);
 577     return kFALSE;
 578   }
 579   if(evtContainer[3] > fVtxRcut)
 580   {
 581     if(fDebug > 2)
 582       printf("%s: %d Event rejected: Vertex position (x-y)\n",(char*)__FILE__,__LINE__);
 583     return kFALSE;
 584   }
 585   if(evtContainer[4] < nVtxContCut)
 586   {
 587     if(fDebug > 2)
 588       printf("%s: %d Event rejected: Vertex contributors \n",(char*)__FILE__, __LINE__);
 589     return kFALSE;
 590   }
 591
 592   return kTRUE;
 593 }
 594
 595 //---------------------------------------------------------------------------------------------------
 596 Int_t AliAnalysisTaskPPJetSpectra::GetListOfTracks(Int_t trackType, TList *trackList) {
 597   if(fDebug)
 598     printf("%s: %d Getting Tracks\n",(char*)__FILE__,__LINE__);
 599   if(trackType == AliAnalysisTaskPPJetSpectra::kNone) {
 600     if(fDebug) printf("%s: %d No track type selected", (char*)__FILE__,__LINE__);
 601     return 0;
 602   }
 603
 604   if(1 == trackType) {
 605     for(Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
 606       Bool_t isOK = kTRUE;
 607       Bool_t isOK16 = kTRUE;
 608       Bool_t isOK256 = kTRUE;
 609       AliAODTrack* track = (AliAODTrack*)fAOD->GetTrack(i);
 610       if(!track) continue;
 611       if(nTrackFilter>0 && !track->TestFilterBit(nTrackFilter) ) isOK = kFALSE;
 612       if(nTrackFilter>0 && !track->TestFilterBit(16) ) isOK16 = kFALSE;
 613       if(nTrackFilter>0 && !track->TestFilterBit(256) ) isOK256 = kFALSE;
 614       if(track->Pt() < trackPtMin || track->Pt() > trackPtMax) isOK = kFALSE;
 615       if(trackEtaAbsMax < TMath::Abs(track->Eta())) isOK = kFALSE;
 616       Double_t tmpContainer[] = {track->Pt(), track->Eta(), track->Phi(), (Double_t)isOK, (Double_t)isOK16, (Double_t)isOK256, (Double_t)track->Charge()/3};
 617       if(isOK) trackList->Add(track);
 618       fhnTracks->Fill(tmpContainer);
 619     }
 620   }
 621   else {
 622     TClonesArray *tca=dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 623     for(Int_t i = 0; i < tca->GetEntries(); i++) {
 624       AliAODMCParticle *particle = (AliAODMCParticle*)tca->At(i);
 625       if(!particle) continue;
 626       if(!particle->IsPhysicalPrimary()) continue;
 627       Double_t charge = particle->Charge();
 628       if( charge > 0 ) charge = 1;
 629       if( charge < 0 ) charge = -1;
 630       if( (trackType == 3) && (charge == 0) ) continue;
 631       trackList->Add(particle);
 632       Double_t tmpContainer[] = {particle->Pt(), particle->Eta(), particle->Phi(), -999, -999, -999, charge };
 633       if(trackType == 2) fhnMC->Fill(tmpContainer);
 634       else fhnMC2->Fill(tmpContainer);
 635     }
 636   }
 637   trackList->Sort();
 638   return trackList->GetEntries();
 639 }
 640
 641 //---------------------------------------------------------------------------------------------------
 642 Int_t AliAnalysisTaskPPJetSpectra::GetTCAJets(TClonesArray* jca, TString branch) {
 643   if( branch.Length() == 0) {
 644     if(fDebug) printf("%s: %d No branch \"%s\" selected\n", (char*)__FILE__,__LINE__, branch.Data() );
 645     return 0;
 646   }
 647
 648   jca = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(branch.Data()));
 649   if(!jca)                jca = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(branch.Data()));
 650   if(!jca)                jca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(branch.Data()));
 651   if(!jca)                jca = dynamic_cast<TClonesArray*>(fAOD->GetList()->FindObject(branch.Data()));
 652   if(fAODExtension&&!jca) jca = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(branch.Data()));
 653   return jca->GetEntries();
 654
 655 }
 656
 657 //---------------------------------------------------------------------------------------------------
 658 Int_t AliAnalysisTaskPPJetSpectra::GetListOfJets(TClonesArray* jca, TList* jetList, Bool_t useCuts) {
 659   if(fDebug)
 660     printf("%s: %d Getting list of jets\n",(char*)__FILE__,__LINE__);
 661
 662   Double_t ptCut = (useCuts   &&  jetPtMin>0 ?jetPtMin :0.1);
 663   Double_t etaCut = (useCuts  &&  jetEtaCut>0?jetEtaCut:0.9);
 664   Double_t zCut  = (useCuts   &&  jetZmax>0  ?jetZmax  :1);
 665
 666   if(!jca){
 667     return -2;
 668   }
 669
 670   for(Int_t i =0 ; i < jca->GetEntries(); i++) {
 671     AliAODJet* jet = (AliAODJet*)jca->At(i);
 672     if(!jet) continue;
 673     Double_t zLeading = -1;
 674 //    if(jet->GetPtLeading() <= 0) {
 675 //      TRefArray* jetRef = (TRefArray*)jet->GetRefTracks();
 676 //      for(int j = 0; j < jetRef->GetEntries(); j++) {
 677 //      AliVTrack *track = (AliVTrack*)jet->GetTrack(j);
 678 //        if(!track) continue;
 679 //        cout<<track->Pt()<<" "<<zLeading<<endl;
 680 //      if(track->Pt() > zLeading) zLeading = track->Pt();
 681 //        cout<<__LINE__<<endl;
 682 //      }
 683 //      if(jet->Pt() > 0) zLeading /= jet->Pt();
 684 //      else zLeading = -1;
 685 //    }
 686     //cuty
 687     if(ptCut > 0  && ptCut > jet->Pt()                ) continue;
 688     if(etaCut > 0 && etaCut < TMath::Abs(jet->Eta())  ) continue;
 689     if(zCut > 0   && zCut < zLeading  && zLeading > 0 ) continue;
 690
 691     jetList->Add(jet);
 692   }
 693
 694   jetList->Sort();
 695   return jetList->GetEntries();
 696 }
 697
 698 //---------------------------------------------------------------------------------------------------
 699 void AliAnalysisTaskPPJetSpectra::FillJetContainer(TList* jets, THnSparseF* container) {
 700   if(!jets || !container) return;
 701
 702   for(Int_t i =0 ; i < jets->GetEntries(); i++) {
 703     AliAODJet* jet = (AliAODJet*)jets->At(i);
 704     if(!jet) continue;
 705   // pt - eta - phi - AreaCH - AreaN -  ptLeadingTrack/ptJet - ptOrig
 706     Double_t tmpContainer[] = {jet->Pt(), jet->Eta(), jet->Phi(), jet->EffectiveAreaCharged(), jet->EffectiveAreaNeutral(), jet->GetPtLeading()/jet->Pt(), (Double_t)((TRefArray*)jet->GetRefTracks())->GetEntries() };
 707     container->Fill(tmpContainer);
 708   }
 709   return;
 710 }
 711
 712 //---------------------------------------------------------------------------------------------------
 713 Double_t AliAnalysisTaskPPJetSpectra::GetUE(TList* jets, TList* tracks, Double_t Rpar,THnSparseF* thnCont) {
 714   if(!jets || !tracks) {
 715     if(fDebug) printf("%s: %d No jets or tracks\n",(char*)__FILE__,__LINE__);
 716     return 0;
 717   }
 718
 719   //find leading jet
 720   AliAODJet *jet = 0;
 721   for(Int_t iTmp = 0; iTmp < jets->GetEntries(); iTmp++) {
 722     AliAODJet* tmp = (AliAODJet*)jets->At(iTmp);
 723     if(!tmp) continue;
 724     if(!jet) jet = tmp;
 725     else
 726       if( jet->Pt() < tmp->Pt() ) jet = tmp;
 727   }
 728   if(!jet) return 0;
 729
 730   //get track in perpendicular cone
 731   //TVector3 for jet and perpendicular cones with same eta axis
 732   TVector3 j(jet->Px(), jet->Py(), jet->Pz());
 733   TVector3 p1(j);
 734   TVector3 p2(j);
 735
 736   p1.RotateZ(TMath::Pi()/2.);
 737   p2.RotateZ(-TMath::Pi()/2.);
 738
 739   Double_t sumAllPt1 = 0;
 740   Double_t sumAllPt2 = 0;
 741
 742   for(int i = 0; i < tracks->GetEntries(); i++)  {
 743     AliVParticle* tr = (AliVParticle*)tracks->At(i);
 744     if(!tr) continue;
 745     TVector3 v(tr->Px(), tr->Py(), tr->Pz());
 746     Double_t dR1 = v.DrEtaPhi(p1);
 747     Double_t dR2 = v.DrEtaPhi(p2);
 748
 749     if(dR1 < Rpar) sumAllPt1+=v.Pt();
 750     if(dR2 < Rpar) sumAllPt2+=v.Pt();
 751   }
 752
 753   if(Rpar != 0) {
 754     sumAllPt1/=(TMath::Pi()*Rpar*Rpar);
 755     sumAllPt2/=(TMath::Pi()*Rpar*Rpar);
 756   }
 757
 758   Double_t container[] = {j.Pt(), j.Eta(), sumAllPt1, sumAllPt2, (sumAllPt1+sumAllPt2)/2.};
 759   thnCont->Fill(container);
 760
 761   return (sumAllPt1+sumAllPt2)/2.;
 762
 763 }
 764
 765 //---------------------------------------------------------------------------------------------------
 766 Double_t AliAnalysisTaskPPJetSpectra::GetBckgUE(TList* jets, Double_t Rpar, Bool_t isSkipped, THnSparseF* cont) {
 767   if(!jets) return 0;
 768
 769   Int_t leading01 = -1;
 770   Int_t leading02 = -1;
 771   Double_t leading01pt = 0;
 772   Double_t leading02pt = 0;
 773
 774   if(!isSkipped) {
 775     for(Int_t i = 0; i < jets->GetEntries(); i++) {
 776       AliAODJet* tmp = (AliAODJet*)jets->At(i);
 777       if(!tmp) continue;
 778       if(tmp->Pt() > leading01pt) {
 779         leading02 = leading01;
 780         leading01 = i;
 781         leading02pt = leading01pt;
 782         leading01pt = tmp->Pt();
 783       } else if(tmp->Pt() > leading02) {
 784         leading02 = i;
 785         leading02pt = tmp->Pt();
 786       }
 787     }
 788   }
 789
 790   Int_t subLeading01 = -1;
 791   Int_t subLeading02 = -1;
 792   Double_t subLeading01pt = 0;
 793   Double_t subLeading02pt = 0;
 794   Double_t subLeading01area = 0;
 795   Double_t subLeading02area = 0;
 796
 797   Double_t avgAll = 0;
 798   Double_t areaAll = 0;
 799   //search for subleading jets assuming leading 2 (skip02) jets are skipped
 800   if(jets->GetEntries() > 0) {
 801     for(Int_t i = 0; i < jets->GetEntries(); i++) {
 802       if(i == leading01 || i == leading02) continue;
 803       AliAODJet* tmp = (AliAODJet*)jets->At(i);
 804       if(!tmp) continue;
 805       if( subLeading01pt < tmp->Pt() ) {
 806         subLeading02 = subLeading01;
 807         subLeading01 = i;
 808         subLeading02pt = subLeading01pt;
 809         subLeading01pt = tmp->Pt();
 810         subLeading02area = subLeading01area;
 811         subLeading01area = tmp->EffectiveAreaCharged();
 812       } else if( subLeading02pt < tmp->Pt() ) {
 813         subLeading02 = i;
 814         subLeading02pt = tmp->Pt();
 815         subLeading02area = tmp->EffectiveAreaCharged();
 816       }
 817     }
 818
 819     //do avg jet pt
 820     Int_t nJets = 0;
 821     for(int i = 0; i < jets->GetEntries();i++) {
 822       if(i == leading01 || i == leading02) continue;
 823       AliAODJet* tmp = (AliAODJet*)jets->At(i);
 824       if(!tmp) continue;
 825       nJets++;
 826       avgAll+= tmp->Pt();
 827       areaAll += tmp->EffectiveAreaCharged();
 828     }
 829   }
 830
 831
 832   if(avgAll!=0) avgAll/=areaAll;
 833   if(subLeading01pt!=0) subLeading01pt/=subLeading01area;
 834   if(subLeading02pt!=0)subLeading02pt/=subLeading02area;
 835   Double_t tmpContainer[]={subLeading01pt,subLeading02pt,avgAll,Rpar*10};
 836
 837   cont->Fill(tmpContainer);
 838
 839   return avgAll;
 840 }
 841
 842 //---------------------------------------------------------------------------------------------------
 843 Int_t  AliAnalysisTaskPPJetSpectra::CorrectForUE(TList* jets, Double_t UE, TList* newList, THnSparseF* container) {
 844   if(!jets) return 0;
 845
 846   for(Int_t i = 0; i < jets->GetEntries(); i++) {
 847     AliAODJet* tmpJet = (AliAODJet*)jets->At(i);
 848     if(!tmpJet) continue;
 849
 850     AliAODJet* jet = (AliAODJet*)tmpJet->Clone(Form("%s%s",tmpJet->GetName(),"_clone"));
 851     Double_t corrUE = UE*(jet->EffectiveAreaCharged() + jet->EffectiveAreaNeutral() );
 852
 853     if(jet->Pt() > corrUE) jet->SetPtEtaPhiM(jet->Pt()- corrUE, jet->Eta(), jet->Phi(), jet->M());
 854     else jet->SetPtEtaPhiM(0, jet->Eta(), jet->Phi(), jet->M());
 855     newList->Add(jet);
 856
 857   // pt - eta - phi - AreaCH - AreaN - UEdensity - ptLeadingTrack/ptJet - ptOrig
 858     Double_t cont[8]={jet->Pt(), jet->Eta(), jet->Phi(), jet->EffectiveAreaCharged(),jet->EffectiveAreaNeutral(), corrUE, (jet->Pt() > 0?jet->GetPtLeading() / jet->Pt():-1),tmpJet->Pt()};
 859     container->Fill(cont);
 860   }
 861   newList->Sort();
 862   return newList->GetEntries();
 863 }
 864
 865 //---------------------------------------------------------------------------------------------------
 866 void AliAnalysisTaskPPJetSpectra::MatchJets(Bool_t doBothWay, TList* recList, TList* genList, Float_t maxDist, THnSparseF* container) {
 867   if(!genList || !recList) return;
 868   Int_t mode = 1;
 869
 870   const Int_t kGenJets = genList->GetEntries();
 871   const Int_t kRecJets = recList->GetEntries();
 872
 873   TArrayI iMatchIndex(kGenJets);
 874   TArrayF fPtFraction(kGenJets);
 875
 876   TArrayI aGenIndex(kRecJets);
 877   TArrayI aRecIndex(kGenJets);
 878
 879   Double_t genPt = -1;
 880   Double_t recPt = -1;
 881   Double_t ptFraction = -1;
 882
 883   //-- Closest jets to generated and checked for E-fraction
 884   if(!doBothWay) {
 885     AliAnalysisHelperJetTasks::GetJetMatching(genList, kGenJets, recList, kRecJets, iMatchIndex, fPtFraction, fDebug, maxDist, mode);
 886     for(Int_t ig = 0; ig < kGenJets; ig++) {
 887       //get gen-jet
 888       AliAODJet* genJet = (AliAODJet*)genList->At(ig);
 889       if(!genJet) continue;
 890       genPt = genJet->Pt();
 891
 892       //get rec-jet
 893       Int_t ir = iMatchIndex[ig];
 894       if(ir<0||ir>=recList->GetEntries()) continue;
 895       AliAODJet* recJet = (AliAODJet*)recList->At(ir);
 896       if(!recJet) continue;
 897       recPt = recJet->Pt();
 898
 899       ptFraction = fPtFraction[ig];
 900
 901       Double_t dR = recJet->DeltaR(genJet);
 902       Double_t dEta = recJet->Eta() - genJet->Eta();
 903       Double_t dPhi2 = dR*dR - dEta*dEta;
 904
 905         Bool_t isSame = kFALSE;
 906         if(CheckPtBin(genPt) == CheckPtBin(recPt)) isSame = kTRUE;
 907       // ptGen - ptRec - ptFraction - maxDist - bothWays - (1 - rec/gen) - dR - dEta  - dPhi
 908       Double_t cont[] = {genPt,recPt,ptFraction,maxDist,(Float_t)doBothWay, (genPt != 0 ? (1. - recPt/genPt):1), dR, dEta, (dPhi2>=0?TMath::Sqrt(dPhi2):-1), (Double_t)isSame};
 909       container->Fill(cont);
 910     }
 911   }
 912   //-- Closest rec jet to gen and vice-versa & check for E-fraction
 913   else {
 914     AliAnalysisHelperJetTasks::GetClosestJets(genList, kGenJets, recList, kRecJets, aGenIndex, aRecIndex,fDebug, maxDist);
 915     for(Int_t ig = 0; ig < kGenJets; ig++) {
 916       //get gen-jet
 917       AliAODJet* genJet = (AliAODJet*)genList->At(ig);
 918       if(!genJet) continue;
 919       genPt = genJet->Pt();
 920
 921       //get rec-jet
 922       Int_t ir = aRecIndex[ig];
 923       if(ir>=0&&ir<recList->GetEntries()) {
 924         AliAODJet *recJet = (AliAODJet*)recList->At(ir);
 925         if(!recJet) continue;
 926         recPt = recJet->Pt();
 927         Double_t dR = recJet->DeltaR(genJet);
 928         Double_t dEta = recJet->Eta() - genJet->Eta();
 929         Double_t dPhi2 = dR*dR - dEta*dEta;
 930
 931         ptFraction = AliAnalysisHelperJetTasks::GetFractionOfJet(recJet,genJet,mode);
 932
 933         Bool_t isSame = kFALSE;
 934         if(CheckPtBin(genPt) == CheckPtBin(recPt)) isSame = kTRUE;
 935         // ptGen - ptRec - ptFraction - maxDist - bothWays - (1 - rec/gen) - dR - dEta  - dPhi
 936         Double_t cont[] = {genPt,recPt,ptFraction,maxDist,(Float_t)doBothWay, (genPt != 0 ? (1. - recPt/genPt):1), dR, dEta, (dPhi2>=0?TMath::Sqrt(dPhi2):-1), (Double_t)isSame  };
 937         container->Fill(cont);
 938       }
 939     }
 940   }
 941
 942   return;
 943 }
 944
 945 //---------------------------------------------------------------------------------------------------
 946 Int_t AliAnalysisTaskPPJetSpectra::CheckPtBin(Double_t pT) {
 947   Double_t bins[]       = {0.,2.,4.,6.,8.,10.,12., 14.,16., 18.,20.,24.,28.,32.,38.,44.,50.,58.,66.,76.,86.,100.,114.,128.,150.,250.};
 948   Double_t nBins        = 25;
 949
 950   for(Int_t i = 0; i <= nBins; i++) {
 951     if(pT < bins[i] ) return i;
 952   }
 953   return nBins+1;
 954 }
 955
 956 //---------------------------------------------------------------------------------------------------
 957 void AliAnalysisTaskPPJetSpectra::DoUEAnalysis(TList* tracks, Double_t PTcut, Double_t ETAcut) {
 958         Double_t PTmax = 0;
 959         int iMax = -1;
 960         for(int i = 0; i < tracks->GetEntries(); i++) {
 961                 AliVParticle *part = (AliVParticle*)tracks->At(i);
 962                 if(!part) continue;
 963
 964                 if(part->Pt() < PTcut) continue;
 965                 if(TMath::Abs(part->Eta()) > ETAcut ) continue;
 966
 967                 if(PTmax < part->Pt()) {
 968                         PTmax = part->Pt();
 969                         iMax = i;
 970                 }
 971         }
 972
 973         if(PTmax < trackPtMin) return;
 974
 975         Double_t ch_part_density = 0;
 976         Double_t ch_part_pt_density = 0;
 977         Double_t avg_track_pt = 0;
 978         Double_t zero_ch_part = 0;
 979         Double_t transMAX_part_density = 0;
 980         Double_t transMIN_part_density = 0;
 981         Double_t transMAX_pt_density = 0;
 982         Double_t transMIN_pt_density = 0;
 983         Double_t transDIF_pt_density = 0;
 984         Double_t transDIF_part_density = 0;
 985
 986         Double_t A1_part_density = 0;
 987         Double_t A2_part_density = 0;
 988         Double_t A1_pt_density = 0;
 989         Double_t A2_pt_density = 0;
 990
 991         AliVParticle *maxParticle = (AliVParticle*)tracks->At(iMax);
 992
 993         Double_t PHImax = maxParticle->Phi();
 994
 995         for(int i = 0; i < tracks->GetEntries(); i++) {
 996                 AliVParticle* part = (AliVParticle*)tracks->At(i);
 997                 if(!part) continue;
 998                 if(part->Pt() < PTcut) continue;
 999
1000                 Bool_t isArea1 = kFALSE;
1001                 Bool_t isArea2 = kFALSE;
1002
1003                 Double_t dPhi = part->Phi() - PHImax;
1004                 if(dPhi > TMath::Pi()) dPhi -= TMath::Pi();
1005                 if(dPhi < -TMath::Pi()) dPhi +=TMath::Pi();
1006
1007                 if(TMath::Abs(dPhi) > TMath::Pi()/3 && TMath::Abs(dPhi) < 2*TMath::Pi()/3 ) {
1008                         if(dPhi>0) isArea1 = kTRUE;
1009                         else isArea2 = kTRUE;
1010                 }
1011
1012                 if(!isArea1 && !isArea2) continue;
1013
1014                 ch_part_density += 1;
1015                 ch_part_pt_density += part->Pt();
1016                 avg_track_pt += part->Pt();
1017
1018                 if(isArea1) {
1019                         A1_part_density += 1;
1020                         A1_pt_density += part->Pt();
1021                 } else if(isArea2) {
1022                         A2_part_density += 1;
1023                         A2_pt_density += part->Pt();
1024                 }
1025         }
1026
1027
1028         if(ch_part_density > 0) avg_track_pt/=ch_part_density;
1029         else zero_ch_part = 1;
1030         ch_part_density*=3./(4.*TMath::Pi()*ETAcut);
1031         ch_part_pt_density*= 3./(4.*TMath::Pi()*ETAcut);
1032
1033         A1_part_density*=3./(4.*TMath::Pi()*ETAcut);
1034         A1_pt_density*= 3./(4.*TMath::Pi()*ETAcut);
1035         A2_part_density*=3./(4.*TMath::Pi()*ETAcut);
1036         A2_pt_density*= 3./(4.*TMath::Pi()*ETAcut);
1037
1038         if(A1_part_density < A2_part_density) {
1039                 transMAX_part_density = A2_part_density;
1040                 transMIN_part_density = A1_part_density;
1041         } else {
1042                 transMAX_part_density = A1_part_density;
1043                 transMIN_part_density = A2_part_density;
1044         }
1045
1046         if(A1_pt_density < A2_pt_density) {
1047                 transMAX_pt_density = A2_pt_density;
1048                 transMIN_pt_density = A1_pt_density;
1049         } else {
1050                 transMAX_pt_density = A1_pt_density;
1051                 transMIN_pt_density = A2_pt_density;
1052         }
1053
1054         transDIF_pt_density     = transMAX_pt_density - transMIN_pt_density;
1055         transDIF_part_density   = transMAX_part_density - transMIN_part_density;
1056
1057         Double_t container[] = {PTmax, ch_part_density, ch_part_pt_density, avg_track_pt, zero_ch_part,
1058                                                         transMAX_part_density, transMIN_part_density, transDIF_part_density,
1059                                                         transMAX_pt_density, transMIN_pt_density, transDIF_pt_density};
1060
1061         fhnTrackUEanal->Fill(container);
1062 }