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