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 #include <TClonesArray.h>
  20 #include <TObjArray.h>
  21
  22 #include "AliAnalysisTaskSE.h"
  23 #include "AliAnalysisManager.h"
  24 #include "AliStack.h"
  25 #include "AliESDtrackCuts.h"
  26 #include "AliESDEvent.h"
  27 #include "AliESDInputHandler.h"
  28 #include "AliAODEvent.h"
  29 #include "AliMCEvent.h"
  30 #include "AliVTrack.h"
  31 #include "AliVCluster.h"
  32 #include "AliEmcalJet.h"
  33 #include "AliEMCALGeometry.h"
  34 #include "AliPicoTrack.h"
  35 #include "Rtypes.h"
  36
  37 ClassImp(AliAnalysisTaskFullpAJets)
  38
  39 //________________________________________________________________________
  40 AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
  41     AliAnalysisTaskSE(),
  42
  43     fOutput(0),
  44     fhTrackPt(0),
  45     fhTrackEta(0),
  46     fhTrackPhi(0),
  47     fhGlobalTrackPt(0),
  48     fhGlobalTrackEta(0),
  49     fhGlobalTrackPhi(0),
  50     fhComplementaryTrackPt(0),
  51     fhComplementaryTrackEta(0),
  52     fhComplementaryTrackPhi(0),
  53     fhClusterPt(0),
  54     fhClusterEta(0),
  55     fhClusterPhi(0),
  56     fhCentrality(0),
  57     fhEMCalCellCounts(0),
  58     fhDeltaRhoN(0),
  59     fhDeltaRhoCMS(0),
  60
  61     fhTrackEtaPhi(0),
  62     fhTrackPhiPt(0),
  63     fhTrackEtaPt(0),
  64     fhGlobalTrackEtaPhi(0),
  65     fhGlobalTrackPhiPt(0),
  66     fhGlobalTrackEtaPt(0),
  67     fhComplementaryTrackEtaPhi(0),
  68     fhComplementaryTrackPhiPt(0),
  69     fhComplementaryTrackEtaPt(0),
  70     fhClusterEtaPhi(0),
  71     fhClusterPhiPt(0),
  72     fhClusterEtaPt(0),
  73     fhJetPtArea(0),
  74     fhJetConstituentPt(0),
  75     fhRhoScale(0),
  76
  77     fhTrackEtaPhiPt(0),
  78     fhGlobalTrackEtaPhiPt(0),
  79     fhComplementaryTrackEtaPhiPt(0),
  80     fhClusterEtaPhiPt(0),
  81
  82     fpEMCalEventMult(0),
  83     fpTPCEventMult(0),
  84     fpRhoScale(0),
  85
  86     fpJetEtaProfile(0),
  87     fpJetAbsEtaProfile(0),
  88     fpChargedJetRProfile(0),
  89     fpJetRProfile(0),
  90
  91     fpTrackPtProfile(0),
  92     fpClusterPtProfile(0),
  93
  94     fpChargedJetEDProfile(0),
  95     fpJetEDProfile(0),
  96
  97     fTPCRawJets(0),
  98     fEMCalRawJets(0),
  99     fRhoFull0(0),
 100     fRhoFull1(0),
 101     fRhoFull2(0),
 102     fRhoFullN(0),
 103     fRhoFullDijet(0),
 104     fRhoFullkT(0),
 105     fRhoFullCMS(0),
 106     fRhoCharged0(0),
 107     fRhoCharged1(0),
 108     fRhoCharged2(0),
 109     fRhoChargedN(0),
 110     fRhoChargedScale(0),
 111     fRhoChargedkT(0),
 112     fRhoChargedkTScale(0),
 113     fRhoChargedCMS(0),
 114     fRhoChargedCMSScale(0),
 115
 116     fTPCJet(0),
 117     fTPCFullJet(0),
 118     fTPCOnlyJet(0),
 119     fTPCkTFullJet(0),
 120     fEMCalJet(0),
 121     fEMCalFullJet(0),
 122     fEMCalPartJet(0),
 123     fEMCalkTFullJet(0),
 124
 125     fIsInitialized(0),
 126     fRJET(4),
 127     fnEvents(0),
 128     fnEventsCharged(0),
 129     fnDiJetEvents(0),
 130     fEMCalPhiMin(1.39626),
 131     fEMCalPhiMax(3.26377),
 132     fEMCalPhiTotal(1.86750),
 133     fEMCalEtaMin(-0.7),
 134     fEMCalEtaMax(0.7),
 135     fEMCalEtaTotal(1.4),
 136     fEMCalArea(2.61450),
 137     fTPCPhiMin(0),
 138     fTPCPhiMax(6.28319),
 139     fTPCPhiTotal(6.28319),
 140     fTPCEtaMin(-0.9),
 141     fTPCEtaMax(0.9),
 142     fTPCEtaTotal(1.8),
 143     fTPCArea(11.30973),
 144     fParticlePtLow(0.0),
 145     fParticlePtUp(200.0),
 146     fParticlePtBins(2000),
 147     fJetR(0.4),
 148     fJetRForRho(0.5),
 149     fJetAreaCutFrac(0.6),
 150     fJetAreaThreshold(0.30159),
 151     fnEMCalCells(12288),
 152     fScaleFactor(1.50),
 153     fNColl(7),
 154     fTrackMinPt(0.15),
 155     fClusterMinPt(0.3),
 156     fCentralityTag("V0A"),
 157     fCentralityBins(10),
 158     fCentralityLow(0),
 159     fCentralityUp(100),
 160     fEventCentrality(0),
 161     fRhoFull(0),
 162     fRhoCharged(0),
 163     fEtaProfileBins(14),
 164     fEtaProfileLow(-0.7),
 165     fEtaProfileUp(0.7),
 166     fEDProfileRBins(50),
 167     fEDProfileRLow(0),
 168     fEDProfileRUp(0.5),
 169     fEDProfilePtBins(100),
 170     fEDProfilePtLow(0),
 171     fEDProfilePtUp(100),
 172     fEDProfileEtaBins(4),
 173     fEDProfileEtaLow(-0.2),
 174     fEDProfileEtaUp(0.2),
 175     fnTracks(0),
 176     fnClusters(0),
 177     fnAKTFullJets(0),
 178     fnAKTChargedJets(0),
 179     fnKTFullJets(0),
 180     fnKTChargedJets(0),
 181     fnBckgClusters(0),
 182     fTPCJetThreshold(5),
 183     fEMCalJetThreshold(5),
 184     fVertexWindow(10),
 185     fVertexMaxR(1),
 186     fTrackName(0),
 187     fClusName(0),
 188     fkTChargedName(0),
 189     fAkTChargedName(0),
 190     fkTFullName(0),
 191     fAkTFullName(0),
 192     fOrgTracks(0),
 193     fOrgClusters(0),
 194     fmyAKTFullJets(0),
 195     fmyAKTChargedJets(0),
 196     fmyKTFullJets(0),
 197     fmyKTChargedJets(0),
 198     fmyTracks(0),
 199     fmyClusters(0),
 200     fEMCalRCBckgFluc(0),
 201     fTPCRCBckgFluc(0),
 202     fEMCalRCBckgFlucSignal(0),
 203     fTPCRCBckgFlucSignal(0),
 204     fEMCalRCBckgFlucNColl(0),
 205     fTPCRCBckgFlucNColl(0)
 206 {
 207     // Dummy constructor ALWAYS needed for I/O.
 208     fpJetEtaProfile = new TProfile *[14];
 209     fpJetAbsEtaProfile = new TProfile *[14];
 210     fpChargedJetRProfile = new TProfile *[8];
 211     fpJetRProfile= new TProfile *[4];
 212     fpChargedJetEDProfile= new TProfile3D *[10];
 213     fpJetEDProfile= new TProfile3D *[10];
 214     fvertex[0]=0.0,fvertex[1]=0.0,fvertex[2]=0.0;
 215 }
 216
 217 //________________________________________________________________________
 218 AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
 219     AliAnalysisTaskSE(name),
 220
 221     fOutput(0),
 222     fhTrackPt(0),
 223     fhTrackEta(0),
 224     fhTrackPhi(0),
 225     fhGlobalTrackPt(0),
 226     fhGlobalTrackEta(0),
 227     fhGlobalTrackPhi(0),
 228     fhComplementaryTrackPt(0),
 229     fhComplementaryTrackEta(0),
 230     fhComplementaryTrackPhi(0),
 231     fhClusterPt(0),
 232     fhClusterEta(0),
 233     fhClusterPhi(0),
 234     fhCentrality(0),
 235     fhEMCalCellCounts(0),
 236     fhDeltaRhoN(0),
 237     fhDeltaRhoCMS(0),
 238
 239     fhTrackEtaPhi(0),
 240     fhTrackPhiPt(0),
 241     fhTrackEtaPt(0),
 242     fhGlobalTrackEtaPhi(0),
 243     fhGlobalTrackPhiPt(0),
 244     fhGlobalTrackEtaPt(0),
 245     fhComplementaryTrackEtaPhi(0),
 246     fhComplementaryTrackPhiPt(0),
 247     fhComplementaryTrackEtaPt(0),
 248     fhClusterEtaPhi(0),
 249     fhClusterPhiPt(0),
 250     fhClusterEtaPt(0),
 251     fhJetPtArea(0),
 252     fhJetConstituentPt(0),
 253     fhRhoScale(0),
 254
 255     fhTrackEtaPhiPt(0),
 256     fhGlobalTrackEtaPhiPt(0),
 257     fhComplementaryTrackEtaPhiPt(0),
 258     fhClusterEtaPhiPt(0),
 259
 260     fpEMCalEventMult(0),
 261     fpTPCEventMult(0),
 262     fpRhoScale(0),
 263
 264     fpJetEtaProfile(0),
 265     fpJetAbsEtaProfile(0),
 266     fpChargedJetRProfile(0),
 267     fpJetRProfile(0),
 268
 269     fpTrackPtProfile(0),
 270     fpClusterPtProfile(0),
 271
 272     fpChargedJetEDProfile(0),
 273     fpJetEDProfile(0),
 274
 275     fTPCRawJets(0),
 276     fEMCalRawJets(0),
 277     fRhoFull0(0),
 278     fRhoFull1(0),
 279     fRhoFull2(0),
 280     fRhoFullN(0),
 281     fRhoFullDijet(0),
 282     fRhoFullkT(0),
 283     fRhoFullCMS(0),
 284     fRhoCharged0(0),
 285     fRhoCharged1(0),
 286     fRhoCharged2(0),
 287     fRhoChargedN(0),
 288     fRhoChargedScale(0),
 289     fRhoChargedkT(0),
 290     fRhoChargedkTScale(0),
 291     fRhoChargedCMS(0),
 292     fRhoChargedCMSScale(0),
 293
 294     fTPCJet(0),
 295     fTPCFullJet(0),
 296     fTPCOnlyJet(0),
 297     fTPCkTFullJet(0),
 298     fEMCalJet(0),
 299     fEMCalFullJet(0),
 300     fEMCalPartJet(0),
 301     fEMCalkTFullJet(0),
 302
 303     fIsInitialized(0),
 304     fRJET(4),
 305     fnEvents(0),
 306     fnEventsCharged(0),
 307     fnDiJetEvents(0),
 308     fEMCalPhiMin(1.39626),
 309     fEMCalPhiMax(3.26377),
 310     fEMCalPhiTotal(1.86750),
 311     fEMCalEtaMin(-0.7),
 312     fEMCalEtaMax(0.7),
 313     fEMCalEtaTotal(1.4),
 314     fEMCalArea(2.61450),
 315     fTPCPhiMin(0),
 316     fTPCPhiMax(6.28319),
 317     fTPCPhiTotal(6.28319),
 318     fTPCEtaMin(-0.9),
 319     fTPCEtaMax(0.9),
 320     fTPCEtaTotal(1.8),
 321     fTPCArea(11.30973),
 322     fParticlePtLow(0.0),
 323     fParticlePtUp(200.0),
 324     fParticlePtBins(2000),
 325     fJetR(0.4),
 326     fJetRForRho(0.5),
 327     fJetAreaCutFrac(0.6),
 328     fJetAreaThreshold(0.30159),
 329     fnEMCalCells(12288),
 330     fScaleFactor(1.50),
 331     fNColl(7),
 332     fTrackMinPt(0.15),
 333     fClusterMinPt(0.3),
 334     fCentralityTag("V0A"),
 335     fCentralityBins(10),
 336     fCentralityLow(0),
 337     fCentralityUp(100),
 338     fEventCentrality(0),
 339     fRhoFull(0),
 340     fRhoCharged(0),
 341     fEtaProfileBins(14),
 342     fEtaProfileLow(-0.7),
 343     fEtaProfileUp(0.7),
 344     fEDProfileRBins(50),
 345     fEDProfileRLow(0),
 346     fEDProfileRUp(0.5),
 347     fEDProfilePtBins(100),
 348     fEDProfilePtLow(0),
 349     fEDProfilePtUp(100),
 350     fEDProfileEtaBins(4),
 351     fEDProfileEtaLow(-0.2),
 352     fEDProfileEtaUp(0.2),
 353     fnTracks(0),
 354     fnClusters(0),
 355     fnAKTFullJets(0),
 356     fnAKTChargedJets(0),
 357     fnKTFullJets(0),
 358     fnKTChargedJets(0),
 359     fnBckgClusters(0),
 360     fTPCJetThreshold(5),
 361     fEMCalJetThreshold(5),
 362     fVertexWindow(10),
 363     fVertexMaxR(1),
 364     fTrackName(0),
 365     fClusName(0),
 366     fkTChargedName(0),
 367     fAkTChargedName(0),
 368     fkTFullName(0),
 369     fAkTFullName(0),
 370     fOrgTracks(0),
 371     fOrgClusters(0),
 372     fmyAKTFullJets(0),
 373     fmyAKTChargedJets(0),
 374     fmyKTFullJets(0),
 375     fmyKTChargedJets(0),
 376     fmyTracks(0),
 377     fmyClusters(0),
 378     fEMCalRCBckgFluc(0),
 379     fTPCRCBckgFluc(0),
 380     fEMCalRCBckgFlucSignal(0),
 381     fTPCRCBckgFlucSignal(0),
 382     fEMCalRCBckgFlucNColl(0),
 383     fTPCRCBckgFlucNColl(0)
 384 {
 385     // Constructor
 386     // Define input and output slots here (never in the dummy constructor)
 387     // Input slot #0 works with a TChain - it is connected to the default input container
 388     // Output slot #1 writes into a TH1 container
 389     fpJetEtaProfile = new TProfile *[14];
 390     fpJetAbsEtaProfile = new TProfile *[14];
 391     fpChargedJetRProfile = new TProfile *[8];
 392     fpJetRProfile = new TProfile *[4];
 393     fpChargedJetEDProfile= new TProfile3D *[10];
 394     fpJetEDProfile= new TProfile3D *[10];
 395     fvertex[0]=0.0,fvertex[1]=0.0,fvertex[2]=0.0;
 396
 397     DefineOutput(1,TList::Class());  // for output list
 398 }
 399
 400 //________________________________________________________________________
 401 AliAnalysisTaskFullpAJets::~AliAnalysisTaskFullpAJets()
 402 {
 403   // Destructor. Clean-up the output list, but not the histograms that are put inside
 404   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 405     if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())
 406     {
 407         delete fOutput;
 408     }
 409 }
 410
 411 //________________________________________________________________________
 412 void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
 413 {
 414     // Create histograms
 415     // Called once (on the worker node)
 416     fIsInitialized=kFALSE;
 417     fOutput = new TList();
 418     fOutput->SetOwner();  // IMPORTANT!
 419
 420     // Initialize Global Variables
 421     fnEvents=0;
 422     fnEventsCharged=0;
 423     fnDiJetEvents=0;
 424
 425     // fRJET=4 -> fJetR=0.4 && fRJET=25 -> fJetR=0.25, but for writing files, should be 4 and 25 respectively
 426     if (fRJET>10)
 427     {
 428         fJetR=(Double_t)fRJET/100.0;
 429     }
 430     else
 431     {
 432         fJetR=(Double_t)fRJET/10.0;
 433     }
 434     fJetRForRho=0.5;
 435
 436     fEMCalPhiMin=(80/(double)360)*2*TMath::Pi();
 437     fEMCalPhiMax=(187/(double)360)*2*TMath::Pi();
 438     fEMCalPhiTotal= fEMCalPhiMax-fEMCalPhiMin;
 439     fEMCalEtaMin=-0.7;
 440     fEMCalEtaMax=0.7;
 441     fEMCalEtaTotal=fEMCalEtaMax-fEMCalEtaMin;
 442     fEMCalArea=fEMCalPhiTotal*fEMCalEtaTotal;
 443
 444     fTPCPhiMin=(0/(double)360)*2*TMath::Pi();
 445     fTPCPhiMax=(360/(double)360)*2*TMath::Pi();
 446     fTPCPhiTotal= fTPCPhiMax-fTPCPhiMin;
 447     fTPCEtaMin=-0.9;
 448     fTPCEtaMax=0.9;
 449     fTPCEtaTotal=fTPCEtaMax-fTPCEtaMin;
 450     fTPCArea=fTPCPhiTotal*fTPCEtaTotal;
 451
 452     fParticlePtLow=0.0;
 453     fParticlePtUp=200.0;
 454     fParticlePtBins=Int_t(fParticlePtUp-fParticlePtLow);
 455
 456     fCentralityBins=10;
 457     fCentralityLow=0.0;
 458     fCentralityUp=100.0;
 459     Int_t CentralityBinMult=10;
 460
 461     fJetAreaCutFrac =0.6; // Fudge factor for selecting on jets with threshold Area or higher
 462     fJetAreaThreshold=fJetAreaCutFrac*TMath::Pi()*fJetR*fJetR;
 463     fTPCJetThreshold=5.0; // Threshold required for an Anti-kT Charged jet to be considered a "true" jet in TPC
 464     fEMCalJetThreshold=5.0; // Threshold required for an Anti-kT Charged+Neutral jet to be considered a "true" jet in EMCal
 465     fVertexWindow=10.0;
 466     fVertexMaxR=1.0;
 467
 468     fnBckgClusters=1;
 469     fEMCalRCBckgFluc = new Double_t[fnBckgClusters];
 470     fTPCRCBckgFluc = new Double_t[fnBckgClusters];
 471     fEMCalRCBckgFlucSignal = new Double_t[fnBckgClusters];
 472     fTPCRCBckgFlucSignal = new Double_t[fnBckgClusters];
 473     fEMCalRCBckgFlucNColl = new Double_t[fnBckgClusters];
 474     fTPCRCBckgFlucNColl = new Double_t[fnBckgClusters];
 475     for (Int_t i=0;i<fnBckgClusters;i++)
 476     {
 477         fEMCalRCBckgFluc[i]=0.0;
 478         fTPCRCBckgFluc[i]=0.0;
 479         fEMCalRCBckgFlucSignal[i]=0.0;
 480         fTPCRCBckgFlucSignal[i]=0.0;
 481         fEMCalRCBckgFlucNColl[i]=0.0;
 482         fTPCRCBckgFlucNColl[i]=0.0;
 483     }
 484
 485     fnEMCalCells=12288;  // sMods 1-10 have 24x48 cells, sMods 11&12 have 8x48 cells...
 486
 487     // Histograms
 488     Int_t JetPtBins = 200;
 489     Double_t JetPtLow = 0.0;
 490     Double_t JetPtUp = 200.0;
 491
 492     Int_t TCBins=100;
 493
 494     // QA Plots
 495     // Hybrid Tracks
 496     fhTrackPt = new TH1D("fhTrackPt","p_{T} distribution of tracks in event",fParticlePtBins,fParticlePtLow,fParticlePtUp);
 497     fhTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 498     fhTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 499     fhTrackPt->Sumw2();
 500
 501     fhTrackPhi = new TH1D("fhTrackPhi","#phi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
 502     fhTrackPhi->GetXaxis()->SetTitle("#phi");
 503     fhTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#phi");
 504     fhTrackPhi->Sumw2();
 505
 506     fhTrackEta = new TH1D("fhTrackEta","#eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
 507     fhTrackEta->GetXaxis()->SetTitle("#eta");
 508     fhTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
 509     fhTrackEta->Sumw2();
 510
 511     fhTrackEtaPhi = new TH2D("fhTrackEtaPhi","#eta-#phi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
 512     fhTrackEtaPhi->GetXaxis()->SetTitle("#eta");
 513     fhTrackEtaPhi->GetYaxis()->SetTitle("#phi");
 514     fhTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
 515     fhTrackEtaPhi->Sumw2();
 516
 517     fhTrackPhiPt = new TH2D("fhTrackPhiPt","#phi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 518     fhTrackPhiPt->GetXaxis()->SetTitle("#phi");
 519     fhTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 520     fhTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#phidp_{T}");
 521     fhTrackPhiPt->Sumw2();
 522
 523     fhTrackEtaPt = new TH2D("fhTrackEtaPt","#eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 524     fhTrackEtaPt->GetXaxis()->SetTitle("#phi");
 525     fhTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 526     fhTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
 527     fhTrackEtaPt->Sumw2();
 528
 529     fhTrackEtaPhiPt = new TH3D("fhTrackEtaPhiPt","#eta-#phi-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 530     fhTrackEtaPhiPt->GetXaxis()->SetTitle("#eta");
 531     fhTrackEtaPhiPt->GetYaxis()->SetTitle("#phi");
 532     fhTrackEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 533     fhTrackEtaPhiPt->Sumw2();
 534
 535     // Global Tracks
 536     fhGlobalTrackPt = new TH1D("fhGlobalTrackPt","Global p_{T} distribution of tracks in event",fParticlePtBins,fParticlePtLow,fParticlePtUp);
 537     fhGlobalTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 538     fhGlobalTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 539     fhGlobalTrackPt->Sumw2();
 540
 541     fhGlobalTrackPhi = new TH1D("fhGlobalTrackPhi","Global #phi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
 542     fhGlobalTrackPhi->GetXaxis()->SetTitle("#phi");
 543     fhGlobalTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#phi");
 544     fhGlobalTrackPhi->Sumw2();
 545
 546     fhGlobalTrackEta = new TH1D("fhGlobalTrackEta","Global #eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
 547     fhGlobalTrackEta->GetXaxis()->SetTitle("#eta");
 548     fhGlobalTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
 549     fhGlobalTrackEta->Sumw2();
 550
 551     fhGlobalTrackEtaPhi = new TH2D("fhGlobalTrackEtaPhi","Global #eta-#phi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
 552     fhGlobalTrackEtaPhi->GetXaxis()->SetTitle("#eta");
 553     fhGlobalTrackEtaPhi->GetYaxis()->SetTitle("#phi");
 554     fhGlobalTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
 555     fhGlobalTrackEtaPhi->Sumw2();
 556
 557     fhGlobalTrackPhiPt = new TH2D("fhGlobalTrackPhiPt","Global #phi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 558     fhGlobalTrackPhiPt->GetXaxis()->SetTitle("#phi");
 559     fhGlobalTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 560     fhGlobalTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#phidp_{T}");
 561     fhGlobalTrackPhiPt->Sumw2();
 562
 563     fhGlobalTrackEtaPt = new TH2D("fhGlobalTrackEtaPt","Global #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 564     fhGlobalTrackEtaPt->GetXaxis()->SetTitle("#phi");
 565     fhGlobalTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 566     fhGlobalTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
 567     fhGlobalTrackEtaPt->Sumw2();
 568
 569     fhGlobalTrackEtaPhiPt = new TH3D("fhGlobalTrackEtaPhiPt","Global #eta-#phi-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 570     fhGlobalTrackEtaPhiPt->GetXaxis()->SetTitle("#eta");
 571     fhGlobalTrackEtaPhiPt->GetYaxis()->SetTitle("#phi");
 572     fhGlobalTrackEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 573     fhGlobalTrackEtaPhiPt->Sumw2();
 574
 575     // Complementary Tracks
 576     fhComplementaryTrackPt = new TH1D("fhComplementaryTrackPt","Complementary p_{T} distribution of tracks in event",fParticlePtBins,fParticlePtLow,fParticlePtUp);
 577     fhComplementaryTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 578     fhComplementaryTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 579     fhComplementaryTrackPt->Sumw2();
 580
 581     fhComplementaryTrackPhi = new TH1D("fhComplementaryTrackPhi","Complementary #phi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
 582     fhComplementaryTrackPhi->GetXaxis()->SetTitle("#phi");
 583     fhComplementaryTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#phi");
 584     fhComplementaryTrackPhi->Sumw2();
 585
 586     fhComplementaryTrackEta = new TH1D("fhComplementaryTrackEta","Complementary #eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
 587     fhComplementaryTrackEta->GetXaxis()->SetTitle("#eta");
 588     fhComplementaryTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
 589     fhComplementaryTrackEta->Sumw2();
 590
 591     fhComplementaryTrackEtaPhi = new TH2D("fhComplementaryTrackEtaPhi","Complementary #eta-#phi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
 592     fhComplementaryTrackEtaPhi->GetXaxis()->SetTitle("#eta");
 593     fhComplementaryTrackEtaPhi->GetYaxis()->SetTitle("#phi");
 594     fhComplementaryTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
 595     fhComplementaryTrackEtaPhi->Sumw2();
 596
 597     fhComplementaryTrackPhiPt = new TH2D("fhComplementaryTrackPhiPt","Complementary #phi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 598     fhComplementaryTrackPhiPt->GetXaxis()->SetTitle("#phi");
 599     fhComplementaryTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 600     fhComplementaryTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#phidp_{T}");
 601     fhComplementaryTrackPhiPt->Sumw2();
 602
 603     fhComplementaryTrackEtaPt = new TH2D("fhComplementaryTrackEtaPt","Complementary #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 604     fhComplementaryTrackEtaPt->GetXaxis()->SetTitle("#phi");
 605     fhComplementaryTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 606     fhComplementaryTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
 607     fhComplementaryTrackEtaPt->Sumw2();
 608
 609     fhComplementaryTrackEtaPhiPt = new TH3D("fhComplementaryTrackEtaPhiPt","Complementary #eta-#phi-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 610     fhComplementaryTrackEtaPhiPt->GetXaxis()->SetTitle("#eta");
 611     fhComplementaryTrackEtaPhiPt->GetYaxis()->SetTitle("#phi");
 612     fhComplementaryTrackEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 613     fhComplementaryTrackEtaPhiPt->Sumw2();
 614
 615     // Corrected Calo Clusters
 616     fhClusterPt = new TH1D("fhClusterPt","p_{T} distribution of clusters in event",10*JetPtBins,JetPtLow,JetPtUp);
 617     fhClusterPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 618     fhClusterPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
 619     fhClusterPt->Sumw2();
 620
 621     fhClusterPhi = new TH1D("fhClusterPhi","#phi distribution of clusters in event",TCBins,fTPCPhiMin,fTPCPhiMax);
 622     fhClusterPhi->GetXaxis()->SetTitle("#phi");
 623     fhClusterPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#phi");
 624     fhClusterPhi->Sumw2();
 625
 626     fhClusterEta = new TH1D("fhClusterEta","#eta distribution of clusters in event",TCBins,fTPCEtaMin,fTPCEtaMax);
 627     fhClusterEta->GetXaxis()->SetTitle("#eta");
 628     fhClusterEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
 629     fhClusterEta->Sumw2();
 630
 631     fhClusterEtaPhi = new TH2D("fhClusterEtaPhi","#eta-#phi distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
 632     fhClusterEtaPhi->GetXaxis()->SetTitle("#eta");
 633     fhClusterEtaPhi->GetYaxis()->SetTitle("#phi");
 634     fhClusterEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
 635     fhClusterEtaPhi->Sumw2();
 636
 637     fhClusterPhiPt = new TH2D("fhClusterPhiPt","#phi-p_{T} distribution of clusters in event",TCBins,fEMCalPhiMin,fEMCalPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 638     fhClusterPhiPt->GetXaxis()->SetTitle("#phi");
 639     fhClusterPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 640     fhClusterPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#phidp_{T}");
 641     fhClusterPhiPt->Sumw2();
 642
 643     fhClusterEtaPt = new TH2D("fhClusterEtaPt","#eta-p_{T} distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 644     fhClusterEtaPt->GetXaxis()->SetTitle("#phi");
 645     fhClusterEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 646     fhClusterEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
 647     fhClusterEtaPt->Sumw2();
 648
 649     fhClusterEtaPhiPt = new TH3D("fhClusterEtaPhiPt","#eta-#phi-p_{T} distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
 650     fhClusterEtaPhiPt->GetXaxis()->SetTitle("#eta");
 651     fhClusterEtaPhiPt->GetYaxis()->SetTitle("#phi");
 652     fhClusterEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 653     fhClusterEtaPhiPt->Sumw2();
 654
 655     fhCentrality = new TH1D("fhCentrality","Event Centrality Distribution",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
 656     fhCentrality->GetXaxis()->SetTitle(fCentralityTag);
 657     fhCentrality->GetYaxis()->SetTitle("1/N_{Events}");
 658     fhCentrality->Sumw2();
 659
 660     fhJetConstituentPt= new TH2D("fhJetConstituentPt","Jet constituents p_{T} distribution",JetPtBins, JetPtLow, JetPtUp,10*JetPtBins, JetPtLow, JetPtUp);
 661     fhJetConstituentPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
 662     fhJetConstituentPt->GetYaxis()->SetTitle("Constituent p_{T} (GeV/c)");
 663     fhJetConstituentPt->Sumw2();
 664
 665     fhEMCalCellCounts = new TH1D("fhEMCalCellCounts","Distribtuion of cluster counts across the EMCal",fnEMCalCells,1,fnEMCalCells);
 666     fhEMCalCellCounts->GetXaxis()->SetTitle("Absoulute Cell Id");
 667     fhEMCalCellCounts->GetYaxis()->SetTitle("Counts per Event");
 668     fhEMCalCellCounts->Sumw2();
 669
 670     // Rho QA Plots
 671     Int_t RhoBins = 1000;
 672     Double_t RhoPtMin = -50.0;
 673     Double_t RhoPtMax = 50.0;
 674
 675     fhDeltaRhoN = new TH1D("fhDeltaRhoN","0-100% #delta#rho_{N} = #rho_{N}^{TPC+EMCal} - #rho_{N}^{TPC+Scale}",RhoBins,RhoPtMin,RhoPtMax);
 676     fhDeltaRhoN->GetXaxis()->SetTitle("#delta#rho (GeV)");
 677     fhDeltaRhoN->GetYaxis()->SetTitle("Counts");
 678     fhDeltaRhoN->Sumw2();
 679
 680     fhDeltaRhoCMS = new TH1D("fhDeltaRhoCMS","0-100% #delta#rho_{CMS} = #rho_{CMS}^{TPC+EMCal} - #rho_{CMS}^{TPC+Scale}",RhoBins,RhoPtMin,RhoPtMax);
 681     fhDeltaRhoCMS->GetXaxis()->SetTitle("#delta#rho (GeV)");
 682     fhDeltaRhoCMS->GetYaxis()->SetTitle("Counts");
 683     fhDeltaRhoCMS->Sumw2();
 684
 685     // Jet Area vs pT Distribution
 686     Int_t JetPtAreaBins=200;
 687     Double_t JetPtAreaLow=0.0;
 688     Double_t JetPtAreaUp=2.0;
 689
 690     Int_t SFBins =100;
 691     Double_t SFLow=0.0;
 692     Double_t SFUp=10.0;
 693
 694     fhJetPtArea = new TH2D("fhJetPtArea","Jet Area Distribution",JetPtBins, JetPtLow,JetPtUp,JetPtAreaBins,JetPtAreaLow,JetPtAreaUp);
 695     fhJetPtArea->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 696     fhJetPtArea->GetYaxis()->SetTitle("A_{jet}");
 697     fhJetPtArea->GetZaxis()->SetTitle("1/N_{Events} dN/dA_{jet}dp_{T}");
 698     fhJetPtArea->Sumw2();
 699
 700     fhRhoScale = new TH2D("fhRhoScale","Scaling Factor",SFBins,SFLow,SFUp,CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
 701     fhRhoScale->GetXaxis()->SetTitle("Scale Factor");
 702     fhRhoScale->GetYaxis()->SetTitle("Centrality");
 703     fhRhoScale->GetZaxis()->SetTitle("Counts");
 704     fhRhoScale->Sumw2();
 705
 706     // Profiles
 707     fpEMCalEventMult = new TProfile("fpEMCalEventMult","EMCal Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
 708     fpEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag);
 709     fpEMCalEventMult->GetYaxis()->SetTitle("Multiplicity");
 710
 711     fpTPCEventMult = new TProfile("fpTPCEventMult","TPC Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
 712     fpTPCEventMult->GetXaxis()->SetTitle(fCentralityTag);
 713     fpTPCEventMult->GetYaxis()->SetTitle("Multiplicity");
 714
 715     fpRhoScale = new TProfile("fpRhoScale","Scaling Factor Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
 716     fpRhoScale->GetXaxis()->SetTitle(fCentralityTag);
 717     fpRhoScale->GetYaxis()->SetTitle("Scale Factor");
 718
 719     // QA::2D Energy Density Profiles for Tracks and Clusters
 720     fpTrackPtProfile = new TProfile2D("fpTrackPtProfile","2D Profile of track pT density throughout the TPC",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
 721     fpTrackPtProfile->GetXaxis()->SetTitle("#eta");
 722     fpTrackPtProfile->GetYaxis()->SetTitle("#phi");
 723     fpTrackPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
 724
 725     fpClusterPtProfile = new TProfile2D("fpClusterPtProfile","2D Profile of cluster pT density throughout the EMCal",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
 726     fpClusterPtProfile->GetXaxis()->SetTitle("#eta");
 727     fpClusterPtProfile->GetYaxis()->SetTitle("#phi");
 728     fpClusterPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
 729
 730     TString temp_name="";
 731     TString title_name="";
 732
 733     fEtaProfileBins=14;
 734     fEtaProfileLow=-0.7;
 735     fEtaProfileUp=0.7;
 736
 737     for (Int_t i=0;i<14;i++)
 738     {
 739         temp_name=Form("fpJetEtaProfile%d",i);
 740         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.);
 741
 742         fpJetEtaProfile[i]= new TProfile(temp_name,title_name,fEtaProfileBins,fEtaProfileLow,fEtaProfileUp);
 743         fpJetEtaProfile[i]->GetXaxis()->SetTitle("#eta");
 744         fpJetEtaProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
 745         fOutput->Add(fpJetEtaProfile[i]);
 746         temp_name="";
 747         title_name="";
 748
 749         temp_name=Form("fpJetAbsEtaProfile%d",i);
 750         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.);
 751
 752         fpJetAbsEtaProfile[i]= new TProfile(temp_name,title_name,fEtaProfileBins,0,2*fEtaProfileUp);
 753         fpJetAbsEtaProfile[i]->GetXaxis()->SetTitle("#Delta#eta");
 754         fpJetAbsEtaProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
 755         fOutput->Add(fpJetAbsEtaProfile[i]);
 756         temp_name="";
 757         title_name="";
 758 }
 759
 760     fEDProfileRBins=50;
 761     fEDProfileRLow=0.0;
 762     fEDProfileRUp=0.5;
 763     fEDProfilePtBins=100;
 764     fEDProfilePtLow=0.0;
 765     fEDProfilePtUp=100.0;
 766     fEDProfileEtaBins=4;
 767     fEDProfileEtaLow=-0.2;
 768     fEDProfileEtaUp=0.2;
 769
 770     for (Int_t i=0;i<8;i++)
 771     {
 772         temp_name=Form("fpChargedJetRProfile%d",i);
 773         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.);
 774
 775         fpChargedJetRProfile[i]= new TProfile(temp_name,title_name,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
 776         fpChargedJetRProfile[i]->GetXaxis()->SetTitle("Radius");
 777         fpChargedJetRProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
 778         fOutput->Add(fpChargedJetRProfile[i]);
 779         temp_name="";
 780         title_name="";
 781     }
 782
 783     for (Int_t i=0;i<4;i++)
 784     {
 785         temp_name=Form("fpJetRProfile%d",i);
 786         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.);
 787
 788         fpJetRProfile[i]= new TProfile(temp_name,title_name,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
 789         fpJetRProfile[i]->GetXaxis()->SetTitle("Radius");
 790         fpJetRProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
 791         fOutput->Add(fpJetRProfile[i]);
 792         temp_name="";
 793         title_name="";
 794     }
 795
 796     for (Int_t i=0;i<10;i++)
 797     {
 798         temp_name=Form("fpChargedJetEDProfile%d0",i);
 799         title_name=Form("Charged Jet Energy Density Profile for %d0-%d0 Centrality",i,i+1);
 800
 801         fpChargedJetEDProfile[i]= new TProfile3D(temp_name,title_name,fEDProfilePtBins,fEDProfilePtLow,fEDProfilePtUp,fEDProfileEtaBins+4,fEDProfileEtaLow-0.2,fEDProfileEtaUp+0.2,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
 802         fpChargedJetEDProfile[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 803         fpChargedJetEDProfile[i]->GetYaxis()->SetTitle("Pseudorapidity");
 804         fpChargedJetEDProfile[i]->GetZaxis()->SetTitle("Radius");
 805         fOutput->Add(fpChargedJetEDProfile[i]);
 806         temp_name="";
 807         title_name="";
 808
 809         temp_name=Form("fpJetEDProfile%d0",i);
 810         title_name=Form("Jet Energy Density Profile for %d0-%d0 Centrality",i,i+1);
 811
 812         fpJetEDProfile[i]= new TProfile3D(temp_name,title_name,fEDProfilePtBins,fEDProfilePtLow,fEDProfilePtUp,fEDProfileEtaBins,fEDProfileEtaLow,fEDProfileEtaUp,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
 813         fpJetEDProfile[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 814         fpJetEDProfile[i]->GetYaxis()->SetTitle("Pseudorapidity");
 815         fpJetEDProfile[i]->GetZaxis()->SetTitle("Radius");
 816         fOutput->Add(fpJetEDProfile[i]);
 817         temp_name="";
 818         title_name="";
 819     }
 820
 821     fTPCRawJets = new AlipAJetHistos("TPCRawJets",fCentralityTag);
 822     fEMCalRawJets = new AlipAJetHistos("EMCalRawJets",fCentralityTag);
 823
 824     fRhoFull0 = new AlipAJetHistos("RhoFull0",fCentralityTag);
 825     fRhoFull1 = new AlipAJetHistos("RhoFull1",fCentralityTag);
 826     fRhoFull2 = new AlipAJetHistos("RhoFull2",fCentralityTag);
 827     fRhoFullN = new AlipAJetHistos("RhoFullN",fCentralityTag);
 828     fRhoFullDijet = new AlipAJetHistos("RhoFullDijet",fCentralityTag);
 829     fRhoFullkT = new AlipAJetHistos("RhoFullkT",fCentralityTag);
 830     fRhoFullCMS = new AlipAJetHistos("RhoFullCMS",fCentralityTag);
 831
 832     fRhoCharged0 = new AlipAJetHistos("RhoCharged0",fCentralityTag);
 833     fRhoCharged1 = new AlipAJetHistos("RhoCharged1",fCentralityTag);
 834     fRhoCharged2 = new AlipAJetHistos("RhoCharged2",fCentralityTag);
 835     fRhoChargedN = new AlipAJetHistos("RhoChargedN",fCentralityTag);
 836     fRhoChargedkT = new AlipAJetHistos("RhoChargedkT",fCentralityTag);
 837     fRhoChargedScale = new AlipAJetHistos("RhoChargedScale",fCentralityTag);
 838     fRhoChargedkTScale = new AlipAJetHistos("RhoChargedkTScale",fCentralityTag);
 839     fRhoChargedCMS = new AlipAJetHistos("RhoChargedCMS",fCentralityTag);
 840     fRhoChargedCMSScale = new AlipAJetHistos("RhoChargedCMSScale",fCentralityTag);
 841
 842     fOutput->Add(fhTrackPt);
 843     fOutput->Add(fhTrackEta);
 844     fOutput->Add(fhTrackPhi);
 845     fOutput->Add(fhTrackEtaPhi);
 846     fOutput->Add(fhTrackPhiPt);
 847     fOutput->Add(fhTrackEtaPt);
 848     fOutput->Add(fhTrackEtaPhiPt);
 849     fOutput->Add(fhGlobalTrackPt);
 850     fOutput->Add(fhGlobalTrackEta);
 851     fOutput->Add(fhGlobalTrackPhi);
 852     fOutput->Add(fhGlobalTrackEtaPhi);
 853     fOutput->Add(fhGlobalTrackPhiPt);
 854     fOutput->Add(fhGlobalTrackEtaPt);
 855     fOutput->Add(fhGlobalTrackEtaPhiPt);
 856     fOutput->Add(fhComplementaryTrackPt);
 857     fOutput->Add(fhComplementaryTrackEta);
 858     fOutput->Add(fhComplementaryTrackPhi);
 859     fOutput->Add(fhComplementaryTrackEtaPhi);
 860     fOutput->Add(fhComplementaryTrackPhiPt);
 861     fOutput->Add(fhComplementaryTrackEtaPt);
 862     fOutput->Add(fhComplementaryTrackEtaPhiPt);
 863     fOutput->Add(fhClusterPt);
 864     fOutput->Add(fhClusterEta);
 865     fOutput->Add(fhClusterPhi);
 866     fOutput->Add(fhClusterEtaPhi);
 867     fOutput->Add(fhClusterPhiPt);
 868     fOutput->Add(fhClusterEtaPt);
 869     fOutput->Add(fhClusterEtaPhiPt);
 870     fOutput->Add(fhCentrality);
 871     fOutput->Add(fhEMCalCellCounts);
 872     fOutput->Add(fhDeltaRhoN);
 873     fOutput->Add(fhDeltaRhoCMS);
 874     fOutput->Add(fhJetPtArea);
 875     fOutput->Add(fhJetConstituentPt);
 876     fOutput->Add(fhRhoScale);
 877
 878     fOutput->Add(fpTPCEventMult);
 879     fOutput->Add(fpEMCalEventMult);
 880     fOutput->Add(fpRhoScale);
 881
 882     fOutput->Add(fpTrackPtProfile);
 883     fOutput->Add(fpClusterPtProfile);
 884
 885     fOutput->Add(fTPCRawJets->GetOutputHistos());
 886     fOutput->Add(fEMCalRawJets->GetOutputHistos());
 887
 888     fOutput->Add(fRhoFull0->GetOutputHistos());
 889     fOutput->Add(fRhoFull1->GetOutputHistos());
 890     fOutput->Add(fRhoFull2->GetOutputHistos());
 891     fOutput->Add(fRhoFullN->GetOutputHistos());
 892     fOutput->Add(fRhoFullDijet->GetOutputHistos());
 893     fOutput->Add(fRhoFullkT->GetOutputHistos());
 894     fOutput->Add(fRhoFullCMS->GetOutputHistos());
 895     fOutput->Add(fRhoCharged0->GetOutputHistos());
 896     fOutput->Add(fRhoCharged1->GetOutputHistos());
 897     fOutput->Add(fRhoCharged2->GetOutputHistos());
 898     fOutput->Add(fRhoChargedN->GetOutputHistos());
 899     fOutput->Add(fRhoChargedkT->GetOutputHistos());
 900     fOutput->Add(fRhoChargedScale->GetOutputHistos());
 901     fOutput->Add(fRhoChargedkTScale->GetOutputHistos());
 902     fOutput->Add(fRhoChargedCMS->GetOutputHistos());
 903     fOutput->Add(fRhoChargedCMSScale->GetOutputHistos());
 904
 905     // Post data for ALL output slots >0 here,
 906     // To get at least an empty histogram
 907     // 1 is the outputnumber of a certain weg of task 1
 908     PostData(1, fOutput);
 909 }
 910
 911 void AliAnalysisTaskFullpAJets::UserExecOnce()
 912 {
 913     // Get the event tracks
 914     fOrgTracks = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fTrackName));
 915
 916     // Get the event caloclusters
 917     fOrgClusters = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fClusName));
 918
 919     // Get charged jets
 920     fmyKTChargedJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fkTChargedName));
 921     fmyAKTChargedJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fAkTChargedName));
 922
 923     // Get the full jets
 924     fmyKTFullJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fkTFullName));
 925     fmyAKTFullJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fAkTFullName));
 926
 927     fIsInitialized=kTRUE;
 928 }
 929 //________________________________________________________________________
 930 void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
 931 {
 932     if (fIsInitialized==kFALSE)
 933     {
 934         UserExecOnce();
 935     }
 936
 937     // Get pointer to reconstructed event
 938     AliVEvent *event = InputEvent();
 939     if (!event)
 940     {
 941         AliError("Pointer == 0, this can not happen!");
 942         return;
 943     }
 944
 945     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 946     AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
 947
 948     if (esd)
 949     {
 950         fEventCentrality=esd->GetCentrality()->GetCentralityPercentile(fCentralityTag);
 951
 952         if (esd->GetPrimaryVertex()->GetNContributors()<1 || (TMath::Abs(esd->GetPrimaryVertex()->GetZ())>fVertexWindow))
 953         {
 954             return;
 955         }
 956         if (TMath::Sqrt(TMath::Power(esd->GetPrimaryVertex()->GetX(),2)+TMath::Power(esd->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
 957         {
 958             return;
 959         }
 960
 961         esd->GetPrimaryVertex()->GetXYZ(fvertex);
 962     }
 963     else if (aod)
 964     {
 965         fEventCentrality=aod->GetCentrality()->GetCentralityPercentile(fCentralityTag);
 966
 967         if (aod->GetPrimaryVertex()->GetNContributors()<1 || (TMath::Abs(aod->GetPrimaryVertex()->GetZ())>fVertexWindow))
 968         {
 969             return;
 970         }
 971         if (TMath::Sqrt(TMath::Power(aod->GetPrimaryVertex()->GetX(),2)+TMath::Power(aod->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
 972         {
 973             return;
 974         }
 975
 976         aod->GetPrimaryVertex()->GetXYZ(fvertex);
 977     }
 978     else
 979     {
 980         AliError("Cannot get AOD/ESD event. Rejecting Event");
 981         return;
 982     }
 983
 984     // Make sure Centrality isn't exactly 100% (to avoid bin filling errors in profile plots. Make it 99.99
 985     if (fEventCentrality>99.99)
 986     {
 987         fEventCentrality=99.99;
 988     }
 989     fhCentrality->Fill(fEventCentrality);
 990
 991     TrackCuts();
 992     // Reject any event that doesn't have any tracks, i.e. TPC is off
 993     if (fnTracks<1)
 994     {
 995         AliWarning("No PicoTracks, Rejecting Event");
 996         return;
 997     }
 998
 999     ClusterCuts();
1000     if (fnClusters<1)
1001     {
1002         AliInfo("No Corrected CaloClusters, using only charged jets");
1003
1004         TrackHisto();
1005         InitChargedJets();
1006         GenerateTPCRandomConesPt();
1007
1008         // Rho's
1009         EstimateChargedRho0();
1010         EstimateChargedRho1();
1011         EstimateChargedRho2();
1012         EstimateChargedRhoN();
1013         EstimateChargedRhokT();
1014         EstimateChargedRhoCMS();
1015
1016         JetPtChargedProfile();
1017         DeleteJetData(kFALSE);
1018
1019         fnEventsCharged++;
1020
1021         PostData(1, fOutput);
1022         return;
1023     }
1024
1025     TrackHisto();
1026     ClusterHisto();
1027
1028     // Prep the jets
1029     InitChargedJets();
1030     InitFullJets();
1031     JetPtArea();
1032     GenerateTPCRandomConesPt();
1033     GenerateEMCalRandomConesPt();
1034
1035     // Rho's
1036     EstimateChargedRho0();
1037     EstimateChargedRho1();
1038     EstimateChargedRho2();
1039     EstimateChargedRhoN();
1040     EstimateChargedRhokT();
1041     EstimateChargedRhoCMS();
1042
1043     EstimateFullRho0();
1044     EstimateFullRho1();
1045     EstimateFullRho2();
1046     EstimateFullRhoN();
1047     EstimateFullRhokT();
1048     EstimateFullRhoCMS();
1049
1050     EstimateChargedRhoScale();
1051     EstimateChargedRhokTScale();
1052     EstimateChargedRhoCMSScale();
1053
1054     // Dijet
1055     if (IsDiJetEvent()==kTRUE)
1056     {
1057         EstimateFullRhoDijet();
1058     }
1059
1060     // Compute Jet Energy Density Profile
1061     JetPtChargedProfile();
1062     JetPtFullProfile();
1063     JetPtEtaProfile();
1064
1065     // Compute differences between TPC+EMCal Rho to TPC&Scaled Rho
1066     if (fRhoChargedScale->GetRho()>0 && fRhoFullN->GetRho()>0)
1067     {
1068         fhDeltaRhoN->Fill(fRhoFullN->GetRho()-fRhoChargedScale->GetRho());
1069     }
1070     if (fRhoChargedCMSScale->GetRho()>0 && fRhoFullCMS->GetRho()>0)
1071     {
1072         fhDeltaRhoCMS->Fill(fRhoFullCMS->GetRho()-fRhoChargedCMSScale->GetRho());
1073     }
1074
1075     // Delete Dynamic Arrays
1076     DeleteJetData(kTRUE);
1077     fnEvents++;
1078
1079     PostData(1, fOutput);
1080 }
1081
1082 //________________________________________________________________________
1083 void AliAnalysisTaskFullpAJets::Terminate(Option_t *) //specify what you want to have done
1084 {
1085     // Called once at the end of the query. Done nothing here.
1086 }
1087
1088 /////////////////////////////////////////////////////////////////////////////////////////
1089 /////////////////     User Defined Sub_Routines   ///////////////////////////////////////
1090 /////////////////////////////////////////////////////////////////////////////////////////
1091
1092 void AliAnalysisTaskFullpAJets::TrackCuts()
1093 {
1094     // Fill a TObjArray with the tracks from a TClonesArray which grabs the picotracks.
1095     Int_t i;
1096
1097     fmyTracks = new TObjArray();
1098     for (i=0;i<fOrgTracks->GetEntries();i++)
1099     {
1100         AliVTrack* vtrack = (AliVTrack*) fOrgTracks->At(i);
1101         if (vtrack->Pt()>=fTrackMinPt)
1102         {
1103             fmyTracks->Add(vtrack);
1104         }
1105     }
1106     fnTracks = fmyTracks->GetEntries();
1107 }
1108
1109 void AliAnalysisTaskFullpAJets::ClusterCuts()
1110 {
1111     // Fill a TObjArray with the clusters from a TClonesArray which grabs the caloclusterscorr.
1112     Int_t i;
1113
1114     fmyClusters = new TObjArray();
1115     if(fOrgClusters) {
1116     for (i=0;i<fOrgClusters->GetEntries();i++)
1117     {
1118         AliVCluster* vcluster = (AliVCluster*) fOrgClusters->At(i);
1119         TLorentzVector *cluster_vec = new TLorentzVector;
1120         vcluster->GetMomentum(*cluster_vec,fvertex);
1121
1122         if (cluster_vec->Pt()>=fClusterMinPt)
1123         {
1124             fmyClusters->Add(vcluster);
1125         }
1126         delete cluster_vec;
1127
1128     }
1129     }
1130     fnClusters = fmyClusters->GetEntries();
1131 }
1132
1133 void AliAnalysisTaskFullpAJets::TrackHisto()
1134 {
1135     // Fill track histograms: Phi,Eta,Pt
1136     Int_t i,j;
1137     Int_t TCBins=100;
1138     TH2D *hdummypT= new TH2D("hdummypT","",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);  //!
1139
1140     for (i=0;i<fnTracks;i++)
1141     {
1142         AliPicoTrack* vtrack =(AliPicoTrack*) fmyTracks->At(i);
1143         fhTrackPt->Fill(vtrack->Pt());
1144         fhTrackEta->Fill(vtrack->Eta());
1145         fhTrackPhi->Fill(vtrack->Phi());
1146         fhTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
1147         fhTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
1148         fhTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
1149         fhTrackEtaPhiPt->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
1150
1151         // Fill Associated Track Distributions for AOD QA Productions
1152         // Global Tracks
1153         if (vtrack->GetTrackType()==0)
1154         {
1155             fhGlobalTrackPt->Fill(vtrack->Pt());
1156             fhGlobalTrackEta->Fill(vtrack->Eta());
1157             fhGlobalTrackPhi->Fill(vtrack->Phi());
1158             fhGlobalTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
1159             fhGlobalTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
1160             fhGlobalTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
1161             fhGlobalTrackEtaPhiPt->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
1162         }
1163         // Complementary Tracks
1164         else if (vtrack->GetTrackType()==1)
1165         {
1166             fhComplementaryTrackPt->Fill(vtrack->Pt());
1167             fhComplementaryTrackEta->Fill(vtrack->Eta());
1168             fhComplementaryTrackPhi->Fill(vtrack->Phi());
1169             fhComplementaryTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
1170             fhComplementaryTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
1171             fhComplementaryTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
1172             fhComplementaryTrackEtaPhiPt->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
1173         }
1174
1175         hdummypT->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
1176     }
1177     for (i=1;i<=TCBins;i++)
1178     {
1179         for (j=1;j<=TCBins;j++)
1180         {
1181             fpTrackPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fTPCArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
1182         }
1183     }
1184     delete hdummypT;
1185 }
1186
1187 void AliAnalysisTaskFullpAJets::ClusterHisto()
1188 {
1189     // Fill cluster histograms: Phi,Eta,Pt
1190     Int_t i,j;
1191     Int_t TCBins=100;
1192     TH2D *hdummypT= new TH2D("hdummypT","",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);  //!
1193     AliEMCALGeometry *myAliEMCGeo = AliEMCALGeometry::GetInstance();
1194     Int_t myCellID=-2;
1195
1196     for (i=0;i<fnClusters;i++)
1197     {
1198         AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1199         TLorentzVector *cluster_vec = new TLorentzVector;
1200         vcluster->GetMomentum(*cluster_vec,fvertex);
1201
1202         fhClusterPt->Fill(cluster_vec->Pt());
1203         fhClusterEta->Fill(cluster_vec->Eta());
1204         fhClusterPhi->Fill(cluster_vec->Phi());
1205         fhClusterEtaPhi->Fill(cluster_vec->Eta(),cluster_vec->Phi());
1206         fhClusterPhiPt->Fill(cluster_vec->Phi(),cluster_vec->Pt());
1207         fhClusterEtaPt->Fill(cluster_vec->Eta(),cluster_vec->Pt());
1208         fhClusterEtaPhiPt->Fill(cluster_vec->Eta(),cluster_vec->Phi(),cluster_vec->Pt());
1209         hdummypT->Fill(cluster_vec->Eta(),cluster_vec->Phi(),cluster_vec->Pt());
1210         myAliEMCGeo->GetAbsCellIdFromEtaPhi(cluster_vec->Eta(),cluster_vec->Phi(),myCellID);
1211         fhEMCalCellCounts->Fill(myCellID);
1212         delete cluster_vec;
1213     }
1214     for (i=1;i<=TCBins;i++)
1215     {
1216         for (j=1;j<=TCBins;j++)
1217         {
1218             fpClusterPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fEMCalArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
1219         }
1220     }
1221     delete hdummypT;
1222 }
1223
1224 void AliAnalysisTaskFullpAJets::InitChargedJets()
1225 {
1226     // Preliminary Jet Placement and Selection Cuts
1227     Int_t i;
1228
1229     fnAKTChargedJets = fmyAKTChargedJets->GetEntries();
1230     fnKTChargedJets = fmyKTChargedJets->GetEntries();
1231
1232     fTPCJet = new AlipAJetData("fTPCJet",kFALSE,fnAKTChargedJets);
1233     fTPCFullJet = new AlipAJetData("fTPCFullJet",kFALSE,fnAKTChargedJets);
1234     fTPCOnlyJet = new AlipAJetData("fTPCOnlyJet",kFALSE,fnAKTChargedJets);
1235
1236     fTPCJet->SetSignalCut(fTPCJetThreshold);
1237     fTPCJet->SetAreaCutFraction(fJetAreaCutFrac);
1238     fTPCJet->SetJetR(fJetR);
1239     fTPCFullJet->SetSignalCut(fTPCJetThreshold);
1240     fTPCFullJet->SetAreaCutFraction(fJetAreaCutFrac);
1241     fTPCFullJet->SetJetR(fJetR);
1242     fTPCOnlyJet->SetSignalCut(fTPCJetThreshold);
1243     fTPCOnlyJet->SetAreaCutFraction(fJetAreaCutFrac);
1244     fTPCOnlyJet->SetJetR(fJetR);
1245
1246     // Initialize Jet Data
1247     for (i=0;i<fnAKTChargedJets;i++)
1248     {
1249         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(i);
1250
1251         fTPCJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kFALSE),i);
1252         fTPCFullJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE),i);
1253         fTPCOnlyJet->SetIsJetInArray(IsInTPCFull(fJetR,myJet->Phi(),myJet->Eta()),i);
1254     }
1255     fTPCJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
1256     fTPCFullJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
1257     fTPCOnlyJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
1258
1259     // kT Jets
1260     fTPCkTFullJet = new AlipAJetData("fTPCkTFullJet",kFALSE,fnKTChargedJets);
1261     fTPCkTFullJet->SetSignalCut(fTPCJetThreshold);
1262     fTPCkTFullJet->SetAreaCutFraction(0.25*fJetAreaCutFrac);
1263     fTPCkTFullJet->SetJetR(fJetR);
1264
1265     for (i=0;i<fnKTChargedJets;i++)
1266     {
1267         AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(i);
1268         fTPCkTFullJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE),i);
1269     }
1270     fTPCkTFullJet->InitializeJetData(fmyKTChargedJets,fnKTChargedJets);
1271
1272     // Raw Charged Jet Spectra
1273     fTPCRawJets->FillBSJS(fEventCentrality,0.0,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
1274 }
1275
1276 void AliAnalysisTaskFullpAJets::InitFullJets()
1277 {
1278     // Preliminary Jet Placement and Selection Cuts
1279     Int_t i;
1280
1281     fnAKTFullJets = fmyAKTFullJets->GetEntries();
1282     fnKTFullJets = fmyKTFullJets->GetEntries();
1283
1284     fEMCalJet = new AlipAJetData("fEMCalJet",kTRUE,fnAKTFullJets);
1285     fEMCalFullJet = new AlipAJetData("fEMCalFullJet",kTRUE,fnAKTFullJets);
1286     fEMCalPartJet = new AlipAJetData("fEMCalPartJet",kTRUE,fnAKTFullJets);
1287
1288     fEMCalJet->SetSignalCut(fEMCalJetThreshold);
1289     fEMCalJet->SetAreaCutFraction(fJetAreaCutFrac);
1290     fEMCalJet->SetJetR(fJetR);
1291     fEMCalFullJet->SetSignalCut(fEMCalJetThreshold);
1292     fEMCalFullJet->SetAreaCutFraction(fJetAreaCutFrac);
1293     fEMCalFullJet->SetJetR(fJetR);
1294     fEMCalPartJet->SetSignalCut(fEMCalJetThreshold);
1295     fEMCalPartJet->SetAreaCutFraction(fJetAreaCutFrac);
1296     fEMCalPartJet->SetJetR(fJetR);
1297
1298     // Initialize Jet Data
1299     for (i=0;i<fnAKTFullJets;i++)
1300     {
1301         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
1302
1303         fEMCalJet->SetIsJetInArray(IsInEMCal(myJet->Phi(),myJet->Eta()),i);
1304         fEMCalFullJet->SetIsJetInArray(IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta()),i);
1305         fEMCalPartJet->SetIsJetInArray(IsInEMCalPart(fJetR,myJet->Phi(),myJet->Eta()),i);
1306     }
1307     fEMCalJet->InitializeJetData(fmyAKTFullJets,fnAKTFullJets);
1308     fEMCalFullJet->InitializeJetData(fmyAKTFullJets,fnAKTFullJets);
1309     fEMCalPartJet->InitializeJetData(fmyAKTFullJets,fnAKTFullJets);
1310
1311     // kT Jets
1312     fEMCalkTFullJet = new AlipAJetData("fEMCalkTFullJet",kTRUE,fnKTFullJets);
1313     fEMCalkTFullJet->SetSignalCut(fEMCalJetThreshold);
1314     fEMCalkTFullJet->SetAreaCutFraction(0.25*fJetAreaCutFrac);
1315     fEMCalkTFullJet->SetJetR(fJetR);
1316
1317     for (i=0;i<fnKTFullJets;i++)
1318     {
1319         AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(i);
1320         fEMCalkTFullJet->SetIsJetInArray(IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta()),i);
1321     }
1322     fEMCalkTFullJet->InitializeJetData(fmyKTFullJets,fnKTFullJets);
1323
1324     // Raw Full Jet Spectra
1325     fEMCalRawJets->FillBSJS(fEventCentrality,0.0,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
1326 }
1327
1328 void AliAnalysisTaskFullpAJets::JetPtArea()
1329 {
1330     Int_t i;
1331
1332     for (i=0;i<fEMCalFullJet->GetTotalJets();i++)
1333     {
1334         AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalFullJet->GetJetIndex(i));
1335         fhJetPtArea->Fill(myJet->Pt(),myJet->Area());
1336     }
1337 }
1338
1339 void AliAnalysisTaskFullpAJets::GenerateTPCRandomConesPt()
1340 {
1341     Int_t i,j;
1342     Double_t E_tracks_total=0.;
1343     TRandom3 u(time(NULL));
1344
1345     Double_t Eta_Center=0.5*(fTPCEtaMin+fTPCEtaMax);
1346     Double_t Phi_Center=0.5*(fTPCPhiMin+fTPCPhiMax);
1347     Int_t event_mult=0;
1348     Int_t clus_mult=0;
1349
1350     for (i=0;i<fnBckgClusters;i++)
1351     {
1352         fTPCRCBckgFluc[i]=0.0;
1353         fTPCRCBckgFlucSignal[i]=0.0;
1354         fTPCRCBckgFlucNColl[i]=0.0;
1355     }
1356
1357     TLorentzVector *dummy= new TLorentzVector;
1358     TLorentzVector *temp_jet= new TLorentzVector;
1359
1360     // First, consider the RC with no spatial restrictions
1361     for (j=0;j<fnBckgClusters;j++)
1362     {
1363         E_tracks_total=0.;
1364
1365         dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1366         //  Loop over all tracks
1367         for (i=0;i<fnTracks;i++)
1368         {
1369             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1370             if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1371             {
1372                 TLorentzVector *track_vec = new TLorentzVector;
1373                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1374                 if (dummy->DeltaR(*track_vec)<fJetR)
1375                 {
1376                     E_tracks_total+=vtrack->Pt();
1377                 }
1378                 delete track_vec;
1379             }
1380         }
1381         fTPCRCBckgFlucSignal[j]=E_tracks_total;
1382     }
1383
1384     // Now, consider the RC where the vertex of RC is at least 2R away from the leading signal
1385     E_tracks_total=0.0;
1386     if (fTPCJet->GetLeadingPt()<0.0)
1387     {
1388         temp_jet->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1389     }
1390     else
1391     {
1392         AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1393         myJet->GetMom(*temp_jet);
1394     }
1395
1396     for (j=0;j<fnBckgClusters;j++)
1397     {
1398         event_mult=0;
1399         clus_mult=0;
1400         E_tracks_total=0.;
1401
1402         dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1403         while (temp_jet->DeltaR(*dummy)<fJetR)
1404         {
1405             dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1406         }
1407         //  Loop over all tracks
1408         for (i=0;i<fnTracks;i++)
1409         {
1410             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1411             if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1412             {
1413                 event_mult++;
1414                 TLorentzVector *track_vec = new TLorentzVector;
1415                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1416                 if (dummy->DeltaR(*track_vec)<fJetR)
1417                 {
1418                     clus_mult++;
1419                     E_tracks_total+=vtrack->Pt();
1420                 }
1421                 delete track_vec;
1422             }
1423         }
1424         fTPCRCBckgFluc[j]=E_tracks_total;
1425     }
1426     fpTPCEventMult->Fill(fEventCentrality,event_mult);
1427     fTPCRawJets->FillDeltaPt(fEventCentrality,0.0,fJetR,fTPCRCBckgFluc,1);
1428
1429     // For the case of partial exclusion, merely allow a superposition of full and no exclusion with probability p=1/Ncoll
1430     Double_t exclusion_prob;
1431     for (j=0;j<fnBckgClusters;j++)
1432     {
1433         exclusion_prob = u.Uniform(0,1);
1434         if (exclusion_prob<(1/fNColl))
1435         {
1436             fTPCRCBckgFlucNColl[j]=fTPCRCBckgFlucSignal[j];
1437         }
1438         else
1439         {
1440             fTPCRCBckgFlucNColl[j]=fTPCRCBckgFluc[j];
1441         }
1442     }
1443
1444     delete dummy;
1445     delete temp_jet;
1446 }
1447
1448 void AliAnalysisTaskFullpAJets::GenerateEMCalRandomConesPt()
1449 {
1450     Int_t i,j;
1451     Double_t E_tracks_total=0.;
1452     Double_t E_caloclusters_total=0.;
1453     TRandom3 u(time(NULL));
1454
1455     Double_t Eta_Center=0.5*(fEMCalEtaMin+fEMCalEtaMax);
1456     Double_t Phi_Center=0.5*(fEMCalPhiMin+fEMCalPhiMax);
1457     Int_t event_mult=0;
1458     Int_t clus_mult=0;
1459
1460     for (i=0;i<fnBckgClusters;i++)
1461     {
1462         fEMCalRCBckgFluc[i]=0.0;
1463         fEMCalRCBckgFlucSignal[i]=0.0;
1464         fEMCalRCBckgFlucNColl[i]=0.0;
1465     }
1466
1467     TLorentzVector *dummy= new TLorentzVector;
1468     TLorentzVector *temp_jet= new TLorentzVector;
1469
1470     // First, consider the RC with no spatial restrictions
1471     for (j=0;j<fnBckgClusters;j++)
1472     {
1473         E_tracks_total=0.;
1474         E_caloclusters_total=0.;
1475
1476         dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1477         //  Loop over all tracks
1478         for (i=0;i<fnTracks;i++)
1479         {
1480             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1481             if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1482             {
1483                 TLorentzVector *track_vec = new TLorentzVector;
1484                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1485                 if (dummy->DeltaR(*track_vec)<fJetR)
1486                 {
1487                     E_tracks_total+=vtrack->Pt();
1488                 }
1489                 delete track_vec;
1490             }
1491         }
1492
1493         //  Loop over all caloclusters
1494         for (i=0;i<fnClusters;i++)
1495         {
1496             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1497             TLorentzVector *cluster_vec = new TLorentzVector;
1498             vcluster->GetMomentum(*cluster_vec,fvertex);
1499             if (dummy->DeltaR(*cluster_vec)<fJetR)
1500             {
1501                 clus_mult++;
1502                 E_caloclusters_total+=vcluster->E();
1503             }
1504             delete cluster_vec;
1505         }
1506         fEMCalRCBckgFlucSignal[j]=E_tracks_total+E_caloclusters_total;
1507     }
1508
1509     // Now, consider the RC where the vertex of RC is at least 2R away from the leading signal
1510     E_tracks_total=0.;
1511     E_caloclusters_total=0.;
1512     if (fEMCalPartJet->GetLeadingPt()<0.0)
1513     {
1514         temp_jet->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1515     }
1516     else
1517     {
1518         AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
1519         myJet->GetMom(*temp_jet);
1520     }
1521
1522     for (j=0;j<fnBckgClusters;j++)
1523     {
1524         event_mult=0;
1525         clus_mult=0;
1526         E_tracks_total=0.;
1527         E_caloclusters_total=0.;
1528
1529         dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1530         while (temp_jet->DeltaR(*dummy)<fJetR)
1531         {
1532             dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1533         }
1534         //  Loop over all tracks
1535         for (i=0;i<fnTracks;i++)
1536         {
1537             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1538             if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1539             {
1540                 event_mult++;
1541                 TLorentzVector *track_vec = new TLorentzVector;
1542                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1543                 if (dummy->DeltaR(*track_vec)<fJetR)
1544                 {
1545                     clus_mult++;
1546                     E_tracks_total+=vtrack->Pt();
1547                 }
1548                 delete track_vec;
1549             }
1550         }
1551
1552         //  Loop over all caloclusters
1553         for (i=0;i<fnClusters;i++)
1554         {
1555             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1556             TLorentzVector *cluster_vec = new TLorentzVector;
1557             vcluster->GetMomentum(*cluster_vec,fvertex);
1558             event_mult++;
1559             if (dummy->DeltaR(*cluster_vec)<fJetR)
1560             {
1561                 clus_mult++;
1562                 E_caloclusters_total+=vcluster->E();
1563             }
1564             delete cluster_vec;
1565         }
1566         fEMCalRCBckgFluc[j]=E_tracks_total+E_caloclusters_total;
1567     }
1568     fpEMCalEventMult->Fill(fEventCentrality,event_mult);
1569     fEMCalRawJets->FillDeltaPt(fEventCentrality,0.0,fJetR,fEMCalRCBckgFluc,1);
1570
1571     // For the case of partial exclusion, merely allow a superposition of full and no exclusion with probability p=1/Ncoll
1572     Double_t exclusion_prob;
1573     for (j=0;j<fnBckgClusters;j++)
1574     {
1575         exclusion_prob = u.Uniform(0,1);
1576         if (exclusion_prob<(1/fNColl))
1577         {
1578             fEMCalRCBckgFlucNColl[j]=fEMCalRCBckgFlucSignal[j];
1579         }
1580         else
1581         {
1582             fEMCalRCBckgFlucNColl[j]=fEMCalRCBckgFluc[j];
1583         }
1584     }
1585
1586     delete dummy;
1587     delete temp_jet;
1588 }
1589
1590 // Charged Rho's
1591 void AliAnalysisTaskFullpAJets::EstimateChargedRho0()
1592 {
1593     Int_t i;
1594     Double_t E_tracks_total=0.0;
1595     Double_t TPC_rho=0.;
1596
1597     //  Loop over all tracks
1598     for (i=0;i<fnTracks;i++)
1599     {
1600         AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1601         if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1602         {
1603             E_tracks_total+=vtrack->Pt();
1604         }
1605     }
1606
1607     //  Calculate the mean Background density
1608     TPC_rho=E_tracks_total/fTPCArea;
1609     fRhoCharged=TPC_rho;
1610
1611     // Fill Histograms
1612     fRhoCharged0->FillRho(fEventCentrality,TPC_rho);
1613     fRhoCharged0->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCJet->GetJets(),fTPCJet->GetTotalJets());
1614     fRhoCharged0->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
1615     fRhoCharged0->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
1616     fRhoCharged0->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucNColl,1);
1617     fRhoCharged0->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
1618     fRhoCharged0->FillLeadingJetPtRho(fTPCJet->GetLeadingPt(),TPC_rho);
1619
1620 }
1621
1622 void AliAnalysisTaskFullpAJets::EstimateChargedRho1()
1623 {
1624     Int_t i;
1625     Double_t E_tracks_total=0.0;
1626     Double_t TPC_rho=0.;
1627
1628     if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
1629     {
1630         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1631         TLorentzVector *temp_jet= new TLorentzVector;
1632         myJet->GetMom(*temp_jet);
1633
1634         //  Loop over all tracks
1635         for (i=0;i<fnTracks;i++)
1636         {
1637             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1638             if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1639             {
1640                 TLorentzVector *track_vec = new TLorentzVector;
1641                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1642                 if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
1643                 {
1644                     E_tracks_total+=vtrack->Pt();
1645                 }
1646                 delete track_vec;
1647             }
1648         }
1649         delete temp_jet;
1650
1651         //  Calculate the mean Background density
1652         TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myJet->Eta()));
1653     }
1654     else  // i.e. No signal jets -> same as total background density
1655     {
1656         //  Loop over all tracks
1657         for (i=0;i<fnTracks;i++)
1658         {
1659             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1660             if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1661             {
1662                 E_tracks_total+=vtrack->Pt();
1663             }
1664         }
1665         //  Calculate the mean Background density
1666         TPC_rho=E_tracks_total/fTPCArea;
1667     }
1668
1669     // Fill histograms
1670     fRhoCharged1->FillRho(fEventCentrality,TPC_rho);
1671     fRhoCharged1->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
1672     fRhoCharged1->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
1673     fRhoCharged1->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
1674     fRhoCharged1->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucNColl,1);
1675     fRhoCharged1->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
1676     fRhoCharged1->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),TPC_rho);
1677 }
1678
1679 void AliAnalysisTaskFullpAJets::EstimateChargedRho2()
1680 {
1681     Int_t i;
1682     Double_t E_tracks_total=0.0;
1683     Double_t TPC_rho=0.;
1684
1685     if ((fTPCJet->GetLeadingPt()>=fTPCJetThreshold) && (fTPCJet->GetSubLeadingPt()>=fTPCJetThreshold))
1686     {
1687         AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1688         TLorentzVector *temp_jet1= new TLorentzVector;
1689         myhJet->GetMom(*temp_jet1);
1690
1691         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSubLeadingIndex());
1692         TLorentzVector *temp_jet2= new TLorentzVector;
1693         myJet->GetMom(*temp_jet2);
1694
1695         //  Loop over all tracks
1696         for (i=0;i<fnTracks;i++)
1697         {
1698             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1699             if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1700             {
1701                 TLorentzVector *track_vec = new TLorentzVector;
1702                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1703                 if ((temp_jet1->DeltaR(*track_vec)>fJetRForRho) && (temp_jet2->DeltaR(*track_vec)>fJetRForRho))
1704                 {
1705                     E_tracks_total+=vtrack->Pt();
1706                 }
1707                 delete track_vec;
1708             }
1709         }
1710         delete temp_jet1;
1711         delete temp_jet2;
1712
1713         //  Calculate the mean Background density
1714         TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myhJet->Eta())-AreaWithinTPC(fJetR,myJet->Eta()));
1715     }
1716     else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
1717     {
1718         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1719         TLorentzVector *temp_jet= new TLorentzVector;
1720         myJet->GetMom(*temp_jet);
1721
1722         //  Loop over all tracks
1723         for (i=0;i<fnTracks;i++)
1724         {
1725             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1726             if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1727             {
1728                 TLorentzVector *track_vec = new TLorentzVector;
1729                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1730                 if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
1731                 {
1732                     E_tracks_total+=vtrack->Pt();
1733                 }
1734                 delete track_vec;
1735             }
1736         }
1737         delete temp_jet;
1738
1739         //  Calculate the mean Background density
1740         TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myJet->Eta()));
1741     }
1742     else  // i.e. No signal jets -> same as total background density
1743     {
1744         //  Loop over all tracks
1745         for (i=0;i<fnTracks;i++)
1746         {
1747             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1748             if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1749             {
1750                 E_tracks_total+=vtrack->Pt();
1751             }
1752         }
1753
1754         //  Calculate the mean Background density
1755         TPC_rho=E_tracks_total/fTPCArea;
1756     }
1757
1758     // Fill histograms
1759     fRhoCharged2->FillRho(fEventCentrality,TPC_rho);
1760     fRhoCharged2->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
1761     fRhoCharged2->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
1762     fRhoCharged2->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
1763     fRhoCharged2->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucNColl,1);
1764     fRhoCharged2->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
1765     fRhoCharged2->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),TPC_rho);
1766 }
1767
1768 void AliAnalysisTaskFullpAJets::EstimateChargedRhoN()
1769 {
1770     Int_t i,j;
1771     Bool_t track_away_from_jet;
1772     Double_t E_tracks_total=0.0;
1773     Double_t TPC_rho=0.0;
1774     Double_t jet_area_total=0.0;
1775
1776     // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
1777     for (i=0;i<fnTracks;i++)
1778     {
1779         // First, check if track is in the EMCal!!
1780         AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
1781         if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1782         {
1783             if (fTPCJet->GetTotalSignalJets()<1)
1784             {
1785                 E_tracks_total+=vtrack->Pt();
1786             }
1787             else
1788             {
1789                 track_away_from_jet=kTRUE;
1790                 j=0;
1791                 TLorentzVector *track_vec = new TLorentzVector;
1792                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1793                 while (track_away_from_jet==kTRUE && j<fTPCJet->GetTotalSignalJets())
1794                 {
1795                     TLorentzVector *jet_vec= new TLorentzVector;
1796                     AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(j));
1797                     myJet->GetMom(*jet_vec);
1798                     if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
1799                     {
1800                         track_away_from_jet=kFALSE;
1801                     }
1802                     delete jet_vec;
1803                     j++;
1804                 }
1805                 if (track_away_from_jet==kTRUE)
1806                 {
1807                     E_tracks_total+=vtrack->Pt();
1808                 }
1809                 delete track_vec;
1810             }
1811         }
1812     }
1813
1814     // Determine area of all Jets that are within the EMCal
1815     if (fTPCJet->GetTotalSignalJets()==0)
1816     {
1817         jet_area_total=0.0;
1818     }
1819     else
1820     {
1821         for (i=0;i<fTPCJet->GetTotalSignalJets();i++)
1822         {
1823             AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(i));
1824             jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
1825         }
1826     }
1827
1828     // Calculate Rho
1829     TPC_rho = E_tracks_total/(fTPCArea-jet_area_total);
1830
1831     // Fill Histogram
1832     fRhoChargedN->FillRho(fEventCentrality,TPC_rho);
1833     fRhoChargedN->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
1834     fRhoChargedN->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
1835     fRhoChargedN->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
1836     fRhoChargedN->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucNColl,1);
1837     fRhoChargedN->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
1838     fRhoChargedN->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),TPC_rho);
1839 }
1840
1841 void AliAnalysisTaskFullpAJets::EstimateChargedRhoScale()
1842 {
1843     Int_t i,j;
1844     Bool_t track_away_from_jet;
1845     Double_t E_tracks_total=0.0;
1846     Double_t TPC_rho=0.0;
1847     Double_t jet_area_total=0.0;
1848
1849     // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
1850     for (i=0;i<fnTracks;i++)
1851     {
1852         // First, check if track is in the EMCal!!
1853         AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
1854         if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1855         {
1856             if (fTPCJet->GetTotalSignalJets()<1)
1857             {
1858                 E_tracks_total+=vtrack->Pt();
1859             }
1860             else
1861             {
1862                 track_away_from_jet=kTRUE;
1863                 j=0;
1864                 TLorentzVector *track_vec = new TLorentzVector;
1865                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1866                 while (track_away_from_jet==kTRUE && j<fTPCJet->GetTotalSignalJets())
1867                 {
1868                     TLorentzVector *jet_vec= new TLorentzVector;
1869                     AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(j));
1870                     myJet->GetMom(*jet_vec);
1871                     if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
1872                     {
1873                         track_away_from_jet=kFALSE;
1874                     }
1875                     delete jet_vec;
1876                     j++;
1877                 }
1878                 if (track_away_from_jet==kTRUE)
1879                 {
1880                     E_tracks_total+=vtrack->Pt();
1881                 }
1882                 delete track_vec;
1883             }
1884         }
1885     }
1886
1887     // Determine area of all Jets that are within the TPC
1888     if (fTPCJet->GetTotalSignalJets()==0)
1889     {
1890         jet_area_total=0.0;
1891     }
1892     else
1893     {
1894         for (i=0;i<fTPCJet->GetTotalSignalJets();i++)
1895         {
1896             AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(i));
1897             jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
1898         }
1899     }
1900
1901     // Calculate Rho
1902     TPC_rho = E_tracks_total/(fTPCArea-jet_area_total);
1903     TPC_rho*=fScaleFactor;
1904
1905     // Fill Histogram
1906     fRhoChargedScale->FillRho(fEventCentrality,TPC_rho);
1907     fRhoChargedScale->FillBSJS(fEventCentrality,TPC_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
1908     fRhoChargedScale->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fEMCalRCBckgFluc,1);
1909     fRhoChargedScale->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fEMCalRCBckgFlucSignal,1);
1910     fRhoChargedScale->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fEMCalRCBckgFlucNColl,1);
1911     fRhoChargedScale->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
1912     fRhoChargedScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),TPC_rho);
1913 }
1914
1915 void AliAnalysisTaskFullpAJets::EstimateChargedRhokT()
1916 {
1917     Int_t i;
1918     Double_t kTRho = 0.0;
1919     Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1920     Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1921
1922     for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1923     {
1924         AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1925         pTArray[i]=myJet->Pt();
1926         RhoArray[i]=myJet->Pt()/myJet->Area();
1927     }
1928
1929     if (fTPCkTFullJet->GetTotalJets()>=2)
1930     {
1931         kTRho=MedianRhokT(pTArray,RhoArray,fTPCkTFullJet->GetTotalJets());
1932
1933         fRhoChargedkT->FillRho(fEventCentrality,kTRho);
1934         fRhoChargedkT->FillBSJS(fEventCentrality,kTRho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
1935         fRhoChargedkT->FillDeltaPt(fEventCentrality,kTRho,fJetR,fTPCRCBckgFluc,1);
1936         fRhoChargedkT->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fTPCRCBckgFlucSignal,1);
1937         fRhoChargedkT->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fTPCRCBckgFlucNColl,1);
1938         fRhoChargedkT->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
1939         fRhoChargedkT->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),kTRho);
1940     }
1941     delete [] RhoArray;
1942     delete [] pTArray;
1943 }
1944
1945 void AliAnalysisTaskFullpAJets::EstimateChargedRhokTScale()
1946 {
1947     Int_t i;
1948     Double_t kTRho = 0.0;
1949     Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1950     Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1951
1952     for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1953     {
1954         AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1955         pTArray[i]=myJet->Pt();
1956         RhoArray[i]=myJet->Pt()/myJet->Area();
1957     }
1958
1959     if (fTPCkTFullJet->GetTotalJets()>=2)
1960     {
1961         kTRho=MedianRhokT(pTArray,RhoArray,fTPCkTFullJet->GetTotalJets());
1962         kTRho*=fScaleFactor;
1963
1964         fRhoChargedkTScale->FillRho(fEventCentrality,kTRho);
1965         fRhoChargedkTScale->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
1966         fRhoChargedkTScale->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
1967         fRhoChargedkTScale->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
1968         fRhoChargedkTScale->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucNColl,1);
1969         fRhoChargedkTScale->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
1970         fRhoChargedkTScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
1971     }
1972     delete [] RhoArray;
1973     delete [] pTArray;
1974 }
1975
1976 void AliAnalysisTaskFullpAJets::EstimateChargedRhoCMS()
1977 {
1978     Int_t i,k;
1979     Double_t kTRho = 0.0;
1980     Double_t CMSTotalkTArea = 0.0;
1981     Double_t CMSTrackArea = 0.0;
1982     Double_t CMSCorrectionFactor = 1.0;
1983     Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1984     Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1985
1986     k=0;
1987     if ((fTPCJet->GetLeadingPt()>=fTPCJetThreshold) && (fTPCJet->GetSubLeadingPt()>=fTPCJetThreshold))
1988     {
1989         AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1990         AliEmcalJet *myJet2 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSubLeadingIndex());
1991
1992         for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1993         {
1994             AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1995
1996             CMSTotalkTArea+=myJet->Area();
1997             if (myJet->GetNumberOfTracks()>0)
1998             {
1999                 CMSTrackArea+=myJet->Area();
2000             }
2001             if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kFALSE)
2002             {
2003                 pTArray[k]=myJet->Pt();
2004                 RhoArray[k]=myJet->Pt()/myJet->Area();
2005                 k++;
2006             }
2007         }
2008         if (k>0)
2009         {
2010             kTRho=MedianRhokT(pTArray,RhoArray,k);
2011         }
2012         else
2013         {
2014             kTRho=0.0;
2015         }
2016     }
2017     else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
2018     {
2019         AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
2020
2021         for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2022         {
2023             AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
2024
2025             CMSTotalkTArea+=myJet->Area();
2026             if (myJet->GetNumberOfTracks()>0)
2027             {
2028                 CMSTrackArea+=myJet->Area();
2029             }
2030             if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE)
2031             {
2032                 pTArray[k]=myJet->Pt();
2033                 RhoArray[k]=myJet->Pt()/myJet->Area();
2034                 k++;
2035             }
2036         }
2037         if (k>0)
2038         {
2039             kTRho=MedianRhokT(pTArray,RhoArray,k);
2040         }
2041         else
2042         {
2043             kTRho=0.0;
2044         }
2045     }
2046     else
2047     {
2048         for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2049         {
2050             AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
2051
2052             CMSTotalkTArea+=myJet->Area();
2053             if (myJet->GetNumberOfTracks()>0)
2054             {
2055                 CMSTrackArea+=myJet->Area();
2056             }
2057             pTArray[k]=myJet->Pt();
2058             RhoArray[k]=myJet->Pt()/myJet->Area();
2059             k++;
2060         }
2061         if (k>0)
2062         {
2063             kTRho=MedianRhokT(pTArray,RhoArray,k);
2064         }
2065         else
2066         {
2067             kTRho=0.0;
2068         }
2069     }
2070     // Scale CMS Rho by Correction factor
2071     if (CMSTotalkTArea==0.0)
2072     {
2073         CMSCorrectionFactor = 1.0;
2074     }
2075     else
2076     {
2077         //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
2078         CMSCorrectionFactor = CMSTrackArea/(fTPCPhiTotal*(fTPCEtaTotal-2*fJetR));  // The total physical area should be reduced by the eta cut due to looping over only fully contained kT jets within the TPC
2079     }
2080     kTRho*=CMSCorrectionFactor;
2081     fRhoChargedCMS->FillRho(fEventCentrality,kTRho);
2082     fRhoChargedCMS->FillBSJS(fEventCentrality,kTRho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
2083     fRhoChargedCMS->FillDeltaPt(fEventCentrality,kTRho,fJetR,fTPCRCBckgFluc,1);
2084     fRhoChargedCMS->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fTPCRCBckgFlucSignal,1);
2085     fRhoChargedCMS->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fTPCRCBckgFlucNColl,1);
2086     fRhoChargedCMS->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
2087     fRhoChargedCMS->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),kTRho);
2088     delete [] RhoArray;
2089     delete [] pTArray;
2090 }
2091
2092 void AliAnalysisTaskFullpAJets::EstimateChargedRhoCMSScale()
2093 {
2094     Int_t i,k;
2095     Double_t kTRho = 0.0;
2096     Double_t CMSTotalkTArea = 0.0;
2097     Double_t CMSTrackArea = 0.0;
2098     Double_t CMSCorrectionFactor = 1.0;
2099     Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2100     Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2101
2102     k=0;
2103     if ((fTPCJet->GetLeadingPt()>=fTPCJetThreshold) && (fTPCJet->GetSubLeadingPt()>=fTPCJetThreshold))
2104     {
2105         AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
2106         AliEmcalJet *myJet2 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSubLeadingIndex());
2107
2108         for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2109         {
2110             AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
2111
2112             CMSTotalkTArea+=myJet->Area();
2113             if (myJet->GetNumberOfTracks()>0)
2114             {
2115                 CMSTrackArea+=myJet->Area();
2116             }
2117             if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kFALSE)
2118             {
2119                 pTArray[k]=myJet->Pt();
2120                 RhoArray[k]=myJet->Pt()/myJet->Area();
2121                 k++;
2122             }
2123         }
2124         if (k>0)
2125         {
2126             kTRho=MedianRhokT(pTArray,RhoArray,k);
2127         }
2128         else
2129         {
2130             kTRho=0.0;
2131         }
2132     }
2133     else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
2134     {
2135         AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
2136
2137         for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2138         {
2139             AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
2140
2141             CMSTotalkTArea+=myJet->Area();
2142             if (myJet->GetNumberOfTracks()>0)
2143             {
2144                 CMSTrackArea+=myJet->Area();
2145             }
2146             if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE)
2147             {
2148                 pTArray[k]=myJet->Pt();
2149                 RhoArray[k]=myJet->Pt()/myJet->Area();
2150                 k++;
2151             }
2152         }
2153         if (k>0)
2154         {
2155             kTRho=MedianRhokT(pTArray,RhoArray,k);
2156         }
2157         else
2158         {
2159             kTRho=0.0;
2160         }
2161     }
2162     else
2163     {
2164         for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2165         {
2166             AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
2167
2168             CMSTotalkTArea+=myJet->Area();
2169             if (myJet->GetNumberOfTracks()>0)
2170             {
2171                 CMSTrackArea+=myJet->Area();
2172             }
2173             pTArray[k]=myJet->Pt();
2174             RhoArray[k]=myJet->Pt()/myJet->Area();
2175             k++;
2176         }
2177         if (k>0)
2178         {
2179             kTRho=MedianRhokT(pTArray,RhoArray,k);
2180         }
2181         else
2182         {
2183             kTRho=0.0;
2184         }
2185     }
2186     kTRho*=fScaleFactor;
2187     // Scale CMS Rho by Correction factor
2188     if (CMSTotalkTArea==0.0)
2189     {
2190         CMSCorrectionFactor = 1.0;
2191     }
2192     else
2193     {
2194         //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
2195         CMSCorrectionFactor = CMSTrackArea/(fTPCPhiTotal*(fTPCEtaTotal-2*fJetR));  // The total physical area should be reduced by the eta cut due to looping over only fully contained kT jets within the TPC
2196     }
2197     kTRho*=CMSCorrectionFactor;
2198
2199     fRhoChargedCMSScale->FillRho(fEventCentrality,kTRho);
2200     fRhoChargedCMSScale->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2201     fRhoChargedCMSScale->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
2202     fRhoChargedCMSScale->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
2203     fRhoChargedCMSScale->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucNColl,1);
2204     fRhoChargedCMSScale->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
2205     fRhoChargedCMSScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
2206     delete [] RhoArray;
2207     delete [] pTArray;
2208 }
2209
2210 // Full Rho's
2211 void AliAnalysisTaskFullpAJets::EstimateFullRho0()
2212 {
2213     Int_t i;
2214     Double_t E_tracks_total=0.0;
2215     Double_t E_caloclusters_total=0.0;
2216     Double_t EMCal_rho=0.0;
2217
2218     //  Loop over all tracks
2219     for (i=0;i<fnTracks;i++)
2220     {
2221         AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2222         if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2223         {
2224             E_tracks_total+=vtrack->Pt();
2225         }
2226     }
2227
2228     //  Loop over all caloclusters
2229     for (i=0;i<fnClusters;i++)
2230     {
2231         AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2232         TLorentzVector *cluster_vec = new TLorentzVector;
2233         vcluster->GetMomentum(*cluster_vec,fvertex);
2234         E_caloclusters_total+=cluster_vec->Pt();
2235         //E_caloclusters_total+=0.5*cluster_vec->Pt();
2236     }
2237
2238     //  Calculate the mean Background density
2239     EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2240     fRhoFull=EMCal_rho;
2241
2242     // Fill Histograms
2243     if (fRhoCharged>0)
2244     {
2245         fpRhoScale->Fill(fEventCentrality,fRhoFull/fRhoCharged);
2246         fhRhoScale->Fill(fRhoFull/fRhoCharged,fEventCentrality);
2247     }
2248
2249     fRhoFull0->FillRho(fEventCentrality,EMCal_rho);
2250     fRhoFull0->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2251     fRhoFull0->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
2252     fRhoFull0->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
2253     fRhoFull0->FillDeltaPtNColl(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucNColl,1);
2254     fRhoFull0->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
2255     fRhoFull0->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
2256 }
2257
2258 void AliAnalysisTaskFullpAJets::EstimateFullRho1()
2259 {
2260     Int_t i;
2261     Double_t E_tracks_total=0.0;
2262     Double_t E_caloclusters_total=0.0;
2263     Double_t EMCal_rho=0.0;
2264
2265     if (fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold)
2266     {
2267         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
2268         TLorentzVector *temp_jet= new TLorentzVector;
2269         myJet->GetMom(*temp_jet);
2270
2271         //  Loop over all tracks
2272         for (i=0;i<fnTracks;i++)
2273         {
2274             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2275             if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2276             {
2277                 TLorentzVector *track_vec = new TLorentzVector;
2278                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2279                 if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
2280                 {
2281                     E_tracks_total+=vtrack->Pt();
2282                 }
2283                 delete track_vec;
2284             }
2285         }
2286
2287         //  Loop over all caloclusters
2288         for (i=0;i<fnClusters;i++)
2289         {
2290             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2291             TLorentzVector *cluster_vec = new TLorentzVector;
2292             vcluster->GetMomentum(*cluster_vec,fvertex);
2293             if (temp_jet->DeltaR(*cluster_vec)>fJetRForRho)
2294             {
2295                 E_caloclusters_total+=vcluster->E();
2296             }
2297             delete cluster_vec;
2298         }
2299         delete temp_jet;
2300         //  Calculate the mean Background density
2301         EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2302     }
2303     else  // i.e. No signal jets -> same as total background density
2304     {
2305         //  Loop over all tracks
2306         for (i=0;i<fnTracks;i++)
2307         {
2308             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2309             if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2310             {
2311                 E_tracks_total+=vtrack->Pt();
2312             }
2313         }
2314
2315         //  Loop over all caloclusters
2316         for (i=0;i<fnClusters;i++)
2317         {
2318             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2319             E_caloclusters_total+=vcluster->E();
2320         }
2321         //  Calculate the mean Background density
2322         EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2323     }
2324
2325     // Fill histograms
2326     fRhoFull1->FillRho(fEventCentrality,EMCal_rho);
2327     fRhoFull1->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2328     fRhoFull1->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
2329     fRhoFull1->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
2330     fRhoFull1->FillDeltaPtNColl(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucNColl,1);
2331     fRhoFull1->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
2332     fRhoFull1->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
2333 }
2334
2335 void AliAnalysisTaskFullpAJets::EstimateFullRho2()
2336 {
2337     Int_t i;
2338     Double_t E_tracks_total=0.0;
2339     Double_t E_caloclusters_total=0.0;
2340     Double_t EMCal_rho=0.0;
2341
2342     if ((fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold) && (fEMCalPartJet->GetSubLeadingPt()>=fEMCalJetThreshold))
2343     {
2344         AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
2345         TLorentzVector *temp_jet1 = new TLorentzVector;
2346         myhJet->GetMom(*temp_jet1);
2347
2348         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSubLeadingIndex());
2349         TLorentzVector *temp_jet2 = new TLorentzVector;
2350         myJet->GetMom(*temp_jet2);
2351
2352         //  Loop over all tracks
2353         for (i=0;i<fnTracks;i++)
2354         {
2355             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2356             if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2357             {
2358                 TLorentzVector *track_vec = new TLorentzVector;
2359                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2360                 if ((temp_jet1->DeltaR(*track_vec)>fJetRForRho) && (temp_jet2->DeltaR(*track_vec)>fJetRForRho))
2361                 {
2362                     E_tracks_total+=vtrack->Pt();
2363                 }
2364                 delete track_vec;
2365             }
2366         }
2367
2368         //  Loop over all caloclusters
2369         for (i=0;i<fnClusters;i++)
2370         {
2371             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2372             TLorentzVector *cluster_vec = new TLorentzVector;
2373             vcluster->GetMomentum(*cluster_vec,fvertex);
2374             if ((temp_jet1->DeltaR(*cluster_vec)>fJetRForRho) && (temp_jet2->DeltaR(*cluster_vec)>fJetRForRho))
2375             {
2376                 E_caloclusters_total+=vcluster->E();
2377             }
2378             delete cluster_vec;
2379         }
2380         delete temp_jet1;
2381         delete temp_jet2;
2382
2383         //  Calculate the mean Background density
2384         EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myhJet->Phi(),myhJet->Eta())-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2385     }
2386     else if (fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold)
2387     {
2388         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
2389         TLorentzVector *temp_jet= new TLorentzVector;
2390         myJet->GetMom(*temp_jet);
2391
2392         //  Loop over all tracks
2393         for (i=0;i<fnTracks;i++)
2394         {
2395             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2396             if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2397             {
2398                 TLorentzVector *track_vec = new TLorentzVector;
2399                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2400                 if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
2401                 {
2402                     E_tracks_total+=vtrack->Pt();
2403                 }
2404                 delete track_vec;
2405             }
2406         }
2407
2408         //  Loop over all caloclusters
2409         for (i=0;i<fnClusters;i++)
2410         {
2411             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2412             TLorentzVector *cluster_vec = new TLorentzVector;
2413             vcluster->GetMomentum(*cluster_vec,fvertex);
2414             if (temp_jet->DeltaR(*cluster_vec)>fJetRForRho)
2415             {
2416                 E_caloclusters_total+=vcluster->E();
2417             }
2418             delete cluster_vec;
2419         }
2420         delete temp_jet;
2421         //  Calculate the mean Background density
2422         EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2423     }
2424     else  // i.e. No signal jets -> same as total background density
2425     {
2426         //  Loop over all tracks
2427         for (i=0;i<fnTracks;i++)
2428         {
2429             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2430             if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2431             {
2432                 E_tracks_total+=vtrack->Pt();
2433             }
2434         }
2435
2436         //  Loop over all caloclusters
2437         for (i=0;i<fnClusters;i++)
2438         {
2439             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2440             E_caloclusters_total+=vcluster->E();
2441         }
2442         //  Calculate the mean Background density
2443         EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2444     }
2445
2446     // Fill histograms
2447     fRhoFull2->FillRho(fEventCentrality,EMCal_rho);
2448     fRhoFull2->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2449     fRhoFull2->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
2450     fRhoFull2->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
2451     fRhoFull2->FillDeltaPtNColl(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucNColl,1);
2452     fRhoFull2->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
2453     fRhoFull2->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
2454 }
2455
2456 void AliAnalysisTaskFullpAJets::EstimateFullRhoN()
2457 {
2458     Int_t i,j;
2459     Bool_t track_away_from_jet;
2460     Bool_t cluster_away_from_jet;
2461     Double_t E_tracks_total=0.0;
2462     Double_t E_caloclusters_total=0.0;
2463     Double_t EMCal_rho=0.0;
2464     Double_t jet_area_total=0.0;
2465
2466     // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
2467     for (i=0;i<fnTracks;i++)
2468     {
2469         // First, check if track is in the EMCal!!
2470         AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
2471         if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2472         {
2473             if (fEMCalPartJet->GetTotalSignalJets()<1)
2474             {
2475                 E_tracks_total+=vtrack->Pt();
2476             }
2477             else
2478             {
2479                 track_away_from_jet=kTRUE;
2480                 j=0;
2481                 TLorentzVector *track_vec = new TLorentzVector;
2482                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2483                 while (track_away_from_jet==kTRUE && j<fEMCalPartJet->GetTotalSignalJets())
2484                 {
2485                     TLorentzVector *jet_vec= new TLorentzVector;
2486                     AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSignalJetIndex(j));
2487                     myJet->GetMom(*jet_vec);
2488                     if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
2489                     {
2490                         track_away_from_jet=kFALSE;
2491                     }
2492                     delete jet_vec;
2493                     j++;
2494                 }
2495                 if (track_away_from_jet==kTRUE)
2496                 {
2497                     E_tracks_total+=vtrack->Pt();
2498                 }
2499                 delete track_vec;
2500             }
2501         }
2502     }
2503
2504     // Next, sum all CaloClusters within the EMCal (obviously all clusters must be within EMCal!!) that are away from jet(s) above Pt Threshold
2505     for (i=0;i<fnClusters;i++)
2506     {
2507         AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2508         if (fEMCalPartJet->GetTotalSignalJets()<1)
2509         {
2510             E_caloclusters_total+=vcluster->E();
2511         }
2512         else
2513         {
2514             cluster_away_from_jet=kTRUE;
2515             j=0;
2516
2517             TLorentzVector *cluster_vec = new TLorentzVector;
2518             vcluster->GetMomentum(*cluster_vec,fvertex);
2519             while (cluster_away_from_jet==kTRUE && j<fEMCalPartJet->GetTotalSignalJets())
2520             {
2521                 TLorentzVector *jet_vec= new TLorentzVector;
2522                 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSignalJetIndex(j));
2523                 myJet->GetMom(*jet_vec);
2524                 if (cluster_vec->DeltaR(*jet_vec)<=fJetRForRho)
2525                 {
2526                     cluster_away_from_jet=kFALSE;
2527                 }
2528                 delete jet_vec;
2529                 j++;
2530             }
2531             if (cluster_away_from_jet==kTRUE)
2532             {
2533                 E_caloclusters_total+=vcluster->E();
2534             }
2535             delete cluster_vec;
2536         }
2537     }
2538
2539     // Determine area of all Jets that are within the EMCal
2540     if (fEMCalPartJet->GetTotalSignalJets()==0)
2541     {
2542         jet_area_total=0.0;
2543     }
2544     else
2545     {
2546         for (i=0;i<fEMCalPartJet->GetTotalSignalJets();i++)
2547         {
2548             AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSignalJetIndex(i));
2549             jet_area_total+=AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta());
2550         }
2551     }
2552
2553     // Calculate Rho
2554     EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-jet_area_total);
2555
2556     // Fill Histogram
2557     fRhoFullN->FillRho(fEventCentrality,EMCal_rho);
2558     fRhoFullN->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2559     fRhoFullN->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
2560     fRhoFullN->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
2561     fRhoFullN->FillDeltaPtNColl(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucNColl,1);
2562     fRhoFullN->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
2563     fRhoFullN->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
2564 }
2565
2566 void AliAnalysisTaskFullpAJets::EstimateFullRhoDijet()
2567 {
2568     Int_t i;
2569     Double_t E_tracks_total=0.0;
2570     Double_t E_caloclusters_total=0.0;
2571     Double_t EMCal_rho=0.0;
2572
2573     //  Loop over all tracks
2574     for (i=0;i<fnTracks;i++)
2575     {
2576         AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2577         if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2578         {
2579             E_tracks_total+=vtrack->Pt();
2580         }
2581     }
2582
2583     //  Loop over all caloclusters
2584     for (i=0;i<fnClusters;i++)
2585     {
2586         AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2587         E_caloclusters_total+=vcluster->E();
2588     }
2589
2590     //  Calculate the mean Background density
2591     EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2592
2593     // Fill Histograms
2594     fRhoFullDijet->FillRho(fEventCentrality,EMCal_rho);
2595     fRhoFullDijet->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2596     fRhoFullDijet->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
2597     fRhoFullDijet->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
2598     fRhoFullDijet->FillDeltaPtNColl(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucNColl,1);
2599     fRhoFullDijet->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
2600     fRhoFullDijet->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
2601 }
2602
2603 void AliAnalysisTaskFullpAJets::EstimateFullRhokT()
2604 {
2605     Int_t i;
2606     Double_t kTRho = 0.0;
2607     Double_t *pTArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2608     Double_t *RhoArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2609
2610     for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2611     {
2612         AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
2613         pTArray[i]=myJet->Pt();
2614         RhoArray[i]=myJet->Pt()/myJet->Area();
2615     }
2616
2617     if (fEMCalkTFullJet->GetTotalJets()>0)
2618     {
2619         kTRho=MedianRhokT(pTArray,RhoArray,fEMCalkTFullJet->GetTotalJets());
2620     }
2621     else
2622     {
2623         kTRho=0.0;
2624     }
2625     fRhoFullkT->FillRho(fEventCentrality,kTRho);
2626     fRhoFullkT->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2627     fRhoFullkT->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
2628     fRhoFullkT->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
2629     fRhoFullkT->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucNColl,1);
2630     fRhoFullkT->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
2631     fRhoFullkT->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
2632     delete [] RhoArray;
2633     delete [] pTArray;
2634 }
2635
2636 void AliAnalysisTaskFullpAJets::EstimateFullRhoCMS()
2637 {
2638     Int_t i,k;
2639     Double_t kTRho = 0.0;
2640     Double_t CMSTotalkTArea = 0.0;
2641     Double_t CMSParticleArea = 0.0;
2642     Double_t CMSCorrectionFactor = 1.0;
2643     Double_t *pTArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2644     Double_t *RhoArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2645
2646     k=0;
2647     if ((fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold) && (fEMCalPartJet->GetSubLeadingPt()>=fEMCalJetThreshold))
2648     {
2649         AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
2650         AliEmcalJet *myJet2 =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSubLeadingIndex());
2651
2652         for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2653         {
2654             AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
2655
2656             CMSTotalkTArea+=myJet->Area();
2657             if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2658             {
2659                 CMSParticleArea+=myJet->Area();
2660             }
2661             if (IsJetOverlap(myJet,myJet1,kTRUE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kTRUE)
2662             {
2663                 pTArray[k]=myJet->Pt();
2664                 RhoArray[k]=myJet->Pt()/myJet->Area();
2665                 k++;
2666             }
2667         }
2668         if (k>0)
2669         {
2670             kTRho=MedianRhokT(pTArray,RhoArray,k);
2671         }
2672         else
2673         {
2674             kTRho=0.0;
2675         }
2676     }
2677     else if (fEMCalJet->GetLeadingPt()>=fEMCalJetThreshold)
2678     {
2679         AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalJet->GetLeadingIndex());
2680
2681         for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2682         {
2683             AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
2684
2685             CMSTotalkTArea+=myJet->Area();
2686             if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2687             {
2688                 CMSParticleArea+=myJet->Area();
2689             }
2690             if (IsJetOverlap(myJet,myJet1,kTRUE)==kFALSE)
2691             {
2692                 pTArray[k]=myJet->Pt();
2693                 RhoArray[k]=myJet->Pt()/myJet->Area();
2694                 k++;
2695             }
2696         }
2697         if (k>0)
2698         {
2699             kTRho=MedianRhokT(pTArray,RhoArray,k);
2700         }
2701         else
2702         {
2703             kTRho=0.0;
2704         }
2705     }
2706     else
2707     {
2708         for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2709         {
2710             AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
2711
2712             CMSTotalkTArea+=myJet->Area();
2713             if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2714             {
2715                 CMSParticleArea+=myJet->Area();
2716             }
2717             pTArray[k]=myJet->Pt();
2718             RhoArray[k]=myJet->Pt()/myJet->Area();
2719             k++;
2720         }
2721         if (k>0)
2722         {
2723             kTRho=MedianRhokT(pTArray,RhoArray,k);
2724         }
2725         else
2726         {
2727             kTRho=0.0;
2728         }
2729     }
2730     // Scale CMS Rho by Correction factor
2731     if (CMSTotalkTArea==0.0)
2732     {
2733         CMSCorrectionFactor = 1.0;
2734     }
2735     else
2736     {
2737         //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
2738         CMSCorrectionFactor = CMSParticleArea/((fEMCalPhiTotal-2*fJetR)*(fEMCalEtaTotal-2*fJetR));  // The total physical area should be reduced by the eta & phi cuts due to looping over only fully contained kT jets within the EMCal
2739     }
2740     kTRho*=CMSCorrectionFactor;
2741
2742     fRhoFullCMS->FillRho(fEventCentrality,kTRho);
2743     fRhoFullCMS->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2744     fRhoFullCMS->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
2745     fRhoFullCMS->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
2746     fRhoFullCMS->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucNColl,1);
2747     fRhoFullCMS->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
2748     fRhoFullCMS->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
2749     delete [] RhoArray;
2750     delete [] pTArray;
2751 }
2752
2753 void AliAnalysisTaskFullpAJets::JetPtChargedProfile()
2754 {
2755     Int_t i,j;
2756     Double_t delta_R;
2757     Double_t ED_pT[fEDProfileRBins];
2758
2759     for (i=0;i<fTPCFullJet->GetTotalSignalJets();i++)
2760     {
2761         AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCFullJet->GetSignalJetIndex(i));
2762         if (InsideRect(myJet->Phi(),fTPCPhiMin,fTPCPhiMax,myJet->Eta(),fTPCEtaMin+fEDProfileRUp,fTPCEtaMax-fEDProfileRUp)==kTRUE)
2763         {
2764             for (j=0;j<fEDProfileRBins;j++)
2765             {
2766                 ED_pT[j]=0;
2767             }
2768             TLorentzVector *jet_vec= new TLorentzVector;
2769             myJet->GetMom(*jet_vec);
2770             // Sum all tracks in concentric rings around jet vertex
2771             for (j=0;j<fnTracks;j++)
2772             {
2773                 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2774                 TLorentzVector *track_vec = new TLorentzVector;
2775                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2776                 delta_R=jet_vec->DeltaR(*track_vec);
2777                 if (delta_R<=fEDProfileRUp)
2778                 {
2779                     ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vtrack->Pt();
2780                 }
2781                 delete track_vec;
2782             }
2783
2784             for (j=0;j<fEDProfileRBins;j++)
2785             {
2786                 ED_pT[j]/=TMath::Pi()*TMath::Power((fEDProfileRUp/fEDProfileRBins),2)*(2*j+1);
2787                 fpChargedJetEDProfile[TMath::FloorNint(fEventCentrality/10.)]->Fill(myJet->Pt(),myJet->Eta(),(fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2788                 if (fEventCentrality<=20)
2789                 {
2790                     fpChargedJetRProfile[4+TMath::FloorNint(myJet->Eta()*10.)]->Fill((fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2791                 }
2792             }
2793             delete jet_vec;
2794         }
2795     }
2796 }
2797
2798 void AliAnalysisTaskFullpAJets::JetPtFullProfile()
2799 {
2800     Int_t i,j;
2801     Double_t delta_R;
2802     Double_t ED_pT[fEDProfileRBins];
2803
2804     for (i=0;i<fEMCalFullJet->GetTotalSignalJets();i++)
2805     {
2806         AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalFullJet->GetSignalJetIndex(i));
2807         if (InsideRect(myJet->Phi(),fEMCalPhiMin+fEDProfileRUp,fEMCalPhiMax-fEDProfileRUp,myJet->Eta(),fEMCalEtaMin+fEDProfileRUp,fEMCalEtaMax-fEDProfileRUp)==kTRUE)
2808         {
2809             for (j=0;j<fEDProfileRBins;j++)
2810             {
2811                 ED_pT[j]=0;
2812             }
2813             TLorentzVector *jet_vec= new TLorentzVector;
2814             myJet->GetMom(*jet_vec);
2815             // Sum all tracks in concentric rings around jet vertex
2816             for (j=0;j<fnTracks;j++)
2817             {
2818                 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2819                 TLorentzVector *track_vec = new TLorentzVector;
2820                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2821                 delta_R=jet_vec->DeltaR(*track_vec);
2822                 if (delta_R<=fEDProfileRUp)
2823                 {
2824                     ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vtrack->Pt();
2825                 }
2826                 delete track_vec;
2827             }
2828
2829             // Sum all clusters in concentric rings around jet vertex
2830             for (j=0;j<fnClusters;j++)
2831             {
2832                 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
2833                 TLorentzVector *cluster_vec = new TLorentzVector;
2834                 vcluster->GetMomentum(*cluster_vec,fvertex);
2835                 delta_R=jet_vec->DeltaR(*cluster_vec);
2836                 if (delta_R<=fEDProfileRUp)
2837                 {
2838                     ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vcluster->E();
2839                 }
2840                 delete cluster_vec;
2841             }
2842
2843             for (j=0;j<fEDProfileRBins;j++)
2844             {
2845                 ED_pT[j]/=TMath::Pi()*TMath::Power((fEDProfileRUp/fEDProfileRBins),2)*(2*j+1);
2846                 fpJetEDProfile[TMath::FloorNint(fEventCentrality/10.)]->Fill(myJet->Pt(),myJet->Eta(),(fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2847                 // Fill profile if a "most" central event (0-20%)
2848                 if (fEventCentrality<=20)
2849                 {
2850                     fpJetRProfile[2+TMath::FloorNint(myJet->Eta()*10.)]->Fill((fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2851                 }
2852             }
2853             delete jet_vec;
2854
2855             // Fill constituent histogram
2856             for (j=0;j<myJet->GetNumberOfTracks();j++)
2857             {
2858                 AliVTrack* vtrack = (AliVTrack*) fOrgTracks->At(myJet->TrackAt(j));
2859                 fhJetConstituentPt->Fill(myJet->Pt(),vtrack->Pt());
2860             }
2861
2862             for (j=0;j<myJet->GetNumberOfClusters();j++)
2863             {
2864                 AliVCluster* vcluster = (AliVCluster*) fOrgClusters->At(myJet->ClusterAt(j));
2865                 TLorentzVector *cluster_vec = new TLorentzVector;
2866                 vcluster->GetMomentum(*cluster_vec,fvertex);
2867                 fhJetConstituentPt->Fill(myJet->Pt(),cluster_vec->Pt());
2868                 delete cluster_vec;
2869             }
2870         }
2871     }
2872 }
2873
2874 void AliAnalysisTaskFullpAJets::JetPtEtaProfile()
2875 {
2876     Int_t i,j;
2877     Double_t eta;
2878     Double_t delta_eta;
2879     Double_t Eta_pT[fEtaProfileBins];
2880     Double_t Eta_abs_pT[Int_t(0.5*fEtaProfileBins)];
2881
2882     for (i=0;i<fEMCalFullJet->GetTotalSignalJets();i++)
2883     {
2884         AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalFullJet->GetSignalJetIndex(i));
2885         if (IsInEMCal(myJet->Phi(),myJet->Eta())==kTRUE)
2886         {
2887             for (j=0;j<fEtaProfileBins;j++)
2888             {
2889                 Eta_pT[j]=0;
2890                 Eta_abs_pT[j]=0;
2891             }
2892
2893             // Sum all tracks in strips of eta away from the jet vertex
2894             for (j=0;j<fnTracks;j++)
2895             {
2896                 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2897                 eta=vtrack->Eta();
2898                 delta_eta=TMath::Abs(vtrack->Eta()-myJet->Eta());
2899                 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2900                 {
2901                     Eta_pT[Int_t(0.5*fEtaProfileBins)+TMath::FloorNint(10*eta)]+=vtrack->Pt();
2902                     Eta_abs_pT[TMath::FloorNint(10*delta_eta)]+=vtrack->Pt();
2903                 }
2904             }
2905
2906             // Sum all clusters in strips of eta away from the jet vertex
2907             for (j=0;j<fnClusters;j++)
2908             {
2909                 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
2910                 TLorentzVector *cluster_vec = new TLorentzVector;
2911                 vcluster->GetMomentum(*cluster_vec,fvertex);
2912                 eta=cluster_vec->Eta();
2913                 delta_eta=TMath::Abs(cluster_vec->Eta()-myJet->Eta());
2914                 Eta_pT[Int_t(0.5*fEtaProfileBins)+TMath::FloorNint(10*eta)]+=vcluster->E();
2915                 Eta_abs_pT[TMath::FloorNint(10*delta_eta)]+=vcluster->E();
2916                 delete cluster_vec;
2917             }
2918
2919             for (j=0;j<fEtaProfileBins;j++)
2920             {
2921                 Eta_pT[j]/=0.1*fEMCalPhiTotal;
2922                 // Fill profile if a "most" central event (0-20%)
2923                 if (j<(10*(fEMCalEtaMax-TMath::Abs(myJet->Eta()))))
2924                 {
2925                     Eta_abs_pT[j]/=0.2*fEMCalPhiTotal;
2926                 }
2927                 else
2928                 {
2929                     Eta_abs_pT[j]/=0.1*fEMCalPhiTotal;
2930                 }
2931                 // Fill profile if a "most" central event (0-20%)
2932                 if (fEventCentrality<=20)
2933                 {
2934                     fpJetAbsEtaProfile[7+TMath::FloorNint(myJet->Eta()*10.)]->Fill(0.1*j,Eta_abs_pT[j]);
2935                     fpJetEtaProfile[7+TMath::FloorNint(myJet->Eta()*10.)]->Fill(0.1*(j-7),Eta_pT[j]);
2936                 }
2937             }
2938         }
2939     }
2940 }
2941
2942 void AliAnalysisTaskFullpAJets::DeleteJetData(Bool_t EMCalOn)
2943 {
2944     delete fmyTracks;
2945     delete fTPCJet;
2946     delete fTPCFullJet;
2947     delete fTPCOnlyJet;
2948     delete fTPCkTFullJet;
2949     if (EMCalOn==kTRUE)
2950     {
2951         delete fmyClusters;
2952         delete fEMCalJet;
2953         delete fEMCalFullJet;
2954         delete fEMCalPartJet;
2955         delete fEMCalkTFullJet;
2956     }
2957 }
2958
2959 /////////////////////////////////////////////////////////////////////////////////////////
2960 /////////////////     User Defined Functions      ///////////////////////////////////////
2961 /////////////////////////////////////////////////////////////////////////////////////////
2962
2963 Bool_t AliAnalysisTaskFullpAJets::IsDiJetEvent()
2964 {
2965     // Determine if event contains a di-jet within the detector. Uses charged jets.
2966     // Requires the delta phi of the jets to be 180 +/- 15 degrees.
2967     // Requires both jets to be outside of the EMCal
2968     // Requires both jets to be signal jets
2969
2970     const Double_t dijet_delta_phi=(180/360.)*2*TMath::Pi();
2971     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
2972     Double_t dummy_phi=0.0;
2973
2974     if (fTPCOnlyJet->GetTotalSignalJets()>1)
2975     {
2976         AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCOnlyJet->GetLeadingIndex());
2977         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCOnlyJet->GetSubLeadingIndex());
2978         dummy_phi=TMath::Min(TMath::Abs(myhJet->Phi()-myJet->Phi()),2*TMath::Pi()-TMath::Abs(myhJet->Phi()-myJet->Phi()));
2979         if (dummy_phi>(dijet_delta_phi-dijet_phi_acceptance))
2980         {
2981             fnDiJetEvents++;
2982             return kTRUE;
2983         }
2984     }
2985     return kFALSE;
2986 }
2987
2988 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)
2989 {
2990     if (phi>phi_min && phi<phi_max)
2991     {
2992         if (eta>eta_min && eta<eta_max)
2993         {
2994             return kTRUE;
2995         }
2996     }
2997     return kFALSE;
2998 }
2999
3000 Bool_t AliAnalysisTaskFullpAJets::IsInEMCal(Double_t phi,Double_t eta)
3001 {
3002     return InsideRect(phi,fEMCalPhiMin,fEMCalPhiMax,eta,fEMCalEtaMin,fEMCalEtaMax);
3003 }
3004
3005 Bool_t AliAnalysisTaskFullpAJets::IsInEMCalFull(Double_t r,Double_t phi,Double_t eta)
3006 {
3007     return InsideRect(phi,fEMCalPhiMin+r,fEMCalPhiMax-r,eta,fEMCalEtaMin+r,fEMCalEtaMax-r);
3008 }
3009
3010 Bool_t AliAnalysisTaskFullpAJets::IsInEMCalPart(Double_t r,Double_t phi,Double_t eta)
3011 {
3012     return InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
3013 }
3014
3015 Bool_t AliAnalysisTaskFullpAJets::IsInTPCFull(Double_t r,Double_t phi,Double_t eta)
3016 {
3017     Bool_t in_EMCal= InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
3018     Bool_t in_TPC= InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
3019
3020     if (in_EMCal==kFALSE && in_TPC==kTRUE)
3021     {
3022         return kTRUE;
3023     }
3024     return kFALSE;
3025 }
3026
3027 Bool_t AliAnalysisTaskFullpAJets::IsInTPC(Double_t r,Double_t phi,Double_t eta,Bool_t Complete)
3028 {
3029     if (Complete==kTRUE)
3030     {
3031         return InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
3032     }
3033     return InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin,fTPCEtaMax);
3034 }
3035
3036 Bool_t AliAnalysisTaskFullpAJets::IsJetOverlap(AliEmcalJet *jet1,AliEmcalJet *jet2,Bool_t EMCalOn)
3037 {
3038     Int_t i,j;
3039     Int_t jetTrack1=0;
3040     Int_t jetTrack2=0;
3041     Int_t jetCluster1=0;
3042     Int_t jetCluster2=0;
3043
3044     for (i=0;i<jet1->GetNumberOfTracks();i++)
3045     {
3046         jetTrack1=jet1->TrackAt(i);
3047         for (j=0;j<jet2->GetNumberOfTracks();j++)
3048         {
3049             jetTrack2=jet2->TrackAt(j);
3050             if (jetTrack1 == jetTrack2)
3051             {
3052                 return kTRUE;
3053             }
3054         }
3055     }
3056     if (EMCalOn == kTRUE)
3057     {
3058         for (i=0;i<jet1->GetNumberOfClusters();i++)
3059         {
3060             jetCluster1=jet1->ClusterAt(i);
3061             for (j=0;j<jet2->GetNumberOfClusters();j++)
3062             {
3063                 jetCluster2=jet2->ClusterAt(j);
3064                 if (jetCluster1 == jetCluster2)
3065                 {
3066                     return kTRUE;
3067                 }
3068             }
3069         }
3070     }
3071     return kFALSE;
3072 }
3073
3074 Double_t AliAnalysisTaskFullpAJets::AreaWithinTPC(Double_t r,Double_t eta)
3075 {
3076     Double_t z;
3077     if (eta<(fTPCEtaMin+r))
3078     {
3079         z=eta-fTPCEtaMin;
3080     }
3081     else if(eta>(fTPCEtaMax-r))
3082     {
3083         z=fTPCEtaMax-eta;
3084     }
3085     else
3086     {
3087         z=r;
3088     }
3089     return r*r*TMath::Pi()-AreaEdge(r,z);
3090 }
3091
3092 Double_t AliAnalysisTaskFullpAJets::AreaWithinEMCal(Double_t r,Double_t phi,Double_t eta)
3093 {
3094     Double_t x,y;
3095
3096     if (phi<(fEMCalPhiMin-r) || phi>(fEMCalPhiMax+r))
3097     {
3098         x=-r;
3099     }
3100     else if (phi<(fEMCalPhiMin+r))
3101     {
3102         x=phi-fEMCalPhiMin;
3103     }
3104     else if (phi>(fEMCalPhiMin+r) && phi<(fEMCalPhiMax-r))
3105     {
3106         x=r;
3107     }
3108     else
3109     {
3110         x=fEMCalPhiMax-phi;
3111     }
3112
3113     if (eta<(fEMCalEtaMin-r) || eta>(fEMCalEtaMax+r))
3114     {
3115         y=-r;
3116     }
3117     else if (eta<(fEMCalEtaMin+r))
3118     {
3119         y=eta-fEMCalEtaMin;
3120     }
3121     else if (eta>(fEMCalEtaMin+r) && eta<(fEMCalEtaMax-r))
3122     {
3123         y=r;
3124     }
3125     else
3126     {
3127         y=fEMCalEtaMax-eta;
3128     }
3129
3130     if (x>=0 && y>=0)
3131     {
3132         if (TMath::Sqrt(x*x+y*y)>=r)
3133         {
3134             return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
3135         }
3136         return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y)+AreaOverlap(r,x,y);
3137     }
3138     else if ((x>=r && y<0) || (y>=r && x<0))
3139     {
3140         return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
3141     }
3142     else if (x>0 && x<r && y<0)
3143     {
3144         Double_t a=TMath::Sqrt(r*r-x*x);
3145         Double_t b=TMath::Sqrt(r*r-y*y);
3146         if ((x-b)>0)
3147         {
3148             return r*r*TMath::ASin(b/r)+y*b;
3149         }
3150         else
3151         {
3152             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);
3153         }
3154     }
3155     else if (y>0 && y<r && x<0)
3156     {
3157         Double_t a=TMath::Sqrt(r*r-x*x);
3158         Double_t b=TMath::Sqrt(r*r-y*y);
3159         if ((y-a)>0)
3160         {
3161             return r*r*TMath::ASin(a/r)+x*a;
3162         }
3163         else
3164         {
3165             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);
3166         }
3167     }
3168     else
3169     {
3170         Double_t a=TMath::Sqrt(r*r-x*x);
3171         Double_t b=TMath::Sqrt(r*r-y*y);
3172         if ((x+b)<0)
3173         {
3174             return 0;
3175         }
3176         else
3177         {
3178             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);
3179         }
3180     }
3181 }
3182
3183 Double_t AliAnalysisTaskFullpAJets::AreaEdge(Double_t r,Double_t z)
3184 {
3185     Double_t a=TMath::Sqrt(r*r-z*z);
3186     return r*r*TMath::ASin(a/r)-a*z;
3187 }
3188
3189 Double_t AliAnalysisTaskFullpAJets::AreaOverlap(Double_t r,Double_t x,Double_t y)
3190 {
3191     Double_t a=TMath::Sqrt(r*r-x*x);
3192     Double_t b=TMath::Sqrt(r*r-y*y);
3193     return x*y-0.5*(x*a+y*b)+0.5*r*r*(TMath::ASin(b/r)-TMath::ASin(x/r));
3194 }
3195
3196 Double_t AliAnalysisTaskFullpAJets::TransverseArea(Double_t r,Double_t psi0,Double_t phi,Double_t eta)
3197 {
3198     Double_t area_left=0;
3199     Double_t area_right=0;
3200     Double_t eta_a=0;
3201     Double_t eta_b=0;
3202     Double_t eta_up=0;
3203     Double_t eta_down=0;
3204
3205     Double_t u=eta-fEMCalEtaMin;
3206     Double_t v=fEMCalEtaMax-eta;
3207
3208     Double_t phi1=phi+u*TMath::Tan(psi0);
3209     Double_t phi2=phi-u*TMath::Tan(psi0);
3210     Double_t phi3=phi+v*TMath::Tan(psi0);
3211     Double_t phi4=phi-v*TMath::Tan(psi0);
3212
3213     //Calculate the Left side area
3214     if (phi1>=fEMCalPhiMax)
3215     {
3216         eta_a=eta-u*((fEMCalPhiMax-phi)/(phi1-phi));
3217     }
3218     if (phi2<=fEMCalPhiMin)
3219     {
3220         eta_b=eta-u*((phi-fEMCalPhiMin)/(phi-phi2));
3221     }
3222
3223     if ((phi1>=fEMCalPhiMax) && (phi2<=fEMCalPhiMin))
3224     {
3225         eta_up=TMath::Max(eta_a,eta_b);
3226         eta_down=TMath::Min(eta_a,eta_b);
3227
3228         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);
3229     }
3230     else if (phi1>=fEMCalPhiMax)
3231     {
3232         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);
3233     }
3234     else if (phi2<=fEMCalPhiMin)
3235     {
3236         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);
3237     }
3238     else
3239     {
3240         area_left=0.5*(phi1-phi2+2*r*TMath::Tan(psi0))*(u-r);
3241     }
3242
3243     // Calculate the Right side area
3244     if (phi3>=fEMCalPhiMax)
3245     {
3246         eta_a=eta+v*((fEMCalPhiMax-phi)/(phi3-phi));
3247     }
3248     if (phi4<=fEMCalPhiMin)
3249     {
3250         eta_b=eta+v*((phi-fEMCalPhiMin)/(phi-phi4));
3251     }
3252
3253     if ((phi3>=fEMCalPhiMax) && (phi4<=fEMCalPhiMin))
3254     {
3255         eta_up=TMath::Max(eta_a,eta_b);
3256         eta_down=TMath::Min(eta_a,eta_b);
3257
3258         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);
3259     }
3260     else if (phi3>=fEMCalPhiMax)
3261     {
3262         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);
3263     }
3264     else if (phi4<=fEMCalPhiMin)
3265     {
3266         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);
3267     }
3268     else
3269     {
3270         area_right=0.5*(phi3-phi4+2*r*TMath::Tan(psi0))*(v-r);
3271     }
3272     return area_left+area_right;
3273 }
3274
3275 Double_t AliAnalysisTaskFullpAJets::MedianRhokT(Double_t *pTkTEntries, Double_t *RhokTEntries, Int_t nEntries)
3276 {
3277     // This function is used to calculate the median Rho kT value. The procedure is:
3278     // - Order the kT cluster array from highest rho value to lowest
3279     // - Exclude highest rho kT cluster
3280     // - Return the median rho value of the remaining subset
3281
3282     // Sort Array
3283     const Double_t rho_min=-9.9999E+99;
3284     Int_t j,k;
3285     Double_t w[nEntries];  // Used for sorting
3286     Double_t smax=rho_min;
3287     Int_t sindex=-1;
3288
3289     Double_t pTmax=0.0;
3290     Int_t pTmaxID=-1;
3291
3292     for (j=0;j<nEntries;j++)
3293     {
3294         w[j]=0.0;
3295     }
3296
3297     for (j=0;j<nEntries;j++)
3298     {
3299         if (pTkTEntries[j]>pTmax)
3300         {
3301             pTmax=pTkTEntries[j];
3302             pTmaxID=j;
3303         }
3304     }
3305
3306     for (j=0;j<nEntries;j++)
3307     {
3308         for (k=0;k<nEntries;k++)
3309         {
3310             if (RhokTEntries[k]>smax)
3311             {
3312                 smax=RhokTEntries[k];
3313                 sindex=k;
3314             }
3315         }
3316         w[j]=smax;
3317         RhokTEntries[sindex]=rho_min;
3318         smax=rho_min;
3319         sindex=-1;
3320     }
3321     return w[nEntries/2];
3322 }
3323
3324
3325 // AlipAJetData Class Member Defs
3326 // Constructors
3327 AliAnalysisTaskFullpAJets::AlipAJetData::AlipAJetData() :
3328
3329     fName(0),
3330     fIsJetsFull(0),
3331     fnTotal(0),
3332     fnJets(0),
3333     fnJetsSC(0),
3334     fJetR(0),
3335     fSignalPt(0),
3336     fAreaCutFrac(0.6),
3337     fPtMaxIndex(0),
3338     fPtMax(0),
3339     fPtSubLeadingIndex(0),
3340     fPtSubLeading(0),
3341     fJetsIndex(0),
3342     fJetsSCIndex(0),
3343     fIsJetInArray(0)
3344 {
3345     fnTotal=0;
3346     // Dummy constructor ALWAYS needed for I/O.
3347 }
3348
3349 AliAnalysisTaskFullpAJets::AlipAJetData::AlipAJetData(const char *name, Bool_t isFull, Int_t nEntries) :
3350
3351     fName(0),
3352     fIsJetsFull(0),
3353     fnTotal(0),
3354     fnJets(0),
3355     fnJetsSC(0),
3356     fJetR(0),
3357     fSignalPt(0),
3358     fAreaCutFrac(0.6),
3359     fPtMaxIndex(0),
3360     fPtMax(0),
3361     fPtSubLeadingIndex(0),
3362     fPtSubLeading(0),
3363     fJetsIndex(0),
3364     fJetsSCIndex(0),
3365     fIsJetInArray(0)
3366 {
3367     SetName(name);
3368     SetIsJetsFull(isFull);
3369     SetTotalEntries(nEntries);
3370     SetLeading(0,-9.99E+099);
3371     SetSubLeading(0,-9.99E+099);
3372     SetSignalCut(0);
3373     SetAreaCutFraction(0.6);
3374     SetJetR(fJetR);
3375 }
3376
3377 // Destructor
3378 AliAnalysisTaskFullpAJets::AlipAJetData::~AlipAJetData()
3379 {
3380     if (fnTotal!=0)
3381     {
3382         SetName("");
3383         SetIsJetsFull(kFALSE);
3384         SetTotalEntries(0);
3385         SetTotalJets(0);
3386         SetTotalSignalJets(0);
3387         SetLeading(0,0);
3388         SetSubLeading(0,0);
3389         SetSignalCut(0);
3390         SetAreaCutFraction(0);
3391         SetJetR(0);
3392
3393         delete [] fJetsIndex;
3394         delete [] fJetsSCIndex;
3395         delete [] fIsJetInArray;
3396     }
3397 }
3398
3399 // User Defined Sub-Routines
3400 void AliAnalysisTaskFullpAJets::AlipAJetData::InitializeJetData(TClonesArray *jetList, Int_t nEntries)
3401 {
3402     Int_t i=0;
3403     Int_t k=0;
3404     Int_t l=0;
3405     Double_t AreaThreshold = fAreaCutFrac*TMath::Pi()*TMath::Power(fJetR,2);
3406
3407     // Initialize Jet Data
3408     for (i=0;i<nEntries;i++)
3409     {
3410         AliEmcalJet *myJet =(AliEmcalJet*) jetList->At(i);
3411
3412         if (fIsJetInArray[i]==kTRUE && myJet->Area()>AreaThreshold)
3413         {
3414             SetJetIndex(i,k);
3415             if (myJet->Pt()>fPtMax)
3416             {
3417                 SetSubLeading(fPtMaxIndex,fPtMax);
3418                 SetLeading(i,myJet->Pt());
3419             }
3420             else if (myJet->Pt()>fPtSubLeading)
3421             {
3422                 SetSubLeading(i,myJet->Pt());
3423             }
3424             if (myJet->Pt()>=fSignalPt)
3425             {
3426                 SetSignalJetIndex(i,l);
3427                 l++;
3428             }
3429             k++;
3430         }
3431     }
3432     SetTotalJets(k);
3433     SetTotalSignalJets(l);
3434 }
3435
3436 // Setters
3437 void AliAnalysisTaskFullpAJets::AlipAJetData::SetName(const char *name)
3438 {
3439     fName = name;
3440 }
3441
3442 void AliAnalysisTaskFullpAJets::AlipAJetData::SetIsJetsFull(Bool_t isFull)
3443 {
3444     fIsJetsFull = isFull;
3445 }
3446
3447 void AliAnalysisTaskFullpAJets::AlipAJetData::SetTotalEntries(Int_t nEntries)
3448 {
3449     fnTotal = nEntries;
3450     fJetsIndex = new Int_t[fnTotal];
3451     fJetsSCIndex = new Int_t[fnTotal];
3452     fIsJetInArray = new Bool_t[fnTotal];
3453 }
3454
3455 void AliAnalysisTaskFullpAJets::AlipAJetData::SetTotalJets(Int_t nJets)
3456 {
3457     fnJets = nJets;
3458 }
3459
3460 void AliAnalysisTaskFullpAJets::AlipAJetData::SetTotalSignalJets(Int_t nSignalJets)
3461 {
3462     fnJetsSC = nSignalJets;
3463 }
3464
3465 void AliAnalysisTaskFullpAJets::AlipAJetData::SetSignalCut(Double_t Pt)
3466 {
3467     fSignalPt = Pt;
3468 }
3469
3470 void AliAnalysisTaskFullpAJets::AlipAJetData::SetLeading(Int_t index, Double_t Pt)
3471 {
3472     fPtMaxIndex = index;
3473     fPtMax = Pt;
3474 }
3475
3476 void AliAnalysisTaskFullpAJets::AlipAJetData::SetSubLeading(Int_t index, Double_t Pt)
3477 {
3478     fPtSubLeadingIndex = index;
3479     fPtSubLeading = Pt;
3480 }
3481
3482 void AliAnalysisTaskFullpAJets::AlipAJetData::SetJetIndex(Int_t index, Int_t At)
3483 {
3484     fJetsIndex[At] = index;
3485 }
3486
3487 void AliAnalysisTaskFullpAJets::AlipAJetData::SetSignalJetIndex(Int_t index, Int_t At)
3488 {
3489     fJetsSCIndex[At] = index;
3490 }
3491
3492 void AliAnalysisTaskFullpAJets::AlipAJetData::SetIsJetInArray(Bool_t isInArray, Int_t At)
3493 {
3494     fIsJetInArray[At] = isInArray;
3495 }
3496
3497 void AliAnalysisTaskFullpAJets::AlipAJetData::SetAreaCutFraction(Double_t areaFraction)
3498 {
3499     fAreaCutFrac = areaFraction;
3500 }
3501
3502 void AliAnalysisTaskFullpAJets::AlipAJetData::SetJetR(Double_t jetR)
3503 {
3504     fJetR = jetR;
3505 }
3506
3507 // Getters
3508 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetTotalEntries()
3509 {
3510     return fnTotal;
3511 }
3512
3513 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetTotalJets()
3514 {
3515     return fnJets;
3516 }
3517
3518 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetTotalSignalJets()
3519 {
3520     return fnJetsSC;
3521 }
3522
3523 Double_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSignalCut()
3524 {
3525     return fSignalPt;
3526 }
3527
3528 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetLeadingIndex()
3529 {
3530     return fPtMaxIndex;
3531 }
3532
3533 Double_t AliAnalysisTaskFullpAJets::AlipAJetData::GetLeadingPt()
3534 {
3535     return fPtMax;
3536 }
3537
3538 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSubLeadingIndex()
3539 {
3540     return fPtSubLeadingIndex;
3541 }
3542
3543 Double_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSubLeadingPt()
3544 {
3545     return fPtSubLeading;
3546 }
3547
3548 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetJetIndex(Int_t At)
3549 {
3550     return fJetsIndex[At];
3551 }
3552
3553 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSignalJetIndex(Int_t At)
3554 {
3555     return fJetsSCIndex[At];
3556 }
3557
3558 Bool_t AliAnalysisTaskFullpAJets::AlipAJetData::GetIsJetInArray(Int_t At)
3559 {
3560     return fIsJetInArray[At];
3561 }
3562
3563
3564 // AlipAJetHistos Class Member Defs
3565 // Constructors
3566 AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos() :
3567
3568     fOutput(0),
3569
3570     fh020Rho(0),
3571     fh80100Rho(0),
3572     fhRho(0),
3573     fhRhoCen(0),
3574     fh020BSPt(0),
3575     fh80100BSPt(0),
3576     fhBSPt(0),
3577     fhBSPtCen(0),
3578     fh020BSPtSignal(0),
3579     fh80100BSPtSignal(0),
3580     fhBSPtSignal(0),
3581     fhBSPtCenSignal(0),
3582     fh020DeltaPt(0),
3583     fh80100DeltaPt(0),
3584     fhDeltaPt(0),
3585     fhDeltaPtCen(0),
3586     fh020DeltaPtSignal(0),
3587     fh80100DeltaPtSignal(0),
3588     fhDeltaPtSignal(0),
3589     fhDeltaPtCenSignal(0),
3590     fh020DeltaPtNColl(0),
3591     fh80100DeltaPtNColl(0),
3592     fhDeltaPtNColl(0),
3593     fhDeltaPtCenNColl(0),
3594     fh020BckgFlucPt(0),
3595     fh80100BckgFlucPt(0),
3596     fhBckgFlucPt(0),
3597     fhBckgFlucPtCen(0),
3598
3599     fpRho(0),
3600     fpLJetRho(0),
3601
3602     fName(0),
3603     fCentralityTag(0),
3604     fCentralityBins(0),
3605     fCentralityLow(0),
3606     fCentralityUp(0),
3607     fPtBins(0),
3608     fPtLow(0),
3609     fPtUp(0),
3610     fRhoPtBins(0),
3611     fRhoPtLow(0),
3612     fRhoPtUp(0),
3613     fDeltaPtBins(0),
3614     fDeltaPtLow(0),
3615     fDeltaPtUp(0),
3616     fBckgFlucPtBins(0),
3617     fBckgFlucPtLow(0),
3618     fBckgFlucPtUp(0),
3619     fLJetPtBins(0),
3620     fLJetPtLow(0),
3621     fLJetPtUp(0),
3622     fRhoValue(0)
3623 {
3624     // Dummy constructor ALWAYS needed for I/O.
3625 }
3626
3627 AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name) :
3628
3629     fOutput(0),
3630
3631     fh020Rho(0),
3632     fh80100Rho(0),
3633     fhRho(0),
3634     fhRhoCen(0),
3635     fh020BSPt(0),
3636     fh80100BSPt(0),
3637     fhBSPt(0),
3638     fhBSPtCen(0),
3639     fh020BSPtSignal(0),
3640     fh80100BSPtSignal(0),
3641     fhBSPtSignal(0),
3642     fhBSPtCenSignal(0),
3643     fh020DeltaPt(0),
3644     fh80100DeltaPt(0),
3645     fhDeltaPt(0),
3646     fhDeltaPtCen(0),
3647     fh020DeltaPtSignal(0),
3648     fh80100DeltaPtSignal(0),
3649     fhDeltaPtSignal(0),
3650     fhDeltaPtCenSignal(0),
3651     fh020DeltaPtNColl(0),
3652     fh80100DeltaPtNColl(0),
3653     fhDeltaPtNColl(0),
3654     fhDeltaPtCenNColl(0),
3655     fh020BckgFlucPt(0),
3656     fh80100BckgFlucPt(0),
3657     fhBckgFlucPt(0),
3658     fhBckgFlucPtCen(0),
3659
3660     fpRho(0),
3661     fpLJetRho(0),
3662
3663     fName(0),
3664     fCentralityTag(0),
3665     fCentralityBins(0),
3666     fCentralityLow(0),
3667     fCentralityUp(0),
3668     fPtBins(0),
3669     fPtLow(0),
3670     fPtUp(0),
3671     fRhoPtBins(0),
3672     fRhoPtLow(0),
3673     fRhoPtUp(0),
3674     fDeltaPtBins(0),
3675     fDeltaPtLow(0),
3676     fDeltaPtUp(0),
3677     fBckgFlucPtBins(0),
3678     fBckgFlucPtLow(0),
3679     fBckgFlucPtUp(0),
3680     fLJetPtBins(0),
3681     fLJetPtLow(0),
3682     fLJetPtUp(0),
3683     fRhoValue(0)
3684 {
3685     SetName(name);
3686     SetCentralityTag("V0A");
3687     SetCentralityRange(100,0,100);
3688     SetPtRange(250,-50,200);
3689     SetRhoPtRange(500,0,50);
3690     SetDeltaPtRange(200,-100,100);
3691     SetBackgroundFluctuationsPtRange(100,0,100);
3692     SetLeadingJetPtRange(200,0,200);
3693
3694     Init();
3695 }
3696
3697 AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, const char *centag) :
3698
3699     fOutput(0),
3700
3701     fh020Rho(0),
3702     fh80100Rho(0),
3703     fhRho(0),
3704     fhRhoCen(0),
3705     fh020BSPt(0),
3706     fh80100BSPt(0),
3707     fhBSPt(0),
3708     fhBSPtCen(0),
3709     fh020BSPtSignal(0),
3710     fh80100BSPtSignal(0),
3711     fhBSPtSignal(0),
3712     fhBSPtCenSignal(0),
3713     fh020DeltaPt(0),
3714     fh80100DeltaPt(0),
3715     fhDeltaPt(0),
3716     fhDeltaPtCen(0),
3717     fh020DeltaPtSignal(0),
3718     fh80100DeltaPtSignal(0),
3719     fhDeltaPtSignal(0),
3720     fhDeltaPtCenSignal(0),
3721     fh020DeltaPtNColl(0),
3722     fh80100DeltaPtNColl(0),
3723     fhDeltaPtNColl(0),
3724     fhDeltaPtCenNColl(0),
3725     fh020BckgFlucPt(0),
3726     fh80100BckgFlucPt(0),
3727     fhBckgFlucPt(0),
3728     fhBckgFlucPtCen(0),
3729
3730     fpRho(0),
3731     fpLJetRho(0),
3732
3733     fName(0),
3734     fCentralityTag(0),
3735     fCentralityBins(0),
3736     fCentralityLow(0),
3737     fCentralityUp(0),
3738     fPtBins(0),
3739     fPtLow(0),
3740     fPtUp(0),
3741     fRhoPtBins(0),
3742     fRhoPtLow(0),
3743     fRhoPtUp(0),
3744     fDeltaPtBins(0),
3745     fDeltaPtLow(0),
3746     fDeltaPtUp(0),
3747     fBckgFlucPtBins(0),
3748     fBckgFlucPtLow(0),
3749     fBckgFlucPtUp(0),
3750     fLJetPtBins(0),
3751     fLJetPtLow(0),
3752     fLJetPtUp(0),
3753     fRhoValue(0)
3754 {
3755     SetName(name);
3756     SetCentralityTag(centag);
3757     SetCentralityRange(100,0,100);
3758     SetPtRange(250,-50,200);
3759     SetRhoPtRange(500,0,50);
3760     SetDeltaPtRange(200,-100,100);
3761     SetBackgroundFluctuationsPtRange(100,0,100);
3762     SetLeadingJetPtRange(200,0,200);
3763
3764     Init();
3765 }
3766
3767 // Destructor
3768 AliAnalysisTaskFullpAJets::AlipAJetHistos::~AlipAJetHistos()
3769 {
3770     if (fOutput)
3771     {
3772         delete fOutput;
3773     }
3774 }
3775
3776 void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
3777 {
3778     fOutput = new TList();
3779     fOutput->SetOwner();
3780     fOutput->SetName(fName);
3781
3782     TString RhoString="";
3783     TString PtString="";
3784     TString DeltaPtString="";
3785     TString BckgFlucPtString="";
3786     TString CentralityString;
3787     CentralityString = Form("Centrality (%s)",fCentralityTag);
3788
3789     // Rho Spectral Plots
3790     RhoString = Form("%d-%d Centrality, Rho Spectrum",0,20);
3791     fh020Rho = new TH1D("fh020Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
3792     fh020Rho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
3793     fh020Rho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
3794     fh020Rho->Sumw2();
3795
3796     RhoString = Form("%d-%d Centrality, Rho Spectrum",80,100);
3797     fh80100Rho = new TH1D("fh80100Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
3798     fh80100Rho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
3799     fh80100Rho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
3800     fh80100Rho->Sumw2();
3801
3802     RhoString = Form("%d-%d Centrality, Rho Spectrum",0,100);
3803     fhRho = new TH1D("fhRho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
3804     fhRho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
3805     fhRho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
3806     fhRho->Sumw2();
3807
3808     RhoString = "Rho Spectrum vs Centrality";
3809     fhRhoCen = new TH2D("fhRhoCen",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3810     fhRhoCen->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
3811     fhRhoCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3812     fhRhoCen->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
3813     fhRhoCen->Sumw2();
3814
3815     // Background Subtracted Plots
3816     PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,20);
3817     fh020BSPt = new TH1D("fh020BSPt",PtString,fPtBins,fPtLow,fPtUp);
3818     fh020BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3819     fh020BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3820     fh020BSPt->Sumw2();
3821
3822     PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",80,100);
3823     fh80100BSPt = new TH1D("fh80100BSPt",PtString,fPtBins,fPtLow,fPtUp);
3824     fh80100BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3825     fh80100BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3826     fh80100BSPt->Sumw2();
3827
3828     PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,100);
3829     fhBSPt = new TH1D("fhBSPt",PtString,fPtBins,fPtLow,fPtUp);
3830     fhBSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3831     fhBSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3832     fhBSPt->Sumw2();
3833
3834     PtString = "Background Subtracted Jet Spectrum vs Centrality";
3835     fhBSPtCen = new TH2D("fhBSPtCen",PtString,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3836     fhBSPtCen->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3837     fhBSPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3838     fhBSPtCen->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3839     fhBSPtCen->Sumw2();
3840
3841     PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,20);
3842     fh020BSPtSignal = new TH1D("fh020BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
3843     fh020BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3844     fh020BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3845     fh020BSPtSignal->Sumw2();
3846
3847     PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",80,100);
3848     fh80100BSPtSignal = new TH1D("fh80100BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
3849     fh80100BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3850     fh80100BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3851     fh80100BSPtSignal->Sumw2();
3852
3853     PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,100);
3854     fhBSPtSignal = new TH1D("fhBSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
3855     fhBSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3856     fhBSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3857     fhBSPtSignal->Sumw2();
3858
3859     PtString = "Background Subtracted Signal Jet Spectrum vs Centrality";
3860     fhBSPtCenSignal = new TH2D("fhBSPtCenSignal",PtString,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3861     fhBSPtCenSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3862     fhBSPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3863     fhBSPtCenSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3864     fhBSPtCenSignal->Sumw2();
3865
3866     // Delta Pt Plots with RC at least 2R away from Leading Signal
3867     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
3868     fh020DeltaPt = new TH1D("fh020DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3869     fh020DeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3870     fh020DeltaPt->GetYaxis()->SetTitle("Probability Density");
3871     fh020DeltaPt->Sumw2();
3872
3873     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
3874     fh80100DeltaPt = new TH1D("fh80100DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3875     fh80100DeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3876     fh80100DeltaPt->GetYaxis()->SetTitle("Probability Density");
3877     fh80100DeltaPt->Sumw2();
3878
3879     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
3880     fhDeltaPt = new TH1D("fhDeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3881     fhDeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3882     fhDeltaPt->GetYaxis()->SetTitle("Probability Density");
3883     fhDeltaPt->Sumw2();
3884
3885     DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
3886     fhDeltaPtCen = new TH2D("fhDeltaPtCen",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3887     fhDeltaPtCen->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3888     fhDeltaPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3889     fhDeltaPtCen->GetZaxis()->SetTitle("Probability Density");
3890     fhDeltaPtCen->Sumw2();
3891
3892     // Delta Pt Plots with no spatial restrictions on RC
3893     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
3894     fh020DeltaPtSignal = new TH1D("fh020DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3895     fh020DeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3896     fh020DeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
3897     fh020DeltaPtSignal->Sumw2();
3898
3899     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
3900     fh80100DeltaPtSignal = new TH1D("fh80100DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3901     fh80100DeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3902     fh80100DeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
3903     fh80100DeltaPtSignal->Sumw2();
3904
3905     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
3906     fhDeltaPtSignal = new TH1D("fhDeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3907     fhDeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3908     fhDeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
3909     fhDeltaPtSignal->Sumw2();
3910
3911     DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
3912     fhDeltaPtCenSignal = new TH2D("fhDeltaPtCenSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3913     fhDeltaPtCenSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3914     fhDeltaPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3915     fhDeltaPtCenSignal->GetZaxis()->SetTitle("Probability Density");
3916     fhDeltaPtCenSignal->Sumw2();
3917
3918     // Delta Pt Plots with NColl restrictions on RC
3919     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
3920     fh020DeltaPtNColl = new TH1D("fh020DeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3921     fh020DeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3922     fh020DeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
3923     fh020DeltaPtNColl->Sumw2();
3924
3925     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
3926     fh80100DeltaPtNColl = new TH1D("fh80100DeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3927     fh80100DeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3928     fh80100DeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
3929     fh80100DeltaPtNColl->Sumw2();
3930
3931     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
3932     fhDeltaPtNColl = new TH1D("fhDeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3933     fhDeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3934     fhDeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
3935     fhDeltaPtNColl->Sumw2();
3936
3937     DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
3938     fhDeltaPtCenNColl = new TH2D("fhDeltaPtCenNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3939     fhDeltaPtCenNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3940     fhDeltaPtCenNColl->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3941     fhDeltaPtCenNColl->GetZaxis()->SetTitle("Probability Density");
3942     fhDeltaPtCenNColl->Sumw2();
3943
3944     // Background Fluctuations Pt Plots
3945     BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,20);
3946     fh020BckgFlucPt = new TH1D("fh020BckgFlucPt",PtString,fPtBins,fPtLow,fPtUp);
3947     fh020BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
3948     fh020BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3949     fh020BckgFlucPt->Sumw2();
3950
3951     BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",80,100);
3952     fh80100BckgFlucPt = new TH1D("fh80100BckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
3953     fh80100BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
3954     fh80100BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3955     fh80100BckgFlucPt->Sumw2();
3956
3957     BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,100);
3958     fhBckgFlucPt = new TH1D("fhBckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
3959     fhBckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
3960     fhBckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3961     fhBckgFlucPt->Sumw2();
3962
3963     BckgFlucPtString = "Background Fluctuation p_{T} Spectrum vs Centrality";
3964     fhBckgFlucPtCen = new TH2D("fhBckgFlucPtCen",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3965     fhBckgFlucPtCen->GetXaxis()->SetTitle("#p_{T} (GeV/c)");
3966     fhBckgFlucPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3967     fhBckgFlucPtCen->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3968     fhBckgFlucPtCen->Sumw2();
3969
3970     // Background Density vs Centrality Profile
3971     RhoString = "Background Density vs Centrality";
3972     fpRho = new TProfile("fpRho",RhoString,fCentralityBins,fCentralityLow,fCentralityUp);
3973     fpRho->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
3974     fpRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
3975
3976     // Background Density vs Leading Jet Profile
3977     fpLJetRho = new TProfile("fpLJetRho","#rho vs Leading Jet p_{T}",fLJetPtBins,fLJetPtLow,fLJetPtUp);
3978     fpLJetRho->GetXaxis()->SetTitle("Leading Jet p_{T}");
3979     fpLJetRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
3980
3981     // Add Histos & Profiles to List
3982     fOutput->Add(fh020Rho);
3983     fOutput->Add(fh80100Rho);
3984     fOutput->Add(fhRho);
3985     fOutput->Add(fhRhoCen);
3986     fOutput->Add(fh020BSPt);
3987     fOutput->Add(fh80100BSPt);
3988     fOutput->Add(fhBSPt);
3989     fOutput->Add(fhBSPtCen);
3990     fOutput->Add(fh020BSPtSignal);
3991     fOutput->Add(fh80100BSPtSignal);
3992     fOutput->Add(fhBSPtSignal);
3993     fOutput->Add(fhBSPtCenSignal);
3994     fOutput->Add(fh020DeltaPt);
3995     fOutput->Add(fh80100DeltaPt);
3996     fOutput->Add(fhDeltaPt);
3997     fOutput->Add(fhDeltaPtCen);
3998     fOutput->Add(fh020DeltaPtSignal);
3999     fOutput->Add(fh80100DeltaPtSignal);
4000     fOutput->Add(fhDeltaPtSignal);
4001     fOutput->Add(fhDeltaPtCenSignal);
4002     fOutput->Add(fh020DeltaPtNColl);
4003     fOutput->Add(fh80100DeltaPtNColl);
4004     fOutput->Add(fhDeltaPtNColl);
4005     fOutput->Add(fhDeltaPtCenNColl);
4006     fOutput->Add(fh020BckgFlucPt);
4007     fOutput->Add(fh80100BckgFlucPt);
4008     fOutput->Add(fhBckgFlucPt);
4009     fOutput->Add(fhBckgFlucPtCen);
4010     fOutput->Add(fpRho);
4011     fOutput->Add(fpLJetRho);
4012 }
4013
4014 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetName(const char *name)
4015 {
4016     fName = name;
4017 }
4018
4019 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetCentralityTag(const char *name)
4020 {
4021     fCentralityTag = name;
4022 }
4023
4024 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetCentralityRange(Int_t bins, Double_t low, Double_t up)
4025 {
4026     fCentralityBins=bins;
4027     fCentralityLow=low;
4028     fCentralityUp=up;
4029 }
4030
4031 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetPtRange(Int_t bins, Double_t low, Double_t up)
4032 {
4033     fPtBins=bins;
4034     fPtLow=low;
4035     fPtUp=up;
4036 }
4037
4038 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetRhoPtRange(Int_t bins, Double_t low, Double_t up)
4039 {
4040     fRhoPtBins=bins;
4041     fRhoPtLow=low;
4042     fRhoPtUp=up;
4043 }
4044
4045 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetDeltaPtRange(Int_t bins, Double_t low, Double_t up)
4046 {
4047     fDeltaPtBins=bins;
4048     fDeltaPtLow=low;
4049     fDeltaPtUp=up;
4050 }
4051
4052 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetBackgroundFluctuationsPtRange(Int_t bins, Double_t low, Double_t up)
4053 {
4054     fBckgFlucPtBins=bins;
4055     fBckgFlucPtLow=low;
4056     fBckgFlucPtUp=up;
4057 }
4058
4059 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetLeadingJetPtRange(Int_t bins, Double_t low, Double_t up)
4060 {
4061     fLJetPtBins=bins;
4062     fLJetPtLow=low;
4063     fLJetPtUp=up;
4064 }
4065
4066 TList* AliAnalysisTaskFullpAJets::AlipAJetHistos::GetOutputHistos()
4067 {
4068     return fOutput;
4069 }
4070
4071 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillRho(Double_t eventCentrality, Double_t rho)
4072 {
4073     fRhoValue = rho;
4074
4075     fhRho->Fill(rho);
4076     fhRhoCen->Fill(rho,eventCentrality);
4077     fpRho->Fill(eventCentrality,rho);
4078
4079     if (eventCentrality<=20)
4080     {
4081         fh020Rho->Fill(rho);
4082     }
4083     else if (eventCentrality>=80)
4084     {
4085         fh80100Rho->Fill(rho);
4086     }
4087 }
4088
4089 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillBSJS(Double_t eventCentrality, Double_t rho, Double_t signalCut, TClonesArray *jetList, Int_t *indexJetList, Int_t nIndexJetList)
4090 {
4091     Int_t i;
4092     Double_t tempPt=0.0;
4093
4094     for (i=0;i<nIndexJetList;i++)
4095     {
4096         AliEmcalJet *myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
4097         tempPt=myJet->Pt()-rho*myJet->Area();
4098
4099         fhBSPt->Fill(tempPt);
4100         fhBSPtCen->Fill(tempPt,eventCentrality);
4101         if (eventCentrality<=20)
4102         {
4103             fh020BSPt->Fill(tempPt);
4104         }
4105         else if (eventCentrality>=80)
4106         {
4107             fh80100BSPt->Fill(tempPt);
4108         }
4109
4110         if (myJet->Pt()>=signalCut)
4111         {
4112             fhBSPtSignal->Fill(tempPt);
4113             fhBSPtCenSignal->Fill(tempPt,eventCentrality);
4114             if (eventCentrality<=20)
4115             {
4116                 fh020BSPtSignal->Fill(tempPt);
4117             }
4118             else if (eventCentrality>=80)
4119             {
4120                 fh80100BSPtSignal->Fill(tempPt);
4121             }
4122         }
4123         tempPt=0.0;
4124     }
4125 }
4126
4127 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillDeltaPt(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
4128 {
4129     Int_t i;
4130     Double_t tempPt=0.0;
4131
4132     for (i=0;i<nRC;i++)
4133     {
4134         tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
4135         fhDeltaPt->Fill(tempPt);
4136         fhDeltaPtCen->Fill(tempPt,eventCentrality);
4137         if (eventCentrality<=20)
4138         {
4139             fh020DeltaPt->Fill(tempPt);
4140         }
4141         else if (eventCentrality>=80)
4142         {
4143             fh80100DeltaPt->Fill(tempPt);
4144         }
4145         tempPt=0.0;
4146     }
4147 }
4148
4149 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillDeltaPtSignal(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
4150 {
4151     Int_t i;
4152     Double_t tempPt=0.0;
4153
4154     for (i=0;i<nRC;i++)
4155     {
4156         tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
4157         fhDeltaPtSignal->Fill(tempPt);
4158         fhDeltaPtCenSignal->Fill(tempPt,eventCentrality);
4159         if (eventCentrality<=20)
4160         {
4161             fh020DeltaPtSignal->Fill(tempPt);
4162         }
4163         else if (eventCentrality>=80)
4164         {
4165             fh80100DeltaPtSignal->Fill(tempPt);
4166         }
4167         tempPt=0.0;
4168     }
4169 }
4170
4171 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillDeltaPtNColl(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
4172 {
4173     Int_t i;
4174     Double_t tempPt=0.0;
4175
4176     for (i=0;i<nRC;i++)
4177     {
4178         tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
4179         fhDeltaPtNColl->Fill(tempPt);
4180         fhDeltaPtCenNColl->Fill(tempPt,eventCentrality);
4181         if (eventCentrality<=20)
4182         {
4183             fh020DeltaPtNColl->Fill(tempPt);
4184         }
4185         else if (eventCentrality>=80)
4186         {
4187             fh80100DeltaPtNColl->Fill(tempPt);
4188         }
4189         tempPt=0.0;
4190     }
4191 }
4192
4193 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillBackgroundFluctuations(Double_t eventCentrality, Double_t rho, Double_t jetRadius)
4194 {
4195     Double_t tempPt=0.0;
4196
4197     tempPt=rho*TMath::Power(jetRadius,2);
4198     fhBckgFlucPt->Fill(tempPt);
4199     fhBckgFlucPtCen->Fill(tempPt,eventCentrality);
4200     if (eventCentrality<=20)
4201     {
4202         fh020BckgFlucPt->Fill(tempPt);
4203     }
4204     else if (eventCentrality>=80)
4205     {
4206         fh80100BckgFlucPt->Fill(tempPt);
4207     }
4208 }
4209
4210 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillLeadingJetPtRho(Double_t jetPt, Double_t rho)
4211 {
4212     fpLJetRho->Fill(jetPt,rho);
4213 }
4214
4215 Double_t AliAnalysisTaskFullpAJets::AlipAJetHistos::GetRho()
4216 {
4217     return fRhoValue;
4218 }
4219
4220