PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullpAJets.cxx

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