PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullpAJets.cxx

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