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