PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullpAJets.cxx

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