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   ,fAODIn(0)
  26   ,fAODOut(0)
  27   ,fAODExt(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   ,fAODIn(0)
  88   ,fAODOut(0)
  89   ,fAODExt(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   if(fAODOut && !tcaRecJets && fRecJetBranch.Length() > 0) tcaRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fRecJetBranch.Data()));
 305   if(fAODOut && !tcaGenJets && fGenJetBranch.Length() > 0) tcaGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fGenJetBranch.Data()));
 306   if(fAODOut && !tcaRecBckg && fRecBckgBranch.Length() > 0) tcaRecBckg = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fRecBckgBranch.Data()));
 307   if(fAODOut && !tcaGenBckg && fGenBckgBranch.Length() > 0) tcaGenBckg = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fGenBckgBranch.Data()));
 308
 309   if(fAODExt) {
 310     if(!tcaRecJets && fRecJetBranch.Length() > 0) tcaRecJets = dynamic_cast<TClonesArray*>(fAODExt->GetAOD()->FindListObject(fRecJetBranch.Data()));
 311     if(!tcaGenJets && fGenJetBranch.Length() > 0) tcaGenJets = dynamic_cast<TClonesArray*>(fAODExt->GetAOD()->FindListObject(fGenJetBranch.Data()));
 312     if(!tcaRecBckg && fRecBckgBranch.Length() > 0) tcaRecBckg = dynamic_cast<TClonesArray*>(fAODExt->GetAOD()->FindListObject(fRecBckgBranch.Data()));
 313     if(!tcaGenBckg && fGenBckgBranch.Length() > 0) tcaGenBckg = dynamic_cast<TClonesArray*>(fAODExt->GetAOD()->FindListObject(fGenBckgBranch.Data()));
 314   }
 315
 316   if(fAODIn) {
 317     if(!tcaRecJets && fRecJetBranch.Length() > 0) tcaRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fRecJetBranch.Data()));
 318     if(!tcaGenJets && fGenJetBranch.Length() > 0) tcaGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fGenJetBranch.Data()));
 319     if(!tcaRecBckg && fRecBckgBranch.Length() > 0) tcaRecBckg = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fRecBckgBranch.Data()));
 320     if(!tcaGenBckg && fGenBckgBranch.Length() > 0) tcaGenBckg = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fGenBckgBranch.Data()));
 321   }
 322
 323   if(!tcaRecJets) rJ = -1;
 324   else rJ = tcaRecJets->GetEntries();
 325   if(!tcaGenJets) gJ = -1;
 326   else gJ = tcaGenJets->GetEntries();
 327   if(!tcaRecBckg) rB = -1;
 328   else rB = tcaRecBckg->GetEntries();
 329   if(!tcaGenBckg) gB = -1;
 330   else gB = tcaGenBckg->GetEntries();
 331
 332   TList jetRecListNoCut;
 333   TList jetGenListNoCut;
 334   Int_t nRecJetsNoCut   = (rJ>0?GetListOfJets(tcaRecJets, &jetRecListNoCut, kFALSE):0);
 335   Int_t nGenJetsNoCut   = (gJ>0?GetListOfJets(tcaGenJets, &jetGenListNoCut, kFALSE):0);
 336
 337   TList jetRecListCut;
 338   TList jetGenListCut;
 339   Int_t nRecJetsCut     = (rJ>0?GetListOfJets(tcaRecJets, &jetRecListCut, kTRUE):0);
 340   Int_t nGenJetsCut     = (gJ>0?GetListOfJets(tcaGenJets, &jetGenListCut, kTRUE):0);
 341
 342   TList bckgRecListCut;
 343   TList bckgGenListCut;
 344   Int_t nRecBckgCut     = (rB>0?GetListOfJets(tcaRecBckg, &bckgRecListCut, kTRUE):0);
 345   Int_t nGenBckgCut     = (gB>0?GetListOfJets(tcaGenBckg, &bckgGenListCut, kTRUE):0);
 346
 347   if(fDebug > 10) std::cout<<nRecJetsNoCut<<" "<<nGenJetsNoCut<<" "<<nRecJetsCut<< " "<< nGenJetsCut<<" "<<nRecBckgCut<<" "<<nGenBckgCut<<std::endl;
 348
 349   evtContainer[6] = (isEventSelected?1:0);
 350   fhnEvent->Fill(evtContainer);
 351
 352   if(!isEventSelected) {
 353     PostData(1,fOutputList);
 354     return;
 355   }
 356
 357   if(kDoUEanalysis) DoUEAnalysis(&trackList, (Double_t)0.15, (Double_t)0.9);
 358
 359   Double_t  recUEdensity = GetUE(&jetRecListCut, &trackList, fRecJetR, fhnTrackUE);
 360   Double_t  genUEdensity = GetUE(&jetGenListCut, &particleList, fGenJetR, fhnParticleUE);
 361
 362   TList  trackUEcorrJetRecListCut;
 363   TList  trackUEcorrJetGenListCut;
 364
 365   CorrectForUE(&jetRecListCut, recUEdensity, &trackUEcorrJetRecListCut, fhnRecJetsTrackUEcor);
 366   CorrectForUE(&jetGenListCut, genUEdensity, &trackUEcorrJetGenListCut, fhnGenJetsTrackUEcor);
 367
 368   FillJetContainer(&jetRecListNoCut, fhnRecJetsNoCut);
 369   FillJetContainer(&jetGenListNoCut, fhnGenJetsNoCut);
 370   FillJetContainer(&jetRecListCut,   fhnRecJetsCut);
 371   FillJetContainer(&jetGenListCut,   fhnGenJetsCut);
 372
 373     if(trackUEcorrJetRecListCut.GetEntries() != 0  && trackUEcorrJetGenListCut.GetEntries() != 0 ) {
 374         MatchJets(kTRUE,  &trackUEcorrJetRecListCut,    &trackUEcorrJetGenListCut, 0.1,  fhnTrackCorrMatching);
 375         MatchJets(kTRUE,  &trackUEcorrJetRecListCut,    &trackUEcorrJetGenListCut, 0.3,  fhnTrackCorrMatching);
 376         MatchJets(kTRUE,  &trackUEcorrJetRecListCut,    &trackUEcorrJetGenListCut, 0.5,  fhnTrackCorrMatching);
 377     }
 378     if(jetRecListCut.GetEntries() != 0  && jetGenListCut.GetEntries() != 0 ) {
 379         MatchJets(kTRUE,  &jetRecListCut,               &jetGenListCut,           0.1,  fhnMatching);
 380         MatchJets(kTRUE,  &jetRecListCut,               &jetGenListCut,           0.3,  fhnMatching);
 381         MatchJets(kTRUE,  &jetRecListCut,               &jetGenListCut,           0.5,  fhnMatching);
 382     }
 383
 384
 385   PostData(1,fOutputList);
 386   return;
 387 }
 388
 389 //---------------------------------------------------------------------------------------------------
 390 void AliAnalysisTaskPPJetSpectra::Terminate(Option_t *)
 391 {
 392   if(fDebug) printf("%s: %d Terminating\n",(char*)__FILE__, __LINE__);
 393   return;
 394 }
 395
 396 //---------------------------------------------------------------------------------------------------
 397 void AliAnalysisTaskPPJetSpectra::SetRecJetBranch(TString s) {
 398   fRecJetBranch = s;
 399   if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fRecJetR = 0.2;
 400   else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fRecJetR = 0.3;
 401   else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fRecJetR = 0.4;
 402   else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fRecJetR = 0.5;
 403   else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fRecJetR = 0.6;
 404   else fRecJetR = 0;
 405   return;
 406 }
 407
 408 //---------------------------------------------------------------------------------------------------
 409 void AliAnalysisTaskPPJetSpectra::SetGenJetBranch(TString s) {
 410   fGenJetBranch = s;
 411   if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fGenJetR = 0.2;
 412   else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fGenJetR = 0.3;
 413   else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fGenJetR = 0.4;
 414   else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fGenJetR = 0.5;
 415   else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fGenJetR = 0.6;
 416   else fGenJetR = 0;
 417   return;
 418 }
 419
 420 //---------------------------------------------------------------------------------------------------
 421 void AliAnalysisTaskPPJetSpectra::SetRecBckgBranch(TString s) {
 422   fRecBckgBranch = s;
 423   if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fRecBckgR = 0.2;
 424   else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fRecBckgR = 0.3;
 425   else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fRecBckgR = 0.4;
 426   else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fRecBckgR = 0.5;
 427   else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fRecBckgR = 0.6;
 428   else fRecBckgR = 0;
 429   return;
 430 }
 431
 432 //---------------------------------------------------------------------------------------------------
 433 void AliAnalysisTaskPPJetSpectra::SetGenBckgBranch(TString s) {
 434   fGenBckgBranch = s;
 435   if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fGenBckgR = 0.2;
 436   else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fGenBckgR = 0.3;
 437   else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fGenBckgR = 0.4;
 438   else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fGenBckgR = 0.5;
 439   else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fGenBckgR = 0.6;
 440   else fGenBckgR = 0;
 441   return;
 442 }
 443
 444 //---------------------------------------------------------------------------------------------------
 445 void AliAnalysisTaskPPJetSpectra::SetVertexCuts(Int_t nCont, Double_t Zcut, Double_t Rcut) {
 446   if(fDebug)
 447     printf("%s: %d Setting Vertex Cuts\n",(char*)__FILE__,__LINE__);
 448   nVtxContCut = nCont;
 449   fVtxZcut = Zcut;
 450   fVtxRcut = Rcut;
 451   return;
 452 }
 453
 454 //---------------------------------------------------------------------------------------------------
 455 void AliAnalysisTaskPPJetSpectra::SetTrackCuts(Double_t ptMin, Double_t ptMax, Double_t etaMax) {
 456   if(fDebug > 1)
 457     printf("%s: %d Track cuts set\n",(char*)__FILE__,__LINE__);
 458   trackPtMin = ptMin;
 459   trackPtMax = ptMax;
 460   trackEtaAbsMax = etaMax;
 461   return;
 462 }
 463
 464 //---------------------------------------------------------------------------------------------------
 465 void AliAnalysisTaskPPJetSpectra::SetJetCuts(Double_t ptMin, Double_t eta, Double_t z) {
 466   jetPtMin = ptMin;
 467   jetEtaCut = eta;
 468   jetZmax = z;
 469   return;
 470 }
 471
 472 //---------------------------------------------------------------------------------------------------
 473 Bool_t AliAnalysisTaskPPJetSpectra::EventSelection(Double_t evtContainer[6]) {
 474   if(fDebug > 1)
 475     printf("%s: %d UserExec analysing event %d.\n", (char*)__FILE__, __LINE__, (Int_t)fEntry);
 476
 477   // Trigger selection
 478   AliInputEventHandler* inputHandler = (AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 479   if(fDebug > 2 )
 480     printf("IsEventSelected: %d eventSelectionMask: %d\n", (Int_t)inputHandler->IsEventSelected(), (Int_t)fEvtSelectionMask);
 481
 482   Bool_t isTriggerSelected = inputHandler->IsEventSelected() & fEvtSelectionMask;
 483   evtContainer[0] = (Int_t)isTriggerSelected;
 484
 485   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 486   if(!fESD && fDebug > 2)
 487     printf("%s: %d No ESD event found\n",(char*)__FILE__, __LINE__);
 488
 489   if(!fESD) fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 490   fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 491
 492   if(!fESD)fAOD = fAODIn;
 493   else fAOD = fAODOut;
 494
 495   if(fNonStdFile.Length()!=0) {
 496     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 497     fAODExt = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 498   }
 499
 500   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 501   Float_t centrality = -1;
 502   if(fEventClass > 0)
 503   {
 504     if(handler->InheritsFrom("AliAODHandler")) centrality = ((AliVAODHeader*)fAODIn->GetHeader())->GetCentrality();
 505     else if(fESD) centrality = fESD->GetCentrality()->GetCentralityPercentile("V0M");
 506     else centrality = AliAnalysisHelperJetTasks::EventClass();
 507   }
 508   evtContainer[1] = centrality;
 509
 510   // Event selection
 511   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
 512   if(!primVtx) {
 513     evtContainer[2]  =  1e10;
 514     evtContainer[3]  =  1e10;
 515     evtContainer[4]  =  1e10;
 516     return kFALSE;
 517   } else {
 518     primVtx->Print();
 519     evtContainer[2] = primVtx->GetZ();
 520     evtContainer[3] = TMath::Sqrt(primVtx->GetX()*primVtx->GetX() + primVtx->GetY()*primVtx->GetY() );
 521     evtContainer[4] = primVtx->GetNContributors();
 522   }
 523
 524   if(fRejectPileUp && AliAnalysisHelperJetTasks::IsPileUp()){
 525     if (fDebug > 1) Printf("%s:%d SPD pileup: event REJECTED...", (char*)__FILE__,__LINE__);
 526     evtContainer[5] = 0.;
 527     return kFALSE;
 528   }
 529   else
 530     evtContainer[5] = 1.;
 531
 532   if(evtContainer[0] == 0)
 533   {
 534     if(fDebug > 2)
 535       printf("%s: %d Event rejected: Trigger\n",(char*)__FILE__, __LINE__);
 536     return kFALSE;
 537   }
 538   if(evtContainer[1] >= 0 && evtContainer[1] < 10*((Int_t)fEventClass/10) && evtContainer[1] > 10*((Int_t)fEventClass/10 + 1) )
 539   {
 540     if(fDebug > 2)
 541       printf("%s: %d Event rejected: Centrality\n",(char*)__FILE__, __LINE__);
 542     return kFALSE;
 543   }
 544   if(TMath::Abs(evtContainer[2]) > fVtxZcut)
 545   {
 546     if(fDebug > 2)
 547       printf("%s: %d Event rejected: Vertex position (z)\n",(char*)__FILE__, __LINE__);
 548     return kFALSE;
 549   }
 550   if(evtContainer[3] > fVtxRcut)
 551   {
 552     if(fDebug > 2)
 553       printf("%s: %d Event rejected: Vertex position (x-y)\n",(char*)__FILE__,__LINE__);
 554     return kFALSE;
 555   }
 556   if(evtContainer[4] < nVtxContCut)
 557   {
 558     if(fDebug > 2)
 559       printf("%s: %d Event rejected: Vertex contributors \n",(char*)__FILE__, __LINE__);
 560     return kFALSE;
 561   }
 562
 563   return kTRUE;
 564 }
 565
 566 //---------------------------------------------------------------------------------------------------
 567 Int_t AliAnalysisTaskPPJetSpectra::GetListOfTracks(Int_t trackType, TList *trackList) {
 568   if(fDebug)
 569     printf("%s: %d Getting Tracks\n",(char*)__FILE__,__LINE__);
 570   if(trackType == AliAnalysisTaskPPJetSpectra::kNone) {
 571     if(fDebug) printf("%s: %d No track type selected", (char*)__FILE__,__LINE__);
 572     return 0;
 573   }
 574
 575   if(1 == trackType) {
 576     for(Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
 577       Bool_t isOK = kTRUE;
 578       Bool_t isOK16 = kTRUE;
 579       Bool_t isOK256 = kTRUE;
 580       AliAODTrack* track = (AliAODTrack*)fAOD->GetTrack(i);
 581       if(!track) continue;
 582       if(nTrackFilter>0 && !track->TestFilterBit(nTrackFilter) ) isOK = kFALSE;
 583       if(nTrackFilter>0 && !track->TestFilterBit(16) ) isOK16 = kFALSE;
 584       if(nTrackFilter>0 && !track->TestFilterBit(256) ) isOK256 = kFALSE;
 585       if(track->Pt() < trackPtMin || track->Pt() > trackPtMax) isOK = kFALSE;
 586       if(trackEtaAbsMax < TMath::Abs(track->Eta())) isOK = kFALSE;
 587       Double_t tmpContainer[] = {track->Pt(), track->Eta(), track->Phi(), (Double_t)isOK, (Double_t)isOK16, (Double_t)isOK256, (Double_t)track->Charge()/3};
 588       if(isOK) trackList->Add(track);
 589       fhnTracks->Fill(tmpContainer);
 590     }
 591   }
 592   else {
 593     TClonesArray *tca=dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 594     if(!tca) {
 595       if(fDebug > 2) printf("no branch %s\n", AliAODMCParticle::StdBranchName());
 596     }
 597     for(Int_t i = 0; i < tca->GetEntries(); i++) {
 598       AliAODMCParticle *particle = (AliAODMCParticle*)tca->At(i);
 599       if(!particle) continue;
 600       if(!particle->IsPhysicalPrimary()) continue;
 601       Double_t charge = particle->Charge();
 602       if( charge > 0 ) charge = 1;
 603       if( charge < 0 ) charge = -1;
 604       if( (trackType == 3) && (charge == 0) ) continue;
 605       trackList->Add(particle);
 606       Double_t tmpContainer[] = {particle->Pt(), particle->Eta(), particle->Phi(), -999, -999, -999, charge };
 607       if(trackType == 2) fhnMC->Fill(tmpContainer);
 608       else fhnMC2->Fill(tmpContainer);
 609     }
 610   }
 611   trackList->Sort();
 612   return trackList->GetEntries();
 613 }
 614
 615 //---------------------------------------------------------------------------------------------------
 616 Int_t AliAnalysisTaskPPJetSpectra::GetListOfJets(TClonesArray* jca, TList* jetList, Bool_t useCuts) {
 617   if(fDebug)
 618     printf("%s: %d Getting list of jets\n",(char*)__FILE__,__LINE__);
 619
 620   Double_t ptCut = (useCuts   &&  jetPtMin>0 ?jetPtMin :0.1);
 621   Double_t etaCut = (useCuts  &&  jetEtaCut>0?jetEtaCut:0.9);
 622   Double_t zCut  = (useCuts   &&  jetZmax>0  ?jetZmax  :1);
 623
 624   if(!jca){
 625     return -2;
 626   }
 627
 628   for(Int_t i =0 ; i < jca->GetEntries(); i++) {
 629     AliAODJet* jet = (AliAODJet*)jca->At(i);
 630     if(!jet) continue;
 631     Double_t zLeading = -1;
 632 //    if(jet->GetPtLeading() <= 0) {
 633 //      TRefArray* jetRef = (TRefArray*)jet->GetRefTracks();
 634 //      for(int j = 0; j < jetRef->GetEntries(); j++) {
 635 //      AliVTrack *track = (AliVTrack*)jet->GetTrack(j);
 636 //        if(!track) continue;
 637 //        cout<<track->Pt()<<" "<<zLeading<<endl;
 638 //      if(track->Pt() > zLeading) zLeading = track->Pt();
 639 //        cout<<__LINE__<<endl;
 640 //      }
 641 //      if(jet->Pt() > 0) zLeading /= jet->Pt();
 642 //      else zLeading = -1;
 643 //    }
 644     //cuty
 645     if(ptCut > 0  && ptCut > jet->Pt()                ) continue;
 646     if(etaCut > 0 && etaCut < TMath::Abs(jet->Eta())  ) continue;
 647     if(zCut > 0   && zCut < zLeading  && zLeading > 0 ) continue;
 648
 649     jetList->Add(jet);
 650   }
 651
 652   jetList->Sort();
 653   return jetList->GetEntries();
 654 }
 655
 656 //---------------------------------------------------------------------------------------------------
 657 void AliAnalysisTaskPPJetSpectra::FillJetContainer(TList* jets, THnSparseF* container) {
 658   if(!jets || !container) return;
 659
 660   for(Int_t i =0 ; i < jets->GetEntries(); i++) {
 661     AliAODJet* jet = (AliAODJet*)jets->At(i);
 662     if(!jet) continue;
 663   // pt - eta - phi - AreaCH - AreaN -  ptLeadingTrack/ptJet - ptOrig
 664     Double_t tmpContainer[] = {jet->Pt(), jet->Eta(), jet->Phi(), jet->EffectiveAreaCharged(), jet->EffectiveAreaNeutral(), jet->GetPtLeading()/jet->Pt(), (Double_t)((TRefArray*)jet->GetRefTracks())->GetEntries() };
 665     container->Fill(tmpContainer);
 666   }
 667   return;
 668 }
 669
 670 //---------------------------------------------------------------------------------------------------
 671 Double_t AliAnalysisTaskPPJetSpectra::GetUE(TList* jets, TList* tracks, Double_t Rpar,THnSparseF* thnCont) {
 672   if(!jets || !tracks) {
 673     if(fDebug) printf("%s: %d No jets or tracks\n",(char*)__FILE__,__LINE__);
 674     return 0;
 675   }
 676
 677   //find leading jet
 678   AliAODJet *jet = 0;
 679   for(Int_t iTmp = 0; iTmp < jets->GetEntries(); iTmp++) {
 680     AliAODJet* tmp = (AliAODJet*)jets->At(iTmp);
 681     if(!tmp) continue;
 682     if(!jet) jet = tmp;
 683     else
 684       if( jet->Pt() < tmp->Pt() ) jet = tmp;
 685   }
 686   if(!jet) return 0;
 687
 688   //get track in perpendicular cone
 689   //TVector3 for jet and perpendicular cones with same eta axis
 690   TVector3 j(jet->Px(), jet->Py(), jet->Pz());
 691   TVector3 p1(j);
 692   TVector3 p2(j);
 693
 694   p1.RotateZ(TMath::Pi()/2.);
 695   p2.RotateZ(-TMath::Pi()/2.);
 696
 697   Double_t sumAllPt1 = 0;
 698   Double_t sumAllPt2 = 0;
 699
 700   for(int i = 0; i < tracks->GetEntries(); i++)  {
 701     AliVParticle* tr = (AliVParticle*)tracks->At(i);
 702     if(!tr) continue;
 703     TVector3 v(tr->Px(), tr->Py(), tr->Pz());
 704     Double_t dR1 = v.DrEtaPhi(p1);
 705     Double_t dR2 = v.DrEtaPhi(p2);
 706
 707     if(dR1 < Rpar) sumAllPt1+=v.Pt();
 708     if(dR2 < Rpar) sumAllPt2+=v.Pt();
 709   }
 710
 711   if(Rpar != 0) {
 712     sumAllPt1/=(TMath::Pi()*Rpar*Rpar);
 713     sumAllPt2/=(TMath::Pi()*Rpar*Rpar);
 714   }
 715
 716   Double_t container[] = {j.Pt(), j.Eta(), sumAllPt1, sumAllPt2, (sumAllPt1+sumAllPt2)/2.};
 717   thnCont->Fill(container);
 718
 719   return (sumAllPt1+sumAllPt2)/2.;
 720
 721 }
 722
 723 //---------------------------------------------------------------------------------------------------
 724 Double_t AliAnalysisTaskPPJetSpectra::GetBckgUE(TList* jets, Double_t Rpar, Bool_t isSkipped, THnSparseF* cont) {
 725   if(!jets) return 0;
 726
 727   Int_t leading01 = -1;
 728   Int_t leading02 = -1;
 729   Double_t leading01pt = 0;
 730   Double_t leading02pt = 0;
 731
 732   if(!isSkipped) {
 733     for(Int_t i = 0; i < jets->GetEntries(); i++) {
 734       AliAODJet* tmp = (AliAODJet*)jets->At(i);
 735       if(!tmp) continue;
 736       if(tmp->Pt() > leading01pt) {
 737         leading02 = leading01;
 738         leading01 = i;
 739         leading02pt = leading01pt;
 740         leading01pt = tmp->Pt();
 741       } else if(tmp->Pt() > leading02) {
 742         leading02 = i;
 743         leading02pt = tmp->Pt();
 744       }
 745     }
 746   }
 747
 748   Int_t subLeading01 = -1;
 749   Int_t subLeading02 = -1;
 750   Double_t subLeading01pt = 0;
 751   Double_t subLeading02pt = 0;
 752   Double_t subLeading01area = 0;
 753   Double_t subLeading02area = 0;
 754
 755   Double_t avgAll = 0;
 756   Double_t areaAll = 0;
 757   //search for subleading jets assuming leading 2 (skip02) jets are skipped
 758   if(jets->GetEntries() > 0) {
 759     for(Int_t i = 0; i < jets->GetEntries(); i++) {
 760       if(i == leading01 || i == leading02) continue;
 761       AliAODJet* tmp = (AliAODJet*)jets->At(i);
 762       if(!tmp) continue;
 763       if( subLeading01pt < tmp->Pt() ) {
 764         subLeading02 = subLeading01;
 765         subLeading01 = i;
 766         subLeading02pt = subLeading01pt;
 767         subLeading01pt = tmp->Pt();
 768         subLeading02area = subLeading01area;
 769         subLeading01area = tmp->EffectiveAreaCharged();
 770       } else if( subLeading02pt < tmp->Pt() ) {
 771         subLeading02 = i;
 772         subLeading02pt = tmp->Pt();
 773         subLeading02area = tmp->EffectiveAreaCharged();
 774       }
 775     }
 776
 777     //do avg jet pt
 778     Int_t nJets = 0;
 779     for(int i = 0; i < jets->GetEntries();i++) {
 780       if(i == leading01 || i == leading02) continue;
 781       AliAODJet* tmp = (AliAODJet*)jets->At(i);
 782       if(!tmp) continue;
 783       nJets++;
 784       avgAll+= tmp->Pt();
 785       areaAll += tmp->EffectiveAreaCharged();
 786     }
 787   }
 788
 789
 790   if(avgAll!=0) avgAll/=areaAll;
 791   if(subLeading01pt!=0) subLeading01pt/=subLeading01area;
 792   if(subLeading02pt!=0)subLeading02pt/=subLeading02area;
 793   Double_t tmpContainer[]={subLeading01pt,subLeading02pt,avgAll,Rpar*10};
 794
 795   cont->Fill(tmpContainer);
 796
 797   return avgAll;
 798 }
 799
 800 //---------------------------------------------------------------------------------------------------
 801 Int_t  AliAnalysisTaskPPJetSpectra::CorrectForUE(TList* jets, Double_t UE, TList* newList, THnSparseF* container) {
 802   if(!jets) return 0;
 803
 804   for(Int_t i = 0; i < jets->GetEntries(); i++) {
 805     AliAODJet* tmpJet = (AliAODJet*)jets->At(i);
 806     if(!tmpJet) continue;
 807
 808     AliAODJet* jet = (AliAODJet*)tmpJet->Clone(Form("%s%s",tmpJet->GetName(),"_clone"));
 809     Double_t corrUE = UE*(jet->EffectiveAreaCharged() + jet->EffectiveAreaNeutral() );
 810
 811     if(jet->Pt() > corrUE) jet->SetPtEtaPhiM(jet->Pt()- corrUE, jet->Eta(), jet->Phi(), jet->M());
 812     else jet->SetPtEtaPhiM(0, jet->Eta(), jet->Phi(), jet->M());
 813     newList->Add(jet);
 814
 815   // pt - eta - phi - AreaCH - AreaN - UEdensity - ptLeadingTrack/ptJet - ptOrig
 816     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()};
 817     container->Fill(cont);
 818   }
 819   newList->Sort();
 820   return newList->GetEntries();
 821 }
 822
 823 //---------------------------------------------------------------------------------------------------
 824 void AliAnalysisTaskPPJetSpectra::MatchJets(Bool_t doBothWay, TList* recList, TList* genList, Float_t maxDist, THnSparseF* container) {
 825   if(!genList || !recList) return;
 826   Int_t mode = 1;
 827
 828   const Int_t kGenJets = genList->GetEntries();
 829   const Int_t kRecJets = recList->GetEntries();
 830
 831   TArrayI iMatchIndex(kGenJets);
 832   TArrayF fPtFraction(kGenJets);
 833
 834   TArrayI aGenIndex(kRecJets);
 835   TArrayI aRecIndex(kGenJets);
 836
 837   Double_t genPt = -1;
 838   Double_t recPt = -1;
 839   Double_t ptFraction = -1;
 840
 841   //-- Closest jets to generated and checked for E-fraction
 842   if(!doBothWay) {
 843     AliAnalysisHelperJetTasks::GetJetMatching(genList, kGenJets, recList, kRecJets, iMatchIndex, fPtFraction, fDebug, maxDist, mode);
 844     for(Int_t ig = 0; ig < kGenJets; ig++) {
 845       //get gen-jet
 846       AliAODJet* genJet = (AliAODJet*)genList->At(ig);
 847       if(!genJet) continue;
 848       genPt = genJet->Pt();
 849
 850       //get rec-jet
 851       Int_t ir = iMatchIndex[ig];
 852       if(ir<0||ir>=recList->GetEntries()) continue;
 853       AliAODJet* recJet = (AliAODJet*)recList->At(ir);
 854       if(!recJet) continue;
 855       recPt = recJet->Pt();
 856
 857       ptFraction = fPtFraction[ig];
 858
 859       Double_t dR = recJet->DeltaR(genJet);
 860       Double_t dEta = recJet->Eta() - genJet->Eta();
 861       Double_t dPhi2 = dR*dR - dEta*dEta;
 862
 863         Bool_t isSame = kFALSE;
 864         if(CheckPtBin(genPt) == CheckPtBin(recPt)) isSame = kTRUE;
 865       // ptGen - ptRec - ptFraction - maxDist - bothWays - (1 - rec/gen) - dR - dEta  - dPhi
 866       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};
 867       container->Fill(cont);
 868     }
 869   }
 870   //-- Closest rec jet to gen and vice-versa & check for E-fraction
 871   else {
 872     AliAnalysisHelperJetTasks::GetClosestJets(genList, kGenJets, recList, kRecJets, aGenIndex, aRecIndex,fDebug, maxDist);
 873     for(Int_t ig = 0; ig < kGenJets; ig++) {
 874       //get gen-jet
 875       AliAODJet* genJet = (AliAODJet*)genList->At(ig);
 876       if(!genJet) continue;
 877       genPt = genJet->Pt();
 878
 879       //get rec-jet
 880       Int_t ir = aRecIndex[ig];
 881       if(ir>=0&&ir<recList->GetEntries()) {
 882         AliAODJet *recJet = (AliAODJet*)recList->At(ir);
 883         if(!recJet) continue;
 884         recPt = recJet->Pt();
 885         Double_t dR = recJet->DeltaR(genJet);
 886         Double_t dEta = recJet->Eta() - genJet->Eta();
 887         Double_t dPhi2 = dR*dR - dEta*dEta;
 888
 889         ptFraction = AliAnalysisHelperJetTasks::GetFractionOfJet(recJet,genJet,mode);
 890
 891         Bool_t isSame = kFALSE;
 892         if(CheckPtBin(genPt) == CheckPtBin(recPt)) isSame = kTRUE;
 893         // ptGen - ptRec - ptFraction - maxDist - bothWays - (1 - rec/gen) - dR - dEta  - dPhi
 894         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  };
 895         container->Fill(cont);
 896       }
 897     }
 898   }
 899
 900   return;
 901 }
 902
 903 //---------------------------------------------------------------------------------------------------
 904 Int_t AliAnalysisTaskPPJetSpectra::CheckPtBin(Double_t pT) {
 905   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.};
 906   Double_t nBins        = 25;
 907
 908   for(Int_t i = 0; i <= nBins; i++) {
 909     if(pT < bins[i] ) return i;
 910   }
 911   return nBins+1;
 912 }
 913
 914 //---------------------------------------------------------------------------------------------------
 915 void AliAnalysisTaskPPJetSpectra::DoUEAnalysis(TList* tracks, Double_t PTcut, Double_t ETAcut) {
 916         Double_t PTmax = 0;
 917         int iMax = -1;
 918         for(int i = 0; i < tracks->GetEntries(); i++) {
 919                 AliVParticle *part = (AliVParticle*)tracks->At(i);
 920                 if(!part) continue;
 921
 922                 if(part->Pt() < PTcut) continue;
 923                 if(TMath::Abs(part->Eta()) > ETAcut ) continue;
 924
 925                 if(PTmax < part->Pt()) {
 926                         PTmax = part->Pt();
 927                         iMax = i;
 928                 }
 929         }
 930
 931         if(PTmax < trackPtMin) return;
 932
 933         Double_t ch_part_density = 0;
 934         Double_t ch_part_pt_density = 0;
 935         Double_t avg_track_pt = 0;
 936         Double_t zero_ch_part = 0;
 937         Double_t transMAX_part_density = 0;
 938         Double_t transMIN_part_density = 0;
 939         Double_t transMAX_pt_density = 0;
 940         Double_t transMIN_pt_density = 0;
 941         Double_t transDIF_pt_density = 0;
 942         Double_t transDIF_part_density = 0;
 943
 944         Double_t A1_part_density = 0;
 945         Double_t A2_part_density = 0;
 946         Double_t A1_pt_density = 0;
 947         Double_t A2_pt_density = 0;
 948
 949         AliVParticle *maxParticle = (AliVParticle*)tracks->At(iMax);
 950
 951         Double_t PHImax = maxParticle->Phi();
 952
 953         for(int i = 0; i < tracks->GetEntries(); i++) {
 954                 AliVParticle* part = (AliVParticle*)tracks->At(i);
 955                 if(!part) continue;
 956                 if(part->Pt() < PTcut) continue;
 957
 958                 Bool_t isArea1 = kFALSE;
 959                 Bool_t isArea2 = kFALSE;
 960
 961                 Double_t dPhi = part->Phi() - PHImax;
 962                 if(dPhi > TMath::Pi()) dPhi -= TMath::Pi();
 963                 if(dPhi < -TMath::Pi()) dPhi +=TMath::Pi();
 964
 965                 if(TMath::Abs(dPhi) > TMath::Pi()/3 && TMath::Abs(dPhi) < 2*TMath::Pi()/3 ) {
 966                         if(dPhi>0) isArea1 = kTRUE;
 967                         else isArea2 = kTRUE;
 968                 }
 969
 970                 if(!isArea1 && !isArea2) continue;
 971
 972                 ch_part_density += 1;
 973                 ch_part_pt_density += part->Pt();
 974                 avg_track_pt += part->Pt();
 975
 976                 if(isArea1) {
 977                         A1_part_density += 1;
 978                         A1_pt_density += part->Pt();
 979                 } else if(isArea2) {
 980                         A2_part_density += 1;
 981                         A2_pt_density += part->Pt();
 982                 }
 983         }
 984
 985
 986         if(ch_part_density > 0) avg_track_pt/=ch_part_density;
 987         else zero_ch_part = 1;
 988         ch_part_density*=3./(4.*TMath::Pi()*ETAcut);
 989         ch_part_pt_density*= 3./(4.*TMath::Pi()*ETAcut);
 990
 991         A1_part_density*=3./(4.*TMath::Pi()*ETAcut);
 992         A1_pt_density*= 3./(4.*TMath::Pi()*ETAcut);
 993         A2_part_density*=3./(4.*TMath::Pi()*ETAcut);
 994         A2_pt_density*= 3./(4.*TMath::Pi()*ETAcut);
 995
 996         if(A1_part_density < A2_part_density) {
 997                 transMAX_part_density = A2_part_density;
 998                 transMIN_part_density = A1_part_density;
 999         } else {
1000                 transMAX_part_density = A1_part_density;
1001                 transMIN_part_density = A2_part_density;
1002         }
1003
1004         if(A1_pt_density < A2_pt_density) {
1005                 transMAX_pt_density = A2_pt_density;
1006                 transMIN_pt_density = A1_pt_density;
1007         } else {
1008                 transMAX_pt_density = A1_pt_density;
1009                 transMIN_pt_density = A2_pt_density;
1010         }
1011
1012         transDIF_pt_density     = transMAX_pt_density - transMIN_pt_density;
1013         transDIF_part_density   = transMAX_part_density - transMIN_part_density;
1014
1015         Double_t container[] = {PTmax, ch_part_density, ch_part_pt_density, avg_track_pt, zero_ch_part,
1016                                                         transMAX_part_density, transMIN_part_density, transDIF_part_density,
1017                                                         transMAX_pt_density, transMIN_pt_density, transDIF_pt_density};
1018
1019         fhnTrackUEanal->Fill(container);
1020 }