PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx

   1 #ifndef ALIANALYSISTASKSE_H
   2 #include <Riostream.h>
   3 #include <TROOT.h>
   4 #include <TFile.h>
   5 #include <TCint.h>
   6 #include <TChain.h>
   7 #include <TTree.h>
   8 #include <TKey.h>
   9 #include <TProfile.h>
  10 #include <TH1F.h>
  11 #include <TH2F.h>
  12 #include <TCanvas.h>
  13 #include <TList.h>
  14 #include <TClonesArray.h>
  15 #include <TObject.h>
  16 #include <TMath.h>
  17 #include <TSystem.h>
  18 #include <TInterpreter.h>
  19 #include <TH1.h>
  20 #include "AliAnalysisTask.h"
  21 #include "AliCentrality.h"
  22 #include "AliStack.h"
  23 #include "AliESDEvent.h"
  24 #include "AliESDInputHandler.h"
  25 #include "AliAODEvent.h"
  26 #include "AliAODHandler.h"
  27 #include "AliAnalysisManager.h"
  28 #include "AliAnalysisTaskSE.h"
  29 #endif
  30
  31 #include <TRandom3.h>
  32 #include "AliGenPythiaEventHeader.h"
  33 #include "AliMCEvent.h"
  34 #include "AliLog.h"
  35 #include <AliEmcalJet.h>
  36 #include "AliVEventHandler.h"
  37 #include "AliVParticle.h"
  38 #include "AliAnalysisUtils.h"
  39
  40
  41 #include "AliAnalysisTaskChargedJetsPA.h"
  42
  43 //TODO: Not accessing the particles when using MC
  44 //TODO: FillHistogram can be done better with virtual TH1(?)
  45 ClassImp(AliAnalysisTaskChargedJetsPA)
  46
  47 // ######################################################################################## DEFINE HISTOGRAMS
  48 void AliAnalysisTaskChargedJetsPA::Init()
  49 {
  50   #ifdef DEBUGMODE
  51     AliInfo("Creating histograms.");
  52   #endif
  53
  54   AddHistogram1D<TH1D>("hNumberEvents", "Number of events (0 = before, 1 = after vertex cuts)", "", 2, 0, 2, "#Delta z(cm)","N^{Events}/cut");
  55
  56   // NOTE: Jet histograms
  57   if (fAnalyzeJets)
  58   {
  59     // ######## Jet spectra
  60     AddHistogram2D<TH2D>("hJetPt", "Jets p_{T} distribution", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  61     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedRC", "Jets p_{T} distribution, RC background subtracted", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  62     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKT", "Jets p_{T} distribution, KT background subtracted", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  63     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedTR", "Jets p_{T} distribution, TR background subtracted", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  64     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedRCPosEta", "Jets p_{T} distribution, RC background subtracted, #eta > 0.2", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  65     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTPosEta", "Jets p_{T} distribution, KT background subtracted, #eta > 0.2", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  66     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedTRPosEta", "Jets p_{T} distribution, TR background subtracted, #eta > 0.2", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  67     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedRCNegEta", "Jets p_{T} distribution, RC background subtracted, #eta < -0.2", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  68     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTNegEta", "Jets p_{T} distribution, KT background subtracted, #eta < -0.2", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  69     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedTRNegEta", "Jets p_{T} distribution, TR background subtracted, #eta < -0.2", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  70
  71     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedRCEtaWindow", "Jets p_{T} distribution, RC background (in #eta window around jet) subtracted", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  72     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTEtaWindow", "Jets p_{T} distribution, KT background (in #eta window around jet) subtracted", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  73     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedTREtaWindow", "Jets p_{T} distribution, TR background (in #eta window around jet) subtracted", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  74
  75     // ######## Jet stuff
  76     AddHistogram1D<TH1D>("hJetCountAll", "Number of Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
  77     AddHistogram1D<TH1D>("hJetCountAccepted", "Number of accepted Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
  78     AddHistogram1D<TH1D>("hLeadingJetPt", "Leading jet p_{T}", "", 500,  0, 100, "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
  79     AddHistogram1D<TH1D>("hSecondLeadingJetPt", "Second Leading jet p_{T}", "", 500,  0, 100, "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
  80     AddHistogram1D<TH1D>("hJetDeltaPhi", "Jets combinatorial #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
  81     AddHistogram1D<TH1D>("hLeadingJetDeltaPhi", "1st and 2nd leading jet #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
  82   }
  83
  84   // NOTE: Jet background histograms
  85   if (fAnalyzeBackground)
  86   {
  87     // ########## Different background estimates
  88     AddHistogram2D<TH2D>("hRCBackground", "RC background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
  89     AddHistogram2D<TH2D>("hKTBackground", "KT background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
  90     AddHistogram2D<TH2D>("hTRBackground", "TR background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
  91
  92     // ########## Delta Pt
  93     AddHistogram2D<TH2D>("hDeltaPtKT", "Background fluctuations #delta p_{T} (KT)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  94     AddHistogram2D<TH2D>("hDeltaPtRC", "Background fluctuations #delta p_{T} (RC)", "", 600, -40., 80., 5, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  95     AddHistogram2D<TH2D>("hDeltaPtTR", "Background fluctuations #delta p_{T} (TR)", "", 600, -40., 80., 5, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  96     AddHistogram2D<TH2D>("hDeltaPtKTNoExcl", "Background fluctuations #delta p_{T} (KT, no leading jet correction)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  97     AddHistogram2D<TH2D>("hDeltaPtRCNoExcl", "Background fluctuations #delta p_{T} (RC, no leading jet correction)", "", 600, -40., 80., 5, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  98     AddHistogram2D<TH2D>("hDeltaPtTRNoExcl", "Background fluctuations #delta p_{T} (TR, no leading jet correction)", "", 600, -40., 80., 5, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  99
 100     // ########## Min bias background in eta bins
 101     AddHistogram2D<TH2D>("hRCBackgroundEtaBins", "RC background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta");
 102     AddHistogram2D<TH2D>("hKTBackgroundEtaBins", "KT background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta");
 103     AddHistogram2D<TH2D>("hTRBackgroundEtaBins", "TR background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta");
 104
 105     // ########## Min bias background in sliding eta bins
 106     AddHistogram2D<TH2D>("hRCBackgroundEtaWindow", "RC background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho d#eta");
 107     AddHistogram2D<TH2D>("hKTBackgroundEtaWindow", "KT background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho d#eta");
 108     AddHistogram2D<TH2D>("hTRBackgroundEtaWindow", "TR background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho d#eta");
 109
 110
 111     // ########## Dijet stuff
 112     AddHistogram1D<TH1D>("hDijetLeadingJetPt", "Dijet leading jet p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 113     AddHistogram1D<TH1D>("hDijetConstituentsPt", "Dijet constituents p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 114     AddHistogram2D<TH2D>("hDijetPtCorrelation", "Dijet constituents p_{T} correlation", "COLZ", 500, 5., 100., 500, 5., 100., "1st leading jet p_{T} (GeV/c)","2nd leading jet p_{T} (GeV/c)","dN^{Dijets}/d^{2}p_{T}");
 115     AddHistogram2D<TH2D>("hDijetBackground", "Background density (dijets excluded)", "", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 116     AddHistogram2D<TH2D>("hDijetBackgroundPerpendicular", "Background density (dijets excluded)", "", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 117   }
 118
 119   // NOTE: Pythia histograms
 120   if (fAnalyzePythia)
 121   {
 122     AddHistogram1D<TH1D>("hPythiaPtHard", "Pythia p_{T} hard distribution", "", 2000, 0, 400, "p_{T} hard","dN^{Events}/dp_{T,hard}");
 123     AddHistogram1D<TProfile>("hPythiaXSection", "Pythia cross section distribution", "", fNumPtHardBins+2, -1, fNumPtHardBins+1, "p_{T} hard bin","dN^{Events}/dp_{T,hard}");
 124     AddHistogram1D<TH1D>("hPythiaNTrials", "Pythia trials (no correction for manual cuts)", "", fNumPtHardBins+2, -1, fNumPtHardBins+1, "p_{T} hard bin", "Trials");
 125   }
 126
 127   // register Histograms
 128   for (Int_t i = 0; i < fHistCount; i++)
 129   {
 130     fOutputList->Add(fHistList->At(i));
 131   }
 132
 133   PostData(1,fOutputList); // important for merging
 134
 135 }
 136
 137 //________________________________________________________________________
 138 AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName) : AliAnalysisTaskSE(name), fOutputList(0), fAnalyzeJets(1), fAnalyzeBackground(1), fAnalyzePythia(0), fHasTracks(0), fHasJets(0), fHasBackgroundJets(0), fIsMC(0), fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fNumPtHardBins(11), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fTRBackgroundConeRadius(0.4), fNumberRandCones(8), fNumberExcludedJets(2), fDijetMaxAngleDeviation(10.0), fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fVertexWindow(10.0), fVertexMaxR(1.0), fMinTrackPt(0.150), fMinJetPt(1.0), fMinJetArea(0.4), fMinBackgroundJetPt(0.15), fMinDijetLeadingPt(10.0), fCentralityType("V0A"), fFirstLeadingJet(0), fSecondLeadingJet(0), fNumberSignalJets(0), fCrossSection(0.0), fTrials(0.0),  fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fHistList(0), fHistCount(0)
 139 {
 140   #ifdef DEBUGMODE
 141     AliInfo("Calling constructor.");
 142   #endif
 143
 144   // Every instance of this task gets his own number
 145   static Int_t instance = 0;
 146   fTaskInstanceCounter = instance;
 147   instance++;
 148
 149   fTrackArrayName = new TString(trackArrayName);
 150   if (fTrackArrayName->Contains("MCParticles")) //TODO: Not working for now
 151     fIsMC = kTRUE;
 152
 153   fJetArrayName = new TString(jetArrayName);
 154   if (strcmp(fJetArrayName->Data(),"") == 0)
 155     fAnalyzeJets = kFALSE;
 156   else
 157     fAnalyzeJets = kTRUE;
 158
 159   fBackgroundJetArrayName = new TString(backgroundJetArrayName);
 160   if (strcmp(fBackgroundJetArrayName->Data(),"") == 0)
 161     fAnalyzeBackground = kFALSE;
 162   else
 163     fAnalyzeBackground = kTRUE;
 164
 165   DefineOutput(1, TList::Class());
 166
 167   fHistList = new TList();
 168
 169   #ifdef DEBUGMODE
 170     AliInfo("Constructor done.");
 171   #endif
 172
 173 }
 174
 175 //________________________________________________________________________
 176 inline Double_t AliAnalysisTaskChargedJetsPA::GetConePt(Double_t eta, Double_t phi, Double_t radius)
 177 {
 178   Double_t tmpConePt = 0.0;
 179
 180   for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
 181   {
 182     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
 183     if (IsTrackInAcceptance(tmpTrack))
 184       if(IsTrackInCone(tmpTrack, eta, phi, radius))
 185         tmpConePt = tmpConePt + tmpTrack->Pt();
 186   }
 187   return tmpConePt;
 188 }
 189
 190
 191 //________________________________________________________________________
 192 inline Double_t AliAnalysisTaskChargedJetsPA::GetPtHard()
 193 {
 194   Double_t tmpPtHard = -1.0;
 195
 196   if (!MCEvent())
 197     AliError("MCEvent not accessible although demanded!");
 198   else
 199   {
 200     AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 201     if (!pythiaHeader)
 202       AliError("Pythia Header not accessible!");
 203     else
 204       tmpPtHard = pythiaHeader->GetPtHard();
 205   }
 206   return tmpPtHard;
 207 }
 208
 209
 210 //________________________________________________________________________
 211 inline Int_t AliAnalysisTaskChargedJetsPA::GetPtHardBin()
 212 {
 213   // ########## PT HARD BIN EDGES
 214   const Int_t kPtHardLowerEdges[] =  { 0, 5,11,21,36,57, 84,117,152,191,234};
 215   const Int_t kPtHardHigherEdges[] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 216
 217   Int_t tmpPtHardBin = 0;
 218   Double_t tmpPtHard = GetPtHard();
 219
 220   for (tmpPtHardBin = 0; tmpPtHardBin <= fNumPtHardBins; tmpPtHardBin++)
 221     if (tmpPtHard >= kPtHardLowerEdges[tmpPtHardBin] && tmpPtHard < kPtHardHigherEdges[tmpPtHardBin])
 222       break;
 223
 224   return tmpPtHardBin;
 225 }
 226
 227
 228 //________________________________________________________________________
 229 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInCone(AliVTrack* track, Double_t eta, Double_t phi, Double_t radius)
 230 {
 231   // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
 232   Double_t trackPhi = 0.0;
 233   if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
 234     trackPhi = track->Phi() - TMath::TwoPi();
 235   else if (track->Phi() < (phi+radius - TMath::TwoPi()))
 236     trackPhi = track->Phi() + TMath::TwoPi();
 237   else
 238     trackPhi = track->Phi();
 239
 240   if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
 241     return kTRUE;
 242
 243   return kFALSE;
 244 }
 245
 246 //________________________________________________________________________
 247 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInAcceptance(AliVParticle* track)
 248 {
 249   if (track != 0)
 250     if (TMath::Abs(track->Eta()) <= fTrackEtaWindow)
 251       if (track->Pt() >= fMinTrackPt)
 252         return kTRUE;
 253
 254   return kFALSE;
 255 }
 256
 257 //________________________________________________________________________
 258 inline Bool_t AliAnalysisTaskChargedJetsPA::IsBackgroundJetInAcceptance(AliEmcalJet *jet)
 259 {
 260   if (jet != 0)
 261     if (TMath::Abs(jet->Eta()) <= fBackgroundJetEtaWindow)
 262       if (jet->Pt() >= fMinBackgroundJetPt)
 263         return kTRUE;
 264
 265   return kFALSE;
 266 }
 267
 268 //________________________________________________________________________
 269 inline Bool_t AliAnalysisTaskChargedJetsPA::IsSignalJetInAcceptance(AliEmcalJet *jet)
 270 {
 271   if (jet != 0)
 272     if (TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow)
 273       if (jet->Pt() >= fMinJetPt)
 274         if (jet->Area() >= fMinJetArea)
 275           return kTRUE;
 276
 277   return kFALSE;
 278 }
 279
 280 //________________________________________________________________________
 281 inline Bool_t AliAnalysisTaskChargedJetsPA::IsDijet(AliEmcalJet *jet1, AliEmcalJet *jet2)
 282 {
 283   // Output from GetDeltaPhi is < pi in any case
 284   if ((jet1 != 0) && (jet2 != 0))
 285     if((TMath::Pi() - GetDeltaPhi(jet1->Phi(),jet2->Phi())) < fDijetMaxAngleDeviation)
 286       if ((jet1->Pt() > fMinDijetLeadingPt) && (jet2->Pt() > fMinDijetLeadingPt)) //TODO: Introduce recoil jet?
 287         return kTRUE;
 288
 289   return kFALSE;
 290 }
 291
 292 //________________________________________________________________________
 293 void AliAnalysisTaskChargedJetsPA::ExecOnce()
 294 {
 295   #ifdef DEBUGMODE
 296     AliInfo("Starting ExecOnce.");
 297   #endif
 298   fInitialized = kTRUE;
 299
 300   // Check for track array
 301   if (strcmp(fTrackArrayName->Data(), "") != 0)
 302   {
 303     fTrackArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrayName->Data()));
 304     fHasTracks = kTRUE;
 305     if (!fTrackArray)
 306     {
 307       AliWarning(Form("%s: Could not retrieve tracks %s! This is OK, if tracks are not demanded.", GetName(), fTrackArrayName->Data()));
 308       fHasTracks = kFALSE;
 309     }
 310     else
 311     {
 312       TClass *cl = fTrackArray->GetClass();
 313       if (!cl->GetBaseClass("AliVParticle"))
 314       {
 315         AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName->Data()));
 316         fTrackArray = 0;
 317         fHasTracks = kFALSE;
 318       }
 319     }
 320   }
 321
 322   // Check for jet array
 323   if (strcmp(fJetArrayName->Data(), "") != 0)
 324   {
 325     fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrayName->Data()));
 326     fHasJets = kTRUE;
 327
 328     if (!fJetArray)
 329     {
 330       AliWarning(Form("%s: Could not retrieve jets %s! This is OK, if jets are not demanded.", GetName(), fJetArrayName->Data()));
 331       fHasJets = kFALSE;
 332     }
 333     else
 334     {
 335       if (!fJetArray->GetClass()->GetBaseClass("AliEmcalJet"))
 336       {
 337         AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName->Data()));
 338         fJetArray = 0;
 339         fHasJets = kFALSE;
 340       }
 341     }
 342   }
 343
 344   // Check for background object
 345   if (strcmp(fBackgroundJetArrayName->Data(), "") != 0)
 346   {
 347     fBackgroundJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fBackgroundJetArrayName->Data()));
 348     fHasBackgroundJets = kTRUE;
 349     if (!fBackgroundJetArray)
 350     {
 351       AliInfo(Form("%s: Could not retrieve background jets %s! This is OK, if background is not demanded.", GetName(), fBackgroundJetArrayName->Data()));
 352       fHasBackgroundJets = kFALSE;
 353     }
 354   }
 355
 356   // Look, if initialization is OK
 357   if (!fHasTracks && fAnalyzeBackground)
 358   {
 359     AliError(Form("%s: Tracks NOT successfully casted although demanded! Deactivating background analysis",GetName()));
 360     fAnalyzeBackground = kFALSE;
 361   }
 362   if ((!fHasJets && fAnalyzeJets) || (!fHasJets && fAnalyzeBackground))
 363   {
 364     AliError(Form("%s: Jets NOT successfully casted although demanded!  Deactivating jet- and background analysis",GetName()));
 365     fAnalyzeJets = kFALSE;
 366     fAnalyzeBackground = kFALSE;
 367   }
 368   if (!fHasBackgroundJets && fAnalyzeBackground)
 369   {
 370     AliError(Form("%s: Background NOT successfully casted although demanded!  Deactivating background analysis",GetName()));
 371     fAnalyzeBackground = kFALSE;
 372   }
 373
 374   // Initialize helper class (for vertex selection)
 375   fHelperClass = new AliAnalysisUtils();
 376
 377   // Histogram init
 378   Init();
 379
 380   #ifdef DEBUGMODE
 381     AliInfo("ExecOnce done.");
 382   #endif
 383
 384 }
 385
 386 //________________________________________________________________________
 387 void AliAnalysisTaskChargedJetsPA::GetSignalJets()
 388 {
 389   // Reset vars
 390   fFirstLeadingJet = NULL;
 391   fSecondLeadingJet = NULL;
 392   fNumberSignalJets = 0;
 393
 394   Float_t maxJetPts[] = {0, 0};
 395   Int_t jetIDArray[]   = {-1, -1};
 396   Int_t jetCount = fJetArray->GetEntries();
 397
 398   // Go through all jets and save signal jets and the two leading ones
 399   for (Int_t i = 0; i < jetCount; i++)
 400   {
 401     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(i));
 402     if (!jet)
 403     {
 404       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 405       continue;
 406     }
 407
 408     if (!IsSignalJetInAcceptance(jet)) continue;
 409
 410     if (jet->Pt() > maxJetPts[0])
 411     {
 412       maxJetPts[1] = maxJetPts[0];
 413       jetIDArray[1] = jetIDArray[0];
 414       maxJetPts[0] = jet->Pt();
 415       jetIDArray[0] = i;
 416     }
 417     else if (jet->Pt() > maxJetPts[1])
 418     {
 419       maxJetPts[1] = jet->Pt();
 420       jetIDArray[1] = i;
 421     }
 422     fSignalJets[fNumberSignalJets] = jet;
 423     fNumberSignalJets++;
 424   }
 425
 426   if (fNumberSignalJets > 0)
 427     fFirstLeadingJet  = static_cast<AliEmcalJet*>(fJetArray->At(jetIDArray[0]));
 428   if (fNumberSignalJets > 1)
 429     fSecondLeadingJet = static_cast<AliEmcalJet*>(fJetArray->At(jetIDArray[1]));
 430
 431 }
 432
 433 //________________________________________________________________________
 434 Int_t AliAnalysisTaskChargedJetsPA::GetLeadingJets(TClonesArray* jetArray, Int_t* jetIDArray, Bool_t isSignalJets)
 435 {
 436 // Writes first two leading jets into already registered array jetIDArray
 437
 438   if (!jetArray)
 439   {
 440     AliError("Could not get the jet array to get leading jets from it!");
 441     return 0;
 442   }
 443
 444   Float_t maxJetPts[] = {0, 0};
 445   jetIDArray[0] = -1;
 446   jetIDArray[1] = -1;
 447
 448   Int_t jetCount = jetArray->GetEntries();
 449   Int_t jetCountAccepted = 0;
 450
 451   for (Int_t i = 0; i < jetCount; i++)
 452   {
 453     AliEmcalJet* jet = static_cast<AliEmcalJet*>(jetArray->At(i));
 454     if (!jet)
 455     {
 456       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 457       continue;
 458     }
 459
 460     if(isSignalJets)
 461     {
 462       if (!IsSignalJetInAcceptance(jet)) continue;
 463     }
 464     else
 465     {
 466       if (!IsBackgroundJetInAcceptance(jet)) continue;
 467     }
 468
 469     if (jet->Pt() > maxJetPts[0])
 470     {
 471       maxJetPts[1] = maxJetPts[0];
 472       jetIDArray[1] = jetIDArray[0];
 473       maxJetPts[0] = jet->Pt();
 474       jetIDArray[0] = i;
 475     }
 476     else if (jet->Pt() > maxJetPts[1])
 477     {
 478       maxJetPts[1] = jet->Pt();
 479       jetIDArray[1] = i;
 480     }
 481     jetCountAccepted++;
 482   }
 483   return jetCountAccepted;
 484 }
 485
 486
 487 //________________________________________________________________________
 488 Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Double_t background)
 489 {
 490   #ifdef DEBUGMODE
 491     AliInfo("Getting corrected jet spectra.");
 492   #endif
 493
 494   if(!jet)
 495   {
 496     AliError("Jet pointer passed to GetCorrectedJet() not valid!");
 497     return -1.0;
 498   }
 499
 500   Double_t correctedPt = -1.0;
 501
 502   // if the passed background is not valid, do not subtract it
 503   if(background < 0)
 504     background = 0;
 505
 506   // Subtract background
 507   correctedPt = jet->Pt() - background * jet->Area();
 508
 509   #ifdef DEBUGMODE
 510     AliInfo("Got corrected jet spectra.");
 511   #endif
 512
 513   return correctedPt;
 514 }
 515
 516
 517
 518 //________________________________________________________________________
 519 void AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t& deltaPt, Double_t rho, Bool_t leadingJetExclusion)
 520 {
 521   #ifdef DEBUGMODE
 522     AliInfo("Getting Delta Pt.");
 523   #endif
 524
 525   // Define an invalid delta pt
 526   deltaPt = -10000.0;
 527
 528   // Define eta range
 529   Double_t etaMin, etaMax;
 530   etaMin = -(fTrackEtaWindow-fRandConeRadius);
 531   etaMax = +(fTrackEtaWindow-fRandConeRadius);
 532
 533   // Define random cone
 534   Bool_t coneValid = kTRUE;
 535   Double_t tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
 536   Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
 537
 538   // if there is a jet, check for overlap if demanded
 539   if(leadingJetExclusion)
 540   {
 541     for (Int_t i = 0; i<fNumberSignalJets; i++)
 542     {
 543       AliEmcalJet* tmpJet = fSignalJets[i];
 544
 545       Double_t excludedJetPhi = tmpJet->Phi();
 546       Double_t excludedJetEta = tmpJet->Eta();
 547       Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, excludedJetPhi);
 548
 549       // Check, if cone has overlap with jet
 550       if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(tmpRandConeEta-excludedJetEta)*TMath::Abs(tmpRandConeEta-excludedJetEta) <= fRandConeRadius*fRandConeRadius)
 551       {
 552         // Define probability to exclude the RC
 553         Double_t probability = 1 - (fNumberSignalJets-1)/fNumberSignalJets;
 554
 555         // Only exclude cone with a given probability
 556         if (fRandom->Rndm()<=probability)
 557         {
 558           coneValid = kFALSE;
 559           break;
 560         }
 561       }
 562     }
 563   }
 564
 565   // Get the cones' pt and calculate delta pt
 566   if (coneValid)
 567     deltaPt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius) - (rho*fRandConeRadius*fRandConeRadius*TMath::Pi());
 568
 569   #ifdef DEBUGMODE
 570     AliInfo("Got Delta Pt.");
 571   #endif
 572 }
 573
 574 //________________________________________________________________________
 575 void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMedian, Double_t& areaMean, Double_t etaMin, Double_t etaMax, Bool_t excludeSignalJets)
 576 {
 577   #ifdef DEBUGMODE
 578     AliInfo("Getting KT background density.");
 579   #endif
 580
 581   // static declaration. Advantage: more speed. Disadvantage: Problematic for events with more than 1024 jets :)
 582   static Double_t tmpRhos[1024];
 583   static Double_t tmpAreas[1024];
 584   Int_t maxJetIds[]   = {-1, -1}; // Indices for excludes jets (up to two)
 585
 586   // Setting invalid values
 587   rhoMedian = -1.0;
 588   areaMean= -1.0;
 589
 590   GetLeadingJets(fBackgroundJetArray, &maxJetIds[0], kFALSE);
 591   // Exclude UP TO numberExcludeLeadingJets
 592   if (fNumberSignalJets < 2)
 593     numberExcludeLeadingJets = fNumberSignalJets;
 594
 595   // Check if etaMin/etaMax is given correctly
 596   if(etaMin < -fBackgroundJetEtaWindow)
 597     etaMin = -fBackgroundJetEtaWindow;
 598   if(etaMax > fBackgroundJetEtaWindow)
 599     etaMax = fBackgroundJetEtaWindow;
 600
 601   if ((etaMin == 0) && (etaMax == 0))
 602   {
 603     etaMin = -fBackgroundJetEtaWindow;
 604     etaMax = +fBackgroundJetEtaWindow;
 605   }
 606
 607   Int_t jetCountAccepted = 0;
 608   for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++)
 609   {
 610     AliEmcalJet* backgroundJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
 611     if (!backgroundJet)
 612     {
 613       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 614       continue;
 615     }
 616
 617     // Exclude signal jets
 618     if (excludeSignalJets)
 619     {
 620       Bool_t isOverlapping = kFALSE;
 621       for(Int_t j=0;j<numberExcludeLeadingJets;j++)
 622       {
 623         AliEmcalJet* signalJet = NULL;
 624
 625         if (j==0)
 626           signalJet = fFirstLeadingJet;
 627         else if (j==1)
 628           signalJet = fSecondLeadingJet;
 629         else
 630           AliFatal("Trying to exclude more than 2 signal jets in KT background -- not implemented.");
 631
 632
 633         Double_t tmpDeltaPhi = GetDeltaPhi(backgroundJet->Phi(), signalJet->Phi());
 634
 635         if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(signalJet->Eta()-backgroundJet->Eta())*TMath::Abs(signalJet->Eta()-backgroundJet->Eta()) <= fSignalJetRadius*fSignalJetRadius)
 636         {
 637           isOverlapping = kTRUE;
 638           break;
 639         }
 640       }
 641       if(isOverlapping)
 642         continue;
 643     }
 644     else  // exclude leading background jets
 645     {
 646       if (numberExcludeLeadingJets > 0)
 647         if (i == maxJetIds[0])
 648           continue;
 649       if (numberExcludeLeadingJets > 1)
 650         if (i == maxJetIds[1])
 651           continue;
 652     }
 653
 654     if (!IsBackgroundJetInAcceptance(backgroundJet))
 655       continue;
 656     if (!((backgroundJet->Eta() >= etaMin) && (backgroundJet->Eta() < etaMax)))
 657       continue;
 658
 659
 660     tmpRhos[jetCountAccepted] = backgroundJet->Pt() / backgroundJet->Area();
 661     tmpAreas[jetCountAccepted] = backgroundJet->Area();
 662     jetCountAccepted++;
 663   }
 664
 665   if (jetCountAccepted > 0)
 666   {
 667     rhoMedian = TMath::Median(jetCountAccepted, tmpRhos);
 668     areaMean   = TMath::Mean(jetCountAccepted, tmpAreas);
 669   }
 670   #ifdef DEBUGMODE
 671     AliInfo("Got KT background density.");
 672   #endif
 673 }
 674
 675
 676 //________________________________________________________________________
 677 Int_t AliAnalysisTaskChargedJetsPA::GetRCBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& rhoMedian, Double_t etaMin, Double_t etaMax, Int_t numberRandCones)
 678 {
 679   #ifdef DEBUGMODE
 680     AliInfo("Getting RC background density.");
 681   #endif
 682
 683   if(numberRandCones == 0)
 684     numberRandCones = fNumberRandCones;
 685
 686   std::vector<AliEmcalJet> tmpCones(numberRandCones);
 687
 688   // Setting invalid values
 689   rhoMean = -1.0;
 690   rhoMedian = -1.0;
 691
 692   // Exclude UP TO numberExcludeLeadingJets
 693   if (fNumberSignalJets < 2)
 694     numberExcludeLeadingJets = fNumberSignalJets;
 695
 696   // Search given amount of RCs
 697   Int_t numAcceptedRCs = 0;
 698   for(Int_t i=0;i<numberRandCones;i++)
 699   {
 700     Double_t tmpRandConeEta = 0.0;
 701     Double_t tmpRandConePhi = 0.0;
 702     Double_t excludedJetEta = 0.0;
 703     Double_t excludedJetPhi = 0.0;
 704
 705     // Search random cone in acceptance with no overlap with already excluded jets (leading jets and random cones)
 706     Bool_t coneValid = kTRUE;
 707
 708     // Check if etaMin/etaMax is given correctly
 709     if(etaMin < -fSignalJetEtaWindow)
 710       etaMin = -fSignalJetEtaWindow;
 711     if(etaMax > fSignalJetEtaWindow)
 712       etaMax = fSignalJetEtaWindow;
 713
 714     // Set the random cone position
 715     if ((etaMin == 0) && (etaMax == 0))
 716       tmpRandConeEta = (fTrackEtaWindow-fRandConeRadius)*(2.0*fRandom->Rndm()-1.0); // full RC is in acceptance
 717     else
 718       tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
 719
 720     tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
 721
 722     // Go through all excluded leading jets and check if there's an overlap
 723
 724     for(Int_t j=0;j<numberExcludeLeadingJets;j++)
 725     {
 726       AliEmcalJet* tmpJet = NULL;
 727
 728       if (j==0)
 729         tmpJet = fFirstLeadingJet;
 730       else if (j==1)
 731         tmpJet = fSecondLeadingJet;
 732       else
 733         AliFatal("Trying to exclude more than 2 jets in RC background -- not implemented.");
 734
 735       excludedJetPhi = tmpJet->Phi();
 736       excludedJetEta = tmpJet->Eta();
 737       Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, excludedJetPhi);
 738
 739       if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(tmpRandConeEta-excludedJetEta)*TMath::Abs(tmpRandConeEta-excludedJetEta) <= fRandConeRadius*fRandConeRadius)
 740       {
 741         coneValid = kFALSE;
 742         break;
 743       }
 744     }
 745
 746     // RC is accepted, so save it
 747     if(coneValid)
 748     {
 749       AliEmcalJet tmpJet(GetConePt(tmpRandConeEta, tmpRandConePhi, fRandConeRadius), tmpRandConeEta, tmpRandConePhi, 0.0);
 750       tmpCones[numAcceptedRCs] = tmpJet;
 751       numAcceptedRCs++;
 752     }
 753   }
 754
 755   // Calculate Rho and the mean from the RCs (no excluded jets are considered!)
 756   if(numAcceptedRCs > 0)
 757   {
 758     std::vector<Double_t> tmpRho(numAcceptedRCs);
 759     for (Int_t i=0; i<numAcceptedRCs;i++)
 760       tmpRho[i] = tmpCones[i].Pt()/(fRandConeRadius*fRandConeRadius*TMath::Pi());
 761
 762     rhoMean = TMath::Mean(tmpRho.begin(), tmpRho.end());
 763     rhoMedian = 0.0; // NOT IMPLEMENTED because TMath::Median is not working with iterators
 764   }
 765
 766   #ifdef DEBUGMODE
 767     AliInfo("Got RC background density.");
 768   #endif
 769   return numAcceptedRCs;
 770 }
 771
 772 //________________________________________________________________________
 773 void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, Double_t etaMin, Double_t etaMax)
 774 {
 775   #ifdef DEBUGMODE
 776     AliInfo("Getting TR background density.");
 777   #endif
 778
 779   Double_t summedTracksPt = 0.0;
 780
 781   // Check if etaMin/etaMax is given correctly
 782   if(etaMin < -fTrackEtaWindow)
 783     etaMin = -fTrackEtaWindow;
 784   if(etaMax > fTrackEtaWindow)
 785     etaMax = fTrackEtaWindow;
 786
 787   if ((etaMin == 0) && (etaMax == 0))
 788   {
 789     etaMin = -fTrackEtaWindow;
 790     etaMax = +fTrackEtaWindow;
 791   }
 792
 793   // Setting invalid values
 794   rhoMean = -1.0;
 795   area = -1.0;
 796   // Exclude UP TO numberExcludeLeadingJets
 797   if (fNumberSignalJets < 2)
 798     numberExcludeLeadingJets = fNumberSignalJets;
 799
 800
 801   Int_t   trackCount = fTrackArray->GetEntries();
 802   Int_t   trackCountAccepted = 0;
 803   for (Int_t i = 0; i < trackCount; i++)
 804   {
 805     Bool_t  trackValid = kTRUE;
 806     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
 807     if (IsTrackInAcceptance(tmpTrack))
 808       if ((tmpTrack->Eta() >= etaMin) && (tmpTrack->Eta() < etaMax))
 809       {
 810         for (Int_t j = 0; j < numberExcludeLeadingJets; j++)
 811         {
 812           AliEmcalJet* tmpJet = NULL;
 813           if (j==0)
 814             tmpJet = fFirstLeadingJet;
 815           else if (j==1)
 816             tmpJet = fSecondLeadingJet;
 817           else
 818             AliFatal("Trying to exclude more than 2 jets in track background -- not implemented.");
 819
 820           if (IsTrackInCone(tmpTrack, tmpJet->Eta(), tmpJet->Phi(), fTRBackgroundConeRadius))
 821           {
 822             trackValid = kFALSE;
 823             break;
 824           }
 825         }
 826         if (trackValid)
 827         {
 828           // Add track pt to array
 829           summedTracksPt = summedTracksPt + tmpTrack->Pt();
 830           trackCountAccepted++;
 831         }
 832       }
 833   }
 834
 835   if (trackCountAccepted > 0)
 836   {
 837     Double_t tmpArea = 0.0;
 838
 839     tmpArea = (2.0*fTrackEtaWindow) * TMath::TwoPi() * (etaMax-etaMin)/(2.0*fTrackEtaWindow); // area of the used eta strip
 840
 841     // Now: exclude the part of the leading jet that is in the strip
 842     if (numberExcludeLeadingJets == 2)
 843       tmpArea = tmpArea*(1.0-MCGetOverlapCircleRectancle(fFirstLeadingJet->Eta(), fFirstLeadingJet->Phi(), fTRBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi()) -MCGetOverlapCircleRectancle(fSecondLeadingJet->Eta(), fSecondLeadingJet->Phi(), fTRBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi()));
 844     else if (numberExcludeLeadingJets == 1)
 845       tmpArea = tmpArea*(1.0-MCGetOverlapCircleRectancle(fFirstLeadingJet->Eta(), fFirstLeadingJet->Phi(), fTRBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi()));
 846
 847     rhoMean = summedTracksPt/tmpArea;
 848     area  = tmpArea;
 849   }
 850
 851   #ifdef DEBUGMODE
 852     AliInfo("Got TR background density.");
 853   #endif
 854 }
 855
 856 //________________________________________________________________________
 857 void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, AliEmcalJet* excludeJet1, AliEmcalJet* excludeJet2, Bool_t doSearchPerpendicular)
 858 {
 859   #ifdef DEBUGMODE
 860     AliInfo("Getting TR background density.");
 861   #endif
 862
 863   // Setting invalid values
 864   Double_t summedTracksPt = 0.0;
 865   rhoMean = -1.0;
 866   area = -1.0;
 867
 868   Double_t tmpRadius = 0.0;
 869   if (doSearchPerpendicular)
 870     tmpRadius = 0.5*TMath::Pi(); // exclude 90 degrees around jets
 871   else
 872     tmpRadius = fSignalJetRadius;
 873
 874   numberExcludeLeadingJets = 2; // dijet is excluded here in any case
 875
 876
 877
 878   if (!fTrackArray || !fJetArray)
 879   {
 880     AliError("Could not get the track/jet array to calculate track rho!");
 881     return;
 882   }
 883
 884   Int_t   trackCount = fTrackArray->GetEntries();
 885   Int_t   trackCountAccepted = 0;
 886   for (Int_t i = 0; i < trackCount; i++)
 887   {
 888     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
 889     if (IsTrackInAcceptance(tmpTrack))
 890     {
 891       if (IsTrackInCone(tmpTrack, excludeJet1->Eta(), excludeJet1->Phi(), tmpRadius))
 892         continue;
 893
 894       if (numberExcludeLeadingJets > 1)
 895         if (IsTrackInCone(tmpTrack, excludeJet2->Eta(), excludeJet2->Phi(), tmpRadius))
 896           continue;
 897
 898         // Add track pt to array
 899         summedTracksPt = summedTracksPt + tmpTrack->Pt();
 900         trackCountAccepted++;
 901     }
 902   }
 903
 904   if (trackCountAccepted > 0)
 905   {
 906     Double_t tmpArea = 2.0*fTrackEtaWindow*TMath::TwoPi()  - 2*(tmpRadius*tmpRadius * TMath::Pi()); //TPC area - excluding jet area
 907     rhoMean = summedTracksPt/tmpArea;
 908     area = tmpArea;
 909   }
 910
 911   #ifdef DEBUGMODE
 912     AliInfo("Got TR background density.");
 913   #endif
 914 }
 915
 916 //________________________________________________________________________
 917 void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
 918 {
 919   #ifdef DEBUGMODE
 920     AliInfo("Starting Calculate().");
 921   #endif
 922   ////////////////////// NOTE: initialization & casting
 923
 924   // Additional cuts
 925   FillHistogram("hNumberEvents", 0.5); // number of events before manual cuts
 926
 927   if(!fHelperClass->IsVertexSelected2013pA(event))
 928     return;
 929
 930   FillHistogram("hNumberEvents", 1.5); // number of events after manual cuts
 931
 932   #ifdef DEBUGMODE
 933     AliInfo("Calculate()::Init done.");
 934   #endif
 935
 936   ////////////////////// NOTE: Get Centrality, (Leading)Signal jets and Background
 937
 938   // Get centrality
 939   AliCentrality* tmpCentrality = NULL;
 940   tmpCentrality = event->GetCentrality();
 941   Double_t centralityPercentile = 0.0;
 942   if (tmpCentrality != NULL)
 943     centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
 944
 945   // Get jets
 946   if (fAnalyzeBackground || fAnalyzeJets)
 947     GetSignalJets();
 948
 949   // Get background estimates
 950   Double_t              backgroundKTMedian = -1.0;
 951   Double_t              backgroundRCMean = -1.0;
 952   Double_t              backgroundRCMedian = -1.0;
 953   Double_t              backgroundTRMean = -1.0;
 954   Double_t              backgroundKTAreaMean = -1.0;
 955   Double_t              backgroundTRAreaMean = -1.0;
 956   Double_t              dijetBackground = -1.0;
 957   Double_t              dijetBackgroundPerpendicular = -1.0;
 958
 959   if (fAnalyzeBackground)
 960   {
 961     GetRCBackgroundDensity    (fNumberExcludedJets, backgroundRCMean, backgroundRCMedian);
 962     GetTRBackgroundDensity    (fNumberExcludedJets, backgroundTRMean, backgroundTRAreaMean);
 963     GetKTBackgroundDensity    (fNumberExcludedJets, backgroundKTMedian, backgroundKTAreaMean);
 964   }
 965
 966   #ifdef DEBUGMODE
 967     AliInfo("Calculate()::Centrality&SignalJets&Background-Calculation done.");
 968   #endif
 969
 970   ////////////////////// NOTE: Jet analysis and calculations
 971
 972   if (fAnalyzeJets)
 973   {
 974     // ### SIGNAL JET ANALYSIS
 975     for (Int_t i = 0; i<fNumberSignalJets; i++)
 976     {
 977       AliEmcalJet* tmpJet = fSignalJets[i];
 978
 979       // Jet spectra
 980       FillHistogram("hJetPt", tmpJet->Pt(), centralityPercentile);
 981       FillHistogram("hJetPtBgrdSubtractedRC", GetCorrectedJetPt(tmpJet, backgroundRCMean), centralityPercentile);
 982       FillHistogram("hJetPtBgrdSubtractedKT", GetCorrectedJetPt(tmpJet, backgroundKTMedian), centralityPercentile);
 983       FillHistogram("hJetPtBgrdSubtractedTR", GetCorrectedJetPt(tmpJet, backgroundTRMean), centralityPercentile);
 984
 985       // Jet spectra  (pos+neg eta)
 986       if(tmpJet->Eta() > 0.2)
 987       {
 988         FillHistogram("hJetPtBgrdSubtractedRCPosEta", GetCorrectedJetPt(tmpJet, backgroundRCMean), centralityPercentile);
 989         FillHistogram("hJetPtBgrdSubtractedKTPosEta", GetCorrectedJetPt(tmpJet, backgroundKTMedian), centralityPercentile);
 990         FillHistogram("hJetPtBgrdSubtractedTRPosEta", GetCorrectedJetPt(tmpJet, backgroundTRMean), centralityPercentile);
 991       }
 992       else if(tmpJet->Eta() < -0.2)
 993       {
 994         FillHistogram("hJetPtBgrdSubtractedRCNegEta", GetCorrectedJetPt(tmpJet, backgroundRCMean), centralityPercentile);
 995         FillHistogram("hJetPtBgrdSubtractedKTNegEta", GetCorrectedJetPt(tmpJet, backgroundKTMedian), centralityPercentile);
 996         FillHistogram("hJetPtBgrdSubtractedTRNegEta", GetCorrectedJetPt(tmpJet, backgroundTRMean), centralityPercentile);
 997       }
 998
 999       // Correct EVERY jet with suitable background
1000       Double_t tmpKTRho = 0.0;
1001       Double_t tmpRCRho = 0.0;
1002       Double_t tmpTRRho = 0.0;
1003       Double_t dummy = 0.0;
1004       // Calculate backgrounds
1005       GetKTBackgroundDensity (fNumberExcludedJets, tmpKTRho, dummy, tmpJet->Eta()-0.1, tmpJet->Eta()+0.1, kTRUE);
1006       GetRCBackgroundDensity (fNumberExcludedJets, tmpRCRho, dummy, tmpJet->Eta()-0.1, tmpJet->Eta()+0.1);
1007       GetTRBackgroundDensity (fNumberExcludedJets, tmpTRRho, dummy, tmpJet->Eta()-0.1, tmpJet->Eta()+0.1);
1008
1009       FillHistogram("hKTBackgroundEtaWindow", tmpKTRho, centralityPercentile);
1010       FillHistogram("hRCBackgroundEtaWindow", tmpRCRho, centralityPercentile);
1011       FillHistogram("hTRBackgroundEtaWindow", tmpTRRho, centralityPercentile);
1012
1013       FillHistogram("hJetPtBgrdSubtractedKTEtaWindow", GetCorrectedJetPt(tmpJet, tmpKTRho), centralityPercentile);
1014       FillHistogram("hJetPtBgrdSubtractedRCEtaWindow", GetCorrectedJetPt(tmpJet, tmpRCRho), centralityPercentile);
1015       FillHistogram("hJetPtBgrdSubtractedTREtaWindow", GetCorrectedJetPt(tmpJet, tmpTRRho), centralityPercentile);
1016
1017
1018       // Signal jet vs. signal jet - "Combinatorial"
1019       for (Int_t j = i+1; j<fNumberSignalJets; j++)
1020         FillHistogram("hJetDeltaPhi", GetDeltaPhi(tmpJet->Phi(), fSignalJets[j]->Phi()));
1021     }
1022
1023     // ### DIJETS
1024     if(fNumberSignalJets >= 2)
1025     {
1026       FillHistogram("hLeadingJetDeltaPhi", GetDeltaPhi(fFirstLeadingJet->Phi(), fSecondLeadingJet->Phi()));
1027
1028       if (IsDijet(fFirstLeadingJet, fSecondLeadingJet))
1029       {
1030         FillHistogram("hDijetConstituentsPt", fFirstLeadingJet->Pt());
1031         FillHistogram("hDijetConstituentsPt", fSecondLeadingJet->Pt());
1032
1033         FillHistogram("hDijetLeadingJetPt", fFirstLeadingJet->Pt());
1034         FillHistogram("hDijetPtCorrelation", fFirstLeadingJet->Pt(), fSecondLeadingJet->Pt());
1035         Double_t dummyArea = 0;
1036         GetTRBackgroundDensity (2, dijetBackground, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kFALSE);
1037         GetTRBackgroundDensity (2, dijetBackgroundPerpendicular, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kTRUE);
1038       }
1039     }
1040
1041     // ### SOME JET PLOTS
1042     FillHistogram("hJetCountAll", fJetArray->GetEntries());
1043     FillHistogram("hJetCountAccepted", fNumberSignalJets);
1044     if (fFirstLeadingJet)
1045       FillHistogram("hLeadingJetPt", fFirstLeadingJet->Pt());
1046     if (fSecondLeadingJet)
1047       FillHistogram("hSecondLeadingJetPt", fSecondLeadingJet->Pt());
1048
1049   } //endif AnalyzeJets
1050
1051   #ifdef DEBUGMODE
1052     AliInfo("Calculate()::Jets done.");
1053   #endif
1054   ////////////////////// NOTE: Background analysis
1055
1056   if (fAnalyzeBackground)
1057   {
1058     // Calculate (non-eta corrected) background in centrality classes
1059     FillHistogram("hRCBackground", backgroundRCMean, centralityPercentile);
1060     FillHistogram("hTRBackground", backgroundTRMean, centralityPercentile);
1061     FillHistogram("hKTBackground", backgroundKTMedian, centralityPercentile);
1062
1063     // Calculate backgrounds in eta bins
1064     for(Int_t i=0;i<5;i++)
1065     {
1066       Double_t dummy = 0.0;
1067       Double_t tmpKTRho = 0.0;
1068       Double_t tmpRCRho = 0.0;
1069       Double_t tmpTRRho = 0.0;
1070
1071       Double_t etaMin = -(fTrackEtaWindow-fSignalJetRadius) + 2*(fTrackEtaWindow-fSignalJetRadius)/5 *  i;
1072       Double_t etaMax = -(fTrackEtaWindow-fSignalJetRadius) + 2*(fTrackEtaWindow-fSignalJetRadius)/5 * (i+1);
1073
1074       // Calculate backgrounds
1075       GetKTBackgroundDensity    (fNumberExcludedJets, tmpKTRho, dummy, etaMin, etaMax);
1076       GetTRBackgroundDensity    (fNumberExcludedJets, tmpTRRho, dummy, etaMin, etaMax);
1077       GetRCBackgroundDensity    (fNumberExcludedJets, tmpRCRho, dummy, etaMin, etaMax);
1078
1079       FillHistogram("hRCBackgroundEtaBins", tmpRCRho, (etaMin+etaMax)/2.0);
1080       FillHistogram("hTRBackgroundEtaBins", tmpTRRho, (etaMin+etaMax)/2.0);
1081       FillHistogram("hKTBackgroundEtaBins", tmpKTRho, (etaMin+etaMax)/2.0);
1082     }
1083
1084     // In case of dijets -> look at the background
1085     if (dijetBackground >= 0)
1086       FillHistogram("hDijetBackground", dijetBackground, centralityPercentile);
1087     if (dijetBackgroundPerpendicular >= 0)
1088       FillHistogram("hDijetBackgroundPerpendicular", dijetBackgroundPerpendicular, centralityPercentile);
1089
1090
1091     // Calculate the delta pt
1092     Double_t tmpDeltaPtKT = 0.0;
1093     Double_t tmpDeltaPtRC = 0.0;
1094     Double_t tmpDeltaPtTR = 0.0;
1095     Double_t tmpDeltaPtKTNoExcl = 0.0;
1096     Double_t tmpDeltaPtRCNoExcl = 0.0;
1097     Double_t tmpDeltaPtTRNoExcl = 0.0;
1098     GetDeltaPt(tmpDeltaPtKT, backgroundKTMedian);
1099     GetDeltaPt(tmpDeltaPtRC, backgroundRCMean);
1100     GetDeltaPt(tmpDeltaPtTR, backgroundTRMean);
1101     GetDeltaPt(tmpDeltaPtKTNoExcl, backgroundKTMedian, kFALSE);
1102     GetDeltaPt(tmpDeltaPtRCNoExcl, backgroundRCMean, kFALSE);
1103     GetDeltaPt(tmpDeltaPtTRNoExcl, backgroundTRMean, kFALSE);
1104
1105     // If valid, fill the delta pt histograms
1106     if(tmpDeltaPtKT > -10000.0)
1107       FillHistogram("hDeltaPtKT", tmpDeltaPtKT, centralityPercentile);
1108     if(tmpDeltaPtRC > -10000.0)
1109       FillHistogram("hDeltaPtRC", tmpDeltaPtRC, centralityPercentile);
1110     if(tmpDeltaPtTR > -10000.0)
1111       FillHistogram("hDeltaPtTR", tmpDeltaPtTR, centralityPercentile);
1112     if(tmpDeltaPtKTNoExcl > -10000.0)
1113       FillHistogram("hDeltaPtKTNoExcl", tmpDeltaPtKTNoExcl, centralityPercentile);
1114     if(tmpDeltaPtRCNoExcl > -10000.0)
1115       FillHistogram("hDeltaPtRCNoExcl", tmpDeltaPtRCNoExcl, centralityPercentile);
1116     if(tmpDeltaPtTRNoExcl > -10000.0)
1117       FillHistogram("hDeltaPtTRNoExcl", tmpDeltaPtTRNoExcl, centralityPercentile);
1118
1119   }
1120
1121   #ifdef DEBUGMODE
1122     AliInfo("Calculate()::Background done.");
1123   #endif
1124
1125   ////////////////////// NOTE: Pythia histograms
1126   if(fAnalyzePythia)
1127   {
1128     FillHistogram("hPythiaPtHard", GetPtHard());
1129     FillHistogram("hPythiaNTrials", GetPtHardBin()-0.1, fTrials);
1130     FillHistogram("hPythiaXSection", GetPtHardBin()-0.1, fCrossSection);
1131
1132     #ifdef DEBUGMODE
1133       AliInfo("Calculate()::Pythia done.");
1134     #endif
1135   }
1136   #ifdef DEBUGMODE
1137     AliInfo("Calculate() done.");
1138   #endif
1139 }
1140
1141 //________________________________________________________________________
1142 Bool_t AliAnalysisTaskChargedJetsPA::Notify()
1143 {
1144   // Implemented Notify() to read the cross sections
1145   // and number of trials from pyxsec.root
1146   //
1147   #ifdef DEBUGMODE
1148     AliInfo("Notify started.");
1149   #endif
1150
1151   if(fAnalyzePythia)
1152   {
1153     TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1154     TFile *currFile = tree->GetCurrentFile();
1155
1156     TString file(currFile->GetName());
1157
1158     if(file.Contains("root_archive.zip#")){
1159       Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1160       Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1161       file.Replace(pos+1,20,"");
1162     }
1163     else {
1164       // not an archive take the basename....
1165       file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1166     }
1167
1168     TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
1169     if(!fxsec){
1170       // next trial fetch the histgram file
1171       fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1172       if(!fxsec){
1173           // not a severe condition but inciate that we have no information
1174         return kFALSE;
1175       }
1176       else{
1177         // find the tlist we want to be independtent of the name so use the Tkey
1178         TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1179         if(!key){
1180           fxsec->Close();
1181           return kFALSE;
1182         }
1183         TList *list = dynamic_cast<TList*>(key->ReadObj());
1184         if(!list){
1185           fxsec->Close();
1186           return kFALSE;
1187         }
1188         fCrossSection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1189         fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1190         fxsec->Close();
1191       }
1192     } // no tree pyxsec.root
1193     else {
1194       TTree *xtree = (TTree*)fxsec->Get("Xsection");
1195       if(!xtree){
1196         fxsec->Close();
1197         return kFALSE;
1198       }
1199       UInt_t   ntrials  = 0;
1200       Double_t  xsection  = 0;
1201       xtree->SetBranchAddress("xsection",&xsection);
1202       xtree->SetBranchAddress("ntrials",&ntrials);
1203       xtree->GetEntry(0);
1204       fTrials = ntrials;
1205       fCrossSection = xsection;
1206       fxsec->Close();
1207     }
1208     #ifdef DEBUGMODE
1209       AliInfo("Notify ended.");
1210     #endif
1211   }
1212   return kTRUE;
1213 }
1214
1215 //________________________________________________________________________
1216 inline Double_t AliAnalysisTaskChargedJetsPA::EtaToTheta(Double_t arg)
1217   {return 2.*atan(exp(-arg));}
1218 //________________________________________________________________________
1219 inline Double_t AliAnalysisTaskChargedJetsPA::ThetaToEta(Double_t arg)
1220 {
1221   if ((arg > TMath::Pi()) || (arg < 0.0))
1222   {
1223     AliError(Form("ThetaToEta got wrong input! (%f)", arg));
1224     return 0.0;
1225   }
1226   return -log(tan(arg/2.));
1227 }
1228 //________________________________________________________________________
1229 inline Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPhi(Double_t phi1, Double_t phi2)
1230   {return min(TMath::Abs(phi1-phi2),TMath::TwoPi() - TMath::Abs(phi1-phi2));}
1231
1232 //________________________________________________________________________
1233 Double_t AliAnalysisTaskChargedJetsPA::MCGetOverlapCircleRectancle(Double_t cPosX, Double_t cPosY, Double_t cRadius, Double_t rPosXmin, Double_t rPosXmax, Double_t rPosYmin, Double_t rPosYmax)
1234 {
1235   const Int_t kTests = 1000;
1236   Int_t hits = 0;
1237   TRandom3 randomGen(0);
1238
1239   // Loop over kTests-many tests
1240   for (Int_t i=0; i<kTests; i++)
1241   {
1242     //Choose random position in rectangle for the tester
1243     Double_t tmpTestX = randomGen.Uniform(rPosXmin, rPosXmax);
1244     Double_t tmpTestY = randomGen.Uniform(rPosYmin, rPosYmax);
1245
1246     //Check, if tester is in circle. If yes, increment circle counter.
1247     Double_t tmpDistance = TMath::Sqrt( (tmpTestX - cPosX)*(tmpTestX - cPosX) + (tmpTestY - cPosY)*(tmpTestY - cPosY) );
1248     if(tmpDistance < cRadius)
1249       hits++;
1250   }
1251
1252   // return ratio
1253   return (static_cast<Double_t>(hits)/static_cast<Double_t>(kTests));
1254 }
1255
1256 //________________________________________________________________________
1257 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x)
1258 {
1259   TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(key)));
1260   if(!tmpHist)
1261   {
1262     AliWarning(Form("Cannot find histogram <%s> ",key)) ;
1263     return;
1264   }
1265
1266   tmpHist->Fill(x);
1267 }
1268
1269 //________________________________________________________________________
1270 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y)
1271 {
1272   TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(key)));
1273   if(!tmpHist)
1274   {
1275     AliWarning(Form("Cannot find histogram <%s> ",key));
1276     return;
1277   }
1278
1279   if (tmpHist->IsA()->GetBaseClass("TH1"))
1280     static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1281   else if (tmpHist->IsA()->GetBaseClass("TH2"))
1282     static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1283 }
1284
1285 //________________________________________________________________________
1286 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y, Double_t add)
1287 {
1288   TH2* tmpHist = static_cast<TH2*>(fOutputList->FindObject(GetHistoName(key)));
1289   if(!tmpHist)
1290   {
1291     AliWarning(Form("Cannot find histogram <%s> ",key));
1292     return;
1293   }
1294
1295   tmpHist->Fill(x,y,add);
1296 }
1297 //________________________________________________________________________
1298 template <class T> T* AliAnalysisTaskChargedJetsPA::AddHistogram1D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, const char* xTitle, const char* yTitle)
1299 {
1300   T* tmpHist = new T(GetHistoName(name), GetHistoName(title), xBins, xMin, xMax);
1301
1302   tmpHist->GetXaxis()->SetTitle(xTitle);
1303   tmpHist->GetYaxis()->SetTitle(yTitle);
1304   tmpHist->SetOption(options);
1305   tmpHist->SetMarkerStyle(kFullCircle);
1306   tmpHist->Sumw2();
1307
1308   fHistList->Add(tmpHist);
1309   fHistCount++;
1310
1311   return tmpHist;
1312 }
1313
1314 //________________________________________________________________________
1315 template <class T> T* AliAnalysisTaskChargedJetsPA::AddHistogram2D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, Int_t yBins, Double_t yMin, Double_t yMax, const char* xTitle, const char* yTitle, const char* zTitle)
1316 {
1317   T* tmpHist = new T(GetHistoName(name), GetHistoName(title), xBins, xMin, xMax, yBins, yMin, yMax);
1318   tmpHist->GetXaxis()->SetTitle(xTitle);
1319   tmpHist->GetYaxis()->SetTitle(yTitle);
1320   tmpHist->GetZaxis()->SetTitle(zTitle);
1321   tmpHist->SetOption(options);
1322   tmpHist->SetMarkerStyle(kFullCircle);
1323   tmpHist->Sumw2();
1324
1325   fHistList->Add(tmpHist);
1326   fHistCount++;
1327
1328   return tmpHist;
1329 }
1330
1331 //________________________________________________________________________
1332 void AliAnalysisTaskChargedJetsPA::Terminate(Option_t *)
1333 {
1334   PostData(1, fOutputList);
1335
1336   // Mandatory
1337   fOutputList = dynamic_cast<TList*> (GetOutputData(1)); // '1' refers to the output slot
1338   if (!fOutputList) {
1339     printf("ERROR: Output list not available\n");
1340     return;
1341   }
1342 }
1343
1344 //________________________________________________________________________
1345 AliAnalysisTaskChargedJetsPA::~AliAnalysisTaskChargedJetsPA()
1346 {
1347   // Destructor. Clean-up the output list, but not the histograms that are put inside
1348   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
1349   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
1350     delete fOutputList;
1351   }
1352 }
1353
1354 //________________________________________________________________________
1355 void AliAnalysisTaskChargedJetsPA::UserCreateOutputObjects()
1356 {
1357   // called once to create user defined output objects like histograms, plots etc.
1358   // and to put it on the output list.
1359   // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects.
1360
1361   fRandom = new TRandom3(0);
1362
1363   fOutputList = new TList();
1364   fOutputList->SetOwner(); // otherwise it produces leaks in merging
1365
1366   PostData(1, fOutputList);
1367 }
1368
1369 //________________________________________________________________________
1370 void AliAnalysisTaskChargedJetsPA::UserExec(Option_t *)
1371 {
1372   #ifdef DEBUGMODE
1373     AliInfo("UserExec() started.");
1374   #endif
1375
1376   if (!InputEvent())
1377   {
1378     AliError("??? Event pointer == 0 ???");
1379     return;
1380   }
1381
1382   if (!fInitialized)
1383     ExecOnce(); // Get tracks, jets, background from arrays if not already given + Init Histos
1384
1385   Calculate(InputEvent());
1386
1387   PostData(1, fOutputList);
1388 }