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