PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullpAJets.cxx

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