PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullpAJets.cxx

   1 #include "AliAnalysisTaskFullpAJets.h"
   2
   3 #include <Riostream.h>
   4 #include <ctime>
   5 #include <TString.h>
   6 #include <TChain.h>
   7 #include <TTree.h>
   8 #include <TH1D.h>
   9 #include <TH2D.h>
  10 #include <TH3D.h>
  11 #include <TCanvas.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14 #include <TProfile.h>
  15 #include <TProfile2D.h>
  16 #include <TProfile3D.h>
  17 #include <TRandom.h>
  18 #include <TRandom3.h>
  19
  20 #include "AliAnalysisTaskSE.h"
  21 #include "AliAnalysisManager.h"
  22 #include "AliStack.h"
  23 #include "AliESDtrackCuts.h"
  24 #include "AliESDEvent.h"
  25 #include "AliESDInputHandler.h"
  26 #include "AliAODEvent.h"
  27 #include "AliMCEvent.h"
  28 #include "AliEmcalJet.h"
  29 #include "AliEMCALGeometry.h"
  30
  31 ClassImp(AliAnalysisTaskFullpAJets)
  32
  33 //________________________________________________________________________
  34 AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
  35     AliAnalysisTaskSE(),
  36
  37     fOutput(0),
  38     fTrackCuts(0),
  39     fhTrackPt(0),
  40     fhTrackEta(0),
  41     fhTrackPhi(0),
  42     fhClusterPt(0),
  43     fhClusterEta(0),
  44     fhClusterPhi(0),
  45     fhCentrality(0),
  46     fhBckgMult(0),
  47     fhBckgFluc(0),
  48     fhChargedJetPt(0),
  49     fhChargedJetPtAreaCut(0),
  50     fhJetPtEMCal(0),
  51     fhJetPtEMCalAreaCut(0),
  52     fhJetPtEMCalAreaCutSignal(0),
  53     fhJetPtTPC(0),
  54     fhJetPtTPCAreaCut(0),
  55     fhJetTPtRhoTotal(0),
  56     fhJetTPtRhoTotalSignal(0),
  57     fhJetTPtRhoNoLeading(0),
  58     fhJetTPtRhoNoLeadingSignal(0),
  59     fhJetTPt1B(0),
  60     fhJetTPt1BSignal(0),
  61     fhEMCalBckg1B(0),
  62     fhJetTPt1C(0),
  63     fhEMCalBckg1C(0),
  64     fhEMCalJet2A(0),
  65     fhJetTPt2B(0),
  66     fhEMCalBckg2B(0),
  67     fhJetTPt3(0),
  68     fhDeltaPtTotal(0),
  69     fhDeltaPtNoLeading(0),
  70     fhDeltaPt1B(0),
  71     fhDeltaRho01(0),
  72     fhEMCalCellCounts(0),
  73     fh020RhoTotal(0),
  74     fh020RhoNoLeading(0),
  75     fh020Rho1B(0),
  76     fh020Rho2B(0),
  77     fh020Rho3(0),
  78     fh020JetPtEMCal(0),
  79     fh020JetPtEMCalAreaCut(0),
  80     fh020JetPtEMCalAreaCutSignal(0),
  81     fh020JetTPtRhoTotal(0),
  82     fh020JetTPtRhoTotalSignal(0),
  83     fh020JetTPtRhoNoLeading(0),
  84     fh020JetTPtRhoNoLeadingSignal(0),
  85     fh020JetTPt1B(0),
  86     fh020JetTPt1BSignal(0),
  87     fh020JetTPt1C(0),
  88     fh020JetTPt2B(0),
  89     fh020JetTPt3(0),
  90     fhDeltaPt2B(0),
  91     fhDeltaPtkT(0),
  92     fh020DiJetAsy(0),
  93     fh020RhokT(0),
  94     fh020EMCalkTClusters(0),
  95     fh020EMCalAkTJets(0),
  96     fh020DiJetDeltaPhi(0),
  97     fhDiJetEMCalLeadingPt(0),
  98     fhDiJetEMCalLeadingDeltaPhi(0),
  99     fh020EMCalJet2A(0),
 100     fhDeltaRho0DiJet(0),
 101     fh020Rho2BCore(0),
 102     fh020Rho3NoJets(0),
 103     fh020Rho3DiJets(0),
 104     fh020Rho3Perp(0),
 105
 106     fhTrackEtaPhi(0),
 107     fhClusterEtaPhi(0),
 108     fhJetPtArea(0),
 109     fhRhoTotal(0),
 110     fhRhoNoLeading(0),
 111     fhRho1B(0),
 112     fhRho1C(0),
 113     fhRho2B(0),
 114     fhRho3(0),
 115     fhJetConstituentPt(0),
 116     fhJetPtCenEMCal(0),
 117     fhJetPtCenEMCalAreaCut(0),
 118     fhJetPtCenEMCalAreaCutSignal(0),
 119     fhJetTPtCenRhoTotal(0),
 120     fhJetTPtCenRhoTotalSignal(0),
 121     fhJetTPtCenRhoNoLeading(0),
 122     fhJetTPtCenRhoNoLeadingSignal(0),
 123     fhJetTPtCen1B(0),
 124     fhJetTPtCen1BSignal(0),
 125     fhJetTPtCen1C(0),
 126     fhJetTPtCen2B(0),
 127     fhJetTPtCen3(0),
 128     fhDiJetCenAsy(0),
 129     fhDiJetCenDeltaPhi(0),
 130     fhEMCalCenJet2A(0),
 131     fhRho2BCore(0),
 132     fhRho3NoJets(0),
 133     fhRho3DiJets(0),
 134     fhRho3Perp(0),
 135
 136     fhJetTrigR1A(0),
 137
 138     fpEventMult(0),
 139     fpRhoTotal(0),
 140     fpRhoNoLeading(0),
 141     fpRho1B(0),
 142     fpRho2B(0),
 143     fpRho3(0),
 144     fpRhoScale(0),
 145     fpRhokT(0),
 146     fpJetPtRhoTotal(0),
 147     fpJetPtRhoNoLeading(0),
 148     fpJetPtRhokT(0),
 149     fpRhoChargedkT(0),
 150     fpRhoScalekT(0),
 151     fpRho2BCore(0),
 152     fpRho3NoJets(0),
 153     fpRho3DiJets(0),
 154     fpRho3Perp(0),
 155
 156     fpTrackPtProfile(0),
 157     fpClusterPtProfile(0),
 158
 159     fIsInitialized(0),
 160     fRJET(0),
 161     fnEvents(0),
 162     fnEventsCharged(0),
 163     fnDiJetEvents(0),
 164     fEMCalPhiMin(0),
 165     fEMCalPhiMax(0),
 166     fEMCalPhiTotal(0),
 167     fEMCalEtaMin(0),
 168     fEMCalEtaMax(0),
 169     fEMCalEtaTotal(0),
 170     fEMCalArea(0),
 171     fTPCPhiMin(0),
 172     fTPCPhiMax(0),
 173     fTPCPhiTotal(0),
 174     fTPCEtaMin(0),
 175     fTPCEtaMax(0),
 176     fTPCEtaTotal(0),
 177     fTPCArea(0),
 178     fJetR(0),
 179     fJetAreaCutFrac(0),
 180     fJetAreaThreshold(0),
 181     fDeltaRho01(0),
 182     fnEMCalCells(0),
 183     fCentralityBins(0),
 184     fCentralityLow(0),
 185     fCentralityUp(0),
 186     fEventCentrality(0),
 187     fRhoTotal(0),
 188     fRhoCharged(0),
 189     fRhokTTotal(0),
 190     fRhokTCharged(0),
 191     fRhoAkTTotal(0),
 192     fEtaProfileBins(0),
 193     fEtaProfileLow(0),
 194     fEtaProfileUp(0),
 195     fEDProfileRBins(0),
 196     fEDProfileRLow(0),
 197     fEDProfileRUp(0),
 198     fEDProfilePtBins(0),
 199     fEDProfilePtLow(0),
 200     fEDProfilePtUp(0),
 201     fEDProfileEtaBins(0),
 202     fEDProfileEtaLow(0),
 203     fEDProfileEtaUp(0),
 204     fnTracks(0),
 205     fnClusters(0),
 206     fnAKTFullJets(0),
 207     fnAKTChargedJets(0),
 208     fnKTFullJets(0),
 209     fnKTChargedJets(0),
 210     fnBckgClusters(0),
 211     fTPCJetThreshold(0),
 212     fEMCalJetThreshold(0),
 213     fVertexWindow(0),
 214     fVertexMaxR(0),
 215     fnJetsPtCut(0),
 216     fnJetsPtTPCCut(0),
 217     fnJetsPtTotalCut(0),
 218     fnJetsChargedPtCut(0),
 219     fnJetskTEMCalFull(0),
 220     fnJetskTTPCFull(0),
 221     fPtMaxID(0),
 222     fPtFullMaxID(0),
 223     fPtTPCMaxID(0),
 224     fPtFullTPCMaxID(0),
 225     fPtTotalMaxID(0),
 226     fPtChargedMaxID(0),
 227     fPtMax(0),
 228     fPtFullMax(0),
 229     fPtTPCMax(0),
 230     fPtFullTPCMax(0),
 231     fPtTotalMax(0),
 232     fPtChargedMax(0),
 233     fChargedBackJetID(0),
 234     fmyTracks(0),
 235     fmyClusters(0),
 236     fmyAKTFullJets(0),
 237     fmyAKTChargedJets(0),
 238     fmyKTFullJets(0),
 239     fmyKTChargedJets(0),
 240     fJetPtCutID(0),
 241     fJetPtTPCCutID(0),
 242     fJetPtTotalCutID(0),
 243     fJetPtChargedCutID(0),
 244     fJetkTEMCalFullID(0),
 245     fJetkTTPCFullID(0),
 246     fInEMCal(0),
 247     fInEMCalFull(0),
 248     fInTPCFull(0),
 249     fInTPCChargedFull(0),
 250     fRCBckgFluc(0)
 251 {
 252     // Dummy constructor ALWAYS needed for I/O.
 253     fpJetEtaProfile = new TProfile *[14];
 254     fpJetAbsEtaProfile = new TProfile *[14];
 255     fpChargedJetRProfile = new TProfile *[8];
 256     fpJetRProfile= new TProfile *[4];
 257     fpChargedJetEDProfile= new TProfile3D *[10];
 258     fpJetEDProfile= new TProfile3D *[10];
 259     fvertex[0]=0.0,fvertex[1]=0.0,fvertex[2]=0.0;
 260 }
 261
 262 //________________________________________________________________________
 263 AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
 264     AliAnalysisTaskSE(name),
 265
 266     fOutput(0),
 267     fTrackCuts(0),
 268     fhTrackPt(0),
 269     fhTrackEta(0),
 270     fhTrackPhi(0),
 271     fhClusterPt(0),
 272     fhClusterEta(0),
 273     fhClusterPhi(0),
 274     fhCentrality(0),
 275     fhBckgMult(0),
 276     fhBckgFluc(0),
 277     fhChargedJetPt(0),
 278     fhChargedJetPtAreaCut(0),
 279     fhJetPtEMCal(0),
 280     fhJetPtEMCalAreaCut(0),
 281     fhJetPtEMCalAreaCutSignal(0),
 282     fhJetPtTPC(0),
 283     fhJetPtTPCAreaCut(0),
 284     fhJetTPtRhoTotal(0),
 285     fhJetTPtRhoTotalSignal(0),
 286     fhJetTPtRhoNoLeading(0),
 287     fhJetTPtRhoNoLeadingSignal(0),
 288     fhJetTPt1B(0),
 289     fhJetTPt1BSignal(0),
 290     fhEMCalBckg1B(0),
 291     fhJetTPt1C(0),
 292     fhEMCalBckg1C(0),
 293     fhEMCalJet2A(0),
 294     fhJetTPt2B(0),
 295     fhEMCalBckg2B(0),
 296     fhJetTPt3(0),
 297     fhDeltaPtTotal(0),
 298     fhDeltaPtNoLeading(0),
 299     fhDeltaPt1B(0),
 300     fhDeltaRho01(0),
 301     fhEMCalCellCounts(0),
 302     fh020RhoTotal(0),
 303     fh020RhoNoLeading(0),
 304     fh020Rho1B(0),
 305     fh020Rho2B(0),
 306     fh020Rho3(0),
 307     fh020JetPtEMCal(0),
 308     fh020JetPtEMCalAreaCut(0),
 309     fh020JetPtEMCalAreaCutSignal(0),
 310     fh020JetTPtRhoTotal(0),
 311     fh020JetTPtRhoTotalSignal(0),
 312     fh020JetTPtRhoNoLeading(0),
 313     fh020JetTPtRhoNoLeadingSignal(0),
 314     fh020JetTPt1B(0),
 315     fh020JetTPt1BSignal(0),
 316     fh020JetTPt1C(0),
 317     fh020JetTPt2B(0),
 318     fh020JetTPt3(0),
 319     fhDeltaPt2B(0),
 320     fhDeltaPtkT(0),
 321     fh020DiJetAsy(0),
 322     fh020RhokT(0),
 323     fh020EMCalkTClusters(0),
 324     fh020EMCalAkTJets(0),
 325     fh020DiJetDeltaPhi(0),
 326     fhDiJetEMCalLeadingPt(0),
 327     fhDiJetEMCalLeadingDeltaPhi(0),
 328     fh020EMCalJet2A(0),
 329     fhDeltaRho0DiJet(0),
 330     fh020Rho2BCore(0),
 331     fh020Rho3NoJets(0),
 332     fh020Rho3DiJets(0),
 333     fh020Rho3Perp(0),
 334
 335     fhTrackEtaPhi(0),
 336     fhClusterEtaPhi(0),
 337     fhJetPtArea(0),
 338     fhRhoTotal(0),
 339     fhRhoNoLeading(0),
 340     fhRho1B(0),
 341     fhRho1C(0),
 342     fhRho2B(0),
 343     fhRho3(0),
 344     fhJetConstituentPt(0),
 345     fhJetPtCenEMCal(0),
 346     fhJetPtCenEMCalAreaCut(0),
 347     fhJetPtCenEMCalAreaCutSignal(0),
 348     fhJetTPtCenRhoTotal(0),
 349     fhJetTPtCenRhoTotalSignal(0),
 350     fhJetTPtCenRhoNoLeading(0),
 351     fhJetTPtCenRhoNoLeadingSignal(0),
 352     fhJetTPtCen1B(0),
 353     fhJetTPtCen1BSignal(0),
 354     fhJetTPtCen1C(0),
 355     fhJetTPtCen2B(0),
 356     fhJetTPtCen3(0),
 357     fhDiJetCenAsy(0),
 358     fhDiJetCenDeltaPhi(0),
 359     fhEMCalCenJet2A(0),
 360     fhRho2BCore(0),
 361     fhRho3NoJets(0),
 362     fhRho3DiJets(0),
 363     fhRho3Perp(0),
 364
 365     fhJetTrigR1A(0),
 366
 367     fpEventMult(0),
 368     fpRhoTotal(0),
 369     fpRhoNoLeading(0),
 370     fpRho1B(0),
 371     fpRho2B(0),
 372     fpRho3(0),
 373     fpRhoScale(0),
 374     fpRhokT(0),
 375     fpJetPtRhoTotal(0),
 376     fpJetPtRhoNoLeading(0),
 377     fpJetPtRhokT(0),
 378     fpRhoChargedkT(0),
 379     fpRhoScalekT(0),
 380     fpRho2BCore(0),
 381     fpRho3NoJets(0),
 382     fpRho3DiJets(0),
 383     fpRho3Perp(0),
 384
 385     fpTrackPtProfile(0),
 386     fpClusterPtProfile(0),
 387
 388     fIsInitialized(0),
 389     fRJET(0),
 390     fnEvents(0),
 391     fnEventsCharged(0),
 392     fnDiJetEvents(0),
 393     fEMCalPhiMin(0),
 394     fEMCalPhiMax(0),
 395     fEMCalPhiTotal(0),
 396     fEMCalEtaMin(0),
 397     fEMCalEtaMax(0),
 398     fEMCalEtaTotal(0),
 399     fEMCalArea(0),
 400     fTPCPhiMin(0),
 401     fTPCPhiMax(0),
 402     fTPCPhiTotal(0),
 403     fTPCEtaMin(0),
 404     fTPCEtaMax(0),
 405     fTPCEtaTotal(0),
 406     fTPCArea(0),
 407     fJetR(0),
 408     fJetAreaCutFrac(0),
 409     fJetAreaThreshold(0),
 410     fDeltaRho01(0),
 411     fnEMCalCells(0),
 412     fCentralityBins(0),
 413     fCentralityLow(0),
 414     fCentralityUp(0),
 415     fEventCentrality(0),
 416     fRhoTotal(0),
 417     fRhoCharged(0),
 418     fRhokTTotal(0),
 419     fRhokTCharged(0),
 420     fRhoAkTTotal(0),
 421     fEtaProfileBins(0),
 422     fEtaProfileLow(0),
 423     fEtaProfileUp(0),
 424     fEDProfileRBins(0),
 425     fEDProfileRLow(0),
 426     fEDProfileRUp(0),
 427     fEDProfilePtBins(0),
 428     fEDProfilePtLow(0),
 429     fEDProfilePtUp(0),
 430     fEDProfileEtaBins(0),
 431     fEDProfileEtaLow(0),
 432     fEDProfileEtaUp(0),
 433     fnTracks(0),
 434     fnClusters(0),
 435     fnAKTFullJets(0),
 436     fnAKTChargedJets(0),
 437     fnKTFullJets(0),
 438     fnKTChargedJets(0),
 439     fnBckgClusters(0),
 440     fTPCJetThreshold(0),
 441     fEMCalJetThreshold(0),
 442     fVertexWindow(0),
 443     fVertexMaxR(0),
 444     fnJetsPtCut(0),
 445     fnJetsPtTPCCut(0),
 446     fnJetsPtTotalCut(0),
 447     fnJetsChargedPtCut(0),
 448     fnJetskTEMCalFull(0),
 449     fnJetskTTPCFull(0),
 450     fPtMaxID(0),
 451     fPtFullMaxID(0),
 452     fPtTPCMaxID(0),
 453     fPtFullTPCMaxID(0),
 454     fPtTotalMaxID(0),
 455     fPtChargedMaxID(0),
 456     fPtMax(0),
 457     fPtFullMax(0),
 458     fPtTPCMax(0),
 459     fPtFullTPCMax(0),
 460     fPtTotalMax(0),
 461     fPtChargedMax(0),
 462     fChargedBackJetID(0),
 463     fmyTracks(0),
 464     fmyClusters(0),
 465     fmyAKTFullJets(0),
 466     fmyAKTChargedJets(0),
 467     fmyKTFullJets(0),
 468     fmyKTChargedJets(0),
 469     fJetPtCutID(0),
 470     fJetPtTPCCutID(0),
 471     fJetPtTotalCutID(0),
 472     fJetPtChargedCutID(0),
 473     fJetkTEMCalFullID(0),
 474     fJetkTTPCFullID(0),
 475     fInEMCal(0),
 476     fInEMCalFull(0),
 477     fInTPCFull(0),
 478     fInTPCChargedFull(0),
 479     fRCBckgFluc(0)
 480 {
 481     // Constructor
 482     // Define input and output slots here (never in the dummy constructor)
 483     // Input slot #0 works with a TChain - it is connected to the default input container
 484     // Output slot #1 writes into a TH1 container
 485     fpJetEtaProfile = new TProfile *[14];
 486     fpJetAbsEtaProfile = new TProfile *[14];
 487     fpChargedJetRProfile = new TProfile *[8];
 488     fpJetRProfile = new TProfile *[4];
 489     fpChargedJetEDProfile= new TProfile3D *[10];
 490     fpJetEDProfile= new TProfile3D *[10];
 491     fvertex[0]=0.0,fvertex[1]=0.0,fvertex[2]=0.0;
 492
 493     DefineOutput(1,TList::Class());  // for output list
 494 }
 495
 496 //________________________________________________________________________
 497 AliAnalysisTaskFullpAJets::~AliAnalysisTaskFullpAJets()
 498 {
 499   // Destructor. Clean-up the output list, but not the histograms that are put inside
 500   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 501     if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())
 502     {
 503         delete fOutput;
 504     }
 505     delete fTrackCuts;
 506 }
 507
 508 //________________________________________________________________________
 509 void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
 510 {
 511     // Create histograms
 512     // Called once (on the worker node)
 513     fIsInitialized=kFALSE;
 514     fOutput = new TList();
 515     fOutput->SetOwner();  // IMPORTANT!
 516
 517     fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 518
 519     // Initialize Global Variables
 520     fnEvents=0;
 521     fnEventsCharged=0;
 522     fnDiJetEvents=0;
 523
 524     // fRJET=4 -> fJetR=0.4 && fRJET=25 -> fJetR=0.25, but for writing files, should be 4 and 25 respectively
 525     if (fRJET>10)
 526     {
 527         fJetR=(Double_t)fRJET/100.0;
 528     }
 529     else
 530     {
 531         fJetR=(Double_t)fRJET/10.0;
 532     }
 533
 534     fEMCalPhiMin=(80/(double)360)*2*TMath::Pi();
 535     fEMCalPhiMax=(187/(double)360)*2*TMath::Pi();
 536     fEMCalPhiTotal= fEMCalPhiMax-fEMCalPhiMin;
 537     fEMCalEtaMin=-0.7;
 538     fEMCalEtaMax=0.7;
 539     fEMCalEtaTotal=fEMCalEtaMax-fEMCalEtaMin;
 540     fEMCalArea=fEMCalPhiTotal*fEMCalEtaTotal;
 541
 542     fTPCPhiMin=(0/(double)360)*2*TMath::Pi();
 543     fTPCPhiMax=(360/(double)360)*2*TMath::Pi();
 544     fTPCPhiTotal= fTPCPhiMax-fTPCPhiMin;
 545     fTPCEtaMin=-0.9;
 546     fTPCEtaMax=0.9;
 547     fTPCEtaTotal=fTPCEtaMax-fTPCEtaMin;
 548     fTPCArea=fTPCPhiTotal*fTPCEtaTotal;
 549
 550     fCentralityBins=10;
 551     fCentralityLow=0.0;
 552     fCentralityUp=100.0;
 553     Int_t CentralityBinMult=10;
 554
 555     fJetAreaCutFrac =0.6;  // Fudge factor for selecting on jets with threshold Area or higher
 556     fJetAreaThreshold=fJetAreaCutFrac*TMath::Pi()*fJetR*fJetR;
 557     fTPCJetThreshold=5.0;  //  Threshold required for an Anti-kt jet to be considered a "true" jet
 558     fEMCalJetThreshold=5.0;
 559     fVertexWindow=10.0;
 560     fVertexMaxR=1.0;
 561
 562     fnBckgClusters=TMath::FloorNint(fEMCalArea/(TMath::Pi()*fJetR*fJetR));  // Select the number of RC that is as close as possible to the area of the EMCal.
 563     fRCBckgFluc = new Double_t[fnBckgClusters];
 564     for (Int_t i=0;i<fnBckgClusters;i++)
 565     {
 566         fRCBckgFluc[i]=0.0;
 567     }
 568
 569     fDeltaRho01=0.0;
 570     fnEMCalCells=12288;  // sMods 1-10 have 24x48 cells, sMods 11&12 have 8x48 cells...
 571
 572     // Histograms
 573     Int_t JetPtBins = 200;
 574     Double_t JetPtLow = 0.0;
 575     Double_t JetPtUp = 200.0;
 576
 577     Int_t TCBins=100;
 578
 579     // QA Plots
 580     fhTrackPt = new TH1D("fhTrackPt","p_{T} distribution of tracks in event",10*JetPtBins,JetPtLow,JetPtUp);
 581     fhTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 582     fhTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 583     fhTrackPt->Sumw2();
 584
 585     fhTrackPhi = new TH1D("fhTrackPhi","#phi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
 586     fhTrackPhi->GetXaxis()->SetTitle("#phi");
 587     fhTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#phi");
 588     fhTrackPhi->Sumw2();
 589
 590     fhTrackEta = new TH1D("fhTrackEta","#eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
 591     fhTrackEta->GetXaxis()->SetTitle("#eta");
 592     fhTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
 593     fhTrackEta->Sumw2();
 594
 595     fhTrackEtaPhi = new TH2D("fhTrackEtaPhi","#eta-#phi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
 596     fhTrackEtaPhi->GetXaxis()->SetTitle("#eta");
 597     fhTrackEtaPhi->GetYaxis()->SetTitle("#phi");
 598     fhTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
 599     fhTrackEtaPhi->Sumw2();
 600
 601     fhClusterPt = new TH1D("fhClusterPt","p_{T} distribution of clusters in event",10*JetPtBins,JetPtLow,JetPtUp);
 602     fhClusterPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 603     fhClusterPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 604     fhClusterPt->Sumw2();
 605
 606     fhClusterPhi = new TH1D("fhClusterPhi","#phi distribution of clusters in event",TCBins,fTPCPhiMin,fTPCPhiMax);
 607     fhClusterPhi->GetXaxis()->SetTitle("#phi");
 608     fhClusterPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#phi");
 609     fhClusterPhi->Sumw2();
 610
 611     fhClusterEta = new TH1D("fhClusterEta","#eta distribution of clusters in event",TCBins,fTPCEtaMin,fTPCEtaMax);
 612     fhClusterEta->GetXaxis()->SetTitle("#eta");
 613     fhClusterEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
 614     fhClusterEta->Sumw2();
 615
 616     fhClusterEtaPhi = new TH2D("fhClusterEtaPhi","#eta-#phi distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
 617     fhClusterEtaPhi->GetXaxis()->SetTitle("#eta");
 618     fhClusterEtaPhi->GetYaxis()->SetTitle("#phi");
 619     fhClusterEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
 620     fhClusterEtaPhi->Sumw2();
 621
 622     fhCentrality = new TH1D("fhCentrality","Event Centrality Distribution",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 623     fhCentrality->GetXaxis()->SetTitle("Centrality (V0M)");
 624     fhCentrality->GetYaxis()->SetTitle("1/N_{Events}");
 625     fhCentrality->Sumw2();
 626
 627     fhBckgFluc = new TH1D("fhBckgFluc",Form("p_{T} distribution of Background Clusters in near central events at center of EMCal with R=%g",fJetR),JetPtBins,JetPtLow,JetPtUp);
 628     fhBckgFluc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 629     fhBckgFluc->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 630     fhBckgFluc->Sumw2();
 631
 632     fhBckgMult = new TH1D("fhBckgMult",Form("Multiplicity distribution of Background Clusters in near central events at center of EMCal with R=%g",fJetR),JetPtBins,JetPtLow,JetPtUp);
 633     fhBckgMult->GetXaxis()->SetTitle("Multiplicity");
 634     fhBckgMult->GetYaxis()->SetTitle("1/N_{Events}");
 635     fhBckgMult->Sumw2();
 636
 637     fhJetConstituentPt= new TH2D("fhJetConstituentPt","Jet constituents p_{T} distribution",JetPtBins, JetPtLow, JetPtUp,10*JetPtBins, JetPtLow, JetPtUp);
 638     fhJetConstituentPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
 639     fhJetConstituentPt->GetYaxis()->SetTitle("Constituent p_{T} (GeV/c)");
 640     fhJetConstituentPt->Sumw2();
 641
 642     fhEMCalCellCounts = new TH1D("fhEMCalCellCounts","Distribtuion of cluster counts across the EMCal",fnEMCalCells,1,fnEMCalCells);
 643     fhEMCalCellCounts->GetXaxis()->SetTitle("Absoulute Cell Id");
 644     fhEMCalCellCounts->GetYaxis()->SetTitle("Counts per Event");
 645     fhEMCalCellCounts->Sumw2();
 646
 647     // Raw Jet Spectra
 648     fhChargedJetPt = new TH1D("fhChargedJetPt","Charged Jet p_{T} distribution for reconstructed Jets",JetPtBins, JetPtLow, JetPtUp);
 649     fhChargedJetPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 650     fhChargedJetPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 651     fhChargedJetPt->Sumw2();
 652
 653     fhChargedJetPtAreaCut = new TH1D("fhChargedJetPtAreaCut"," Charged Jet p_{T} distribution for reconstructed Jets with Standard Area Cut",JetPtBins, JetPtLow, JetPtUp);
 654     fhChargedJetPtAreaCut->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 655     fhChargedJetPtAreaCut->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 656     fhChargedJetPtAreaCut->Sumw2();
 657
 658     fhJetPtTPC = new TH1D("fhJetPtTPC","Jet p_{T} distribution for reconstructed Jets",JetPtBins, JetPtLow, JetPtUp);
 659     fhJetPtTPC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 660     fhJetPtTPC->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 661     fhJetPtTPC->Sumw2();
 662
 663     fhJetPtTPCAreaCut = new TH1D("fhJetPtTPCAreaCut","Jet p_{T} distribution for reconstructed Jets",JetPtBins, JetPtLow, JetPtUp);
 664     fhJetPtTPCAreaCut->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 665     fhJetPtTPCAreaCut->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 666     fhJetPtTPCAreaCut->Sumw2();
 667
 668     fhJetPtEMCal = new TH1D("fhJetPtEMCal","Jet p_{T} distribution for reconstructed Jets within the EMCal",JetPtBins, JetPtLow, JetPtUp);
 669     fhJetPtEMCal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 670     fhJetPtEMCal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 671     fhJetPtEMCal->Sumw2();
 672
 673     fhJetPtEMCalAreaCut = new TH1D("fhJetPtEMCalAreaCut","Jet p_{T} distribution for reconstructed Jets within the EMCal with Standard Area Cut",JetPtBins, JetPtLow, JetPtUp);
 674     fhJetPtEMCalAreaCut->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 675     fhJetPtEMCalAreaCut->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 676     fhJetPtEMCalAreaCut->Sumw2();
 677
 678     fhJetPtEMCalAreaCutSignal = new TH1D("fhJetPtEMCalAreaCutSignal","Jet p_{T} distribution for reconstructed Jets within the EMCal with Standard Area Cut and Signal Cut",JetPtBins, JetPtLow, JetPtUp);
 679     fhJetPtEMCalAreaCutSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 680     fhJetPtEMCalAreaCutSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 681     fhJetPtEMCalAreaCutSignal->Sumw2();
 682
 683     fh020JetPtEMCal = new TH1D("fh020JetPtEMCal","0-20% Jet p_{T} distribution for reconstructed Jets within the EMCal",JetPtBins, JetPtLow, JetPtUp);
 684     fh020JetPtEMCal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 685     fh020JetPtEMCal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 686     fh020JetPtEMCal->Sumw2();
 687
 688     fh020JetPtEMCalAreaCut = new TH1D("fh020JetPtEMCalAreaCut","0-20% Jet p_{T} distribution for reconstructed Jets within the EMCal with Standard Area Cut",JetPtBins, JetPtLow, JetPtUp);
 689     fh020JetPtEMCalAreaCut->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 690     fh020JetPtEMCalAreaCut->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 691     fh020JetPtEMCalAreaCut->Sumw2();
 692
 693     fh020JetPtEMCalAreaCutSignal = new TH1D("fh020JetPtEMCalAreaCutSignal","0-20% Jet p_{T} distribution for reconstructed Jets within the EMCal with Standard Area Cut and Signal Cut",JetPtBins, JetPtLow, JetPtUp);
 694     fh020JetPtEMCalAreaCutSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 695     fh020JetPtEMCalAreaCutSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 696     fh020JetPtEMCalAreaCutSignal->Sumw2();
 697
 698
 699     fhJetPtCenEMCal = new TH2D("fhJetPtCenEMCal","Jet p_{T} distribution for reconstructed Jets within the EMCal vs Centrality",JetPtBins, JetPtLow, JetPtUp, fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 700     fhJetPtCenEMCal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 701     fhJetPtCenEMCal->GetYaxis()->SetTitle("Centrality (V0M)");
 702     fhJetPtCenEMCal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 703     fhJetPtCenEMCal->Sumw2();
 704
 705     fhJetPtCenEMCalAreaCut = new TH2D("fhJetPtCenEMCalAreaCut","Jet p_{T} distribution for reconstructed Jets within the EMCal with Area Cut vs Centrality",JetPtBins, JetPtLow, JetPtUp, fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 706     fhJetPtCenEMCalAreaCut->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 707     fhJetPtCenEMCalAreaCut->GetYaxis()->SetTitle("Centrality (V0M)");
 708     fhJetPtCenEMCalAreaCut->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 709     fhJetPtCenEMCalAreaCut->Sumw2();
 710
 711     fhJetPtCenEMCalAreaCutSignal = new TH2D("fhJetPtCenEMCalAreaCutSignal","Jet p_{T} distribution for reconstructed Jets within the EMCal with Area and Signal Cut vs Centrality",JetPtBins, JetPtLow, JetPtUp, fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 712     fhJetPtCenEMCalAreaCutSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 713     fhJetPtCenEMCalAreaCutSignal->GetYaxis()->SetTitle("Centrality (V0M)");
 714     fhJetPtCenEMCalAreaCutSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 715     fhJetPtCenEMCalAreaCutSignal->Sumw2();
 716
 717     // Jet Area vs pT Distribution
 718     Int_t JetPtAreaBins=200;
 719     Double_t JetPtAreaLow=0.0;
 720     Double_t JetPtAreaUp=2.0;
 721
 722     fhJetPtArea = new TH2D("fhJetPtArea","Jet Area Distribution",JetPtBins, JetPtLow,JetPtUp,JetPtAreaBins,JetPtAreaLow,JetPtAreaUp);
 723     fhJetPtArea->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 724     fhJetPtArea->GetYaxis()->SetTitle("A_{jet}");
 725     fhJetPtArea->GetZaxis()->SetTitle("1/N_{Events} dN/dA_{jet}dp_{T}");
 726     fhJetPtArea->Sumw2();
 727
 728     Int_t A1_PtBins = 100;
 729     Double_t A1_PtLow=0.0;
 730     Double_t A1_PtUp=100.0;
 731     Int_t A1_TrigBins = 100;
 732     Double_t A1_TrigLow=0.0;
 733     Double_t A1_TrigUp=100.0;
 734     Int_t A1_RBins = 16;
 735     Double_t A1_RLow=0.4;
 736     Double_t A1_RUp=2.0;
 737
 738     fhJetTrigR1A = new TH3D("fhJetTrigR1A","Jet p_{T} distribution for reconstructed Jets vs Trigger Jet p_{T} vs #DeltaR",A1_PtBins, A1_PtLow, A1_PtUp,A1_TrigBins, A1_TrigLow, A1_TrigUp,A1_RBins, A1_RLow, A1_RUp);
 739     fhJetTrigR1A->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 740     fhJetTrigR1A->GetYaxis()->SetTitle("Trigger Jet p_{T} (GeV/c)");
 741     fhJetTrigR1A->GetZaxis()->SetTitle("R");
 742     fhJetTrigR1A->Sumw2();
 743
 744     // Corrected Jet Spectra
 745     Int_t JetTPtBins = 250;
 746     Double_t JetTPtLow=-50.0;
 747     Double_t JetTPtUp=200.0;
 748
 749     fhJetTPtRhoTotal = new TH1D("fhJetTPtRhoTotal","True Jet p_{T} distribution for reconstructed Jets in the EMCal",JetTPtBins, JetTPtLow, JetTPtUp);
 750     fhJetTPtRhoTotal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 751     fhJetTPtRhoTotal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 752     fhJetTPtRhoTotal->Sumw2();
 753
 754     fhJetTPtRhoTotalSignal = new TH1D("fhJetTPtRhoTotalSignal","True Jet p_{T} distribution for reconstructed Jets in the EMCal with Signal cut",JetTPtBins, JetTPtLow, JetTPtUp);
 755     fhJetTPtRhoTotalSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 756     fhJetTPtRhoTotalSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 757     fhJetTPtRhoTotalSignal->Sumw2();
 758
 759     fhJetTPtRhoNoLeading = new TH1D("fhJetTPtRhoNoLeading","True Jet p_{T} distribution for reconstructed Jets in the EMCal",JetTPtBins, JetTPtLow, JetTPtUp);
 760     fhJetTPtRhoNoLeading->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 761     fhJetTPtRhoNoLeading->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 762     fhJetTPtRhoNoLeading->Sumw2();
 763
 764     fhJetTPtRhoNoLeadingSignal = new TH1D("fhJetTPtRhoNoLeadingSignal","True Jet p_{T} distribution for reconstructed Jets in the EMCal with Signal cut",JetTPtBins, JetTPtLow, JetTPtUp);
 765     fhJetTPtRhoNoLeadingSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 766     fhJetTPtRhoNoLeadingSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 767     fhJetTPtRhoNoLeadingSignal->Sumw2();
 768
 769     fhJetTPt1B = new TH1D("fhJetTPt1B","True Jet p_{T} distribution for reconstructed Jets in the EMCal",JetTPtBins, JetTPtLow, JetTPtUp);
 770     fhJetTPt1B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 771     fhJetTPt1B->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 772     fhJetTPt1B->Sumw2();
 773
 774     fhJetTPt1BSignal = new TH1D("fhJetTPt1BSignal","True Jet p_{T} distribution for reconstructed Jets in the EMCal",JetTPtBins, JetTPtLow, JetTPtUp);
 775     fhJetTPt1BSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 776     fhJetTPt1BSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 777     fhJetTPt1BSignal->Sumw2();
 778
 779     fhJetTPt1C = new TH1D("fhJetTPt1C","True Jet p_{T} distribution for reconstructed Jets in the EMCal",JetPtBins,JetPtLow,JetPtUp);
 780     fhJetTPt1C->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 781     fhJetTPt1C->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 782     fhJetTPt1C->Sumw2();
 783
 784     fhJetTPt2B = new TH1D("fhJetTPt2B","True Jet p_{T} distribution for reconstructed Jets in the EMCal with dijet Trigger in TPC",JetTPtBins, JetTPtLow, JetTPtUp);
 785     fhJetTPt2B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 786     fhJetTPt2B->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 787     fhJetTPt2B->Sumw2();
 788
 789     fhJetTPt3 = new TH1D("fhJetTPt3","True Charged jet p_{T} distribution for reconstructed Jets in the TPC",JetTPtBins, JetTPtLow, JetTPtUp);
 790     fhJetTPt3->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 791     fhJetTPt3->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 792     fhJetTPt3->Sumw2();
 793
 794     // 0-20% Centrality Corrected Spectra (fh020JetTPt...)
 795     fh020JetTPtRhoTotal = new TH1D("fh020JetTPtRhoTotal","0-20% True Jet p_{T} distribution for reconstructed Jets in the EMCal",JetTPtBins, JetTPtLow, JetTPtUp);
 796     fh020JetTPtRhoTotal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 797     fh020JetTPtRhoTotal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 798     fh020JetTPtRhoTotal->Sumw2();
 799
 800     fh020JetTPtRhoTotalSignal = new TH1D("fh020JetTPtRhoTotalSignal","0-20% True Jet p_{T} distribution for reconstructed Jets in the EMCal with Signal cut",JetTPtBins, JetTPtLow, JetTPtUp);
 801     fh020JetTPtRhoTotalSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 802     fh020JetTPtRhoTotalSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 803     fh020JetTPtRhoTotalSignal->Sumw2();
 804
 805     fh020JetTPtRhoNoLeading = new TH1D("fh020JetTPtRhoNoLeading","0-20% True Jet p_{T} distribution for reconstructed Jets in the EMCal",JetTPtBins, JetTPtLow, JetTPtUp);
 806     fh020JetTPtRhoNoLeading->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 807     fh020JetTPtRhoNoLeading->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 808     fh020JetTPtRhoNoLeading->Sumw2();
 809
 810     fh020JetTPtRhoNoLeadingSignal = new TH1D("fh020JetTPtRhoNoLeadingSignal","0-20% True Jet p_{T} distribution for reconstructed Jets in the EMCal with Signal cut",JetTPtBins, JetTPtLow, JetTPtUp);
 811     fh020JetTPtRhoNoLeadingSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 812     fh020JetTPtRhoNoLeadingSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 813     fh020JetTPtRhoNoLeadingSignal->Sumw2();
 814
 815     fh020JetTPt1B = new TH1D("fh020JetTPt1B","0-20% True Jet p_{T} distribution for reconstructed Jets in the EMCal",JetTPtBins, JetTPtLow, JetTPtUp);
 816     fh020JetTPt1B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 817     fh020JetTPt1B->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 818     fh020JetTPt1B->Sumw2();
 819
 820     fh020JetTPt1BSignal = new TH1D("fh020JetTPt1BSignal","0-20% True Jet p_{T} distribution for reconstructed Jets in the EMCal",JetTPtBins, JetTPtLow, JetTPtUp);
 821     fh020JetTPt1BSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 822     fh020JetTPt1BSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 823     fh020JetTPt1BSignal->Sumw2();
 824
 825     fh020JetTPt1C = new TH1D("fh020JetTPt1C","0-20% True Jet p_{T} distribution for reconstructed Jets in the EMCal",JetPtBins,JetPtLow,JetPtUp);
 826     fh020JetTPt1C->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 827     fh020JetTPt1C->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 828     fh020JetTPt1C->Sumw2();
 829
 830     fh020JetTPt2B = new TH1D("fh020JetTPt2B","0-20% True Jet p_{T} distribution for reconstructed Jets in the EMCal with dijet Trigger in TPC",JetTPtBins, JetTPtLow, JetTPtUp);
 831     fh020JetTPt2B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 832     fh020JetTPt2B->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 833     fh020JetTPt2B->Sumw2();
 834
 835     fh020JetTPt3 = new TH1D("fh020JetTPt3","0-20% True Charged jet p_{T} distribution for reconstructed Jets in the TPC",JetTPtBins, JetTPtLow, JetTPtUp);
 836     fh020JetTPt3->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 837     fh020JetTPt3->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 838     fh020JetTPt3->Sumw2();
 839
 840     // 2D Corrected Spectra (fhJetTPtCen...)
 841     fhJetTPtCenRhoTotal = new TH2D("fhJetTPtCenRhoTotal","True Jet p_{T} distribution for reconstructed Jets in the EMCal vs Centrality",JetTPtBins, JetTPtLow, JetTPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 842     fhJetTPtCenRhoTotal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 843     fhJetTPtCenRhoTotal->GetYaxis()->SetTitle("Centrality (V0M)");
 844     fhJetTPtCenRhoTotal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 845     fhJetTPtCenRhoTotal->Sumw2();
 846
 847     fhJetTPtCenRhoTotalSignal = new TH2D("fhJetTPtCenRhoTotalSignal","True Jet p_{T} distribution for reconstructed Jets in the EMCal with Signal cut vs Centrality",JetTPtBins, JetTPtLow, JetTPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 848     fhJetTPtCenRhoTotalSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 849     fhJetTPtCenRhoTotalSignal->GetYaxis()->SetTitle("Centrality (V0M)");
 850     fhJetTPtCenRhoTotalSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 851     fhJetTPtCenRhoTotalSignal->Sumw2();
 852
 853     fhJetTPtCenRhoNoLeading = new TH2D("fhJetTPtCenRhoNoLeading","True Jet p_{T} distribution for reconstructed Jets in the EMCal vs Centrality",JetTPtBins, JetTPtLow, JetTPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 854     fhJetTPtCenRhoNoLeading->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 855     fhJetTPtCenRhoNoLeading->GetYaxis()->SetTitle("Centrality (V0M)");
 856     fhJetTPtCenRhoNoLeading->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 857     fhJetTPtCenRhoNoLeading->Sumw2();
 858
 859     fhJetTPtCenRhoNoLeadingSignal = new TH2D("fhJetTPtCenRhoNoLeadingSignal","True Jet p_{T} distribution for reconstructed Jets in the EMCal with Signal cut vs Centrality",JetTPtBins, JetTPtLow, JetTPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 860     fhJetTPtCenRhoNoLeadingSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 861     fhJetTPtCenRhoNoLeadingSignal->GetYaxis()->SetTitle("Centrality (V0M)");
 862     fhJetTPtCenRhoNoLeadingSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 863     fhJetTPtCenRhoNoLeadingSignal->Sumw2();
 864
 865     fhJetTPtCen1B = new TH2D("fhJetTPtCen1B","True Jet p_{T} distribution for reconstructed Jets in the EMCal vs Centrality",JetTPtBins, JetTPtLow, JetTPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 866     fhJetTPtCen1B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 867     fhJetTPtCen1B->GetYaxis()->SetTitle("Centrality (V0M)");
 868     fhJetTPtCen1B->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 869     fhJetTPtCen1B->Sumw2();
 870
 871     fhJetTPtCen1BSignal = new TH2D("fhJetTPtCen1BSignal","True Jet p_{T} distribution for reconstructed Jets in the EMCal vs Centrality",JetTPtBins, JetTPtLow, JetTPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 872     fhJetTPtCen1BSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 873     fhJetTPtCen1BSignal->GetYaxis()->SetTitle("Centrality (V0M)");
 874     fhJetTPtCen1BSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 875     fhJetTPtCen1BSignal->Sumw2();
 876
 877     fhJetTPtCen1C = new TH2D("fhJetTPtCen1C","True Jet p_{T} distribution for reconstructed Jets in the EMCal vs Centrality",JetPtBins,JetPtLow,JetPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 878     fhJetTPtCen1C->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 879     fhJetTPtCen1C->GetYaxis()->SetTitle("Centrality (V0M)");
 880     fhJetTPtCen1C->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 881     fhJetTPtCen1C->Sumw2();
 882
 883     fhJetTPtCen2B = new TH2D("fhJetTPtCen2B","True Jet p_{T} distribution for reconstructed Jets in the EMCal with dijet Trigger in TPC vs Centrality",JetTPtBins, JetTPtLow, JetTPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 884     fhJetTPtCen2B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 885     fhJetTPtCen2B->GetYaxis()->SetTitle("Centrality (V0M)");
 886     fhJetTPtCen2B->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 887     fhJetTPtCen2B->Sumw2();
 888
 889     fhJetTPtCen3 = new TH2D("fhJetTPtCen3","True Charged jet p_{T} distribution for reconstructed Jets in the TPC vs Centrality",JetTPtBins, JetTPtLow, JetTPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 890     fhJetTPtCen3->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 891     fhJetTPtCen3->GetYaxis()->SetTitle("Centrality (V0M)");
 892     fhJetTPtCen3->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 893     fhJetTPtCen3->Sumw2();
 894
 895     // Cluster Plots
 896     fhEMCalBckg1B = new TH1D("fhEMCalBckg1B","Cluster p_{T} distribution for reconstructed Tracks and Calocluster at least 0.5 away from all jets in the EMCal",JetPtBins, JetPtLow,JetPtUp);
 897     fhEMCalBckg1B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 898     fhEMCalBckg1B->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 899     fhEMCalBckg1B->Sumw2();
 900
 901     fhEMCalBckg1C = new TH1D("fhEMCalBckg1C","Cluster p_{T} distribution for reconstructed Tracks and Calocluster in Transverse area with R=0.4",JetPtBins, JetPtLow, JetPtUp);
 902     fhEMCalBckg1C->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 903     fhEMCalBckg1C->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 904     fhEMCalBckg1C->Sumw2();
 905
 906     fhEMCalJet2A = new TH1D("fhEMCalJet2A","Cluster p_{T} distribution for jets within EMCal from di-jets in TPC",A1_PtBins,A1_PtLow,A1_PtUp);
 907     fhEMCalJet2A->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 908     fhEMCalJet2A->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 909     fhEMCalJet2A->Sumw2();
 910
 911     fh020EMCalJet2A = new TH1D("fh020EMCalJet2A","0-20% Centrality, Cluster p_{T} distribution for jets within EMCal from di-jets in TPC",A1_PtBins,A1_PtLow,A1_PtUp);
 912     fh020EMCalJet2A->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 913     fh020EMCalJet2A->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 914     fh020EMCalJet2A->Sumw2();
 915
 916     fhEMCalCenJet2A = new TH2D("fhEMCalCenJet2A","Cluster p_{T} distribution for jets within EMCal from di-jets in TPC",A1_PtBins,A1_PtLow,A1_PtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 917     fhEMCalCenJet2A->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 918     fhEMCalCenJet2A->GetYaxis()->SetTitle("Centrality (V0M)");
 919     fhEMCalCenJet2A->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 920     fhEMCalCenJet2A->Sumw2();
 921
 922     fhEMCalBckg2B = new TH1D("fhEMCalBckg2B","Cluster p_{T} distribution for reconstructed Tracks and Calocluster with dijet Trigger in TPC with R=0.4",JetPtBins, JetPtLow, JetPtUp);
 923     fhEMCalBckg2B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 924     fhEMCalBckg2B->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 925     fhEMCalBckg2B->Sumw2();
 926
 927
 928     Int_t EMCalClusterBins=100;
 929     Int_t EMCalClusterLow=0;
 930     Int_t EMCalClusterUp=100;
 931
 932     fh020EMCalkTClusters = new TH1D("fh020EMCalkTClusters","0-20 % Centrality, Number of k_{T} clusters per event",EMCalClusterBins,EMCalClusterLow,EMCalClusterUp);
 933     fh020EMCalkTClusters->GetXaxis()->SetTitle("# of Clusters");
 934     fh020EMCalkTClusters->GetYaxis()->SetTitle("1/N_{Events}");
 935     fh020EMCalkTClusters->Sumw2();
 936
 937     fh020EMCalAkTJets = new TH1D("fh020EMCalAkTJets","0-20 % Centrality, Number of anti-k_{T} jets per event",EMCalClusterBins,EMCalClusterLow,EMCalClusterUp);
 938     fh020EMCalAkTJets->GetXaxis()->SetTitle("# of jets");
 939     fh020EMCalAkTJets->GetYaxis()->SetTitle("1/N_{Events}");
 940     fh020EMCalAkTJets->Sumw2();
 941
 942     // Background Density Plots
 943     // 2D Rho plots (fhRho...)
 944     Int_t RhoPtBins = 500;
 945     Double_t RhoPtLow=0.0;
 946     Double_t RhoPtUp=50.0;
 947
 948     Int_t DeltaRhoPtBins=100;
 949     Double_t DeltaRhoPtLow=-5.0;
 950     Double_t DeltaRhoPtUp=5.0;
 951
 952     fhRhoTotal= new TH2D("fhRhoTotal","Background Density #rho_{0}",RhoPtBins,RhoPtLow,RhoPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 953     fhRhoTotal->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
 954     fhRhoTotal->GetYaxis()->SetTitle("Centrality (V0M)");
 955     fhRhoTotal->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
 956     fhRhoTotal->Sumw2();
 957
 958     fhRhoNoLeading= new TH2D("fhRhoNoLeading","Background Density #rho_{1}",RhoPtBins,RhoPtLow,RhoPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 959     fhRhoNoLeading->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
 960     fhRhoNoLeading->GetYaxis()->SetTitle("Centrality (V0M)");
 961     fhRhoNoLeading->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
 962     fhRhoNoLeading->Sumw2();
 963
 964     fhRho1B = new TH2D("fhRho1B","Background Density #rho_{n} ",RhoPtBins, RhoPtLow, RhoPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 965     fhRho1B->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
 966     fhRho1B->GetYaxis()->SetTitle("Centrality (V0M)");
 967     fhRho1B->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
 968     fhRho1B->Sumw2();
 969
 970     fhRho1C = new TH2D("fhRho1C","Background Density #rho (Method 1C)",RhoPtBins, RhoPtLow, RhoPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 971     fhRho1C->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
 972     fhRho1C->GetYaxis()->SetTitle("Centrality (V0M)");
 973     fhRho1C->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
 974     fhRho1C->Sumw2();
 975
 976     fhRho2B = new TH2D("fhRho2B","Background Density #rho_{dijet}",RhoPtBins, RhoPtLow, RhoPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 977     fhRho2B->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
 978     fhRho2B->GetYaxis()->SetTitle("Centrality (V0M)");
 979     fhRho2B->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
 980     fhRho2B->Sumw2();
 981
 982     fhRho2BCore = new TH2D("fhRho2BCore","Background Density #rho_{dijet}",RhoPtBins, RhoPtLow, RhoPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 983     fhRho2BCore->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
 984     fhRho2BCore->GetYaxis()->SetTitle("Centrality (V0M)");
 985     fhRho2BCore->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
 986     fhRho2BCore->Sumw2();
 987
 988     fhRho3 = new TH2D("fhRho3","Charged Background Density #rho_{char}",RhoPtBins, RhoPtLow, RhoPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 989     fhRho3->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
 990     fhRho3->GetYaxis()->SetTitle("Centrality (V0M)");
 991     fhRho3->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
 992     fhRho3->Sumw2();
 993
 994     fhRho3NoJets = new TH2D("fhRho3NoJets","Charged Background Density #rho_{char} for Events with No Signal Jets",RhoPtBins, RhoPtLow, RhoPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 995     fhRho3NoJets->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
 996     fhRho3NoJets->GetYaxis()->SetTitle("Centrality (V0M)");
 997     fhRho3NoJets->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
 998     fhRho3NoJets->Sumw2();
 999
1000     fhRho3DiJets = new TH2D("fhRho3DiJets","Charged Background Density #rho_{char} for Events with a DiJet",RhoPtBins, RhoPtLow, RhoPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
1001     fhRho3DiJets->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1002     fhRho3DiJets->GetYaxis()->SetTitle("Centrality (V0M)");
1003     fhRho3DiJets->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
1004     fhRho3DiJets->Sumw2();
1005
1006     fhRho3Perp = new TH2D("fhRho3Perp","Charged Background Density #rho_{char} for Events with a DiJet at Median Angle Between Dijets",RhoPtBins, RhoPtLow, RhoPtUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
1007     fhRho3Perp->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1008     fhRho3Perp->GetYaxis()->SetTitle("Centrality (V0M)");
1009     fhRho3Perp->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
1010     fhRho3Perp->Sumw2();
1011
1012     // 0-20% Centrality Plots (rh020Rho...)
1013     fh020RhoTotal = new TH1D("fh020RhoTotal","0-20% Background Density #rho_{0}",RhoPtBins,RhoPtLow,RhoPtUp);
1014     fh020RhoTotal->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1015     fh020RhoTotal->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
1016     fh020RhoTotal->Sumw2();
1017
1018     fh020RhoNoLeading = new TH1D("fh020RhoNoLeading","0-20% Background Density #rho_{1}",RhoPtBins,RhoPtLow,RhoPtUp);
1019     fh020RhoNoLeading->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1020     fh020RhoNoLeading->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
1021     fh020RhoNoLeading->Sumw2();
1022
1023     fh020Rho1B = new TH1D("fh020Rho1B","0-20% Background Density #rho_{n}",RhoPtBins,RhoPtLow,RhoPtUp);
1024     fh020Rho1B->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1025     fh020Rho1B->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
1026     fh020Rho1B->Sumw2();
1027
1028     fh020Rho2B = new TH1D("fh020Rho2B","0-20% Background Density #rho_{dijet}",RhoPtBins,RhoPtLow,RhoPtUp);
1029     fh020Rho2B->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1030     fh020Rho2B->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
1031     fh020Rho2B->Sumw2();
1032
1033     fh020Rho2BCore = new TH1D("fh020Rho2BCore","0-20% Background Density #rho_{dijet}",RhoPtBins,RhoPtLow,RhoPtUp);
1034     fh020Rho2BCore->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1035     fh020Rho2BCore->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
1036     fh020Rho2BCore->Sumw2();
1037
1038     fh020Rho3 = new TH1D("fh020Rho3","0-20% Charged Background Density #rho_{char}",RhoPtBins,RhoPtLow,RhoPtUp);
1039     fh020Rho3->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1040     fh020Rho3->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
1041     fh020Rho3->Sumw2();
1042
1043     fh020Rho3NoJets = new TH1D("fh020Rho3NoJets","0-20% Charged Background Density #rho_{char} for Events with No Signal Jets",RhoPtBins,RhoPtLow,RhoPtUp);
1044     fh020Rho3NoJets->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1045     fh020Rho3NoJets->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
1046     fh020Rho3NoJets->Sumw2();
1047
1048     fh020Rho3DiJets = new TH1D("fh020Rho3DiJets","0-20% Charged Background Density #rho_{char} for Events with a DiJet",RhoPtBins,RhoPtLow,RhoPtUp);
1049     fh020Rho3DiJets->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1050     fh020Rho3DiJets->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
1051     fh020Rho3DiJets->Sumw2();
1052
1053     fh020Rho3Perp = new TH1D("fh020Rho3Perp","0-20% Charged Background Density #rho_{char} for Events with a DiJet at Median Angle Between Dijets",RhoPtBins,RhoPtLow,RhoPtUp);
1054     fh020Rho3Perp->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1055     fh020Rho3Perp->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
1056     fh020Rho3Perp->Sumw2();
1057
1058     fh020RhokT = new TH1D("fh020RhokT","0-20% Background Density #rho_{k_{T}}",RhoPtBins,RhoPtLow,RhoPtUp);
1059     fh020RhokT->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
1060     fh020RhokT->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
1061     fh020RhokT->Sumw2();
1062
1063     // Delta pT Plots
1064     Int_t DeltaPtBins=150;
1065     Double_t DeltaPtLow=-50.0;
1066     Double_t DeltaPtUp=100.0;
1067
1068     fhDeltaPtTotal = new TH1D("fhDeltaPtTotal","#deltap_{T} distribution for Random Cones using total #rho",DeltaPtBins,DeltaPtLow,DeltaPtUp);
1069     fhDeltaPtTotal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
1070     fhDeltaPtTotal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
1071     fhDeltaPtTotal->Sumw2();
1072
1073     fhDeltaPtNoLeading = new TH1D("fhDeltaPtNoLeading","#deltap_{T} distribution for Random Cones using leading jet #rho",DeltaPtBins,DeltaPtLow,DeltaPtUp);
1074     fhDeltaPtNoLeading->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
1075     fhDeltaPtNoLeading->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
1076     fhDeltaPtNoLeading->Sumw2();
1077
1078     fhDeltaPt1B = new TH1D("fhDeltaPt1B","#deltap_{T} distribution for Random Cones using method 1B #rho",DeltaPtBins,DeltaPtLow,DeltaPtUp);
1079     fhDeltaPt1B->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
1080     fhDeltaPt1B->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
1081     fhDeltaPt1B->Sumw2();
1082
1083     fhDeltaPt2B = new TH1D("fhDeltaPt2B","#deltap_{T} distribution for Random Cones using method 2B #rho",DeltaPtBins,DeltaPtLow,DeltaPtUp);
1084     fhDeltaPt2B->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
1085     fhDeltaPt2B->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
1086     fhDeltaPt2B->Sumw2();
1087
1088     fhDeltaPtkT = new TH1D("fhDeltaPtkT","#deltap_{T} distribution for Random Cones using kT #rho",DeltaPtBins,DeltaPtLow,DeltaPtUp);
1089     fhDeltaPtkT->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
1090     fhDeltaPtkT->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
1091     fhDeltaPtkT->Sumw2();
1092
1093     fhDeltaRho01 = new TH1D("fhDeltaRho01","Event by Event Differential between #rho_{0} and #rho_{1}",RhoPtBins,RhoPtLow,RhoPtUp);
1094     fhDeltaRho01->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1095     fhDeltaRho01->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
1096     fhDeltaRho01->Sumw2();
1097
1098     fhDeltaRho0DiJet = new TH1D("fhDeltaRho0DiJet","Event by Event Differential between #rho_{0} and #rho_{dijet}",DeltaRhoPtBins,DeltaRhoPtLow,DeltaRhoPtUp);
1099     fhDeltaRho0DiJet->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1100     fhDeltaRho0DiJet->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
1101     fhDeltaRho0DiJet->Sumw2();
1102
1103     // DiJet Plots
1104     Int_t DiJetBins=20;
1105     Double_t DiJetLow=0.0;
1106     Double_t DiJetUp=1.0;
1107
1108     fh020DiJetAsy = new TH1D("fh020DiJetAsy","0-20% Centrality, Di-Jet Asymmetry",DiJetBins,DiJetLow,DiJetUp);
1109     fh020DiJetAsy->GetXaxis()->SetTitle("A_{j}");
1110     fh020DiJetAsy->GetYaxis()->SetTitle("1/N_{Events} dN/dA_{j}");
1111     fh020DiJetAsy->Sumw2();
1112
1113     fhDiJetCenAsy = new TH2D("fhDiJetCenAsy","Di-Jet Asymmetry",DiJetBins,DiJetLow,DiJetUp,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
1114     fhDiJetCenAsy->GetXaxis()->SetTitle("A_{j}");
1115     fhDiJetCenAsy->GetYaxis()->SetTitle("Centrality (V0M)");
1116     fhDiJetCenAsy->GetZaxis()->SetTitle("1/N_{Events} dN/dA_{j}");
1117     fhDiJetCenAsy->Sumw2();
1118
1119     fh020DiJetDeltaPhi = new TH1D("fh020DiJetDeltaPhi","0-20% Centrality, Di-Jet #Delta#phi",TCBins,fTPCPhiMin,0.5*fTPCPhiMax);
1120     fh020DiJetDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
1121     fh020DiJetDeltaPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#Delta#phi");
1122     fh020DiJetDeltaPhi->Sumw2();
1123
1124     fhDiJetCenDeltaPhi = new TH2D("fhDiJetCenDeltaPhi","Di-Jet #Delta#phi",TCBins,fTPCPhiMin,0.5*fTPCPhiMax,fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
1125     fhDiJetCenDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
1126     fhDiJetCenDeltaPhi->GetYaxis()->SetTitle("Centrality (V0M)");
1127     fhDiJetCenDeltaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#Delta#phi");
1128     fhDiJetCenDeltaPhi->Sumw2();
1129
1130     fhDiJetEMCalLeadingPt = new TH1D("fhDiJetEMCalLeadingPt","Leading EMCal anti-k_{T} cluster p_{T} in a dijet event",EMCalClusterBins,EMCalClusterLow,EMCalClusterUp);
1131     fhDiJetEMCalLeadingPt->GetXaxis()->SetTitle("p_{T}");
1132     fhDiJetEMCalLeadingPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
1133     fhDiJetEMCalLeadingPt->Sumw2();
1134
1135     fhDiJetEMCalLeadingDeltaPhi = new TH1D("fhDiJetEMCalLeadingDeltaPhi","Leading EMCal anti-k_{T} cluster #Delta#phi in a dijet event",TCBins,fTPCPhiMin,0.5*fTPCPhiMax);
1136     fhDiJetEMCalLeadingDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
1137     fhDiJetEMCalLeadingDeltaPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#Delta#phi");
1138     fhDiJetEMCalLeadingDeltaPhi->Sumw2();
1139
1140     // Profiles
1141     // Background Density vs Centrality (fpRho...)
1142     fpRhoTotal = new TProfile("fpRhoTotal","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1143     fpRhoTotal->GetXaxis()->SetTitle("Centrality (V0M)");
1144     fpRhoTotal->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1145
1146     fpRhoNoLeading = new TProfile("fpRhoNoLeading","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1147     fpRhoNoLeading->GetXaxis()->SetTitle("Centrality (V0M)");
1148     fpRhoNoLeading->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1149
1150     fpRho1B = new TProfile("fpRho1B","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1151     fpRho1B->GetXaxis()->SetTitle("Centrality (V0M)");
1152     fpRho1B->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1153
1154     fpRho2B = new TProfile("fpRho2B","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1155     fpRho2B->GetXaxis()->SetTitle("Centrality (V0M)");
1156     fpRho2B->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1157
1158     fpRho2BCore = new TProfile("fpRho2BCore","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1159     fpRho2BCore->GetXaxis()->SetTitle("Centrality (V0M)");
1160     fpRho2BCore->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1161
1162     fpRho3 = new TProfile("fpRho3","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1163     fpRho3->GetXaxis()->SetTitle("Centrality (V0M)");
1164     fpRho3->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1165
1166     fpRho3NoJets = new TProfile("fpRho3NoJets","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1167     fpRho3NoJets->GetXaxis()->SetTitle("Centrality (V0M)");
1168     fpRho3NoJets->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1169
1170     fpRho3DiJets = new TProfile("fpRho3DiJets","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1171     fpRho3DiJets->GetXaxis()->SetTitle("Centrality (V0M)");
1172     fpRho3DiJets->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1173
1174     fpRho3Perp = new TProfile("fpRho3Perp","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1175     fpRho3Perp->GetXaxis()->SetTitle("Centrality (V0M)");
1176     fpRho3Perp->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1177
1178     fpRhoScale = new TProfile("fpRhoScale","Scaling Factor Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1179     fpRhoScale->GetXaxis()->SetTitle("Centrality (V0M)");
1180     fpRhoScale->GetYaxis()->SetTitle("Scale Factor");
1181
1182     fpRhokT = new TProfile("fpRhokT","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1183     fpRhokT->GetXaxis()->SetTitle("Centrality (V0M)");
1184     fpRhokT->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1185
1186     fpRhoChargedkT = new TProfile("fpRhoChargedkT","Background Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1187     fpRhoChargedkT->GetXaxis()->SetTitle("Centrality (V0M)");
1188     fpRhoChargedkT->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1189
1190     fpRhoScalekT = new TProfile("fpRhoScalekT","Scaling Factor Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1191     fpRhoScalekT->GetXaxis()->SetTitle("Centrality (V0M)");
1192     fpRhoScalekT->GetYaxis()->SetTitle("Scale Factor");
1193
1194     // Others
1195     fpEventMult = new TProfile("fpEventMult","Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
1196     fpEventMult->GetXaxis()->SetTitle("Centrality (V0M)");
1197     fpEventMult->GetYaxis()->SetTitle("Multiplicity");
1198
1199     fpJetPtRhoTotal = new TProfile("fpJetPtRhoTotal","#rho_{0} Profile vs Leading jet p_{T}",JetPtBins,JetPtLow,JetPtUp);
1200     fpJetPtRhoTotal->GetXaxis()->SetTitle("Leading jet p_{T} (GeV/c)");
1201     fpJetPtRhoTotal->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1202
1203     fpJetPtRhoNoLeading = new TProfile("fpJetPtRhoNoLeading","#rho_{1} Profile vs Leading jet p_{T}",JetPtBins,JetPtLow,JetPtUp);
1204     fpJetPtRhoNoLeading->GetXaxis()->SetTitle("Leading jet p_{T} (GeV/c)");
1205     fpJetPtRhoNoLeading->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1206
1207     fpJetPtRhokT = new TProfile("fpJetPtRhokT","#rho_{k_{T}} Profile vs Leading jet p_{T}",JetPtBins,JetPtLow,JetPtUp);
1208     fpJetPtRhokT->GetXaxis()->SetTitle("Leading jet p_{T} (GeV/c)");
1209     fpJetPtRhokT->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1210
1211     // QA::2D Energy Density Profiles for Tracks and CLusters
1212     fpTrackPtProfile = new TProfile2D("fpTrackPtProfile","2D Profile of track pT density throughout the TPC",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
1213     fpTrackPtProfile->GetXaxis()->SetTitle("#eta");
1214     fpTrackPtProfile->GetYaxis()->SetTitle("#phi");
1215     fpTrackPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
1216
1217     fpClusterPtProfile = new TProfile2D("fpClusterPtProfile","2D Profile of cluster pT density throughout the EMCal",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
1218     fpClusterPtProfile->GetXaxis()->SetTitle("#eta");
1219     fpClusterPtProfile->GetYaxis()->SetTitle("#phi");
1220     fpClusterPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
1221
1222     TString temp_name="";
1223     TString title_name="";
1224
1225     fEtaProfileBins=14;
1226     fEtaProfileLow=-0.7;
1227     fEtaProfileUp=0.7;
1228
1229     for (Int_t i=0;i<14;i++)
1230     {
1231         temp_name=Form("fpJetEtaProfile%d",i);
1232         title_name=Form("Jet Energy Density #eta Profile for ALL p_{T}, 0-20 Centrality, and eta=%g to %g",(i-7)/10.,(i-6)/10.);
1233
1234         fpJetEtaProfile[i]= new TProfile(temp_name,title_name,fEtaProfileBins,fEtaProfileLow,fEtaProfileUp);
1235         fpJetEtaProfile[i]->GetXaxis()->SetTitle("#eta");
1236         fpJetEtaProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1237         fOutput->Add(fpJetEtaProfile[i]);
1238         temp_name="";
1239         title_name="";
1240
1241         temp_name=Form("fpJetAbsEtaProfile%d",i);
1242         title_name=Form("Jet Energy Density #Delta #eta Profile for ALL p_{T}, 0-20 Centrality, and eta=%g to %g",(i-7)/10.,(i-6)/10.);
1243
1244         fpJetAbsEtaProfile[i]= new TProfile(temp_name,title_name,fEtaProfileBins,0,2*fEtaProfileUp);
1245         fpJetAbsEtaProfile[i]->GetXaxis()->SetTitle("#Delta#eta");
1246         fpJetAbsEtaProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1247         fOutput->Add(fpJetAbsEtaProfile[i]);
1248         temp_name="";
1249         title_name="";
1250 }
1251
1252     fEDProfileRBins=50;
1253     fEDProfileRLow=0.0;
1254     fEDProfileRUp=0.5;
1255     fEDProfilePtBins=100;
1256     fEDProfilePtLow=0.0;
1257     fEDProfilePtUp=100.0;
1258     fEDProfileEtaBins=4;
1259     fEDProfileEtaLow=-0.2;
1260     fEDProfileEtaUp=0.2;
1261
1262     for (Int_t i=0;i<8;i++)
1263     {
1264         temp_name=Form("fpChargedJetRProfile%d",i);
1265         title_name=Form("Charged Jet Energy Density Radius Profile for ALL p_{T}, 0-20 Centrality, and eta=%g to %g",(i-4)/10.,(i-3)/10.);
1266
1267         fpChargedJetRProfile[i]= new TProfile(temp_name,title_name,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
1268         fpChargedJetRProfile[i]->GetXaxis()->SetTitle("Radius");
1269         fpChargedJetRProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1270         fOutput->Add(fpChargedJetRProfile[i]);
1271         temp_name="";
1272         title_name="";
1273     }
1274
1275     for (Int_t i=0;i<4;i++)
1276     {
1277         temp_name=Form("fpJetRProfile%d",i);
1278         title_name=Form("Jet Energy Density Radius Profile for ALL p_{T}, 0-20 Centrality, and eta=%g to %g",(i-2)/10.,(i-1)/10.);
1279
1280         fpJetRProfile[i]= new TProfile(temp_name,title_name,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
1281         fpJetRProfile[i]->GetXaxis()->SetTitle("Radius");
1282         fpJetRProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
1283         fOutput->Add(fpJetRProfile[i]);
1284         temp_name="";
1285         title_name="";
1286     }
1287
1288     for (Int_t i=0;i<10;i++)
1289     {
1290         temp_name=Form("fpChargedJetEDProfile%d0",i);
1291         title_name=Form("Charged Jet Energy Density Profile for %d0-%d0 Centrality",i,i+1);
1292
1293         fpChargedJetEDProfile[i]= new TProfile3D(temp_name,title_name,fEDProfilePtBins,fEDProfilePtLow,fEDProfilePtUp,fEDProfileEtaBins+4,fEDProfileEtaLow-0.2,fEDProfileEtaUp+0.2,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
1294         fpChargedJetEDProfile[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1295         fpChargedJetEDProfile[i]->GetYaxis()->SetTitle("Pseudorapidity");
1296         fpChargedJetEDProfile[i]->GetZaxis()->SetTitle("Radius");
1297         fOutput->Add(fpChargedJetEDProfile[i]);
1298         temp_name="";
1299         title_name="";
1300
1301         temp_name=Form("fpJetEDProfile%d0",i);
1302         title_name=Form("Jet Energy Density Profile for %d0-%d0 Centrality",i,i+1);
1303
1304         fpJetEDProfile[i]= new TProfile3D(temp_name,title_name,fEDProfilePtBins,fEDProfilePtLow,fEDProfilePtUp,fEDProfileEtaBins,fEDProfileEtaLow,fEDProfileEtaUp,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
1305         fpJetEDProfile[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1306         fpJetEDProfile[i]->GetYaxis()->SetTitle("Pseudorapidity");
1307         fpJetEDProfile[i]->GetZaxis()->SetTitle("Radius");
1308         fOutput->Add(fpJetEDProfile[i]);
1309         temp_name="";
1310         title_name="";
1311     }
1312
1313     fOutput->Add(fhTrackPt);
1314     fOutput->Add(fhTrackEta);
1315     fOutput->Add(fhTrackPhi);
1316     fOutput->Add(fhTrackEtaPhi);
1317     fOutput->Add(fhClusterPt);
1318     fOutput->Add(fhClusterEta);
1319     fOutput->Add(fhClusterPhi);
1320     fOutput->Add(fhClusterEtaPhi);
1321     fOutput->Add(fhCentrality);
1322     fOutput->Add(fhBckgMult);
1323     fOutput->Add(fhBckgFluc);
1324     fOutput->Add(fhChargedJetPt);
1325     fOutput->Add(fhChargedJetPtAreaCut);
1326     fOutput->Add(fhJetPtEMCal);
1327     fOutput->Add(fhJetPtEMCalAreaCut);
1328     fOutput->Add(fhJetPtEMCalAreaCutSignal);
1329     fOutput->Add(fhJetPtTPC);
1330     fOutput->Add(fhJetPtTPCAreaCut);
1331     fOutput->Add(fhJetPtArea);
1332     fOutput->Add(fhRhoTotal);
1333     fOutput->Add(fhRhoNoLeading);
1334     fOutput->Add(fhJetTPtRhoTotal);
1335     fOutput->Add(fhJetTPtRhoTotalSignal);
1336     fOutput->Add(fhJetTPtRhoNoLeading);
1337     fOutput->Add(fhJetTPtRhoNoLeadingSignal);
1338     fOutput->Add(fhJetTrigR1A);
1339     fOutput->Add(fhJetTPt1B);
1340     fOutput->Add(fhJetTPt1BSignal);
1341     fOutput->Add(fhEMCalBckg1B);
1342     fOutput->Add(fhRho1B);
1343     fOutput->Add(fhJetTPt1C);
1344     fOutput->Add(fhRho1C);
1345     fOutput->Add(fhEMCalBckg1C);
1346     fOutput->Add(fhEMCalJet2A);
1347     fOutput->Add(fh020EMCalJet2A);
1348     fOutput->Add(fhJetTPt2B);
1349     fOutput->Add(fhEMCalBckg2B);
1350     fOutput->Add(fhRho2B);
1351     fOutput->Add(fhRho3);
1352     fOutput->Add(fhJetTPt3);
1353     fOutput->Add(fhDeltaPtTotal);
1354     fOutput->Add(fhDeltaPtNoLeading);
1355     fOutput->Add(fhDeltaPt1B);
1356     fOutput->Add(fhJetConstituentPt);
1357     fOutput->Add(fhDeltaRho01);
1358     fOutput->Add(fhEMCalCellCounts);
1359     fOutput->Add(fh020RhoTotal);
1360     fOutput->Add(fh020RhoNoLeading);
1361     fOutput->Add(fh020Rho1B);
1362     fOutput->Add(fh020Rho2B);
1363     fOutput->Add(fh020Rho3);
1364     fOutput->Add(fh020JetPtEMCal);
1365     fOutput->Add(fh020JetPtEMCalAreaCut);
1366     fOutput->Add(fh020JetPtEMCalAreaCutSignal);
1367     fOutput->Add(fh020JetTPtRhoTotal);
1368     fOutput->Add(fh020JetTPtRhoTotalSignal);
1369     fOutput->Add(fh020JetTPtRhoNoLeading);
1370     fOutput->Add(fh020JetTPtRhoNoLeadingSignal);
1371     fOutput->Add(fh020JetTPt1B);
1372     fOutput->Add(fh020JetTPt1BSignal);
1373     fOutput->Add(fh020JetTPt1C);
1374     fOutput->Add(fh020JetTPt2B);
1375     fOutput->Add(fh020JetTPt3);
1376     fOutput->Add(fhDeltaPt2B);
1377     fOutput->Add(fhDeltaPtkT);
1378     fOutput->Add(fh020DiJetAsy);
1379     fOutput->Add(fh020RhokT);
1380     fOutput->Add(fh020EMCalkTClusters);
1381     fOutput->Add(fh020EMCalAkTJets);
1382     fOutput->Add(fh020DiJetDeltaPhi);
1383     fOutput->Add(fhDiJetEMCalLeadingPt);
1384     fOutput->Add(fhDiJetEMCalLeadingDeltaPhi);
1385     fOutput->Add(fhJetPtCenEMCal);
1386     fOutput->Add(fhJetPtCenEMCalAreaCut);
1387     fOutput->Add(fhJetPtCenEMCalAreaCutSignal);
1388     fOutput->Add(fhJetTPtCenRhoTotal);
1389     fOutput->Add(fhJetTPtCenRhoTotalSignal);
1390     fOutput->Add(fhJetTPtCenRhoNoLeading);
1391     fOutput->Add(fhJetTPtCenRhoNoLeadingSignal);
1392     fOutput->Add(fhJetTPtCen1B);
1393     fOutput->Add(fhJetTPtCen1BSignal);
1394     fOutput->Add(fhJetTPtCen1C);
1395     fOutput->Add(fhJetTPtCen2B);
1396     fOutput->Add(fhJetTPtCen3);
1397     fOutput->Add(fhDiJetCenAsy);
1398     fOutput->Add(fhDiJetCenDeltaPhi);
1399     fOutput->Add(fhEMCalCenJet2A);
1400     fOutput->Add(fhDeltaRho0DiJet);
1401     fOutput->Add(fh020Rho2BCore);
1402     fOutput->Add(fhRho2BCore);
1403     fOutput->Add(fh020Rho3NoJets);
1404     fOutput->Add(fh020Rho3DiJets);
1405     fOutput->Add(fhRho3NoJets);
1406     fOutput->Add(fhRho3DiJets);
1407     fOutput->Add(fh020Rho3Perp);
1408     fOutput->Add(fhRho3Perp);
1409
1410     fOutput->Add(fpEventMult);
1411     fOutput->Add(fpRhoTotal);
1412     fOutput->Add(fpRhoNoLeading);
1413     fOutput->Add(fpRho1B);
1414     fOutput->Add(fpRho2B);
1415     fOutput->Add(fpRho2BCore);
1416     fOutput->Add(fpRho3);
1417     fOutput->Add(fpRho3NoJets);
1418     fOutput->Add(fpRho3DiJets);
1419     fOutput->Add(fpRho3Perp);
1420     fOutput->Add(fpRhoScale);
1421     fOutput->Add(fpRhokT);
1422     fOutput->Add(fpRhoChargedkT);
1423     fOutput->Add(fpRhoScalekT);
1424     fOutput->Add(fpJetPtRhoTotal);
1425     fOutput->Add(fpJetPtRhoNoLeading);
1426     fOutput->Add(fpJetPtRhokT);
1427
1428     fOutput->Add(fpTrackPtProfile);
1429     fOutput->Add(fpClusterPtProfile);
1430
1431     // Post data for ALL output slots >0 here,
1432     // To get at least an empty histogram
1433     // 1 is the outputnumber of a certain weg of task 1
1434     PostData(1, fOutput);
1435 }
1436
1437 void AliAnalysisTaskFullpAJets::UserExecOnce()
1438 {
1439     // Get the event tracks from PicoTracks
1440     TString track_name="PicoTracks";
1441     fmyTracks =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(track_name));
1442
1443     // Get the event caloclusters from CaloClustersCorr
1444     TString cluster_name="CaloClustersCorr";
1445     fmyClusters =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(cluster_name));
1446
1447     // Get charged jets
1448     TString jet_algorithm=Form("Jet_AKTChargedR0%d0_PicoTracks_pT0150",fRJET);
1449     fmyAKTChargedJets =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
1450
1451     jet_algorithm=Form("Jet_KTChargedR0%d0_PicoTracks_pT0150",fRJET);
1452     fmyKTChargedJets =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
1453
1454     // Get the full jets
1455     jet_algorithm=Form("Jet_AKTFullR0%d0_PicoTracks_pT0150_CaloClustersCorr_ET0300",fRJET);
1456     fmyAKTFullJets =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
1457
1458     jet_algorithm=Form("Jet_KTFullR0%d0_PicoTracks_pT0150_CaloClustersCorr_ET0300",fRJET);
1459     fmyKTFullJets =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
1460
1461     jet_algorithm="";
1462     fIsInitialized=kTRUE;
1463 }
1464 //________________________________________________________________________
1465 void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
1466 {
1467     if (fIsInitialized==kFALSE)
1468     {
1469         UserExecOnce();
1470     }
1471
1472     // Get pointer to reconstructed event
1473     AliVEvent *event = InputEvent();
1474     if (!event)
1475     {
1476         AliError("Pointer == 0, this can not happen!");
1477         return;
1478     }
1479
1480     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
1481     AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1482
1483     if (esd)
1484     {
1485         fEventCentrality=esd->GetCentrality()->GetCentralityPercentile("V0M");
1486
1487         if (esd->GetPrimaryVertex()->GetNContributors()<1 || (TMath::Abs(esd->GetPrimaryVertex()->GetZ())>fVertexWindow))
1488         {
1489             return;
1490         }
1491         if (TMath::Sqrt(TMath::Power(esd->GetPrimaryVertex()->GetX(),2)+TMath::Power(esd->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
1492         {
1493             return;
1494         }
1495
1496         esd->GetPrimaryVertex()->GetXYZ(fvertex);
1497     }
1498     else if (aod)
1499     {
1500         //cout<<"Hello AOD"<<endl;
1501
1502         fEventCentrality=aod->GetCentrality()->GetCentralityPercentile("V0M");
1503
1504         if (aod->GetPrimaryVertex()->GetNContributors()<1 || (TMath::Abs(aod->GetPrimaryVertex()->GetZ())>fVertexWindow))
1505         {
1506             return;
1507         }
1508         if (TMath::Sqrt(TMath::Power(aod->GetPrimaryVertex()->GetX(),2)+TMath::Power(aod->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
1509         {
1510             return;
1511         }
1512
1513         aod->GetPrimaryVertex()->GetXYZ(fvertex);
1514     }
1515     else
1516     {
1517         AliError("Cannot get AOD/ESD event. Rejecting Event");
1518         return;
1519     }
1520
1521     // Make sure Centrality isn't exactly 100% (to avoid bin filling errors in profile plots. Make it 99.99
1522     if (fEventCentrality>99.99)
1523     {
1524         fEventCentrality=99.99;
1525     }
1526     fhCentrality->Fill(fEventCentrality);
1527
1528     fnTracks = fmyTracks->GetEntries();
1529     //cout<<"n Tracks:"<<fnTracks<<endl;
1530     // Reject any event that doesn't have any tracks, i.e. TPC is off
1531     if (fnTracks<1)
1532     {
1533         AliWarning("No PicoTracks, Rejecting Event");
1534         return;
1535     }
1536
1537     fnClusters = fmyClusters->GetEntries();
1538     //cout<<"n Cluster:"<<fnClusters<<endl;
1539     if (fnClusters<1)
1540     {
1541         AliInfo("No Corrected CaloClusters, using only charged jets");
1542
1543         TrackHisto();
1544         InitChargedJets();
1545         Method3(kFALSE);
1546         JetPtChargedProfile();
1547         DeleteArrays(kFALSE);
1548
1549         fnEventsCharged++;
1550         return;
1551     }
1552
1553     TrackHisto();
1554     ClusterHisto();
1555
1556     // Prep the jets
1557     InitChargedJets();
1558     InitFullJets();
1559     EventHistos();
1560
1561     EstimateTotalBackground();
1562     EstimateBackgoundMinusLJet();
1563
1564     Method1A();
1565     Method1B();
1566     Method1C();
1567
1568     // Method 2
1569     if (IsDiJetEvent()==kTRUE)
1570     {
1571         Method2A();
1572         Method2B();
1573         Method3DiJet();
1574         Method3Perp();
1575     }
1576
1577     Method3(kTRUE);
1578
1579     // Compute Jet Energy Density Profile
1580     JetPtChargedProfile();
1581     JetPtFullProfile();
1582     JetPtEtaProfile();
1583
1584     // Delete Dynamic Arrays
1585     DeleteArrays(kTRUE);
1586     fnEvents++;
1587
1588     PostData(1, fOutput);
1589 }
1590
1591 //________________________________________________________________________
1592 void AliAnalysisTaskFullpAJets::Terminate(Option_t *) //specify what you want to have done
1593 {
1594     // Called once at the end of the query. Done nothing here.
1595 }
1596
1597 /////////////////////////////////////////////////////////////////////////////////////////
1598 /////////////////     User Defined Sub_Routines   ///////////////////////////////////////
1599 /////////////////////////////////////////////////////////////////////////////////////////
1600
1601 void AliAnalysisTaskFullpAJets::TrackHisto()
1602 {
1603     // Fill track histograms: Phi,Eta,Pt
1604     Int_t i,j;
1605     Int_t TCBins=100;
1606     TH2D *hdummypT= new TH2D("hdummypT","",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);  //!
1607
1608     for (i=0;i<fnTracks;i++)
1609     {
1610         AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1611         fhTrackPt->Fill(vtrack->Pt());
1612         fhTrackEta->Fill(vtrack->Eta());
1613         fhTrackPhi->Fill(vtrack->Phi());
1614         fhTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
1615         hdummypT->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
1616     }
1617     for (i=1;i<=TCBins;i++)
1618     {
1619         for (j=1;j<=TCBins;j++)
1620         {
1621             fpTrackPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fTPCArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
1622         }
1623     }
1624     delete hdummypT;
1625 }
1626
1627 void AliAnalysisTaskFullpAJets::ClusterHisto()
1628 {
1629     // Fill cluster histograms: Phi,Eta,Pt
1630     Int_t i,j;
1631     Int_t TCBins=100;
1632     TH2D *hdummypT= new TH2D("hdummypT","",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);  //!
1633     AliEMCALGeometry *myAliEMCGeo = AliEMCALGeometry::GetInstance();
1634     Int_t myCellID=-2;
1635
1636     for (i=0;i<fnClusters;i++)
1637     {
1638         AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1639         TLorentzVector *cluster_vec = new TLorentzVector;
1640         vcluster->GetMomentum(*cluster_vec,fvertex);
1641
1642         fhClusterPt->Fill(cluster_vec->Pt());
1643         fhClusterEta->Fill(cluster_vec->Eta());
1644         fhClusterPhi->Fill(cluster_vec->Phi());
1645         fhClusterEtaPhi->Fill(cluster_vec->Eta(),cluster_vec->Phi());
1646         hdummypT->Fill(cluster_vec->Eta(),cluster_vec->Phi(),cluster_vec->Pt());
1647         myAliEMCGeo->GetAbsCellIdFromEtaPhi(cluster_vec->Eta(),cluster_vec->Phi(),myCellID);
1648         fhEMCalCellCounts->Fill(myCellID);
1649         //cout<<"Cluster ID:"<<i<<"  (Eta,Phi): ("<<cluster_vec->Eta()<<","<<cluster_vec->Phi()<<")  Cell ID:"<<myCellID<<endl;
1650         delete cluster_vec;
1651     }
1652     for (i=1;i<=TCBins;i++)
1653     {
1654         for (j=1;j<=TCBins;j++)
1655         {
1656             fpClusterPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fEMCalArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
1657         }
1658     }
1659     //myAliEMCGeo->GetAbsCellIdFromEtaPhi(0.38,2.88,myCellID);
1660     //cout<<"Cell ID Test:"<<myCellID<<endl;
1661     delete hdummypT;
1662 }
1663
1664 void AliAnalysisTaskFullpAJets::EventHistos()
1665 {
1666     Int_t i,j;
1667     Double_t E_tracks_total=0.;
1668     Double_t E_caloclusters_total=0.;
1669     TRandom3 u(time(NULL));
1670
1671     Double_t Eta_Center=0.5*(fEMCalEtaMin+fEMCalEtaMax);
1672     Double_t Phi_Center=0.5*(fEMCalPhiMin+fEMCalPhiMax);
1673     Int_t event_mult=0;
1674     Int_t clus_mult=0;
1675
1676     TLorentzVector *dummy= new TLorentzVector;
1677
1678     for (j=0;j<fnBckgClusters;j++)
1679     {
1680         event_mult=0;
1681         clus_mult=0;
1682         E_tracks_total=0.;
1683         E_caloclusters_total=0.;
1684         dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1685
1686         //  Loop over all tracks
1687         for (i=0;i<fnTracks;i++)
1688         {
1689             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1690             if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1691             {
1692                 event_mult++;
1693                 TLorentzVector *track_vec = new TLorentzVector;
1694                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1695                 if (dummy->DeltaR(*track_vec)<fJetR)
1696                 {
1697                     clus_mult++;
1698                     E_tracks_total+=vtrack->Pt();
1699                 }
1700                 delete track_vec;
1701             }
1702         }
1703
1704         //  Loop over all caloclusters
1705         for (i=0;i<fnClusters;i++)
1706         {
1707             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1708             TLorentzVector *cluster_vec = new TLorentzVector;
1709             vcluster->GetMomentum(*cluster_vec,fvertex);
1710             event_mult++;
1711             if (dummy->DeltaR(*cluster_vec)<fJetR)
1712             {
1713                 clus_mult++;
1714                 E_caloclusters_total+=vcluster->E();
1715             }
1716             delete cluster_vec;
1717         }
1718         //  Fill Histograms
1719         if (fEventCentrality<=20)
1720         {
1721             fhBckgMult->Fill(clus_mult);
1722             fhBckgFluc->Fill(E_tracks_total+E_caloclusters_total);
1723             fRCBckgFluc[j]=E_tracks_total+E_caloclusters_total;
1724         }
1725     }
1726
1727     fpEventMult->Fill(fEventCentrality,event_mult);
1728     delete dummy;
1729 }
1730
1731 void AliAnalysisTaskFullpAJets::InitChargedJets()
1732 {
1733     // Preliminary Jet Placement and Selection Cuts
1734     Int_t i,j;
1735     Double_t kTRho=0.0;
1736     Double_t delta_phi=0.0;
1737
1738     fnAKTChargedJets = fmyAKTChargedJets->GetEntries();
1739     fnKTChargedJets = fmyKTChargedJets->GetEntries();
1740     fJetPtChargedCutID = new Int_t[fnAKTChargedJets+1];
1741     fJetkTTPCFullID = new Int_t[fnKTChargedJets+1];
1742     fnJetsChargedPtCut=0;
1743     fnJetskTTPCFull=0;
1744     fPtChargedMaxID=-1;  // Initialize to not have any jet(s) fully contained within
1745     fPtChargedMax=0.0;
1746
1747     fInTPCChargedFull = new Bool_t[fnAKTChargedJets+1];
1748
1749     for (i=0;i<fnAKTChargedJets;i++)
1750     {
1751         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(i);
1752
1753         // Determine where the jet is
1754         fInTPCChargedFull[i]=IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE);
1755
1756         // Fill Histograms with appropriate Jet Kinematics
1757         if (fInTPCChargedFull[i]==kTRUE)
1758         {
1759             fhChargedJetPt->Fill(myJet->Pt());
1760
1761             if (myJet->Pt()>=fPtChargedMax)
1762             {
1763                 fPtChargedMax=myJet->Pt();
1764                 fPtChargedMaxID=i;
1765             }
1766             //  Now determine if the jet is above the EMCal Jet Pt Threshold
1767             if (myJet->Area()>=fJetAreaThreshold)
1768             {
1769                 fhChargedJetPtAreaCut->Fill(myJet->Pt());
1770             }
1771             if (myJet->Pt()>=fTPCJetThreshold && myJet->Area()>=fJetAreaThreshold)
1772             {
1773                 fJetPtChargedCutID[fnJetsChargedPtCut]=i;
1774                 fnJetsChargedPtCut++;
1775             }
1776         }
1777     }
1778
1779     // Fill dijet delta phi plots
1780     if (fnJetsChargedPtCut>1)
1781     {
1782         AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fPtChargedMaxID);
1783         j=0;
1784         while (j<fnJetsChargedPtCut)
1785         {
1786             if (fJetPtChargedCutID[j]==fPtChargedMaxID)
1787             {
1788                 j++;
1789             }
1790             else
1791             {
1792                 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fJetPtChargedCutID[j]);
1793
1794                 delta_phi=TMath::Min(TMath::Abs(myhJet->Phi()-myJet->Phi()),2*TMath::Pi()-TMath::Abs(myhJet->Phi()-myJet->Phi()));
1795                 fhDiJetCenDeltaPhi->Fill(delta_phi,fEventCentrality);
1796                 if (fEventCentrality<=20)
1797                 {
1798                     fh020DiJetDeltaPhi->Fill(delta_phi);
1799                 }
1800                 j++;
1801             }
1802         }
1803     }
1804
1805     fRhokTCharged=0.0;
1806     // kT jets. Used for calculating rho
1807     Int_t nkT_mid=-1;
1808     for (i=0;i<fnKTChargedJets;i++)
1809     {
1810         AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(i);
1811
1812         if (IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE)==kTRUE)
1813         {
1814             fJetkTTPCFullID[fnJetskTTPCFull]=i;
1815             fnJetskTTPCFull++;
1816         }
1817     }
1818
1819     if (fnJetskTTPCFull>=2)
1820     {
1821         nkT_mid=fnJetskTTPCFull/2;
1822     }
1823     else if (fnJetskTTPCFull==1)
1824     {
1825         nkT_mid=0;
1826     }
1827
1828     if (nkT_mid>=0)
1829     {
1830         AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fJetkTTPCFullID[nkT_mid]);
1831         kTRho=myJet->Pt()/myJet->Area();
1832         fRhokTCharged=kTRho;
1833         fpRhoChargedkT->Fill(fEventCentrality,kTRho);
1834     }
1835 }
1836
1837 void AliAnalysisTaskFullpAJets::InitFullJets()
1838 {
1839     // Preliminary Jet Placement and Selection Cuts
1840     Int_t i;
1841     Int_t EMCalFullCount=0;
1842     Double_t kTRho=0.0;
1843
1844     fnAKTFullJets = fmyAKTFullJets->GetEntries();
1845     fnKTFullJets = fmyKTFullJets->GetEntries();
1846
1847     fJetPtCutID = new Int_t[fnAKTFullJets+1];
1848     fJetPtTPCCutID= new Int_t[fnAKTFullJets+1];
1849     fJetPtTotalCutID= new Int_t[fnAKTFullJets+1];
1850     fJetkTEMCalFullID= new Int_t[fnKTFullJets+1];
1851
1852     fnJetsPtCut=0;
1853     fnJetsPtTPCCut=0;
1854     fnJetsPtTotalCut=0;
1855     fnJetskTEMCalFull=0;
1856
1857     fPtMaxID=-1;  // Initialize to not have any jet(s) in EMCal
1858     fPtFullMaxID=-1;  // Initialize to not have any jet(s) fully contained within EMCal
1859     fPtTPCMaxID=-1;  // Initialize to not have any jet(s) in TPC
1860     fPtFullTPCMaxID=-1;  // Initialize to not have any jet(s) fully contained within TPC
1861     fPtTotalMaxID=-1;  // Initialize to not have any jet(s) in Total Acceptance
1862
1863     fPtMax=0.;
1864     fPtFullMax=0.;
1865     fPtTPCMax=0.;
1866     fPtFullTPCMax=0.;
1867     fPtTotalMax=0.;
1868
1869     fInEMCal = new Bool_t[fnAKTFullJets+1];
1870     fInEMCalFull = new Bool_t[fnAKTFullJets+1];
1871     fInTPCFull = new Bool_t[fnAKTFullJets+1];
1872
1873     for (i=0;i<fnAKTFullJets;i++)
1874     {
1875         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
1876
1877         // Area distribution of the jet
1878         fhJetPtArea->Fill(myJet->Pt(),myJet->Area());
1879
1880         // Determine where the jet is
1881         fInEMCal[i]=IsInEMCalPart(fJetR,myJet->Phi(),myJet->Eta());
1882         fInEMCalFull[i]=IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta());
1883         fInTPCFull[i]=IsInTPCFull(fJetR,myJet->Phi(),myJet->Eta());
1884
1885         // Fill Histograms with appropriate Jet Kinematics
1886         if (fInEMCal[i]==kTRUE)
1887         {
1888             // Method 1A
1889             if (myJet->Pt()>fPtMax)
1890             {
1891                 fPtMax=myJet->Pt();
1892                 fPtMaxID=i;
1893             }
1894
1895             // Method 1B
1896             if (fInEMCalFull[i]==kTRUE)
1897             {
1898                 // Fill Jet Pt Distribution
1899                 fhJetPtEMCal->Fill(myJet->Pt());
1900                 fhJetPtCenEMCal->Fill(myJet->Pt(),fEventCentrality);
1901                 if (fEventCentrality<=20)
1902                 {
1903                     fh020JetPtEMCal->Fill(myJet->Pt());
1904                 }
1905                 if (myJet->Area()>=fJetAreaThreshold)
1906                 {
1907                     fhJetPtEMCalAreaCut->Fill(myJet->Pt());
1908                     fhJetPtCenEMCalAreaCut->Fill(myJet->Pt(),fEventCentrality);
1909                     if (fEventCentrality<=20)
1910                     {
1911                         fh020JetPtEMCalAreaCut->Fill(myJet->Pt());
1912                         EMCalFullCount++;
1913                     }
1914                     if (myJet->Pt()>=fEMCalJetThreshold)
1915                     {
1916                         fhJetPtEMCalAreaCutSignal->Fill(myJet->Pt());
1917                         fhJetPtCenEMCalAreaCutSignal->Fill(myJet->Pt(),fEventCentrality);
1918                         if (fEventCentrality<=20)
1919                         {
1920                             fh020JetPtEMCalAreaCutSignal->Fill(myJet->Pt());
1921                         }
1922                     }
1923                 }
1924
1925                 if (myJet->Pt()>=fPtFullMax)
1926                 {
1927                     fPtFullMax=myJet->Pt();
1928                     fPtFullMaxID=i;
1929                 }
1930
1931                 //  Now determine if the jet is above the EMCal Jet Pt Threshold
1932                 if (myJet->Pt()>=fEMCalJetThreshold)
1933                 {
1934                     fJetPtCutID[fnJetsPtCut]=i;
1935                     fnJetsPtCut++;
1936                 }
1937             }
1938         }
1939         else
1940         {
1941             // Method 2A
1942             if (myJet->Pt()>fPtTPCMax)
1943             {
1944                 fPtTPCMax=myJet->Pt();
1945                 fPtTPCMaxID=i;
1946             }
1947             if (fInTPCFull[i]==kTRUE)
1948             {
1949                 // Jet Pt distribution
1950                 fhJetPtTPC->Fill(myJet->Pt());
1951                 if (myJet->Area()>=fJetAreaThreshold)
1952                 {
1953                     fhJetPtTPCAreaCut->Fill(myJet->Pt());
1954                 }
1955
1956                 if (myJet->Pt()>fPtFullTPCMax)
1957                 {
1958                     fPtFullTPCMax=myJet->Pt();
1959                     fPtFullTPCMaxID=i;
1960                 }
1961             }
1962
1963             //  Now determine if the jet is above the TPC Jet Pt Threshold
1964             if (myJet->Pt()>=fTPCJetThreshold)
1965             {
1966                 fJetPtTPCCutID[fnJetsPtTPCCut]=i;
1967                 fnJetsPtTPCCut++;
1968             }
1969         }
1970         // Find all jet(s) above the threshold within the Detector (TPC+EMCal; No Fudicial cut) for Plan2
1971         if (myJet->Pt()>fPtTotalMax)
1972         {
1973             fPtTotalMax=myJet->Pt();
1974             fPtTotalMaxID=i;
1975         }
1976         // And if they are above the threshold?
1977         if (myJet->Pt()>=fTPCJetThreshold)
1978         {
1979             fJetPtTotalCutID[fnJetsPtTotalCut]=i;
1980             fnJetsPtTotalCut++;
1981         }
1982     }
1983     fh020EMCalAkTJets->Fill(EMCalFullCount);
1984     EMCalFullCount=0;
1985
1986     fRhokTTotal=0.0;
1987     // kT jets. Used for calculating rho
1988     Int_t nkT_mid=-1;
1989     for (i=0;i<fnKTFullJets;i++)
1990     {
1991         AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(i);
1992
1993         if (IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta())==kTRUE)
1994         {
1995             fJetkTEMCalFullID[fnJetskTEMCalFull]=i;
1996             fnJetskTEMCalFull++;
1997         }
1998     }
1999     if (fEventCentrality<=20)
2000     {
2001         fh020EMCalkTClusters->Fill(fnJetskTEMCalFull);
2002     }
2003
2004     if (fnJetskTEMCalFull>=2)
2005     {
2006         nkT_mid=fnJetskTEMCalFull/2;
2007     }
2008     else if (fnJetskTEMCalFull==1)
2009     {
2010         nkT_mid=0;
2011     }
2012
2013     if (nkT_mid>=0)
2014     {
2015         AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fJetkTEMCalFullID[nkT_mid]);
2016         kTRho=myJet->Pt()/myJet->Area();
2017         fRhokTTotal=kTRho;
2018         fpRhokT->Fill(fEventCentrality,kTRho);
2019         fpJetPtRhokT->Fill(fPtFullMax,kTRho);
2020         if (fRhokTCharged!=0.0)
2021         {
2022             fpRhoScalekT->Fill(fEventCentrality,fRhokTTotal/fRhokTCharged);
2023         }
2024         if (fEventCentrality<=20)
2025         {
2026             FillBckgFlucDeltaPt(fhDeltaPtkT,kTRho);
2027             fh020RhokT->Fill(kTRho);
2028         }
2029     }
2030 }
2031
2032 void AliAnalysisTaskFullpAJets::EstimateTotalBackground()
2033 {
2034     Int_t i;
2035     Double_t E_tracks_total=0.;
2036     Double_t E_caloclusters_total=0.;
2037     Double_t EMCal_rho=0.;
2038     fDeltaRho01=0.0;
2039
2040     //  Loop over all tracks
2041     for (i=0;i<fnTracks;i++)
2042     {
2043         AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2044         if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2045         {
2046             E_tracks_total+=vtrack->Pt();
2047         }
2048     }
2049
2050     //  Loop over all caloclusters
2051     for (i=0;i<fnClusters;i++)
2052     {
2053         AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2054         E_caloclusters_total+=vcluster->E();
2055     }
2056
2057     //  Calculate the mean Background density
2058     EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2059     fRhoAkTTotal=EMCal_rho;
2060
2061     //  Fill histograms
2062     fhRhoTotal->Fill(EMCal_rho,fEventCentrality);
2063     fpRhoTotal->Fill(fEventCentrality,EMCal_rho);
2064     fpJetPtRhoTotal->Fill(fPtFullMax,EMCal_rho);
2065     FillFullCorrJetPt(fhJetTPtRhoTotal,EMCal_rho,kFALSE);
2066     FillFullCorrJetPt(fhJetTPtRhoTotalSignal,EMCal_rho,kTRUE);
2067     fDeltaRho01=EMCal_rho;
2068
2069     FillFullCorrJetPt(fhJetTPtCenRhoTotal,EMCal_rho,kFALSE);
2070     FillFullCorrJetPt(fhJetTPtCenRhoTotalSignal,EMCal_rho,kTRUE);
2071     // Fill Background fluctuation spectrum F(A)
2072     if (fEventCentrality<=20)
2073     {
2074         FillFullCorrJetPt(fh020JetTPtRhoTotal,EMCal_rho,kFALSE);
2075         FillFullCorrJetPt(fh020JetTPtRhoTotalSignal,EMCal_rho,kTRUE);
2076         FillBckgFlucDeltaPt(fhDeltaPtTotal,EMCal_rho);
2077         fh020RhoTotal->Fill(EMCal_rho);
2078     }
2079 }
2080
2081 void AliAnalysisTaskFullpAJets::EstimateBackgoundMinusLJet()
2082 {
2083     Int_t i;
2084     Double_t E_tracks_total=0.;
2085     Double_t E_caloclusters_total=0.;
2086     Double_t EMCal_rho=0.;
2087
2088     if (fPtFullMax>=fEMCalJetThreshold)
2089     {
2090         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fPtMaxID);
2091         TLorentzVector *temp_jet= new TLorentzVector;
2092         myJet->GetMom(*temp_jet);
2093
2094         //  Loop over all tracks
2095         for (i=0;i<fnTracks;i++)
2096         {
2097             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2098             if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2099             {
2100                 TLorentzVector *track_vec = new TLorentzVector;
2101                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2102                 if (temp_jet->DeltaR(*track_vec)>fJetR)
2103                 {
2104                     E_tracks_total+=vtrack->Pt();
2105                 }
2106                 delete track_vec;
2107             }
2108         }
2109
2110         //  Loop over all caloclusters
2111         for (i=0;i<fnClusters;i++)
2112         {
2113             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2114             TLorentzVector *cluster_vec = new TLorentzVector;
2115             vcluster->GetMomentum(*cluster_vec,fvertex);
2116             if (temp_jet->DeltaR(*cluster_vec)>fJetR)
2117             {
2118                 E_caloclusters_total+=vcluster->E();
2119             }
2120             delete cluster_vec;
2121         }
2122         delete temp_jet;
2123         //  Calculate the mean Background density
2124         EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-TMath::Pi()*TMath::Power(fJetR,2));
2125     }
2126     else  // i.e. No signal jets-> same as total background density
2127     {
2128         //  Loop over all tracks
2129         for (i=0;i<fnTracks;i++)
2130         {
2131             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2132             if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2133             {
2134                 E_tracks_total+=vtrack->Pt();
2135             }
2136         }
2137
2138         //  Loop over all caloclusters
2139         for (i=0;i<fnClusters;i++)
2140         {
2141             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2142             E_caloclusters_total+=vcluster->E();
2143         }
2144         //  Calculate the mean Background density
2145         EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2146     }
2147
2148     //  Fill histograms
2149     fhRhoNoLeading->Fill(EMCal_rho,fEventCentrality);
2150     fpRhoNoLeading->Fill(fEventCentrality,EMCal_rho);
2151     fpJetPtRhoNoLeading->Fill(fPtFullMax,EMCal_rho);
2152     FillFullCorrJetPt(fhJetTPtRhoNoLeading,EMCal_rho,kFALSE);
2153     FillFullCorrJetPt(fhJetTPtRhoNoLeadingSignal,EMCal_rho,kTRUE);
2154     fDeltaRho01-=EMCal_rho;
2155     FillFullDeltaRho(fhDeltaRho01,fDeltaRho01,kTRUE);
2156     fDeltaRho01=0.0;
2157
2158     FillFullCorrJetPt(fhJetTPtCenRhoNoLeading,EMCal_rho,kFALSE);
2159     FillFullCorrJetPt(fhJetTPtCenRhoNoLeadingSignal,EMCal_rho,kTRUE);
2160     // Fill Background fluctuation spectrum F(A)
2161     if (fEventCentrality<=20)
2162     {
2163         FillFullCorrJetPt(fh020JetTPtRhoNoLeading,EMCal_rho,kFALSE);
2164         FillFullCorrJetPt(fh020JetTPtRhoNoLeadingSignal,EMCal_rho,kTRUE);
2165         FillBckgFlucDeltaPt(fhDeltaPtNoLeading,EMCal_rho);
2166         fh020RhoNoLeading->Fill(EMCal_rho);
2167     }
2168 }
2169
2170 void AliAnalysisTaskFullpAJets::Method1A()
2171 {
2172     Int_t i;
2173     Double_t delta_R;
2174
2175     if (fPtMax>=fEMCalJetThreshold && fnAKTFullJets>1)
2176     {
2177         TLorentzVector *high_jet= new TLorentzVector;
2178         TLorentzVector *temp_jet= new TLorentzVector;
2179
2180         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fPtMaxID);
2181         myJet->GetMom(*high_jet);
2182         //cout<<"HJ Phi="<<high_jet->Phi()<<"   HJ Eta="<<high_jet->Eta()<<endl;
2183         for(i=0;i<fnAKTFullJets;i++)
2184         {
2185             if (i!=fPtMaxID && fInEMCalFull[i]==kTRUE)
2186             {
2187                 AliEmcalJet *myBckg =(AliEmcalJet*) fmyAKTFullJets->At(i);
2188                 myBckg->GetMom(*temp_jet);
2189                 delta_R=temp_jet->DeltaR(*high_jet);
2190                 //cout<<"TJ Phi="<<temp_jet->Phi()<<"   TJ Eta="<<temp_jet->Eta()<<endl;
2191                 //cout<<"Delta R="<<delta_R<<endl;
2192                 if (delta_R>=(2*fJetR))
2193                 {
2194                     fhJetTrigR1A->Fill(myBckg->Pt(),fPtMax,delta_R);
2195                 }
2196             }
2197         }
2198         delete high_jet;
2199         delete temp_jet;
2200     }
2201 }
2202
2203 void AliAnalysisTaskFullpAJets::Method1B()
2204 {
2205     Int_t i,j;
2206     Bool_t track_away_from_jet;
2207     Bool_t cluster_away_from_jet;
2208     Double_t E_tracks_total=0.;
2209     Double_t E_caloclusters_total=0.;
2210     Double_t EMCal_rho=0.;
2211     Double_t jet_area_total=0.;
2212
2213     // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
2214     fRhoTotal=0;
2215     for (i=0;i<fnTracks;i++)
2216     {
2217         // First, check if track is in the EMCal!!
2218         AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i); // pointer to reconstructed to track
2219         if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2220         {
2221             if (fnJetsPtCut<1)
2222             {
2223                 E_tracks_total+=vtrack->Pt();
2224             }
2225             else
2226             {
2227                 track_away_from_jet=kTRUE;
2228                 j=0;
2229                 TLorentzVector *track_vec = new TLorentzVector;
2230                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2231                 //cout<<endl<<endl<<endl<<"Event # :"<<fnEvents+1<<"  njets="<<fnAKTFullJets<<"  Threshold EMCal jets="<<fnJetsPtCut<<"  tracks # :"<<i<<","<<fnTracks<<endl;
2232                 while (track_away_from_jet==kTRUE && j<fnJetsPtCut)
2233                 {
2234                     //cout<<"j="<<j<<endl<<endl<<endl;
2235
2236                     TLorentzVector *jet_vec= new TLorentzVector;
2237                     AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fJetPtCutID[j]);
2238                     myJet->GetMom(*jet_vec);
2239                     if (track_vec->DeltaR(*jet_vec)<=(fJetR))
2240                     {
2241                         track_away_from_jet=kFALSE;
2242                     }
2243                     delete jet_vec;
2244                     j++;
2245                 }
2246                 if (track_away_from_jet==kTRUE)
2247                 {
2248                     E_tracks_total+=vtrack->Pt();
2249                 }
2250                 delete track_vec;
2251             }
2252         }
2253     }
2254
2255     // Next, sum all CaloClusters within the EMCal (obviously all clusters must be within EMCal!!) that are away from jet(s) above Pt Threshold
2256
2257     for (i=0;i<fnClusters;i++)
2258     {
2259         AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i); // pointer to reconstructed to cluster
2260         if (fnJetsPtCut<1)
2261         {
2262             E_caloclusters_total+=vcluster->E();
2263         }
2264         else
2265         {
2266             cluster_away_from_jet=kTRUE;
2267             j=0;
2268
2269             TLorentzVector *cluster_vec = new TLorentzVector;
2270             vcluster->GetMomentum(*cluster_vec,fvertex);
2271             //cout<<endl<<endl<<endl<<"Event # :"<<fnEvents+1<<"  njets="<<fnAKTFullJets<<"  Threshold EMCal jets="<<fnJetsPtCut<<"  cluster # :"<<i<<","<<fnClusters<<endl;
2272
2273             while (cluster_away_from_jet==kTRUE && j<fnJetsPtCut)
2274             {
2275                 TLorentzVector *jet_vec= new TLorentzVector;
2276                 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fJetPtCutID[j]);
2277                 myJet->GetMom(*jet_vec);
2278                 if (cluster_vec->DeltaR(*jet_vec)<=(fJetR))
2279                 {
2280                     cluster_away_from_jet=kFALSE;
2281                 }
2282                 delete jet_vec;
2283                 j++;
2284             }
2285             if (cluster_away_from_jet==kTRUE)
2286             {
2287                 E_caloclusters_total+=vcluster->E();
2288             }
2289             delete cluster_vec;
2290         }
2291     }
2292
2293     // Determine area of all Jets that are within the EMCal
2294     if (fnJetsPtCut==0)
2295     {
2296         jet_area_total=0.;
2297     }
2298     else
2299     {
2300         for (i=0;i<fnJetsPtCut;i++)
2301         {
2302             AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fJetPtCutID[i]);
2303             jet_area_total+=AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta());
2304         }
2305     }
2306
2307     // Calculate Rho
2308     EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-jet_area_total);
2309     fRhoTotal=EMCal_rho;
2310
2311     // Fill Background Histogram
2312     fhEMCalBckg1B->Fill(EMCal_rho*TMath::Pi()*TMath::Power(fJetR,2));
2313     fhRho1B->Fill(EMCal_rho,fEventCentrality);
2314     fpRho1B->Fill(fEventCentrality,EMCal_rho);
2315     FillFullCorrJetPt(fhJetTPt1B,EMCal_rho,kFALSE);
2316     FillFullCorrJetPt(fhJetTPt1BSignal,EMCal_rho,kTRUE);
2317
2318     FillFullCorrJetPt(fhJetTPtCen1B,EMCal_rho,kFALSE);
2319     FillFullCorrJetPt(fhJetTPtCen1BSignal,EMCal_rho,kTRUE);
2320     // Fill delta pT
2321     if (fEventCentrality<=20)
2322     {
2323         FillFullCorrJetPt(fh020JetTPt1B,EMCal_rho,kFALSE);
2324         FillFullCorrJetPt(fh020JetTPt1BSignal,EMCal_rho,kTRUE);
2325         FillBckgFlucDeltaPt(fhDeltaPt1B,EMCal_rho);
2326         fh020Rho1B->Fill(EMCal_rho);
2327     }
2328 }
2329
2330 void AliAnalysisTaskFullpAJets::Method1C()
2331 {
2332     const Double_t psi_ref=0.5*(45/360.)*2*TMath::Pi(); //Input the total acceptance within the paraenthesis to be +/- psi_ref
2333     Int_t i;
2334     Double_t E_tracks_total=0.;
2335     Double_t E_caloclusters_total=0.;
2336     Double_t EMCal_rho=0.;
2337     Double_t psi;
2338
2339     if (fPtFullMaxID !=-1)
2340     {
2341         AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fPtFullMaxID);
2342
2343         if (myJet->Area()>=fJetAreaThreshold)
2344         {
2345             //  Loop over all tracks
2346             for (i=0;i<fnTracks;i++)
2347             {
2348                 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2349                 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2350                 {
2351                     if ((vtrack->Eta()>=(myJet->Eta()+fJetR)) || (vtrack->Eta()<=(myJet->Eta()-fJetR)))
2352                     {
2353                         psi=TMath::ATan((vtrack->Phi()-myJet->Phi())/(vtrack->Eta()-myJet->Eta()));
2354                         if ((psi>=(-1*psi_ref)) && (psi<=psi_ref))
2355                         {
2356                             //  The Track made the Cut!!
2357                             E_tracks_total+=vtrack->Pt();
2358                         }
2359                     }
2360                 }
2361             }
2362
2363             //  Loop over all caloclusters
2364             for (i=0;i<fnClusters;i++)
2365             {
2366                 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2367                 TLorentzVector *cluster_vec = new TLorentzVector;
2368                 vcluster->GetMomentum(*cluster_vec,fvertex);
2369
2370                 if ((cluster_vec->Eta()>=(myJet->Eta()+fJetR)) || (cluster_vec->Eta()<=(myJet->Eta()-fJetR)))
2371                 {
2372                     psi=TMath::ATan((cluster_vec->Phi()-myJet->Phi())/(cluster_vec->Eta()-myJet->Eta()));
2373                     if ((psi>=(-1*psi_ref)) && (psi<=psi_ref))
2374                     {
2375                         //  The CaloCluster made the Cut!!
2376                         E_caloclusters_total+=vcluster->E();
2377                     }
2378                 }
2379             }
2380
2381             //  Calculate the mean Background density
2382             EMCal_rho=(E_tracks_total+E_caloclusters_total)/TransverseArea(fJetR,psi_ref,myJet->Phi(),myJet->Eta());
2383
2384             //  Fill histograms
2385             fhEMCalBckg1C->Fill(EMCal_rho*TMath::Pi()*fJetR*fJetR);
2386             fhRho1C->Fill(EMCal_rho,fEventCentrality);
2387             fhJetTPt1C->Fill(myJet->Pt()-EMCal_rho*myJet->Area());
2388             FillFullCorrJetPt(fhJetTPt1C,EMCal_rho,kFALSE);
2389             FillFullCorrJetPt(fhJetTPtCen1C,EMCal_rho,kFALSE);
2390
2391
2392             if (fEventCentrality<=20)
2393             {
2394                 FillFullCorrJetPt(fh020JetTPt1C,EMCal_rho,kFALSE);
2395             }
2396         }
2397     }
2398 }
2399
2400 void AliAnalysisTaskFullpAJets::Method2A()
2401 {
2402     Int_t i;
2403
2404     for (i=0;i<fnAKTFullJets;i++)
2405     {
2406         if (fInEMCalFull[i]==kTRUE)
2407         {
2408             AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
2409             if (myJet->Area()>=fJetAreaThreshold)
2410             {
2411                 fhEMCalCenJet2A->Fill(myJet->Pt(),fEventCentrality);
2412                 fhEMCalJet2A->Fill(myJet->Pt());
2413                 if (fEventCentrality<=20)
2414                 {
2415                     fh020EMCalJet2A->Fill(myJet->Pt());
2416                 }
2417             }
2418         }
2419     }
2420 }
2421
2422 void AliAnalysisTaskFullpAJets::Method2B()
2423 {
2424     Int_t i;
2425
2426     Double_t E_tracks_total=0.0;
2427     Double_t E_caloclusters_total=0.0;
2428     Double_t EMCal_rho=0.0;
2429     Double_t E_tracks_core_total=0.0;
2430     Double_t E_caloclusters_core_total=0.0;
2431     Double_t EMCal_core_rho=0.0;
2432     Double_t RCore=0.4;
2433
2434     //  Loop over all tracks
2435     for (i=0;i<fnTracks;i++)
2436     {
2437         AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2438         if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2439         {
2440             E_tracks_total+=vtrack->Pt();
2441             if (InsideRect(vtrack->Phi(),fEMCalPhiMin+RCore,fEMCalPhiMax-RCore,vtrack->Eta(),fEMCalEtaMin+RCore,fEMCalEtaMax-RCore)==kTRUE)
2442             {
2443                 E_tracks_core_total+=vtrack->Pt();
2444             }
2445         }
2446     }
2447
2448     //  Loop over all caloclusters
2449     for (i=0;i<fnClusters;i++)
2450     {
2451         AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2452         E_caloclusters_total+=vcluster->E();
2453         TLorentzVector *cluster_vec = new TLorentzVector;
2454         vcluster->GetMomentum(*cluster_vec,fvertex);
2455         if (InsideRect(cluster_vec->Phi(),fEMCalPhiMin+RCore,fEMCalPhiMax-RCore,cluster_vec->Eta(),fEMCalEtaMin+RCore,fEMCalEtaMax-RCore)==kTRUE)
2456         {
2457             E_caloclusters_core_total+=vcluster->E();
2458         }
2459         delete cluster_vec;
2460     }
2461
2462     //  Calculate the mean Background density
2463     EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2464     EMCal_core_rho=(E_tracks_core_total+E_caloclusters_core_total)/((fEMCalPhiTotal-2*RCore)*(fEMCalEtaTotal-2*RCore));
2465
2466     //Fill Background Cluster Histogram
2467     fhEMCalBckg2B->Fill(EMCal_rho*TMath::Pi()*TMath::Power(fJetR,2));
2468     fhRho2B->Fill(EMCal_rho,fEventCentrality);
2469     fpRho2B->Fill(fEventCentrality,EMCal_rho);
2470     FillFullCorrJetPt(fhJetTPt2B,EMCal_rho,kFALSE);
2471     FillFullCorrJetPt(fhJetTPtCen2B,EMCal_rho,kFALSE);
2472     fhDeltaRho0DiJet->Fill(fRhoAkTTotal-EMCal_rho);
2473     fhRho2BCore->Fill(EMCal_core_rho,fEventCentrality);
2474     fpRho2BCore->Fill(fEventCentrality,EMCal_core_rho);
2475
2476     if (fEventCentrality<=20)
2477     {
2478         FillFullCorrJetPt(fh020JetTPt2B,EMCal_rho,kFALSE);
2479         FillBckgFlucDeltaPt(fhDeltaPt2B,EMCal_rho);
2480         fh020Rho2B->Fill(EMCal_rho);
2481         fh020Rho2BCore->Fill(EMCal_core_rho);
2482     }
2483 }
2484
2485 void AliAnalysisTaskFullpAJets::Method3(Bool_t EMCalOn)
2486 {
2487     Int_t i,j;
2488     Bool_t track_away_from_jet;
2489     Double_t E_tracks_total=0.0;
2490     Double_t TPC_rho=0.0;
2491     Double_t jet_area_total=0.0;
2492
2493     // Calculate charged background density in events with no signal jets
2494     if (fnJetsChargedPtCut==0)
2495     {
2496         //  Loop over all tracks
2497         for (i=0;i<fnTracks;i++)
2498         {
2499             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2500             E_tracks_total+=vtrack->Pt();
2501         }
2502         TPC_rho=E_tracks_total/fTPCArea;
2503
2504         fhRho3NoJets->Fill(TPC_rho,fEventCentrality);
2505         fpRho3NoJets->Fill(fEventCentrality,TPC_rho);
2506         if (fEventCentrality<=20)
2507         {
2508             fh020Rho3NoJets->Fill(TPC_rho);
2509         }
2510         E_tracks_total=0.0;
2511         TPC_rho=0.0;
2512     }
2513
2514     // First, sum all tracks within the TPC that are away from jet(s) above Pt Threshold
2515     fRhoCharged=0;
2516     for (i=0;i<fnTracks;i++)
2517     {
2518         AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
2519         if (fnJetsChargedPtCut<1)
2520         {
2521             E_tracks_total+=vtrack->Pt();
2522         }
2523         else
2524         {
2525             track_away_from_jet=kTRUE;
2526             j=0;
2527             TLorentzVector *track_vec = new TLorentzVector;
2528             track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2529             while (track_away_from_jet==kTRUE && j<fnJetsChargedPtCut)
2530             {
2531                 TLorentzVector *jet_vec= new TLorentzVector;
2532                 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fJetPtChargedCutID[j]);
2533                 myJet->GetMom(*jet_vec);
2534                 if (track_vec->DeltaR(*jet_vec)<=fJetR)
2535                 {
2536                     track_away_from_jet=kFALSE;
2537                 }
2538                 delete jet_vec;
2539                 j++;
2540             }
2541             if (track_away_from_jet==kTRUE)
2542             {
2543                 E_tracks_total+=vtrack->Pt();
2544             }
2545             delete track_vec;
2546         }
2547     }
2548
2549     // Determine area of all Jets that are within the EMCal
2550     if (fnJetsChargedPtCut==0)
2551     {
2552         jet_area_total=0.;
2553     }
2554     else
2555     {
2556         for (i=0;i<fnJetsChargedPtCut;i++)
2557         {
2558             AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fJetPtChargedCutID[i]);
2559             jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
2560         }
2561     }
2562
2563     //Calculate Rho
2564     TPC_rho=(E_tracks_total)/(fTPCArea-jet_area_total);
2565     fRhoCharged=TPC_rho;
2566
2567     //Fill Background Histogram
2568     fhRho3->Fill(TPC_rho,fEventCentrality);
2569     fpRho3->Fill(fEventCentrality,TPC_rho);
2570
2571     //Fill "True" Jet Pt Spectrum
2572     for (i=0;i<fnAKTChargedJets;i++)
2573     {
2574         if (fInTPCChargedFull[i]==kTRUE)
2575         {
2576             AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(i);
2577             if (myJet->Area()>=fJetAreaThreshold && myJet->Pt()>=fTPCJetThreshold)
2578             {
2579                 fhJetTPt3->Fill(myJet->Pt()-(TPC_rho*myJet->Area()));
2580                 fhJetTPtCen3->Fill(myJet->Pt()-(TPC_rho*myJet->Area()),fEventCentrality);
2581                 if (fEventCentrality<=20)
2582                 {
2583                     fh020JetTPt3->Fill(myJet->Pt()-(TPC_rho*myJet->Area()));
2584                 }
2585             }
2586         }
2587     }
2588
2589     // Fill Background Scaling factor profile.
2590     if (EMCalOn==kTRUE && fRhoCharged!=0)
2591     {
2592         fpRhoScale->Fill(fEventCentrality,(fRhoTotal/fRhoCharged));
2593     }
2594
2595     if (fEventCentrality<=20)
2596     {
2597         fh020Rho3->Fill(TPC_rho);
2598     }
2599 }
2600
2601 void AliAnalysisTaskFullpAJets::Method3DiJet()
2602 {
2603     Int_t i;
2604     Double_t E_tracks_total=0.0;
2605     Double_t TPC_rho=0.0;
2606     Double_t jet_area_total=0.0;
2607     Double_t jet_pT_total=0.0;
2608
2609     //  Loop over all tracks
2610     for (i=0;i<fnTracks;i++)
2611     {
2612         AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2613         E_tracks_total+=vtrack->Pt();
2614     }
2615
2616     AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fPtChargedMaxID);
2617     AliEmcalJet *mybJet =(AliEmcalJet*) fmyAKTChargedJets->At(fChargedBackJetID);
2618
2619     jet_area_total=myhJet->Area()+mybJet->Area();
2620     jet_pT_total=myhJet->Pt()+mybJet->Pt();
2621     TPC_rho=(E_tracks_total-jet_pT_total)/(fTPCArea-jet_area_total);
2622
2623     fhRho3DiJets->Fill(TPC_rho,fEventCentrality);
2624     fpRho3DiJets->Fill(fEventCentrality,TPC_rho);
2625     if (fEventCentrality<=20)
2626     {
2627         fh020Rho3DiJets->Fill(TPC_rho);
2628     }
2629 }
2630
2631 void AliAnalysisTaskFullpAJets::Method3Perp()
2632 {
2633     Int_t i;
2634     Double_t E_tracks_total=0.0;
2635     Double_t TPC_rho=0.0;
2636     Double_t phi_perp=0.0;  // By construction, this angle must be between 90 and 270
2637     const Double_t delta_phi_prep=(30/360.0)*2*TMath::Pi();  // Half angle...
2638
2639     AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fPtChargedMaxID);
2640     AliEmcalJet *mybJet =(AliEmcalJet*) fmyAKTChargedJets->At(fChargedBackJetID);
2641
2642     phi_perp=(0.5*(myhJet->Phi()+mybJet->Phi())/360.0)*2*TMath::Pi();
2643
2644     //  Loop over all tracks
2645     for (i=0;i<fnTracks;i++)
2646     {
2647         AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2648         // Select only those tracks that are within the "strip"
2649         if (vtrack->Phi() >= (phi_perp-delta_phi_prep) && vtrack->Phi() <= (phi_perp+delta_phi_prep))
2650         {
2651             E_tracks_total+=vtrack->Pt();
2652         }
2653     }
2654
2655     TPC_rho=(E_tracks_total)/(fTPCEtaMax*2*delta_phi_prep);
2656
2657     fhRho3Perp->Fill(TPC_rho,fEventCentrality);
2658     fpRho3Perp->Fill(fEventCentrality,TPC_rho);
2659     if (fEventCentrality<=20)
2660     {
2661         fh020Rho3Perp->Fill(TPC_rho);
2662     }
2663 }
2664
2665 void AliAnalysisTaskFullpAJets::JetPtChargedProfile()
2666 {
2667     Int_t i,j;
2668     Double_t delta_R;
2669     Double_t ED_pT[fEDProfileRBins];
2670
2671     for (i=0;i<fnJetsChargedPtCut;i++)
2672     {
2673         AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fJetPtChargedCutID[i]);
2674         if (InsideRect(myJet->Phi(),fTPCPhiMin,fTPCPhiMax,myJet->Eta(),fTPCEtaMin+fEDProfileRUp,fTPCEtaMax-fEDProfileRUp)==kTRUE)
2675         {
2676             for (j=0;j<fEDProfileRBins;j++)
2677             {
2678                 ED_pT[j]=0;
2679             }
2680             TLorentzVector *jet_vec= new TLorentzVector;
2681             myJet->GetMom(*jet_vec);
2682             // Sum all tracks in concentric rings around jet vertex
2683             for (j=0;j<fnTracks;j++)
2684             {
2685                 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2686                 TLorentzVector *track_vec = new TLorentzVector;
2687                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2688                 delta_R=jet_vec->DeltaR(*track_vec);
2689                 if (delta_R<=fEDProfileRUp)
2690                 {
2691                     ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vtrack->Pt();
2692                 }
2693                 delete track_vec;
2694             }
2695
2696             //cout<<"Event Centrality:"<<fEventCentrality<<endl;
2697             //cout<<endl<<endl<<"Event Centrality bin:"<<TMath::FloorNint(fEventCentrality/10.)<<endl;
2698             for (j=0;j<fEDProfileRBins;j++)
2699             {
2700                 ED_pT[j]/=TMath::Pi()*TMath::Power((fEDProfileRUp/fEDProfileRBins),2)*(2*j+1);
2701                 //cout<<"Strip "<<j<<"  ED="<<ED_pT[j]<<endl;
2702                 fpChargedJetEDProfile[TMath::FloorNint(fEventCentrality/10.)]->Fill(myJet->Pt(),myJet->Eta(),(fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2703                 if (fEventCentrality<=20)
2704                 {
2705                     fpChargedJetRProfile[4+TMath::FloorNint(myJet->Eta()*10.)]->Fill((fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2706                 }
2707             }
2708             delete jet_vec;
2709         }
2710     }
2711 }
2712
2713 void AliAnalysisTaskFullpAJets::JetPtFullProfile()
2714 {
2715     Int_t i,j;
2716     Double_t delta_R;
2717     Double_t ED_pT[fEDProfileRBins];
2718
2719     for (i=0;i<fnJetsPtCut;i++)
2720     {
2721         AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fJetPtCutID[i]);
2722         if (InsideRect(myJet->Phi(),fEMCalPhiMin+fEDProfileRUp,fEMCalPhiMax-fEDProfileRUp,myJet->Eta(),fEMCalEtaMin+fEDProfileRUp,fEMCalEtaMax-fEDProfileRUp)==kTRUE)
2723         {
2724             for (j=0;j<fEDProfileRBins;j++)
2725             {
2726                 ED_pT[j]=0;
2727             }
2728             TLorentzVector *jet_vec= new TLorentzVector;
2729             myJet->GetMom(*jet_vec);
2730             // Sum all tracks in concentric rings around jet vertex
2731             for (j=0;j<fnTracks;j++)
2732             {
2733                 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2734                 TLorentzVector *track_vec = new TLorentzVector;
2735                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2736                 delta_R=jet_vec->DeltaR(*track_vec);
2737                 if (delta_R<=fEDProfileRUp)
2738                 {
2739                     ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vtrack->Pt();
2740                 }
2741                 delete track_vec;
2742             }
2743
2744             // Sum all clusters in concentric rings around jet vertex
2745             for (j=0;j<fnClusters;j++)
2746             {
2747                 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
2748                 TLorentzVector *cluster_vec = new TLorentzVector;
2749                 vcluster->GetMomentum(*cluster_vec,fvertex);
2750                 delta_R=jet_vec->DeltaR(*cluster_vec);
2751                 if (delta_R<=fEDProfileRUp)
2752                 {
2753                     ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vcluster->E();
2754                 }
2755                 delete cluster_vec;
2756             }
2757
2758             for (j=0;j<fEDProfileRBins;j++)
2759             {
2760                 ED_pT[j]/=TMath::Pi()*TMath::Power((fEDProfileRUp/fEDProfileRBins),2)*(2*j+1);
2761                 fpJetEDProfile[TMath::FloorNint(fEventCentrality/10.)]->Fill(myJet->Pt(),myJet->Eta(),(fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2762                 // Fill profile if a "most" central event (0-20%)
2763                 if (fEventCentrality<=20)
2764                 {
2765                     fpJetRProfile[2+TMath::FloorNint(myJet->Eta()*10.)]->Fill((fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2766                 }
2767             }
2768             delete jet_vec;
2769
2770             // Fill constituent histogram
2771             for (j=0;j<myJet->GetNumberOfTracks();j++)
2772             {
2773                 AliVParticle* vparticle = (AliVParticle*) myJet->TrackAt(j,fmyTracks);
2774                 fhJetConstituentPt->Fill(myJet->Pt(),vparticle->Pt());
2775             }
2776
2777             for (j=0;j<myJet->GetNumberOfClusters();j++)
2778             {
2779                 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(myJet->ClusterAt(j));
2780                 fhJetConstituentPt->Fill(myJet->Pt(),vcluster->E());
2781             }
2782         }
2783     }
2784 }
2785
2786 void AliAnalysisTaskFullpAJets::JetPtEtaProfile()
2787 {
2788     Int_t i,j;
2789     Double_t eta;
2790     Double_t delta_eta;
2791     Double_t Eta_pT[fEtaProfileBins];
2792     Double_t Eta_abs_pT[Int_t(0.5*fEtaProfileBins)];
2793
2794     for (i=0;i<fnJetsPtCut;i++)
2795     {
2796         AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fJetPtCutID[i]);
2797         if (IsInEMCal(myJet->Phi(),myJet->Eta())==kTRUE)
2798         {
2799             for (j=0;j<fEtaProfileBins;j++)
2800             {
2801                 Eta_pT[j]=0;
2802                 Eta_abs_pT[j]=0;
2803             }
2804
2805             // Sum all tracks in strips of eta away from the jet vertex
2806             for (j=0;j<fnTracks;j++)
2807             {
2808                 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2809                 eta=vtrack->Eta();
2810                 delta_eta=TMath::Abs(vtrack->Eta()-myJet->Eta());
2811                 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2812                 {
2813                     Eta_pT[Int_t(0.5*fEtaProfileBins)+TMath::FloorNint(10*eta)]+=vtrack->Pt();
2814                     Eta_abs_pT[TMath::FloorNint(10*delta_eta)]+=vtrack->Pt();
2815                 }
2816             }
2817
2818             // Sum all clusters in strips of eta away from the jet vertex
2819             for (j=0;j<fnClusters;j++)
2820             {
2821                 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
2822                 TLorentzVector *cluster_vec = new TLorentzVector;
2823                 vcluster->GetMomentum(*cluster_vec,fvertex);
2824                 eta=cluster_vec->Eta();
2825                 delta_eta=TMath::Abs(cluster_vec->Eta()-myJet->Eta());
2826                 Eta_pT[Int_t(0.5*fEtaProfileBins)+TMath::FloorNint(10*eta)]+=vcluster->E();
2827                 Eta_abs_pT[TMath::FloorNint(10*delta_eta)]+=vcluster->E();
2828                 delete cluster_vec;
2829             }
2830
2831             for (j=0;j<fEtaProfileBins;j++)
2832             {
2833                 Eta_pT[j]/=0.1*fEMCalPhiTotal;
2834                 // Fill profile if a "most" central event (0-20%)
2835                 if (j<(10*(fEMCalEtaMax-TMath::Abs(myJet->Eta()))))
2836                 {
2837                     Eta_abs_pT[j]/=0.2*fEMCalPhiTotal;
2838                 }
2839                 else
2840                 {
2841                     Eta_abs_pT[j]/=0.1*fEMCalPhiTotal;
2842                 }
2843                 // Fill profile if a "most" central event (0-20%)
2844                 if (fEventCentrality<=20)
2845                 {
2846                     fpJetAbsEtaProfile[7+TMath::FloorNint(myJet->Eta()*10.)]->Fill(0.1*j,Eta_abs_pT[j]);
2847                     fpJetEtaProfile[7+TMath::FloorNint(myJet->Eta()*10.)]->Fill(0.1*(j-7),Eta_pT[j]);
2848                 }
2849             }
2850         }
2851     }
2852 }
2853
2854 void AliAnalysisTaskFullpAJets::FillFullCorrJetPt(TH1D *myHisto,Double_t rho, Bool_t signal_cut)
2855 {
2856     Int_t i;
2857     // Fill "True" Jet Pt Spectrum. Always demand that jet area is greater then area threshold.
2858     for (i=0;i<fnAKTFullJets;i++)
2859     {
2860         if (fInEMCalFull[i]==kTRUE)
2861         {
2862             AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(i);
2863             if (myJet->Area()>=fJetAreaThreshold)
2864             {
2865                 if (signal_cut==kTRUE && myJet->Pt()>=fEMCalJetThreshold)
2866                 {
2867                     myHisto->Fill(myJet->Pt()-(rho*myJet->Area()));
2868                 }
2869                 else if (signal_cut==kFALSE)
2870                 {
2871                     myHisto->Fill(myJet->Pt()-(rho*myJet->Area()));
2872                 }
2873             }
2874         }
2875     }
2876 }
2877
2878 void AliAnalysisTaskFullpAJets::FillFullCorrJetPt(TH2D *myHisto,Double_t rho, Bool_t signal_cut)
2879 {
2880     Int_t i;
2881     // Fill "True" Jet Pt Spectrum. Always demand that jet area is greater then area threshold.
2882     for (i=0;i<fnAKTFullJets;i++)
2883     {
2884         if (fInEMCalFull[i]==kTRUE)
2885         {
2886             AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(i);
2887             if (myJet->Area()>=fJetAreaThreshold)
2888             {
2889                 if (signal_cut==kTRUE && myJet->Pt()>=fEMCalJetThreshold)
2890                 {
2891                     myHisto->Fill(myJet->Pt()-(rho*myJet->Area()),fEventCentrality);
2892                 }
2893                 else if (signal_cut==kFALSE)
2894                 {
2895                     myHisto->Fill(myJet->Pt()-(rho*myJet->Area()),fEventCentrality);
2896                 }
2897             }
2898         }
2899     }
2900 }
2901
2902 void AliAnalysisTaskFullpAJets::FillFullDeltaRho(TH1D *myHisto,Double_t delta_rho,Bool_t signal_cut)
2903 {
2904     Int_t i;
2905     // Fill "True" Jet Pt Spectrum. Always demand that jet area is greater then area threshold.
2906     for (i=0;i<fnAKTFullJets;i++)
2907     {
2908         if (fInEMCalFull[i]==kTRUE)
2909         {
2910             AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(i);
2911             if (myJet->Area()>=fJetAreaThreshold)
2912             {
2913                 if (signal_cut==kTRUE && myJet->Pt()>=fEMCalJetThreshold)
2914                 {
2915                     myHisto->Fill(delta_rho);
2916                 }
2917                 else if (signal_cut==kFALSE)
2918                 {
2919                     myHisto->Fill(delta_rho);
2920                 }
2921             }
2922         }
2923     }
2924 }
2925
2926 void AliAnalysisTaskFullpAJets::FillBckgFlucDeltaPt(TH1D *myHisto, Double_t rho)
2927 {
2928     Int_t i;
2929
2930     for (i=0;i<fnBckgClusters;i++)
2931     {
2932         myHisto->Fill(fRCBckgFluc[i]-rho*TMath::Pi()*fJetR*fJetR);
2933     }
2934 }
2935
2936
2937 void AliAnalysisTaskFullpAJets::DeleteArrays(Bool_t EMCalOn)
2938 {
2939     delete [] fJetPtChargedCutID;
2940     delete [] fInTPCChargedFull;
2941     delete [] fJetkTTPCFullID;
2942     if (EMCalOn==kTRUE)
2943     {
2944         delete [] fJetPtCutID;
2945         delete [] fJetPtTPCCutID;
2946         delete [] fJetPtTotalCutID;
2947         delete [] fJetkTEMCalFullID;
2948         delete [] fInEMCal;
2949         delete [] fInEMCalFull;
2950         delete [] fInTPCFull;
2951     }
2952 }
2953
2954 /////////////////////////////////////////////////////////////////////////////////////////
2955 /////////////////     User Defined Functions      ///////////////////////////////////////
2956 /////////////////////////////////////////////////////////////////////////////////////////
2957
2958 Bool_t AliAnalysisTaskFullpAJets::IsDiJetEvent()
2959 {
2960     // Determine if event contains a di-jet within the detector. Uses charged jets.
2961     // Requires the delta phi of the jets to be 180 +/- 15 degrees.
2962     // Requires both jets to be outside of the EMCal
2963     // Requires both jets to be signal jets
2964
2965     Int_t i;
2966     const Double_t dijet_delta_phi=(180/360.)*2*TMath::Pi();
2967     const Double_t dijet_phi_acceptance=0.5*(30/360.)*2*TMath::Pi(); //Input the total acceptance within the paraenthesis to be +/- dijet_phi_acceptance
2968     Double_t dummy_phi=0.0;
2969     Double_t dijet_asymmetry=0.0;
2970     Double_t delta_phi=0.0;
2971     fChargedBackJetID=-1;
2972
2973     if (fnJetsChargedPtCut>1)
2974     {
2975         AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fPtChargedMaxID);
2976         i=0;
2977         if (IsInTPCFull(fJetR,myhJet->Phi(),myhJet->Eta())==kFALSE)
2978         {
2979             return kFALSE;
2980         }
2981         while (i<fnJetsChargedPtCut)
2982         {
2983             AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fJetPtChargedCutID[i]);
2984             dummy_phi=TMath::Min(TMath::Abs(myhJet->Phi()-myJet->Phi()),2*TMath::Pi()-TMath::Abs(myhJet->Phi()-myJet->Phi()));
2985             if ((dummy_phi>(dijet_delta_phi-dijet_phi_acceptance)) && (IsInTPCFull(fJetR,myJet->Phi(),myJet->Eta())==kTRUE))
2986             {
2987                 fChargedBackJetID=fJetPtChargedCutID[i];
2988                 fnDiJetEvents++;
2989                 dijet_asymmetry=(myhJet->Pt()-myJet->Pt())/(myhJet->Pt()+myJet->Pt());
2990                 fhDiJetCenAsy->Fill(dijet_asymmetry,fEventCentrality);
2991                 // Add Anti-kT Plots here...
2992                 if (fPtFullMaxID!=-1)
2993                 {
2994                     AliEmcalJet *myTestJet = (AliEmcalJet*) fmyAKTFullJets->At(fPtFullMaxID);
2995
2996                     fhDiJetEMCalLeadingPt->Fill(myTestJet->Pt());
2997                     delta_phi=TMath::Min(TMath::Abs(myhJet->Phi()-myTestJet->Phi()),2*TMath::Pi()-TMath::Abs(myhJet->Phi()-myTestJet->Phi()));
2998                     fhDiJetEMCalLeadingDeltaPhi->Fill(delta_phi);
2999                 }
3000                 if (fEventCentrality<=20)
3001                 {
3002                     fh020DiJetAsy->Fill(dijet_asymmetry);
3003                 }
3004                 return kTRUE;
3005             }
3006             i++;
3007         }
3008     }
3009     return kFALSE;
3010 }
3011
3012 Bool_t AliAnalysisTaskFullpAJets::InsideRect(Double_t phi,Double_t phi_min,Double_t phi_max,Double_t eta,Double_t eta_min,Double_t eta_max)
3013 {
3014     if (phi>phi_min && phi<phi_max)
3015     {
3016         if (eta>eta_min && eta<eta_max)
3017         {
3018             return kTRUE;
3019         }
3020     }
3021     return kFALSE;
3022 }
3023
3024 Bool_t AliAnalysisTaskFullpAJets::IsInEMCal(Double_t phi,Double_t eta)
3025 {
3026     return InsideRect(phi,fEMCalPhiMin,fEMCalPhiMax,eta,fEMCalEtaMin,fEMCalEtaMax);
3027 }
3028
3029 Bool_t AliAnalysisTaskFullpAJets::IsInEMCalFull(Double_t r,Double_t phi,Double_t eta)
3030 {
3031     return InsideRect(phi,fEMCalPhiMin+r,fEMCalPhiMax-r,eta,fEMCalEtaMin+r,fEMCalEtaMax-r);
3032 }
3033
3034 Bool_t AliAnalysisTaskFullpAJets::IsInEMCalPart(Double_t r,Double_t phi,Double_t eta)
3035 {
3036     return InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
3037 }
3038
3039 Bool_t AliAnalysisTaskFullpAJets::IsInTPCFull(Double_t r,Double_t phi,Double_t eta)
3040 {
3041     Bool_t in_EMCal= InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
3042     Bool_t in_TPC= InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
3043
3044     if (in_EMCal==kFALSE && in_TPC==kTRUE)
3045     {
3046         return kTRUE;
3047     }
3048     return kFALSE;
3049 }
3050
3051 Bool_t AliAnalysisTaskFullpAJets::IsInTPC(Double_t r,Double_t phi,Double_t eta,Bool_t Complete)
3052 {
3053     if (Complete==kTRUE)
3054     {
3055         return InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
3056     }
3057     return InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin,fTPCEtaMax);
3058 }
3059
3060 Double_t AliAnalysisTaskFullpAJets::AreaWithinTPC(Double_t r,Double_t eta)
3061 {
3062     Double_t z;
3063     if (eta<(fTPCEtaMin+r))
3064     {
3065         z=eta-fTPCEtaMin;
3066     }
3067     else if(eta>(fTPCEtaMax-r))
3068     {
3069         z=fTPCEtaMax-eta;
3070     }
3071     else
3072     {
3073         z=r;
3074     }
3075     return r*r*TMath::Pi()-AreaEdge(r,z);
3076 }
3077
3078 Double_t AliAnalysisTaskFullpAJets::AreaWithinEMCal(Double_t r,Double_t phi,Double_t eta)
3079 {
3080     Double_t x,y;
3081
3082     if (phi<(fEMCalPhiMin-r) || phi>(fEMCalPhiMax+r))
3083     {
3084         x=-r;
3085     }
3086     else if (phi<(fEMCalPhiMin+r))
3087     {
3088         x=phi-fEMCalPhiMin;
3089     }
3090     else if (phi>(fEMCalPhiMin+r) && phi<(fEMCalPhiMax-r))
3091     {
3092         x=r;
3093     }
3094     else
3095     {
3096         x=fEMCalPhiMax-phi;
3097     }
3098
3099     if (eta<(fEMCalEtaMin-r) || eta>(fEMCalEtaMax+r))
3100     {
3101         y=-r;
3102     }
3103     else if (eta<(fEMCalEtaMin+r))
3104     {
3105         y=eta-fEMCalEtaMin;
3106     }
3107     else if (eta>(fEMCalEtaMin+r) && eta<(fEMCalEtaMax-r))
3108     {
3109         y=r;
3110     }
3111     else
3112     {
3113         y=fEMCalEtaMax-eta;
3114     }
3115
3116     if (x>=0 && y>=0)
3117     {
3118         if (TMath::Sqrt(x*x+y*y)>=r)
3119         {
3120             return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
3121         }
3122         return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y)+AreaOverlap(r,x,y);
3123     }
3124     else if ((x>=r && y<0) || (y>=r && x<0))
3125     {
3126         return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
3127     }
3128     else if (x>0 && x<r && y<0)
3129     {
3130         Double_t a=TMath::Sqrt(r*r-x*x);
3131         Double_t b=TMath::Sqrt(r*r-y*y);
3132         if ((x-b)>0)
3133         {
3134             return r*r*TMath::ASin(b/r)+y*b;
3135         }
3136         else
3137         {
3138             return 0.5*x*a+0.5*r*r*TMath::ASin(x/r)+0.5*y*b+x*y+0.5*r*r*TMath::ASin(b/r);
3139         }
3140     }
3141     else if (y>0 && y<r && x<0)
3142     {
3143         Double_t a=TMath::Sqrt(r*r-x*x);
3144         Double_t b=TMath::Sqrt(r*r-y*y);
3145         if ((y-a)>0)
3146         {
3147             return r*r*TMath::ASin(a/r)+x*a;
3148         }
3149         else
3150         {
3151             return 0.5*y*b+0.5*r*r*TMath::ASin(y/r)+0.5*x*a+x*y+0.5*r*r*TMath::ASin(a/r);
3152         }
3153     }
3154     else
3155     {
3156         Double_t a=TMath::Sqrt(r*r-x*x);
3157         Double_t b=TMath::Sqrt(r*r-y*y);
3158         if ((x+b)<0)
3159         {
3160             return 0;
3161         }
3162         else
3163         {
3164             return 0.5*x*a+0.5*r*r*TMath::ASin(x/r)+0.5*y*b+x*y+0.5*r*r*TMath::ASin(b/r);
3165         }
3166     }
3167 }
3168
3169 Double_t AliAnalysisTaskFullpAJets::AreaEdge(Double_t r,Double_t z)
3170 {
3171     Double_t a=TMath::Sqrt(r*r-z*z);
3172     return r*r*TMath::ASin(a/r)-a*z;
3173 }
3174
3175 Double_t AliAnalysisTaskFullpAJets::AreaOverlap(Double_t r,Double_t x,Double_t y)
3176 {
3177     Double_t a=TMath::Sqrt(r*r-x*x);
3178     Double_t b=TMath::Sqrt(r*r-y*y);
3179     return x*y-0.5*(x*a+y*b)+0.5*r*r*(TMath::ASin(b/r)-TMath::ASin(x/r));
3180 }
3181
3182 Double_t AliAnalysisTaskFullpAJets::TransverseArea(Double_t r,Double_t psi0,Double_t phi,Double_t eta)
3183 {
3184     Double_t area_left,area_right;
3185     Double_t eta_a,eta_b,eta_up,eta_down;
3186
3187     Double_t u=eta-fEMCalEtaMin;
3188     Double_t v=fEMCalEtaMax-eta;
3189
3190     Double_t phi1=phi+u*TMath::Tan(psi0);
3191     Double_t phi2=phi-u*TMath::Tan(psi0);
3192     Double_t phi3=phi+v*TMath::Tan(psi0);
3193     Double_t phi4=phi-v*TMath::Tan(psi0);
3194
3195     //Calculate the Left side area
3196     if (phi1>=fEMCalPhiMax)
3197     {
3198         eta_a=eta-u*((fEMCalPhiMax-phi)/(phi1-phi));
3199     }
3200     if (phi2<=fEMCalPhiMin)
3201     {
3202         eta_b=eta-u*((phi-fEMCalPhiMin)/(phi-phi2));
3203     }
3204
3205     if ((phi1>=fEMCalPhiMax) && (phi2<=fEMCalPhiMin))
3206     {
3207         eta_up=TMath::Max(eta_a,eta_b);
3208         eta_down=TMath::Min(eta_a,eta_b);
3209
3210         area_left=(eta_down-fEMCalEtaMin)*fEMCalPhiTotal + 0.5*(fEMCalPhiTotal+2*(eta-eta_up)*TMath::Tan(psi0))*(eta_up-eta_down) + (eta-eta_up+r)*TMath::Tan(psi0)*(eta-eta_up-r);
3211     }
3212     else if (phi1>=fEMCalPhiMax)
3213     {
3214         area_left=0.5*(fEMCalPhiMax-phi2+2*(eta-eta_a)*TMath::Tan(psi0))*(eta_a-fEMCalEtaMin) + (eta-eta_a+r)*TMath::Tan(psi0)*(eta-eta_a-r);
3215     }
3216     else if (phi2<=fEMCalPhiMin)
3217     {
3218         area_left=0.5*(phi1-fEMCalPhiMin+2*(eta-eta_b)*TMath::Tan(psi0))*(eta_b-fEMCalEtaMin) + (eta-eta_b+r)*TMath::Tan(psi0)*(eta-eta_b-r);
3219     }
3220     else
3221     {
3222         area_left=0.5*(phi1-phi2+2*r*TMath::Tan(psi0))*(u-r);
3223     }
3224
3225     //Calculate the Right side area
3226     if (phi3>=fEMCalPhiMax)
3227     {
3228         eta_a=eta+v*((fEMCalPhiMax-phi)/(phi3-phi));
3229     }
3230     if (phi4<=fEMCalPhiMin)
3231     {
3232         eta_b=eta+v*((phi-fEMCalPhiMin)/(phi-phi4));
3233     }
3234
3235     if ((phi3>=fEMCalPhiMax) && (phi4<=fEMCalPhiMin))
3236     {
3237         eta_up=TMath::Max(eta_a,eta_b);
3238         eta_down=TMath::Min(eta_a,eta_b);
3239
3240         area_right=(fEMCalEtaMax-eta_up)*fEMCalPhiTotal + 0.5*(fEMCalPhiTotal+2*(eta_down-eta)*TMath::Tan(psi0))*(eta_down-eta_up) + (eta_down-eta+r)*TMath::Tan(psi0)*(eta_up-eta-r);
3241     }
3242     else if (phi3>=fEMCalPhiMax)
3243     {
3244         area_right=0.5*(fEMCalPhiMax-phi4+2*(eta_a-eta)*TMath::Tan(psi0))*(fEMCalEtaMax-eta_a) + (eta_a-eta+r)*TMath::Tan(psi0)*(eta_a-eta-r);
3245     }
3246     else if (phi4<=fEMCalPhiMin)
3247     {
3248         area_right=0.5*(phi3-fEMCalPhiMin+2*(eta_b-eta)*TMath::Tan(psi0))*(fEMCalEtaMax-eta_b) + (eta_b-eta+r)*TMath::Tan(psi0)*(eta_b-eta-r);
3249     }
3250     else
3251     {
3252         area_right=0.5*(phi3-phi4+2*r*TMath::Tan(psi0))*(v-r);
3253     }
3254     return area_left+area_right;
3255 }