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