1 #include "AliAnalysisTaskFullpAJets.h"
13 #include <TLorentzVector.h>
15 #include <TProfile2D.h>
16 #include <TProfile3D.h>
19 #include <TClonesArray.h>
20 #include <TObjArray.h>
22 #include "AliAnalysisTaskSE.h"
23 #include "AliAnalysisManager.h"
25 #include "AliESDtrackCuts.h"
26 #include "AliESDEvent.h"
27 #include "AliESDInputHandler.h"
28 #include "AliAODEvent.h"
29 #include "AliMCEvent.h"
30 #include "AliVTrack.h"
31 #include "AliVCluster.h"
32 #include "AliEmcalJet.h"
33 #include "AliEMCALGeometry.h"
36 ClassImp(AliAnalysisTaskFullpAJets)
38 //________________________________________________________________________
39 AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
55 fhJetConstituentPt(0),
63 fpClusterPtProfile(0),
80 fRhoChargedkTScale(0),
82 fRhoChargedCMSScale(0),
98 fEMCalPhiMin(1.39626),
99 fEMCalPhiMax(3.26377),
100 fEMCalPhiTotal(1.86750),
107 fTPCPhiTotal(6.28319),
114 fJetAreaCutFrac(0.6),
115 fJetAreaThreshold(0.30159),
127 fEtaProfileLow(-0.7),
132 fEDProfilePtBins(100),
135 fEDProfileEtaBins(4),
136 fEDProfileEtaLow(-0.2),
137 fEDProfileEtaUp(0.2),
146 fEMCalJetThreshold(5),
152 fmyAKTChargedJets(0),
159 fEMCalRCBckgFlucSignal(0),
160 fTPCRCBckgFlucSignal(0)
162 // Dummy constructor ALWAYS needed for I/O.
163 fpJetEtaProfile = new TProfile *[14];
164 fpJetAbsEtaProfile = new TProfile *[14];
165 fpChargedJetRProfile = new TProfile *[8];
166 fpJetRProfile= new TProfile *[4];
167 fpChargedJetEDProfile= new TProfile3D *[10];
168 fpJetEDProfile= new TProfile3D *[10];
169 fvertex[0]=0.0,fvertex[1]=0.0,fvertex[2]=0.0;
172 //________________________________________________________________________
173 AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
174 AliAnalysisTaskSE(name),
184 fhEMCalCellCounts(0),
189 fhJetConstituentPt(0),
197 fpClusterPtProfile(0),
214 fRhoChargedkTScale(0),
216 fRhoChargedCMSScale(0),
232 fEMCalPhiMin(1.39626),
233 fEMCalPhiMax(3.26377),
234 fEMCalPhiTotal(1.86750),
241 fTPCPhiTotal(6.28319),
248 fJetAreaCutFrac(0.6),
249 fJetAreaThreshold(0.30159),
261 fEtaProfileLow(-0.7),
266 fEDProfilePtBins(100),
269 fEDProfileEtaBins(4),
270 fEDProfileEtaLow(-0.2),
271 fEDProfileEtaUp(0.2),
280 fEMCalJetThreshold(5),
286 fmyAKTChargedJets(0),
293 fEMCalRCBckgFlucSignal(0),
294 fTPCRCBckgFlucSignal(0)
297 // Define input and output slots here (never in the dummy constructor)
298 // Input slot #0 works with a TChain - it is connected to the default input container
299 // Output slot #1 writes into a TH1 container
300 fpJetEtaProfile = new TProfile *[14];
301 fpJetAbsEtaProfile = new TProfile *[14];
302 fpChargedJetRProfile = new TProfile *[8];
303 fpJetRProfile = new TProfile *[4];
304 fpChargedJetEDProfile= new TProfile3D *[10];
305 fpJetEDProfile= new TProfile3D *[10];
306 fvertex[0]=0.0,fvertex[1]=0.0,fvertex[2]=0.0;
308 DefineOutput(1,TList::Class()); // for output list
311 //________________________________________________________________________
312 AliAnalysisTaskFullpAJets::~AliAnalysisTaskFullpAJets()
314 // Destructor. Clean-up the output list, but not the histograms that are put inside
315 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
316 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())
322 //________________________________________________________________________
323 void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
326 // Called once (on the worker node)
327 fIsInitialized=kFALSE;
328 fOutput = new TList();
329 fOutput->SetOwner(); // IMPORTANT!
331 // Initialize Global Variables
336 // fRJET=4 -> fJetR=0.4 && fRJET=25 -> fJetR=0.25, but for writing files, should be 4 and 25 respectively
339 fJetR=(Double_t)fRJET/100.0;
343 fJetR=(Double_t)fRJET/10.0;
347 fEMCalPhiMin=(80/(double)360)*2*TMath::Pi();
348 fEMCalPhiMax=(187/(double)360)*2*TMath::Pi();
349 fEMCalPhiTotal= fEMCalPhiMax-fEMCalPhiMin;
352 fEMCalEtaTotal=fEMCalEtaMax-fEMCalEtaMin;
353 fEMCalArea=fEMCalPhiTotal*fEMCalEtaTotal;
355 fTPCPhiMin=(0/(double)360)*2*TMath::Pi();
356 fTPCPhiMax=(360/(double)360)*2*TMath::Pi();
357 fTPCPhiTotal= fTPCPhiMax-fTPCPhiMin;
360 fTPCEtaTotal=fTPCEtaMax-fTPCEtaMin;
361 fTPCArea=fTPCPhiTotal*fTPCEtaTotal;
366 Int_t CentralityBinMult=10;
368 fJetAreaCutFrac =0.6; // Fudge factor for selecting on jets with threshold Area or higher
369 fJetAreaThreshold=fJetAreaCutFrac*TMath::Pi()*fJetR*fJetR;
370 fTPCJetThreshold=5.0; // Threshold required for an Anti-kT Charged jet to be considered a "true" jet in TPC
371 fEMCalJetThreshold=5.0; // Threshold required for an Anti-kT Charged+Neutral jet to be considered a "true" jet in EMCal
376 fEMCalRCBckgFluc = new Double_t[fnBckgClusters];
377 fTPCRCBckgFluc = new Double_t[fnBckgClusters];
378 fEMCalRCBckgFlucSignal = new Double_t[fnBckgClusters];
379 fTPCRCBckgFlucSignal = new Double_t[fnBckgClusters];
380 for (Int_t i=0;i<fnBckgClusters;i++)
382 fEMCalRCBckgFluc[i]=0.0;
383 fTPCRCBckgFluc[i]=0.0;
384 fEMCalRCBckgFlucSignal[i]=0.0;
385 fTPCRCBckgFlucSignal[i]=0.0;
388 fnEMCalCells=12288; // sMods 1-10 have 24x48 cells, sMods 11&12 have 8x48 cells...
391 Int_t JetPtBins = 200;
392 Double_t JetPtLow = 0.0;
393 Double_t JetPtUp = 200.0;
398 fhTrackPt = new TH1D("fhTrackPt","p_{T} distribution of tracks in event",10*JetPtBins,JetPtLow,JetPtUp);
399 fhTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
400 fhTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
403 fhTrackPhi = new TH1D("fhTrackPhi","#phi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
404 fhTrackPhi->GetXaxis()->SetTitle("#phi");
405 fhTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#phi");
408 fhTrackEta = new TH1D("fhTrackEta","#eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
409 fhTrackEta->GetXaxis()->SetTitle("#eta");
410 fhTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
413 fhTrackEtaPhi = new TH2D("fhTrackEtaPhi","#eta-#phi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
414 fhTrackEtaPhi->GetXaxis()->SetTitle("#eta");
415 fhTrackEtaPhi->GetYaxis()->SetTitle("#phi");
416 fhTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
417 fhTrackEtaPhi->Sumw2();
419 fhClusterPt = new TH1D("fhClusterPt","p_{T} distribution of clusters in event",10*JetPtBins,JetPtLow,JetPtUp);
420 fhClusterPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
421 fhClusterPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
422 fhClusterPt->Sumw2();
424 fhClusterPhi = new TH1D("fhClusterPhi","#phi distribution of clusters in event",TCBins,fTPCPhiMin,fTPCPhiMax);
425 fhClusterPhi->GetXaxis()->SetTitle("#phi");
426 fhClusterPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#phi");
427 fhClusterPhi->Sumw2();
429 fhClusterEta = new TH1D("fhClusterEta","#eta distribution of clusters in event",TCBins,fTPCEtaMin,fTPCEtaMax);
430 fhClusterEta->GetXaxis()->SetTitle("#eta");
431 fhClusterEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
432 fhClusterEta->Sumw2();
434 fhClusterEtaPhi = new TH2D("fhClusterEtaPhi","#eta-#phi distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
435 fhClusterEtaPhi->GetXaxis()->SetTitle("#eta");
436 fhClusterEtaPhi->GetYaxis()->SetTitle("#phi");
437 fhClusterEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
438 fhClusterEtaPhi->Sumw2();
440 fhCentrality = new TH1D("fhCentrality","Event Centrality Distribution",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
441 fhCentrality->GetXaxis()->SetTitle(fCentralityTag);
442 fhCentrality->GetYaxis()->SetTitle("1/N_{Events}");
443 fhCentrality->Sumw2();
445 fhJetConstituentPt= new TH2D("fhJetConstituentPt","Jet constituents p_{T} distribution",JetPtBins, JetPtLow, JetPtUp,10*JetPtBins, JetPtLow, JetPtUp);
446 fhJetConstituentPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
447 fhJetConstituentPt->GetYaxis()->SetTitle("Constituent p_{T} (GeV/c)");
448 fhJetConstituentPt->Sumw2();
450 fhEMCalCellCounts = new TH1D("fhEMCalCellCounts","Distribtuion of cluster counts across the EMCal",fnEMCalCells,1,fnEMCalCells);
451 fhEMCalCellCounts->GetXaxis()->SetTitle("Absoulute Cell Id");
452 fhEMCalCellCounts->GetYaxis()->SetTitle("Counts per Event");
453 fhEMCalCellCounts->Sumw2();
455 // Jet Area vs pT Distribution
456 Int_t JetPtAreaBins=200;
457 Double_t JetPtAreaLow=0.0;
458 Double_t JetPtAreaUp=2.0;
464 fhJetPtArea = new TH2D("fhJetPtArea","Jet Area Distribution",JetPtBins, JetPtLow,JetPtUp,JetPtAreaBins,JetPtAreaLow,JetPtAreaUp);
465 fhJetPtArea->GetXaxis()->SetTitle("p_{T} (GeV/c)");
466 fhJetPtArea->GetYaxis()->SetTitle("A_{jet}");
467 fhJetPtArea->GetZaxis()->SetTitle("1/N_{Events} dN/dA_{jet}dp_{T}");
468 fhJetPtArea->Sumw2();
470 fhRhoScale = new TH2D("fhRhoScale","Scaling Factor",SFBins,SFLow,SFUp,CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
471 fhRhoScale->GetXaxis()->SetTitle("Scale Factor");
472 fhRhoScale->GetYaxis()->SetTitle("Centrality");
473 fhRhoScale->GetZaxis()->SetTitle("Counts");
477 fpEMCalEventMult = new TProfile("fpEMCalEventMult","EMCal Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
478 fpEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag);
479 fpEMCalEventMult->GetYaxis()->SetTitle("Multiplicity");
481 fpTPCEventMult = new TProfile("fpTPCEventMult","TPC Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
482 fpTPCEventMult->GetXaxis()->SetTitle(fCentralityTag);
483 fpTPCEventMult->GetYaxis()->SetTitle("Multiplicity");
485 fpRhoScale = new TProfile("fpRhoScale","Scaling Factor Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
486 fpRhoScale->GetXaxis()->SetTitle(fCentralityTag);
487 fpRhoScale->GetYaxis()->SetTitle("Scale Factor");
489 // QA::2D Energy Density Profiles for Tracks and Clusters
490 fpTrackPtProfile = new TProfile2D("fpTrackPtProfile","2D Profile of track pT density throughout the TPC",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
491 fpTrackPtProfile->GetXaxis()->SetTitle("#eta");
492 fpTrackPtProfile->GetYaxis()->SetTitle("#phi");
493 fpTrackPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
495 fpClusterPtProfile = new TProfile2D("fpClusterPtProfile","2D Profile of cluster pT density throughout the EMCal",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
496 fpClusterPtProfile->GetXaxis()->SetTitle("#eta");
497 fpClusterPtProfile->GetYaxis()->SetTitle("#phi");
498 fpClusterPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
500 TString temp_name="";
501 TString title_name="";
507 for (Int_t i=0;i<14;i++)
509 temp_name=Form("fpJetEtaProfile%d",i);
510 title_name=Form("Jet Energy Density #eta Profile for ALL p_{T}, 0-20 Centrality, and eta=%g to %g",(i-7)/10.,(i-6)/10.);
512 fpJetEtaProfile[i]= new TProfile(temp_name,title_name,fEtaProfileBins,fEtaProfileLow,fEtaProfileUp);
513 fpJetEtaProfile[i]->GetXaxis()->SetTitle("#eta");
514 fpJetEtaProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
515 fOutput->Add(fpJetEtaProfile[i]);
519 temp_name=Form("fpJetAbsEtaProfile%d",i);
520 title_name=Form("Jet Energy Density #Delta #eta Profile for ALL p_{T}, 0-20 Centrality, and eta=%g to %g",(i-7)/10.,(i-6)/10.);
522 fpJetAbsEtaProfile[i]= new TProfile(temp_name,title_name,fEtaProfileBins,0,2*fEtaProfileUp);
523 fpJetAbsEtaProfile[i]->GetXaxis()->SetTitle("#Delta#eta");
524 fpJetAbsEtaProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
525 fOutput->Add(fpJetAbsEtaProfile[i]);
533 fEDProfilePtBins=100;
535 fEDProfilePtUp=100.0;
537 fEDProfileEtaLow=-0.2;
540 for (Int_t i=0;i<8;i++)
542 temp_name=Form("fpChargedJetRProfile%d",i);
543 title_name=Form("Charged Jet Energy Density Radius Profile for ALL p_{T}, 0-20 Centrality, and eta=%g to %g",(i-4)/10.,(i-3)/10.);
545 fpChargedJetRProfile[i]= new TProfile(temp_name,title_name,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
546 fpChargedJetRProfile[i]->GetXaxis()->SetTitle("Radius");
547 fpChargedJetRProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
548 fOutput->Add(fpChargedJetRProfile[i]);
553 for (Int_t i=0;i<4;i++)
555 temp_name=Form("fpJetRProfile%d",i);
556 title_name=Form("Jet Energy Density Radius Profile for ALL p_{T}, 0-20 Centrality, and eta=%g to %g",(i-2)/10.,(i-1)/10.);
558 fpJetRProfile[i]= new TProfile(temp_name,title_name,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
559 fpJetRProfile[i]->GetXaxis()->SetTitle("Radius");
560 fpJetRProfile[i]->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
561 fOutput->Add(fpJetRProfile[i]);
566 for (Int_t i=0;i<10;i++)
568 temp_name=Form("fpChargedJetEDProfile%d0",i);
569 title_name=Form("Charged Jet Energy Density Profile for %d0-%d0 Centrality",i,i+1);
571 fpChargedJetEDProfile[i]= new TProfile3D(temp_name,title_name,fEDProfilePtBins,fEDProfilePtLow,fEDProfilePtUp,fEDProfileEtaBins+4,fEDProfileEtaLow-0.2,fEDProfileEtaUp+0.2,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
572 fpChargedJetEDProfile[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
573 fpChargedJetEDProfile[i]->GetYaxis()->SetTitle("Pseudorapidity");
574 fpChargedJetEDProfile[i]->GetZaxis()->SetTitle("Radius");
575 fOutput->Add(fpChargedJetEDProfile[i]);
579 temp_name=Form("fpJetEDProfile%d0",i);
580 title_name=Form("Jet Energy Density Profile for %d0-%d0 Centrality",i,i+1);
582 fpJetEDProfile[i]= new TProfile3D(temp_name,title_name,fEDProfilePtBins,fEDProfilePtLow,fEDProfilePtUp,fEDProfileEtaBins,fEDProfileEtaLow,fEDProfileEtaUp,fEDProfileRBins,fEDProfileRLow,fEDProfileRUp);
583 fpJetEDProfile[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
584 fpJetEDProfile[i]->GetYaxis()->SetTitle("Pseudorapidity");
585 fpJetEDProfile[i]->GetZaxis()->SetTitle("Radius");
586 fOutput->Add(fpJetEDProfile[i]);
591 fTPCRawJets = new AlipAJetHistos("TPCRawJets",fCentralityTag);
592 fEMCalRawJets = new AlipAJetHistos("EMCalRawJets",fCentralityTag);
594 fRhoFull0 = new AlipAJetHistos("RhoFull0",fCentralityTag);
595 fRhoFull1 = new AlipAJetHistos("RhoFull1",fCentralityTag);
596 fRhoFull2 = new AlipAJetHistos("RhoFull2",fCentralityTag);
597 fRhoFullN = new AlipAJetHistos("RhoFullN",fCentralityTag);
598 fRhoFullDijet = new AlipAJetHistos("RhoFullDijet",fCentralityTag);
599 fRhoFullkT = new AlipAJetHistos("RhoFullkT",fCentralityTag);
600 fRhoFullCMS = new AlipAJetHistos("RhoFullCMS",fCentralityTag);
602 fRhoCharged0 = new AlipAJetHistos("RhoCharged0",fCentralityTag);
603 fRhoCharged1 = new AlipAJetHistos("RhoCharged1",fCentralityTag);
604 fRhoCharged2 = new AlipAJetHistos("RhoCharged2",fCentralityTag);
605 fRhoChargedN = new AlipAJetHistos("RhoChargedN",fCentralityTag);
606 fRhoChargedkT = new AlipAJetHistos("RhoChargedkT",fCentralityTag);
607 fRhoChargedScale = new AlipAJetHistos("RhoChargedScale",fCentralityTag);
608 fRhoChargedkTScale = new AlipAJetHistos("RhoChargedkTScale",fCentralityTag);
609 fRhoChargedCMS = new AlipAJetHistos("RhoChargedCMS",fCentralityTag);
610 fRhoChargedCMSScale = new AlipAJetHistos("RhoChargedCMSScale",fCentralityTag);
612 fOutput->Add(fhTrackPt);
613 fOutput->Add(fhTrackEta);
614 fOutput->Add(fhTrackPhi);
615 fOutput->Add(fhTrackEtaPhi);
616 fOutput->Add(fhClusterPt);
617 fOutput->Add(fhClusterEta);
618 fOutput->Add(fhClusterPhi);
619 fOutput->Add(fhClusterEtaPhi);
620 fOutput->Add(fhCentrality);
621 fOutput->Add(fhEMCalCellCounts);
622 fOutput->Add(fhJetPtArea);
623 fOutput->Add(fhJetConstituentPt);
624 fOutput->Add(fhRhoScale);
626 fOutput->Add(fpTPCEventMult);
627 fOutput->Add(fpEMCalEventMult);
628 fOutput->Add(fpRhoScale);
630 fOutput->Add(fpTrackPtProfile);
631 fOutput->Add(fpClusterPtProfile);
633 fOutput->Add(fTPCRawJets->GetOutputHistos());
634 fOutput->Add(fEMCalRawJets->GetOutputHistos());
636 fOutput->Add(fRhoFull0->GetOutputHistos());
637 fOutput->Add(fRhoFull1->GetOutputHistos());
638 fOutput->Add(fRhoFull2->GetOutputHistos());
639 fOutput->Add(fRhoFullN->GetOutputHistos());
640 fOutput->Add(fRhoFullDijet->GetOutputHistos());
641 fOutput->Add(fRhoFullkT->GetOutputHistos());
642 fOutput->Add(fRhoFullCMS->GetOutputHistos());
643 fOutput->Add(fRhoCharged0->GetOutputHistos());
644 fOutput->Add(fRhoCharged1->GetOutputHistos());
645 fOutput->Add(fRhoCharged2->GetOutputHistos());
646 fOutput->Add(fRhoChargedN->GetOutputHistos());
647 fOutput->Add(fRhoChargedkT->GetOutputHistos());
648 fOutput->Add(fRhoChargedScale->GetOutputHistos());
649 fOutput->Add(fRhoChargedkTScale->GetOutputHistos());
650 fOutput->Add(fRhoChargedCMS->GetOutputHistos());
651 fOutput->Add(fRhoChargedCMSScale->GetOutputHistos());
653 // Post data for ALL output slots >0 here,
654 // To get at least an empty histogram
655 // 1 is the outputnumber of a certain weg of task 1
656 PostData(1, fOutput);
659 void AliAnalysisTaskFullpAJets::UserExecOnce()
661 // Get the event tracks from PicoTracks
662 TString track_name="PicoTracks";
663 fOrgTracks = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(track_name));
665 // Get the event caloclusters from CaloClustersCorr
666 TString cluster_name="caloClustersCorr";
667 //TString cluster_name="EmcalClusters";
668 fOrgClusters = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(cluster_name));
671 TString jet_algorithm=Form("Jet_AKTChargedR0%d0_PicoTracks_pT0150",fRJET);
672 fmyAKTChargedJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
674 jet_algorithm=Form("Jet_KTChargedR0%d0_PicoTracks_pT0150",fRJET);
675 fmyKTChargedJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
678 jet_algorithm=Form("Jet_AKTFullR0%d0_PicoTracks_pT0150_caloClustersCorr_ET0300",fRJET);
679 fmyAKTFullJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
681 jet_algorithm=Form("Jet_KTFullR0%d0_PicoTracks_pT0150_caloClustersCorr_ET0300",fRJET);
682 fmyKTFullJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
685 fIsInitialized=kTRUE;
687 //________________________________________________________________________
688 void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
690 if (fIsInitialized==kFALSE)
695 // Get pointer to reconstructed event
696 AliVEvent *event = InputEvent();
699 AliError("Pointer == 0, this can not happen!");
703 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
704 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
708 fEventCentrality=esd->GetCentrality()->GetCentralityPercentile(fCentralityTag);
710 if (esd->GetPrimaryVertex()->GetNContributors()<1 || (TMath::Abs(esd->GetPrimaryVertex()->GetZ())>fVertexWindow))
714 if (TMath::Sqrt(TMath::Power(esd->GetPrimaryVertex()->GetX(),2)+TMath::Power(esd->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
719 esd->GetPrimaryVertex()->GetXYZ(fvertex);
723 fEventCentrality=aod->GetCentrality()->GetCentralityPercentile(fCentralityTag);
725 if (aod->GetPrimaryVertex()->GetNContributors()<1 || (TMath::Abs(aod->GetPrimaryVertex()->GetZ())>fVertexWindow))
729 if (TMath::Sqrt(TMath::Power(aod->GetPrimaryVertex()->GetX(),2)+TMath::Power(aod->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
734 aod->GetPrimaryVertex()->GetXYZ(fvertex);
738 AliError("Cannot get AOD/ESD event. Rejecting Event");
742 // Make sure Centrality isn't exactly 100% (to avoid bin filling errors in profile plots. Make it 99.99
743 if (fEventCentrality>99.99)
745 fEventCentrality=99.99;
747 fhCentrality->Fill(fEventCentrality);
751 // Reject any event that doesn't have any tracks, i.e. TPC is off
754 AliWarning("No PicoTracks, Rejecting Event");
762 AliInfo("No Corrected CaloClusters, using only charged jets");
766 GenerateTPCRandomConesPt();
769 EstimateChargedRho0();
770 EstimateChargedRho1();
771 EstimateChargedRho2();
772 EstimateChargedRhoN();
773 EstimateChargedRhokT();
774 EstimateChargedRhoCMS();
776 JetPtChargedProfile();
777 DeleteJetData(kFALSE);
781 PostData(1, fOutput);
792 GenerateTPCRandomConesPt();
793 GenerateEMCalRandomConesPt();
796 EstimateChargedRho0();
797 EstimateChargedRho1();
798 EstimateChargedRho2();
799 EstimateChargedRhoN();
800 EstimateChargedRhokT();
801 EstimateChargedRhoCMS();
808 EstimateFullRhoCMS();
810 EstimateChargedRhoScale();
811 EstimateChargedRhokTScale();
812 EstimateChargedRhoCMSScale();
815 if (IsDiJetEvent()==kTRUE)
817 EstimateFullRhoDijet();
820 // Compute Jet Energy Density Profile
821 JetPtChargedProfile();
825 // Delete Dynamic Arrays
826 DeleteJetData(kTRUE);
829 PostData(1, fOutput);
832 //________________________________________________________________________
833 void AliAnalysisTaskFullpAJets::Terminate(Option_t *) //specify what you want to have done
835 // Called once at the end of the query. Done nothing here.
838 /////////////////////////////////////////////////////////////////////////////////////////
839 ///////////////// User Defined Sub_Routines ///////////////////////////////////////
840 /////////////////////////////////////////////////////////////////////////////////////////
842 void AliAnalysisTaskFullpAJets::TrackCuts()
844 // Fill a TObjArray with the tracks from a TClonesArray which grabs the picotracks.
847 fmyTracks = new TObjArray();
848 for (i=0;i<fOrgTracks->GetEntries();i++)
850 AliVTrack* vtrack = (AliVTrack*) fOrgTracks->At(i);
851 if (vtrack->Pt()>=fTrackMinPt)
853 fmyTracks->Add(vtrack);
856 fnTracks = fmyTracks->GetEntries();
859 void AliAnalysisTaskFullpAJets::ClusterCuts()
861 // Fill a TObjArray with the clusters from a TClonesArray which grabs the caloclusterscorr.
864 fmyClusters = new TObjArray();
866 for (i=0;i<fOrgClusters->GetEntries();i++)
868 AliVCluster* vcluster = (AliVCluster*) fOrgClusters->At(i);
869 TLorentzVector *cluster_vec = new TLorentzVector;
870 vcluster->GetMomentum(*cluster_vec,fvertex);
872 if (cluster_vec->Pt()>=fClusterMinPt)
874 fmyClusters->Add(vcluster);
880 fnClusters = fmyClusters->GetEntries();
883 void AliAnalysisTaskFullpAJets::TrackHisto()
885 // Fill track histograms: Phi,Eta,Pt
888 TH2D *hdummypT= new TH2D("hdummypT","",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax); //!
890 for (i=0;i<fnTracks;i++)
892 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
893 fhTrackPt->Fill(vtrack->Pt());
894 fhTrackEta->Fill(vtrack->Eta());
895 fhTrackPhi->Fill(vtrack->Phi());
896 fhTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
897 hdummypT->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
899 for (i=1;i<=TCBins;i++)
901 for (j=1;j<=TCBins;j++)
903 fpTrackPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fTPCArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
909 void AliAnalysisTaskFullpAJets::ClusterHisto()
911 // Fill cluster histograms: Phi,Eta,Pt
914 TH2D *hdummypT= new TH2D("hdummypT","",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax); //!
915 AliEMCALGeometry *myAliEMCGeo = AliEMCALGeometry::GetInstance();
918 for (i=0;i<fnClusters;i++)
920 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
921 TLorentzVector *cluster_vec = new TLorentzVector;
922 vcluster->GetMomentum(*cluster_vec,fvertex);
924 fhClusterPt->Fill(cluster_vec->Pt());
925 fhClusterEta->Fill(cluster_vec->Eta());
926 fhClusterPhi->Fill(cluster_vec->Phi());
927 fhClusterEtaPhi->Fill(cluster_vec->Eta(),cluster_vec->Phi());
928 hdummypT->Fill(cluster_vec->Eta(),cluster_vec->Phi(),cluster_vec->Pt());
929 myAliEMCGeo->GetAbsCellIdFromEtaPhi(cluster_vec->Eta(),cluster_vec->Phi(),myCellID);
930 fhEMCalCellCounts->Fill(myCellID);
933 for (i=1;i<=TCBins;i++)
935 for (j=1;j<=TCBins;j++)
937 fpClusterPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fEMCalArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
943 void AliAnalysisTaskFullpAJets::InitChargedJets()
945 // Preliminary Jet Placement and Selection Cuts
948 fnAKTChargedJets = fmyAKTChargedJets->GetEntries();
949 fnKTChargedJets = fmyKTChargedJets->GetEntries();
951 fTPCJet = new AlipAJetData("fTPCJet",kFALSE,fnAKTChargedJets);
952 fTPCFullJet = new AlipAJetData("fTPCFullJet",kFALSE,fnAKTChargedJets);
953 fTPCOnlyJet = new AlipAJetData("fTPCOnlyJet",kFALSE,fnAKTChargedJets);
955 fTPCJet->SetSignalCut(fTPCJetThreshold);
956 fTPCJet->SetAreaCutFraction(fJetAreaCutFrac);
957 fTPCJet->SetJetR(fJetR);
958 fTPCFullJet->SetSignalCut(fTPCJetThreshold);
959 fTPCFullJet->SetAreaCutFraction(fJetAreaCutFrac);
960 fTPCFullJet->SetJetR(fJetR);
961 fTPCOnlyJet->SetSignalCut(fTPCJetThreshold);
962 fTPCOnlyJet->SetAreaCutFraction(fJetAreaCutFrac);
963 fTPCOnlyJet->SetJetR(fJetR);
965 // Initialize Jet Data
966 for (i=0;i<fnAKTChargedJets;i++)
968 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(i);
970 fTPCJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kFALSE),i);
971 fTPCFullJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE),i);
972 fTPCOnlyJet->SetIsJetInArray(IsInTPCFull(fJetR,myJet->Phi(),myJet->Eta()),i);
974 fTPCJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
975 fTPCFullJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
976 fTPCOnlyJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
979 fTPCkTFullJet = new AlipAJetData("fTPCkTFullJet",kFALSE,fnKTChargedJets);
980 fTPCkTFullJet->SetSignalCut(fTPCJetThreshold);
981 fTPCkTFullJet->SetAreaCutFraction(0.25*fJetAreaCutFrac);
982 fTPCkTFullJet->SetJetR(fJetR);
984 for (i=0;i<fnKTChargedJets;i++)
986 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(i);
987 fTPCkTFullJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE),i);
989 fTPCkTFullJet->InitializeJetData(fmyKTChargedJets,fnKTChargedJets);
991 // Raw Charged Jet Spectra
992 fTPCRawJets->FillBSJS(fEventCentrality,0.0,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
995 void AliAnalysisTaskFullpAJets::InitFullJets()
997 // Preliminary Jet Placement and Selection Cuts
1000 fnAKTFullJets = fmyAKTFullJets->GetEntries();
1001 fnKTFullJets = fmyKTFullJets->GetEntries();
1003 fEMCalJet = new AlipAJetData("fEMCalJet",kTRUE,fnAKTFullJets);
1004 fEMCalFullJet = new AlipAJetData("fEMCalFullJet",kTRUE,fnAKTFullJets);
1005 fEMCalPartJet = new AlipAJetData("fEMCalPartJet",kTRUE,fnAKTFullJets);
1007 fEMCalJet->SetSignalCut(fEMCalJetThreshold);
1008 fEMCalJet->SetAreaCutFraction(fJetAreaCutFrac);
1009 fEMCalJet->SetJetR(fJetR);
1010 fEMCalFullJet->SetSignalCut(fEMCalJetThreshold);
1011 fEMCalFullJet->SetAreaCutFraction(fJetAreaCutFrac);
1012 fEMCalFullJet->SetJetR(fJetR);
1013 fEMCalPartJet->SetSignalCut(fEMCalJetThreshold);
1014 fEMCalPartJet->SetAreaCutFraction(fJetAreaCutFrac);
1015 fEMCalPartJet->SetJetR(fJetR);
1017 // Initialize Jet Data
1018 for (i=0;i<fnAKTFullJets;i++)
1020 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
1022 fEMCalJet->SetIsJetInArray(IsInEMCal(myJet->Phi(),myJet->Eta()),i);
1023 fEMCalFullJet->SetIsJetInArray(IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta()),i);
1024 fEMCalPartJet->SetIsJetInArray(IsInEMCalPart(fJetR,myJet->Phi(),myJet->Eta()),i);
1026 fEMCalJet->InitializeJetData(fmyAKTFullJets,fnAKTFullJets);
1027 fEMCalFullJet->InitializeJetData(fmyAKTFullJets,fnAKTFullJets);
1028 fEMCalPartJet->InitializeJetData(fmyAKTFullJets,fnAKTFullJets);
1031 fEMCalkTFullJet = new AlipAJetData("fEMCalkTFullJet",kTRUE,fnKTFullJets);
1032 fEMCalkTFullJet->SetSignalCut(fEMCalJetThreshold);
1033 fEMCalkTFullJet->SetAreaCutFraction(0.25*fJetAreaCutFrac);
1034 fEMCalkTFullJet->SetJetR(fJetR);
1036 for (i=0;i<fnKTFullJets;i++)
1038 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(i);
1039 fEMCalkTFullJet->SetIsJetInArray(IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta()),i);
1041 fEMCalkTFullJet->InitializeJetData(fmyKTFullJets,fnKTFullJets);
1043 // Raw Full Jet Spectra
1044 fEMCalRawJets->FillBSJS(fEventCentrality,0.0,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
1047 void AliAnalysisTaskFullpAJets::JetPtArea()
1051 for (i=0;i<fEMCalFullJet->GetTotalJets();i++)
1053 AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalFullJet->GetJetIndex(i));
1054 fhJetPtArea->Fill(myJet->Pt(),myJet->Area());
1058 void AliAnalysisTaskFullpAJets::GenerateTPCRandomConesPt()
1061 Double_t E_tracks_total=0.;
1062 TRandom3 u(time(NULL));
1064 Double_t Eta_Center=0.5*(fTPCEtaMin+fTPCEtaMax);
1065 Double_t Phi_Center=0.5*(fTPCPhiMin+fTPCPhiMax);
1069 for (i=0;i<fnBckgClusters;i++)
1071 fTPCRCBckgFluc[i]=0.0;
1072 fTPCRCBckgFlucSignal[i]=0.0;
1075 TLorentzVector *dummy= new TLorentzVector;
1076 TLorentzVector *temp_jet= new TLorentzVector;
1078 // First, consider the RC with no spatial restrictions
1079 for (j=0;j<fnBckgClusters;j++)
1083 dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1084 // Loop over all tracks
1085 for (i=0;i<fnTracks;i++)
1087 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1088 if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1090 TLorentzVector *track_vec = new TLorentzVector;
1091 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1092 if (dummy->DeltaR(*track_vec)<fJetR)
1094 E_tracks_total+=vtrack->Pt();
1099 fTPCRCBckgFlucSignal[j]=E_tracks_total;
1102 // Now, consider the RC where the vertex of RC is at least 2R away from the leading signal
1104 if (fTPCJet->GetLeadingPt()<0.0)
1106 temp_jet->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1110 AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1111 myJet->GetMom(*temp_jet);
1114 for (j=0;j<fnBckgClusters;j++)
1120 dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1121 while (temp_jet->DeltaR(*dummy)<fJetR)
1123 dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1125 // Loop over all tracks
1126 for (i=0;i<fnTracks;i++)
1128 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1129 if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1132 TLorentzVector *track_vec = new TLorentzVector;
1133 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1134 if (dummy->DeltaR(*track_vec)<fJetR)
1137 E_tracks_total+=vtrack->Pt();
1142 fTPCRCBckgFluc[j]=E_tracks_total;
1144 fpTPCEventMult->Fill(fEventCentrality,event_mult);
1145 fTPCRawJets->FillDeltaPt(fEventCentrality,0.0,fJetR,fTPCRCBckgFluc,1);
1151 void AliAnalysisTaskFullpAJets::GenerateEMCalRandomConesPt()
1154 Double_t E_tracks_total=0.;
1155 Double_t E_caloclusters_total=0.;
1156 TRandom3 u(time(NULL));
1158 Double_t Eta_Center=0.5*(fEMCalEtaMin+fEMCalEtaMax);
1159 Double_t Phi_Center=0.5*(fEMCalPhiMin+fEMCalPhiMax);
1163 for (i=0;i<fnBckgClusters;i++)
1165 fEMCalRCBckgFluc[i]=0.0;
1166 fEMCalRCBckgFlucSignal[i]=0.0;
1169 TLorentzVector *dummy= new TLorentzVector;
1170 TLorentzVector *temp_jet= new TLorentzVector;
1172 // First, consider the RC with no spatial restrictions
1173 for (j=0;j<fnBckgClusters;j++)
1176 E_caloclusters_total=0.;
1178 dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1179 // Loop over all tracks
1180 for (i=0;i<fnTracks;i++)
1182 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1183 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1185 TLorentzVector *track_vec = new TLorentzVector;
1186 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1187 if (dummy->DeltaR(*track_vec)<fJetR)
1189 E_tracks_total+=vtrack->Pt();
1195 // Loop over all caloclusters
1196 for (i=0;i<fnClusters;i++)
1198 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1199 TLorentzVector *cluster_vec = new TLorentzVector;
1200 vcluster->GetMomentum(*cluster_vec,fvertex);
1201 if (dummy->DeltaR(*cluster_vec)<fJetR)
1204 E_caloclusters_total+=vcluster->E();
1208 fEMCalRCBckgFlucSignal[j]=E_tracks_total+E_caloclusters_total;
1211 // Now, consider the RC where the vertex of RC is at least 2R away from the leading signal
1213 E_caloclusters_total=0.;
1214 if (fEMCalPartJet->GetLeadingPt()<0.0)
1216 temp_jet->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1220 AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
1221 myJet->GetMom(*temp_jet);
1224 for (j=0;j<fnBckgClusters;j++)
1229 E_caloclusters_total=0.;
1231 dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1232 while (temp_jet->DeltaR(*dummy)<fJetR)
1234 dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1236 // Loop over all tracks
1237 for (i=0;i<fnTracks;i++)
1239 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1240 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1243 TLorentzVector *track_vec = new TLorentzVector;
1244 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1245 if (dummy->DeltaR(*track_vec)<fJetR)
1248 E_tracks_total+=vtrack->Pt();
1254 // Loop over all caloclusters
1255 for (i=0;i<fnClusters;i++)
1257 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1258 TLorentzVector *cluster_vec = new TLorentzVector;
1259 vcluster->GetMomentum(*cluster_vec,fvertex);
1261 if (dummy->DeltaR(*cluster_vec)<fJetR)
1264 E_caloclusters_total+=vcluster->E();
1268 fEMCalRCBckgFluc[j]=E_tracks_total+E_caloclusters_total;
1270 fpEMCalEventMult->Fill(fEventCentrality,event_mult);
1271 fEMCalRawJets->FillDeltaPt(fEventCentrality,0.0,fJetR,fEMCalRCBckgFluc,1);
1278 void AliAnalysisTaskFullpAJets::EstimateChargedRho0()
1281 Double_t E_tracks_total=0.0;
1282 Double_t TPC_rho=0.;
1284 // Loop over all tracks
1285 for (i=0;i<fnTracks;i++)
1287 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1288 if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1290 E_tracks_total+=vtrack->Pt();
1294 // Calculate the mean Background density
1295 TPC_rho=E_tracks_total/fTPCArea;
1296 fRhoCharged=TPC_rho;
1299 fRhoCharged0->FillRho(fEventCentrality,TPC_rho);
1300 fRhoCharged0->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCJet->GetJets(),fTPCJet->GetTotalJets());
1301 fRhoCharged0->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
1302 fRhoCharged0->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
1303 fRhoCharged0->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
1304 fRhoCharged0->FillLeadingJetPtRho(fTPCJet->GetLeadingPt(),TPC_rho);
1308 void AliAnalysisTaskFullpAJets::EstimateChargedRho1()
1311 Double_t E_tracks_total=0.0;
1312 Double_t TPC_rho=0.;
1314 if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
1316 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1317 TLorentzVector *temp_jet= new TLorentzVector;
1318 myJet->GetMom(*temp_jet);
1320 // Loop over all tracks
1321 for (i=0;i<fnTracks;i++)
1323 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1324 if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1326 TLorentzVector *track_vec = new TLorentzVector;
1327 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1328 if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
1330 E_tracks_total+=vtrack->Pt();
1337 // Calculate the mean Background density
1338 TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myJet->Eta()));
1340 else // i.e. No signal jets -> same as total background density
1342 // Loop over all tracks
1343 for (i=0;i<fnTracks;i++)
1345 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1346 if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1348 E_tracks_total+=vtrack->Pt();
1351 // Calculate the mean Background density
1352 TPC_rho=E_tracks_total/fTPCArea;
1356 fRhoCharged1->FillRho(fEventCentrality,TPC_rho);
1357 fRhoCharged1->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
1358 fRhoCharged1->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
1359 fRhoCharged1->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
1360 fRhoCharged1->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
1361 fRhoCharged1->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),TPC_rho);
1364 void AliAnalysisTaskFullpAJets::EstimateChargedRho2()
1367 Double_t E_tracks_total=0.0;
1368 Double_t TPC_rho=0.;
1370 if ((fTPCJet->GetLeadingPt()>=fTPCJetThreshold) && (fTPCJet->GetSubLeadingPt()>=fTPCJetThreshold))
1372 AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1373 TLorentzVector *temp_jet1= new TLorentzVector;
1374 myhJet->GetMom(*temp_jet1);
1376 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSubLeadingIndex());
1377 TLorentzVector *temp_jet2= new TLorentzVector;
1378 myJet->GetMom(*temp_jet2);
1380 // Loop over all tracks
1381 for (i=0;i<fnTracks;i++)
1383 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1384 if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1386 TLorentzVector *track_vec = new TLorentzVector;
1387 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1388 if ((temp_jet1->DeltaR(*track_vec)>fJetRForRho) && (temp_jet2->DeltaR(*track_vec)>fJetRForRho))
1390 E_tracks_total+=vtrack->Pt();
1398 // Calculate the mean Background density
1399 TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myhJet->Eta())-AreaWithinTPC(fJetR,myJet->Eta()));
1401 else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
1403 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1404 TLorentzVector *temp_jet= new TLorentzVector;
1405 myJet->GetMom(*temp_jet);
1407 // Loop over all tracks
1408 for (i=0;i<fnTracks;i++)
1410 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1411 if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1413 TLorentzVector *track_vec = new TLorentzVector;
1414 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1415 if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
1417 E_tracks_total+=vtrack->Pt();
1424 // Calculate the mean Background density
1425 TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myJet->Eta()));
1427 else // i.e. No signal jets -> same as total background density
1429 // Loop over all tracks
1430 for (i=0;i<fnTracks;i++)
1432 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1433 if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1435 E_tracks_total+=vtrack->Pt();
1439 // Calculate the mean Background density
1440 TPC_rho=E_tracks_total/fTPCArea;
1444 fRhoCharged2->FillRho(fEventCentrality,TPC_rho);
1445 fRhoCharged2->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
1446 fRhoCharged2->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
1447 fRhoCharged2->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
1448 fRhoCharged2->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
1449 fRhoCharged2->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),TPC_rho);
1452 void AliAnalysisTaskFullpAJets::EstimateChargedRhoN()
1455 Bool_t track_away_from_jet;
1456 Double_t E_tracks_total=0.0;
1457 Double_t TPC_rho=0.0;
1458 Double_t jet_area_total=0.0;
1460 // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
1461 for (i=0;i<fnTracks;i++)
1463 // First, check if track is in the EMCal!!
1464 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
1465 if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1467 if (fTPCJet->GetTotalSignalJets()<1)
1469 E_tracks_total+=vtrack->Pt();
1473 track_away_from_jet=kTRUE;
1475 TLorentzVector *track_vec = new TLorentzVector;
1476 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1477 while (track_away_from_jet==kTRUE && j<fTPCJet->GetTotalSignalJets())
1479 TLorentzVector *jet_vec= new TLorentzVector;
1480 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(j));
1481 myJet->GetMom(*jet_vec);
1482 if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
1484 track_away_from_jet=kFALSE;
1489 if (track_away_from_jet==kTRUE)
1491 E_tracks_total+=vtrack->Pt();
1498 // Determine area of all Jets that are within the EMCal
1499 if (fTPCJet->GetTotalSignalJets()==0)
1505 for (i=0;i<fTPCJet->GetTotalSignalJets();i++)
1507 AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(i));
1508 jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
1513 TPC_rho = E_tracks_total/(fTPCArea-jet_area_total);
1516 fRhoChargedN->FillRho(fEventCentrality,TPC_rho);
1517 fRhoChargedN->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
1518 fRhoChargedN->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
1519 fRhoChargedN->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
1520 fRhoChargedN->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
1521 fRhoChargedN->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),TPC_rho);
1524 void AliAnalysisTaskFullpAJets::EstimateChargedRhoScale()
1527 Bool_t track_away_from_jet;
1528 Double_t E_tracks_total=0.0;
1529 Double_t TPC_rho=0.0;
1530 Double_t jet_area_total=0.0;
1532 // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
1533 for (i=0;i<fnTracks;i++)
1535 // First, check if track is in the EMCal!!
1536 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
1537 if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1539 if (fTPCJet->GetTotalSignalJets()<1)
1541 E_tracks_total+=vtrack->Pt();
1545 track_away_from_jet=kTRUE;
1547 TLorentzVector *track_vec = new TLorentzVector;
1548 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1549 while (track_away_from_jet==kTRUE && j<fTPCJet->GetTotalSignalJets())
1551 TLorentzVector *jet_vec= new TLorentzVector;
1552 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(j));
1553 myJet->GetMom(*jet_vec);
1554 if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
1556 track_away_from_jet=kFALSE;
1561 if (track_away_from_jet==kTRUE)
1563 E_tracks_total+=vtrack->Pt();
1570 // Determine area of all Jets that are within the EMCal
1571 if (fTPCJet->GetTotalSignalJets()==0)
1577 for (i=0;i<fTPCJet->GetTotalSignalJets();i++)
1579 AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(i));
1580 jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
1585 TPC_rho = E_tracks_total/(fTPCArea-jet_area_total);
1586 TPC_rho*=fScaleFactor;
1589 fRhoChargedScale->FillRho(fEventCentrality,TPC_rho);
1590 fRhoChargedScale->FillBSJS(fEventCentrality,TPC_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
1591 fRhoChargedScale->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fEMCalRCBckgFluc,1);
1592 fRhoChargedScale->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fEMCalRCBckgFlucSignal,1);
1593 fRhoChargedScale->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
1594 fRhoChargedScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),TPC_rho);
1597 void AliAnalysisTaskFullpAJets::EstimateChargedRhokT()
1600 Double_t kTRho = 0.0;
1601 Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1602 Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1604 for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1606 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1607 pTArray[i]=myJet->Pt();
1608 RhoArray[i]=myJet->Pt()/myJet->Area();
1611 if (fTPCkTFullJet->GetTotalJets()>=2)
1613 kTRho=MedianRhokT(pTArray,RhoArray,fTPCkTFullJet->GetTotalJets());
1615 fRhoChargedkT->FillRho(fEventCentrality,kTRho);
1616 fRhoChargedkT->FillBSJS(fEventCentrality,kTRho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
1617 fRhoChargedkT->FillDeltaPt(fEventCentrality,kTRho,fJetR,fTPCRCBckgFluc,1);
1618 fRhoChargedkT->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fTPCRCBckgFlucSignal,1);
1619 fRhoChargedkT->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
1620 fRhoChargedkT->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),kTRho);
1626 void AliAnalysisTaskFullpAJets::EstimateChargedRhokTScale()
1629 Double_t kTRho = 0.0;
1630 Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1631 Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1633 for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1635 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1636 pTArray[i]=myJet->Pt();
1637 RhoArray[i]=myJet->Pt()/myJet->Area();
1640 if (fTPCkTFullJet->GetTotalJets()>=2)
1642 kTRho=MedianRhokT(pTArray,RhoArray,fTPCkTFullJet->GetTotalJets());
1643 kTRho*=fScaleFactor;
1645 fRhoChargedkTScale->FillRho(fEventCentrality,kTRho);
1646 fRhoChargedkTScale->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
1647 fRhoChargedkTScale->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
1648 fRhoChargedkTScale->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
1649 fRhoChargedkTScale->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
1650 fRhoChargedkTScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
1656 void AliAnalysisTaskFullpAJets::EstimateChargedRhoCMS()
1659 Double_t kTRho = 0.0;
1660 Double_t CMSTotalkTArea = 0.0;
1661 Double_t CMSTrackArea = 0.0;
1662 Double_t CMSCorrectionFactor = 1.0;
1663 Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1664 Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1667 if ((fTPCJet->GetLeadingPt()>=fTPCJetThreshold) && (fTPCJet->GetSubLeadingPt()>=fTPCJetThreshold))
1669 AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1670 AliEmcalJet *myJet2 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSubLeadingIndex());
1672 for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1674 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1676 CMSTotalkTArea+=myJet->Area();
1677 if (myJet->GetNumberOfTracks()>0)
1679 CMSTrackArea+=myJet->Area();
1681 if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kFALSE)
1683 pTArray[k]=myJet->Pt();
1684 RhoArray[k]=myJet->Pt()/myJet->Area();
1690 kTRho=MedianRhokT(pTArray,RhoArray,k);
1697 else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
1699 AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1701 for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1703 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1705 CMSTotalkTArea+=myJet->Area();
1706 if (myJet->GetNumberOfTracks()>0)
1708 CMSTrackArea+=myJet->Area();
1710 if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE)
1712 pTArray[k]=myJet->Pt();
1713 RhoArray[k]=myJet->Pt()/myJet->Area();
1719 kTRho=MedianRhokT(pTArray,RhoArray,k);
1728 for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1730 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1732 CMSTotalkTArea+=myJet->Area();
1733 if (myJet->GetNumberOfTracks()>0)
1735 CMSTrackArea+=myJet->Area();
1737 pTArray[k]=myJet->Pt();
1738 RhoArray[k]=myJet->Pt()/myJet->Area();
1743 kTRho=MedianRhokT(pTArray,RhoArray,k);
1750 // Scale CMS Rho by Correction factor
1751 if (CMSTotalkTArea==0.0)
1753 CMSCorrectionFactor = 1.0;
1757 //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
1758 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
1760 kTRho*=CMSCorrectionFactor;
1761 fRhoChargedCMS->FillRho(fEventCentrality,kTRho);
1762 fRhoChargedCMS->FillBSJS(fEventCentrality,kTRho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
1763 fRhoChargedCMS->FillDeltaPt(fEventCentrality,kTRho,fJetR,fTPCRCBckgFluc,1);
1764 fRhoChargedCMS->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fTPCRCBckgFlucSignal,1);
1765 fRhoChargedCMS->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
1766 fRhoChargedCMS->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),kTRho);
1771 void AliAnalysisTaskFullpAJets::EstimateChargedRhoCMSScale()
1774 Double_t kTRho = 0.0;
1775 Double_t CMSTotalkTArea = 0.0;
1776 Double_t CMSTrackArea = 0.0;
1777 Double_t CMSCorrectionFactor = 1.0;
1778 Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1779 Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1782 if ((fTPCJet->GetLeadingPt()>=fTPCJetThreshold) && (fTPCJet->GetSubLeadingPt()>=fTPCJetThreshold))
1784 AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1785 AliEmcalJet *myJet2 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSubLeadingIndex());
1787 for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1789 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1791 CMSTotalkTArea+=myJet->Area();
1792 if (myJet->GetNumberOfTracks()>0)
1794 CMSTrackArea+=myJet->Area();
1796 if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kFALSE)
1798 pTArray[k]=myJet->Pt();
1799 RhoArray[k]=myJet->Pt()/myJet->Area();
1805 kTRho=MedianRhokT(pTArray,RhoArray,k);
1812 else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
1814 AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
1816 for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1818 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1820 CMSTotalkTArea+=myJet->Area();
1821 if (myJet->GetNumberOfTracks()>0)
1823 CMSTrackArea+=myJet->Area();
1825 if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE)
1827 pTArray[k]=myJet->Pt();
1828 RhoArray[k]=myJet->Pt()/myJet->Area();
1834 kTRho=MedianRhokT(pTArray,RhoArray,k);
1843 for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1845 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
1847 CMSTotalkTArea+=myJet->Area();
1848 if (myJet->GetNumberOfTracks()>0)
1850 CMSTrackArea+=myJet->Area();
1852 pTArray[k]=myJet->Pt();
1853 RhoArray[k]=myJet->Pt()/myJet->Area();
1858 kTRho=MedianRhokT(pTArray,RhoArray,k);
1865 kTRho*=fScaleFactor;
1866 // Scale CMS Rho by Correction factor
1867 if (CMSTotalkTArea==0.0)
1869 CMSCorrectionFactor = 1.0;
1873 //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
1874 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
1876 kTRho*=CMSCorrectionFactor;
1878 fRhoChargedCMSScale->FillRho(fEventCentrality,kTRho);
1879 fRhoChargedCMSScale->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
1880 fRhoChargedCMSScale->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
1881 fRhoChargedCMSScale->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
1882 fRhoChargedCMSScale->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
1883 fRhoChargedCMSScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
1889 void AliAnalysisTaskFullpAJets::EstimateFullRho0()
1892 Double_t E_tracks_total=0.0;
1893 Double_t E_caloclusters_total=0.0;
1894 Double_t EMCal_rho=0.0;
1896 // Loop over all tracks
1897 for (i=0;i<fnTracks;i++)
1899 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1900 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1902 E_tracks_total+=vtrack->Pt();
1906 // Loop over all caloclusters
1907 for (i=0;i<fnClusters;i++)
1909 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1910 TLorentzVector *cluster_vec = new TLorentzVector;
1911 vcluster->GetMomentum(*cluster_vec,fvertex);
1912 E_caloclusters_total+=cluster_vec->Pt();
1913 //E_caloclusters_total+=0.5*cluster_vec->Pt();
1916 // Calculate the mean Background density
1917 EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
1923 fpRhoScale->Fill(fEventCentrality,fRhoFull/fRhoCharged);
1924 fhRhoScale->Fill(fRhoFull/fRhoCharged,fEventCentrality);
1927 fRhoFull0->FillRho(fEventCentrality,EMCal_rho);
1928 fRhoFull0->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
1929 fRhoFull0->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
1930 fRhoFull0->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
1931 fRhoFull0->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
1932 fRhoFull0->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
1935 void AliAnalysisTaskFullpAJets::EstimateFullRho1()
1938 Double_t E_tracks_total=0.0;
1939 Double_t E_caloclusters_total=0.0;
1940 Double_t EMCal_rho=0.0;
1942 if (fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold)
1944 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
1945 TLorentzVector *temp_jet= new TLorentzVector;
1946 myJet->GetMom(*temp_jet);
1948 // Loop over all tracks
1949 for (i=0;i<fnTracks;i++)
1951 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1952 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1954 TLorentzVector *track_vec = new TLorentzVector;
1955 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1956 if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
1958 E_tracks_total+=vtrack->Pt();
1964 // Loop over all caloclusters
1965 for (i=0;i<fnClusters;i++)
1967 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1968 TLorentzVector *cluster_vec = new TLorentzVector;
1969 vcluster->GetMomentum(*cluster_vec,fvertex);
1970 if (temp_jet->DeltaR(*cluster_vec)>fJetRForRho)
1972 E_caloclusters_total+=vcluster->E();
1977 // Calculate the mean Background density
1978 EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
1980 else // i.e. No signal jets -> same as total background density
1982 // Loop over all tracks
1983 for (i=0;i<fnTracks;i++)
1985 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1986 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1988 E_tracks_total+=vtrack->Pt();
1992 // Loop over all caloclusters
1993 for (i=0;i<fnClusters;i++)
1995 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1996 E_caloclusters_total+=vcluster->E();
1998 // Calculate the mean Background density
1999 EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2003 fRhoFull1->FillRho(fEventCentrality,EMCal_rho);
2004 fRhoFull1->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2005 fRhoFull1->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
2006 fRhoFull1->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
2007 fRhoFull1->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
2008 fRhoFull1->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
2011 void AliAnalysisTaskFullpAJets::EstimateFullRho2()
2014 Double_t E_tracks_total=0.0;
2015 Double_t E_caloclusters_total=0.0;
2016 Double_t EMCal_rho=0.0;
2018 if ((fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold) && (fEMCalPartJet->GetSubLeadingPt()>=fEMCalJetThreshold))
2020 AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
2021 TLorentzVector *temp_jet1 = new TLorentzVector;
2022 myhJet->GetMom(*temp_jet1);
2024 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSubLeadingIndex());
2025 TLorentzVector *temp_jet2 = new TLorentzVector;
2026 myJet->GetMom(*temp_jet2);
2028 // Loop over all tracks
2029 for (i=0;i<fnTracks;i++)
2031 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2032 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2034 TLorentzVector *track_vec = new TLorentzVector;
2035 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2036 if ((temp_jet1->DeltaR(*track_vec)>fJetRForRho) && (temp_jet2->DeltaR(*track_vec)>fJetRForRho))
2038 E_tracks_total+=vtrack->Pt();
2044 // Loop over all caloclusters
2045 for (i=0;i<fnClusters;i++)
2047 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2048 TLorentzVector *cluster_vec = new TLorentzVector;
2049 vcluster->GetMomentum(*cluster_vec,fvertex);
2050 if ((temp_jet1->DeltaR(*cluster_vec)>fJetRForRho) && (temp_jet2->DeltaR(*cluster_vec)>fJetRForRho))
2052 E_caloclusters_total+=vcluster->E();
2059 // Calculate the mean Background density
2060 EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myhJet->Phi(),myhJet->Eta())-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2062 else if (fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold)
2064 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
2065 TLorentzVector *temp_jet= new TLorentzVector;
2066 myJet->GetMom(*temp_jet);
2068 // Loop over all tracks
2069 for (i=0;i<fnTracks;i++)
2071 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2072 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2074 TLorentzVector *track_vec = new TLorentzVector;
2075 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2076 if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
2078 E_tracks_total+=vtrack->Pt();
2084 // Loop over all caloclusters
2085 for (i=0;i<fnClusters;i++)
2087 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2088 TLorentzVector *cluster_vec = new TLorentzVector;
2089 vcluster->GetMomentum(*cluster_vec,fvertex);
2090 if (temp_jet->DeltaR(*cluster_vec)>fJetRForRho)
2092 E_caloclusters_total+=vcluster->E();
2097 // Calculate the mean Background density
2098 EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2100 else // i.e. No signal jets -> same as total background density
2102 // Loop over all tracks
2103 for (i=0;i<fnTracks;i++)
2105 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2106 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2108 E_tracks_total+=vtrack->Pt();
2112 // Loop over all caloclusters
2113 for (i=0;i<fnClusters;i++)
2115 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2116 E_caloclusters_total+=vcluster->E();
2118 // Calculate the mean Background density
2119 EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2123 fRhoFull2->FillRho(fEventCentrality,EMCal_rho);
2124 fRhoFull2->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2125 fRhoFull2->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
2126 fRhoFull2->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
2127 fRhoFull2->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
2128 fRhoFull2->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
2131 void AliAnalysisTaskFullpAJets::EstimateFullRhoN()
2134 Bool_t track_away_from_jet;
2135 Bool_t cluster_away_from_jet;
2136 Double_t E_tracks_total=0.0;
2137 Double_t E_caloclusters_total=0.0;
2138 Double_t EMCal_rho=0.0;
2139 Double_t jet_area_total=0.0;
2141 // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
2142 for (i=0;i<fnTracks;i++)
2144 // First, check if track is in the EMCal!!
2145 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
2146 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2148 if (fEMCalPartJet->GetTotalSignalJets()<1)
2150 E_tracks_total+=vtrack->Pt();
2154 track_away_from_jet=kTRUE;
2156 TLorentzVector *track_vec = new TLorentzVector;
2157 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2158 while (track_away_from_jet==kTRUE && j<fEMCalPartJet->GetTotalSignalJets())
2160 TLorentzVector *jet_vec= new TLorentzVector;
2161 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSignalJetIndex(j));
2162 myJet->GetMom(*jet_vec);
2163 if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
2165 track_away_from_jet=kFALSE;
2170 if (track_away_from_jet==kTRUE)
2172 E_tracks_total+=vtrack->Pt();
2179 // Next, sum all CaloClusters within the EMCal (obviously all clusters must be within EMCal!!) that are away from jet(s) above Pt Threshold
2180 for (i=0;i<fnClusters;i++)
2182 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2183 if (fEMCalPartJet->GetTotalSignalJets()<1)
2185 E_caloclusters_total+=vcluster->E();
2189 cluster_away_from_jet=kTRUE;
2192 TLorentzVector *cluster_vec = new TLorentzVector;
2193 vcluster->GetMomentum(*cluster_vec,fvertex);
2194 while (cluster_away_from_jet==kTRUE && j<fEMCalPartJet->GetTotalSignalJets())
2196 TLorentzVector *jet_vec= new TLorentzVector;
2197 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSignalJetIndex(j));
2198 myJet->GetMom(*jet_vec);
2199 if (cluster_vec->DeltaR(*jet_vec)<=fJetRForRho)
2201 cluster_away_from_jet=kFALSE;
2206 if (cluster_away_from_jet==kTRUE)
2208 E_caloclusters_total+=vcluster->E();
2214 // Determine area of all Jets that are within the EMCal
2215 if (fEMCalPartJet->GetTotalSignalJets()==0)
2221 for (i=0;i<fEMCalPartJet->GetTotalSignalJets();i++)
2223 AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSignalJetIndex(i));
2224 jet_area_total+=AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta());
2229 EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-jet_area_total);
2232 fRhoFullN->FillRho(fEventCentrality,EMCal_rho);
2233 fRhoFullN->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2234 fRhoFullN->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
2235 fRhoFullN->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
2236 fRhoFullN->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
2237 fRhoFullN->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
2240 void AliAnalysisTaskFullpAJets::EstimateFullRhoDijet()
2243 Double_t E_tracks_total=0.0;
2244 Double_t E_caloclusters_total=0.0;
2245 Double_t EMCal_rho=0.0;
2247 // Loop over all tracks
2248 for (i=0;i<fnTracks;i++)
2250 AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2251 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2253 E_tracks_total+=vtrack->Pt();
2257 // Loop over all caloclusters
2258 for (i=0;i<fnClusters;i++)
2260 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2261 E_caloclusters_total+=vcluster->E();
2264 // Calculate the mean Background density
2265 EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2268 fRhoFullDijet->FillRho(fEventCentrality,EMCal_rho);
2269 fRhoFullDijet->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2270 fRhoFullDijet->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
2271 fRhoFullDijet->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
2272 fRhoFullDijet->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
2273 fRhoFullDijet->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
2276 void AliAnalysisTaskFullpAJets::EstimateFullRhokT()
2279 Double_t kTRho = 0.0;
2280 Double_t *pTArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2281 Double_t *RhoArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2283 for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2285 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
2286 pTArray[i]=myJet->Pt();
2287 RhoArray[i]=myJet->Pt()/myJet->Area();
2290 if (fEMCalkTFullJet->GetTotalJets()>0)
2292 kTRho=MedianRhokT(pTArray,RhoArray,fEMCalkTFullJet->GetTotalJets());
2298 fRhoFullkT->FillRho(fEventCentrality,kTRho);
2299 fRhoFullkT->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2300 fRhoFullkT->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
2301 fRhoFullkT->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
2302 fRhoFullkT->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
2303 fRhoFullkT->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
2308 void AliAnalysisTaskFullpAJets::EstimateFullRhoCMS()
2311 Double_t kTRho = 0.0;
2312 Double_t CMSTotalkTArea = 0.0;
2313 Double_t CMSParticleArea = 0.0;
2314 Double_t CMSCorrectionFactor = 1.0;
2315 Double_t *pTArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2316 Double_t *RhoArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2319 if ((fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold) && (fEMCalPartJet->GetSubLeadingPt()>=fEMCalJetThreshold))
2321 AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
2322 AliEmcalJet *myJet2 =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSubLeadingIndex());
2324 for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2326 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
2328 CMSTotalkTArea+=myJet->Area();
2329 if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2331 CMSParticleArea+=myJet->Area();
2333 if (IsJetOverlap(myJet,myJet1,kTRUE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kTRUE)
2335 pTArray[k]=myJet->Pt();
2336 RhoArray[k]=myJet->Pt()/myJet->Area();
2342 kTRho=MedianRhokT(pTArray,RhoArray,k);
2349 else if (fEMCalJet->GetLeadingPt()>=fEMCalJetThreshold)
2351 AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalJet->GetLeadingIndex());
2353 for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2355 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
2357 CMSTotalkTArea+=myJet->Area();
2358 if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2360 CMSParticleArea+=myJet->Area();
2362 if (IsJetOverlap(myJet,myJet1,kTRUE)==kFALSE)
2364 pTArray[k]=myJet->Pt();
2365 RhoArray[k]=myJet->Pt()/myJet->Area();
2371 kTRho=MedianRhokT(pTArray,RhoArray,k);
2380 for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2382 AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
2384 CMSTotalkTArea+=myJet->Area();
2385 if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2387 CMSParticleArea+=myJet->Area();
2389 pTArray[k]=myJet->Pt();
2390 RhoArray[k]=myJet->Pt()/myJet->Area();
2395 kTRho=MedianRhokT(pTArray,RhoArray,k);
2402 // Scale CMS Rho by Correction factor
2403 if (CMSTotalkTArea==0.0)
2405 CMSCorrectionFactor = 1.0;
2409 //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
2410 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
2412 kTRho*=CMSCorrectionFactor;
2414 fRhoFullCMS->FillRho(fEventCentrality,kTRho);
2415 fRhoFullCMS->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
2416 fRhoFullCMS->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
2417 fRhoFullCMS->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
2418 fRhoFullCMS->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
2419 fRhoFullCMS->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
2424 void AliAnalysisTaskFullpAJets::JetPtChargedProfile()
2428 Double_t ED_pT[fEDProfileRBins];
2430 for (i=0;i<fTPCFullJet->GetTotalSignalJets();i++)
2432 AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCFullJet->GetSignalJetIndex(i));
2433 if (InsideRect(myJet->Phi(),fTPCPhiMin,fTPCPhiMax,myJet->Eta(),fTPCEtaMin+fEDProfileRUp,fTPCEtaMax-fEDProfileRUp)==kTRUE)
2435 for (j=0;j<fEDProfileRBins;j++)
2439 TLorentzVector *jet_vec= new TLorentzVector;
2440 myJet->GetMom(*jet_vec);
2441 // Sum all tracks in concentric rings around jet vertex
2442 for (j=0;j<fnTracks;j++)
2444 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2445 TLorentzVector *track_vec = new TLorentzVector;
2446 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2447 delta_R=jet_vec->DeltaR(*track_vec);
2448 if (delta_R<=fEDProfileRUp)
2450 ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vtrack->Pt();
2455 for (j=0;j<fEDProfileRBins;j++)
2457 ED_pT[j]/=TMath::Pi()*TMath::Power((fEDProfileRUp/fEDProfileRBins),2)*(2*j+1);
2458 fpChargedJetEDProfile[TMath::FloorNint(fEventCentrality/10.)]->Fill(myJet->Pt(),myJet->Eta(),(fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2459 if (fEventCentrality<=20)
2461 fpChargedJetRProfile[4+TMath::FloorNint(myJet->Eta()*10.)]->Fill((fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2469 void AliAnalysisTaskFullpAJets::JetPtFullProfile()
2473 Double_t ED_pT[fEDProfileRBins];
2475 for (i=0;i<fEMCalFullJet->GetTotalSignalJets();i++)
2477 AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalFullJet->GetSignalJetIndex(i));
2478 if (InsideRect(myJet->Phi(),fEMCalPhiMin+fEDProfileRUp,fEMCalPhiMax-fEDProfileRUp,myJet->Eta(),fEMCalEtaMin+fEDProfileRUp,fEMCalEtaMax-fEDProfileRUp)==kTRUE)
2480 for (j=0;j<fEDProfileRBins;j++)
2484 TLorentzVector *jet_vec= new TLorentzVector;
2485 myJet->GetMom(*jet_vec);
2486 // Sum all tracks in concentric rings around jet vertex
2487 for (j=0;j<fnTracks;j++)
2489 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2490 TLorentzVector *track_vec = new TLorentzVector;
2491 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2492 delta_R=jet_vec->DeltaR(*track_vec);
2493 if (delta_R<=fEDProfileRUp)
2495 ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vtrack->Pt();
2500 // Sum all clusters in concentric rings around jet vertex
2501 for (j=0;j<fnClusters;j++)
2503 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
2504 TLorentzVector *cluster_vec = new TLorentzVector;
2505 vcluster->GetMomentum(*cluster_vec,fvertex);
2506 delta_R=jet_vec->DeltaR(*cluster_vec);
2507 if (delta_R<=fEDProfileRUp)
2509 ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vcluster->E();
2514 for (j=0;j<fEDProfileRBins;j++)
2516 ED_pT[j]/=TMath::Pi()*TMath::Power((fEDProfileRUp/fEDProfileRBins),2)*(2*j+1);
2517 fpJetEDProfile[TMath::FloorNint(fEventCentrality/10.)]->Fill(myJet->Pt(),myJet->Eta(),(fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2518 // Fill profile if a "most" central event (0-20%)
2519 if (fEventCentrality<=20)
2521 fpJetRProfile[2+TMath::FloorNint(myJet->Eta()*10.)]->Fill((fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
2526 // Fill constituent histogram
2527 for (j=0;j<myJet->GetNumberOfTracks();j++)
2529 AliVTrack* vtrack = (AliVTrack*) fOrgTracks->At(myJet->TrackAt(j));
2530 fhJetConstituentPt->Fill(myJet->Pt(),vtrack->Pt());
2533 for (j=0;j<myJet->GetNumberOfClusters();j++)
2535 AliVCluster* vcluster = (AliVCluster*) fOrgClusters->At(myJet->ClusterAt(j));
2536 TLorentzVector *cluster_vec = new TLorentzVector;
2537 vcluster->GetMomentum(*cluster_vec,fvertex);
2538 fhJetConstituentPt->Fill(myJet->Pt(),cluster_vec->Pt());
2545 void AliAnalysisTaskFullpAJets::JetPtEtaProfile()
2550 Double_t Eta_pT[fEtaProfileBins];
2551 Double_t Eta_abs_pT[Int_t(0.5*fEtaProfileBins)];
2553 for (i=0;i<fEMCalFullJet->GetTotalSignalJets();i++)
2555 AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalFullJet->GetSignalJetIndex(i));
2556 if (IsInEMCal(myJet->Phi(),myJet->Eta())==kTRUE)
2558 for (j=0;j<fEtaProfileBins;j++)
2564 // Sum all tracks in strips of eta away from the jet vertex
2565 for (j=0;j<fnTracks;j++)
2567 AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2569 delta_eta=TMath::Abs(vtrack->Eta()-myJet->Eta());
2570 if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2572 Eta_pT[Int_t(0.5*fEtaProfileBins)+TMath::FloorNint(10*eta)]+=vtrack->Pt();
2573 Eta_abs_pT[TMath::FloorNint(10*delta_eta)]+=vtrack->Pt();
2577 // Sum all clusters in strips of eta away from the jet vertex
2578 for (j=0;j<fnClusters;j++)
2580 AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
2581 TLorentzVector *cluster_vec = new TLorentzVector;
2582 vcluster->GetMomentum(*cluster_vec,fvertex);
2583 eta=cluster_vec->Eta();
2584 delta_eta=TMath::Abs(cluster_vec->Eta()-myJet->Eta());
2585 Eta_pT[Int_t(0.5*fEtaProfileBins)+TMath::FloorNint(10*eta)]+=vcluster->E();
2586 Eta_abs_pT[TMath::FloorNint(10*delta_eta)]+=vcluster->E();
2590 for (j=0;j<fEtaProfileBins;j++)
2592 Eta_pT[j]/=0.1*fEMCalPhiTotal;
2593 // Fill profile if a "most" central event (0-20%)
2594 if (j<(10*(fEMCalEtaMax-TMath::Abs(myJet->Eta()))))
2596 Eta_abs_pT[j]/=0.2*fEMCalPhiTotal;
2600 Eta_abs_pT[j]/=0.1*fEMCalPhiTotal;
2602 // Fill profile if a "most" central event (0-20%)
2603 if (fEventCentrality<=20)
2605 fpJetAbsEtaProfile[7+TMath::FloorNint(myJet->Eta()*10.)]->Fill(0.1*j,Eta_abs_pT[j]);
2606 fpJetEtaProfile[7+TMath::FloorNint(myJet->Eta()*10.)]->Fill(0.1*(j-7),Eta_pT[j]);
2613 void AliAnalysisTaskFullpAJets::DeleteJetData(Bool_t EMCalOn)
2619 delete fTPCkTFullJet;
2624 delete fEMCalFullJet;
2625 delete fEMCalPartJet;
2626 delete fEMCalkTFullJet;
2630 /////////////////////////////////////////////////////////////////////////////////////////
2631 ///////////////// User Defined Functions ///////////////////////////////////////
2632 /////////////////////////////////////////////////////////////////////////////////////////
2634 Bool_t AliAnalysisTaskFullpAJets::IsDiJetEvent()
2636 // Determine if event contains a di-jet within the detector. Uses charged jets.
2637 // Requires the delta phi of the jets to be 180 +/- 15 degrees.
2638 // Requires both jets to be outside of the EMCal
2639 // Requires both jets to be signal jets
2641 const Double_t dijet_delta_phi=(180/360.)*2*TMath::Pi();
2642 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
2643 Double_t dummy_phi=0.0;
2645 if (fTPCOnlyJet->GetTotalSignalJets()>1)
2647 AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCOnlyJet->GetLeadingIndex());
2648 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCOnlyJet->GetSubLeadingIndex());
2649 dummy_phi=TMath::Min(TMath::Abs(myhJet->Phi()-myJet->Phi()),2*TMath::Pi()-TMath::Abs(myhJet->Phi()-myJet->Phi()));
2650 if (dummy_phi>(dijet_delta_phi-dijet_phi_acceptance))
2659 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)
2661 if (phi>phi_min && phi<phi_max)
2663 if (eta>eta_min && eta<eta_max)
2671 Bool_t AliAnalysisTaskFullpAJets::IsInEMCal(Double_t phi,Double_t eta)
2673 return InsideRect(phi,fEMCalPhiMin,fEMCalPhiMax,eta,fEMCalEtaMin,fEMCalEtaMax);
2676 Bool_t AliAnalysisTaskFullpAJets::IsInEMCalFull(Double_t r,Double_t phi,Double_t eta)
2678 return InsideRect(phi,fEMCalPhiMin+r,fEMCalPhiMax-r,eta,fEMCalEtaMin+r,fEMCalEtaMax-r);
2681 Bool_t AliAnalysisTaskFullpAJets::IsInEMCalPart(Double_t r,Double_t phi,Double_t eta)
2683 return InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
2686 Bool_t AliAnalysisTaskFullpAJets::IsInTPCFull(Double_t r,Double_t phi,Double_t eta)
2688 Bool_t in_EMCal= InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
2689 Bool_t in_TPC= InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
2691 if (in_EMCal==kFALSE && in_TPC==kTRUE)
2698 Bool_t AliAnalysisTaskFullpAJets::IsInTPC(Double_t r,Double_t phi,Double_t eta,Bool_t Complete)
2700 if (Complete==kTRUE)
2702 return InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
2704 return InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin,fTPCEtaMax);
2707 Bool_t AliAnalysisTaskFullpAJets::IsJetOverlap(AliEmcalJet *jet1,AliEmcalJet *jet2,Bool_t EMCalOn)
2712 Int_t jetCluster1=0;
2713 Int_t jetCluster2=0;
2715 for (i=0;i<jet1->GetNumberOfTracks();i++)
2717 jetTrack1=jet1->TrackAt(i);
2718 for (j=0;j<jet2->GetNumberOfTracks();j++)
2720 jetTrack2=jet2->TrackAt(j);
2721 if (jetTrack1 == jetTrack2)
2727 if (EMCalOn == kTRUE)
2729 for (i=0;i<jet1->GetNumberOfClusters();i++)
2731 jetCluster1=jet1->ClusterAt(i);
2732 for (j=0;j<jet2->GetNumberOfClusters();j++)
2734 jetCluster2=jet2->ClusterAt(j);
2735 if (jetCluster1 == jetCluster2)
2745 Double_t AliAnalysisTaskFullpAJets::AreaWithinTPC(Double_t r,Double_t eta)
2748 if (eta<(fTPCEtaMin+r))
2752 else if(eta>(fTPCEtaMax-r))
2760 return r*r*TMath::Pi()-AreaEdge(r,z);
2763 Double_t AliAnalysisTaskFullpAJets::AreaWithinEMCal(Double_t r,Double_t phi,Double_t eta)
2767 if (phi<(fEMCalPhiMin-r) || phi>(fEMCalPhiMax+r))
2771 else if (phi<(fEMCalPhiMin+r))
2775 else if (phi>(fEMCalPhiMin+r) && phi<(fEMCalPhiMax-r))
2784 if (eta<(fEMCalEtaMin-r) || eta>(fEMCalEtaMax+r))
2788 else if (eta<(fEMCalEtaMin+r))
2792 else if (eta>(fEMCalEtaMin+r) && eta<(fEMCalEtaMax-r))
2803 if (TMath::Sqrt(x*x+y*y)>=r)
2805 return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
2807 return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y)+AreaOverlap(r,x,y);
2809 else if ((x>=r && y<0) || (y>=r && x<0))
2811 return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
2813 else if (x>0 && x<r && y<0)
2815 Double_t a=TMath::Sqrt(r*r-x*x);
2816 Double_t b=TMath::Sqrt(r*r-y*y);
2819 return r*r*TMath::ASin(b/r)+y*b;
2823 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);
2826 else if (y>0 && y<r && x<0)
2828 Double_t a=TMath::Sqrt(r*r-x*x);
2829 Double_t b=TMath::Sqrt(r*r-y*y);
2832 return r*r*TMath::ASin(a/r)+x*a;
2836 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);
2841 Double_t a=TMath::Sqrt(r*r-x*x);
2842 Double_t b=TMath::Sqrt(r*r-y*y);
2849 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);
2854 Double_t AliAnalysisTaskFullpAJets::AreaEdge(Double_t r,Double_t z)
2856 Double_t a=TMath::Sqrt(r*r-z*z);
2857 return r*r*TMath::ASin(a/r)-a*z;
2860 Double_t AliAnalysisTaskFullpAJets::AreaOverlap(Double_t r,Double_t x,Double_t y)
2862 Double_t a=TMath::Sqrt(r*r-x*x);
2863 Double_t b=TMath::Sqrt(r*r-y*y);
2864 return x*y-0.5*(x*a+y*b)+0.5*r*r*(TMath::ASin(b/r)-TMath::ASin(x/r));
2867 Double_t AliAnalysisTaskFullpAJets::TransverseArea(Double_t r,Double_t psi0,Double_t phi,Double_t eta)
2869 Double_t area_left=0;
2870 Double_t area_right=0;
2874 Double_t eta_down=0;
2876 Double_t u=eta-fEMCalEtaMin;
2877 Double_t v=fEMCalEtaMax-eta;
2879 Double_t phi1=phi+u*TMath::Tan(psi0);
2880 Double_t phi2=phi-u*TMath::Tan(psi0);
2881 Double_t phi3=phi+v*TMath::Tan(psi0);
2882 Double_t phi4=phi-v*TMath::Tan(psi0);
2884 //Calculate the Left side area
2885 if (phi1>=fEMCalPhiMax)
2887 eta_a=eta-u*((fEMCalPhiMax-phi)/(phi1-phi));
2889 if (phi2<=fEMCalPhiMin)
2891 eta_b=eta-u*((phi-fEMCalPhiMin)/(phi-phi2));
2894 if ((phi1>=fEMCalPhiMax) && (phi2<=fEMCalPhiMin))
2896 eta_up=TMath::Max(eta_a,eta_b);
2897 eta_down=TMath::Min(eta_a,eta_b);
2899 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);
2901 else if (phi1>=fEMCalPhiMax)
2903 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);
2905 else if (phi2<=fEMCalPhiMin)
2907 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);
2911 area_left=0.5*(phi1-phi2+2*r*TMath::Tan(psi0))*(u-r);
2914 // Calculate the Right side area
2915 if (phi3>=fEMCalPhiMax)
2917 eta_a=eta+v*((fEMCalPhiMax-phi)/(phi3-phi));
2919 if (phi4<=fEMCalPhiMin)
2921 eta_b=eta+v*((phi-fEMCalPhiMin)/(phi-phi4));
2924 if ((phi3>=fEMCalPhiMax) && (phi4<=fEMCalPhiMin))
2926 eta_up=TMath::Max(eta_a,eta_b);
2927 eta_down=TMath::Min(eta_a,eta_b);
2929 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);
2931 else if (phi3>=fEMCalPhiMax)
2933 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);
2935 else if (phi4<=fEMCalPhiMin)
2937 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);
2941 area_right=0.5*(phi3-phi4+2*r*TMath::Tan(psi0))*(v-r);
2943 return area_left+area_right;
2946 Double_t AliAnalysisTaskFullpAJets::MedianRhokT(Double_t *pTkTEntries, Double_t *RhokTEntries, Int_t nEntries)
2948 // This function is used to calculate the median Rho kT value. The procedure is:
2949 // - Order the kT cluster array from highest rho value to lowest
2950 // - Exclude highest rho kT cluster
2951 // - Return the median rho value of the remaining subset
2954 const Double_t rho_min=-9.9999E+99;
2956 Double_t w[nEntries]; // Used for sorting
2957 Double_t smax=rho_min;
2963 for (j=0;j<nEntries;j++)
2968 for (j=0;j<nEntries;j++)
2970 if (pTkTEntries[j]>pTmax)
2972 pTmax=pTkTEntries[j];
2977 for (j=0;j<nEntries;j++)
2979 for (k=0;k<nEntries;k++)
2981 if (RhokTEntries[k]>smax)
2983 smax=RhokTEntries[k];
2988 RhokTEntries[sindex]=rho_min;
2992 return w[nEntries/2];
2996 // AlipAJetData Class Member Defs
2998 AliAnalysisTaskFullpAJets::AlipAJetData::AlipAJetData()
3001 // Dummy constructor ALWAYS needed for I/O.
3004 AliAnalysisTaskFullpAJets::AlipAJetData::AlipAJetData(const char *name, Bool_t isFull, Int_t nEntries)
3007 SetIsJetsFull(isFull);
3008 SetTotalEntries(nEntries);
3009 SetLeading(0,-9.99E+099);
3010 SetSubLeading(0,-9.99E+099);
3012 SetAreaCutFraction(0.6);
3017 AliAnalysisTaskFullpAJets::AlipAJetData::~AlipAJetData()
3022 SetIsJetsFull(kFALSE);
3025 SetTotalSignalJets(0);
3029 SetAreaCutFraction(0);
3032 delete [] fJetsIndex;
3033 delete [] fJetsSCIndex;
3034 delete [] fIsJetInArray;
3038 // User Defined Sub-Routines
3039 void AliAnalysisTaskFullpAJets::AlipAJetData::InitializeJetData(TClonesArray *jetList, Int_t nEntries)
3044 Double_t AreaThreshold = fAreaCutFrac*TMath::Pi()*TMath::Power(fJetR,2);
3046 // Initialize Jet Data
3047 for (i=0;i<nEntries;i++)
3049 AliEmcalJet *myJet =(AliEmcalJet*) jetList->At(i);
3051 if (fIsJetInArray[i]==kTRUE && myJet->Area()>AreaThreshold)
3054 if (myJet->Pt()>fPtMax)
3056 SetSubLeading(fPtMaxIndex,fPtMax);
3057 SetLeading(i,myJet->Pt());
3059 else if (myJet->Pt()>fPtSubLeading)
3061 SetSubLeading(i,myJet->Pt());
3063 if (myJet->Pt()>=fSignalPt)
3065 SetSignalJetIndex(i,l);
3072 SetTotalSignalJets(l);
3076 void AliAnalysisTaskFullpAJets::AlipAJetData::SetName(const char *name)
3081 void AliAnalysisTaskFullpAJets::AlipAJetData::SetIsJetsFull(Bool_t isFull)
3083 fIsJetsFull = isFull;
3086 void AliAnalysisTaskFullpAJets::AlipAJetData::SetTotalEntries(Int_t nEntries)
3089 fJetsIndex = new Int_t[fnTotal];
3090 fJetsSCIndex = new Int_t[fnTotal];
3091 fIsJetInArray = new Bool_t[fnTotal];
3094 void AliAnalysisTaskFullpAJets::AlipAJetData::SetTotalJets(Int_t nJets)
3099 void AliAnalysisTaskFullpAJets::AlipAJetData::SetTotalSignalJets(Int_t nSignalJets)
3101 fnJetsSC = nSignalJets;
3104 void AliAnalysisTaskFullpAJets::AlipAJetData::SetSignalCut(Double_t Pt)
3109 void AliAnalysisTaskFullpAJets::AlipAJetData::SetLeading(Int_t index, Double_t Pt)
3111 fPtMaxIndex = index;
3115 void AliAnalysisTaskFullpAJets::AlipAJetData::SetSubLeading(Int_t index, Double_t Pt)
3117 fPtSubLeadingIndex = index;
3121 void AliAnalysisTaskFullpAJets::AlipAJetData::SetJetIndex(Int_t index, Int_t At)
3123 fJetsIndex[At] = index;
3126 void AliAnalysisTaskFullpAJets::AlipAJetData::SetSignalJetIndex(Int_t index, Int_t At)
3128 fJetsSCIndex[At] = index;
3131 void AliAnalysisTaskFullpAJets::AlipAJetData::SetIsJetInArray(Bool_t isInArray, Int_t At)
3133 fIsJetInArray[At] = isInArray;
3136 void AliAnalysisTaskFullpAJets::AlipAJetData::SetAreaCutFraction(Double_t areaFraction)
3138 fAreaCutFrac = areaFraction;
3141 void AliAnalysisTaskFullpAJets::AlipAJetData::SetJetR(Double_t jetR)
3147 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetTotalEntries()
3152 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetTotalJets()
3157 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetTotalSignalJets()
3162 Double_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSignalCut()
3167 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetLeadingIndex()
3172 Double_t AliAnalysisTaskFullpAJets::AlipAJetData::GetLeadingPt()
3177 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSubLeadingIndex()
3179 return fPtSubLeadingIndex;
3182 Double_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSubLeadingPt()
3184 return fPtSubLeading;
3187 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetJetIndex(Int_t At)
3189 return fJetsIndex[At];
3192 Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSignalJetIndex(Int_t At)
3194 return fJetsSCIndex[At];
3197 Bool_t AliAnalysisTaskFullpAJets::AlipAJetData::GetIsJetInArray(Int_t At)
3199 return fIsJetInArray[At];
3203 // AlipAJetHistos Class Member Defs
3205 AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos()
3207 // Dummy constructor ALWAYS needed for I/O.
3210 AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name)
3213 SetCentralityTag("V0A");
3214 SetCentralityRange(100,0,100);
3215 SetPtRange(250,-50,200);
3216 SetRhoPtRange(500,0,50);
3217 SetDeltaPtRange(200,-100,100);
3218 SetBackgroundFluctuationsPtRange(100,0,100);
3219 SetLeadingJetPtRange(200,0,200);
3224 AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, const char *centag)
3227 SetCentralityTag(centag);
3228 SetCentralityRange(100,0,100);
3229 SetPtRange(250,-50,200);
3230 SetRhoPtRange(500,0,50);
3231 SetDeltaPtRange(200,-100,100);
3232 SetBackgroundFluctuationsPtRange(100,0,100);
3233 SetLeadingJetPtRange(200,0,200);
3239 AliAnalysisTaskFullpAJets::AlipAJetHistos::~AlipAJetHistos()
3247 void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
3249 fOutput = new TList();
3250 fOutput->SetOwner();
3251 fOutput->SetName(fName);
3253 TString RhoString="";
3254 TString PtString="";
3255 TString DeltaPtString="";
3256 TString BckgFlucPtString="";
3257 TString CentralityString;
3258 CentralityString = Form("Centrality (%s)",fCentralityTag);
3260 // Rho Spectral Plots
3261 RhoString = Form("%d-%d Centrality, Rho Spectrum",0,20);
3262 fh020Rho = new TH1D("fh020Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
3263 fh020Rho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
3264 fh020Rho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
3267 RhoString = Form("%d-%d Centrality, Rho Spectrum",80,100);
3268 fh80100Rho = new TH1D("fh80100Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
3269 fh80100Rho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
3270 fh80100Rho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
3271 fh80100Rho->Sumw2();
3273 RhoString = Form("%d-%d Centrality, Rho Spectrum",0,100);
3274 fhRho = new TH1D("fhRho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
3275 fhRho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
3276 fhRho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
3279 RhoString = "Rho Spectrum vs Centrality";
3280 fhRhoCen = new TH2D("fhRhoCen",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3281 fhRhoCen->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
3282 fhRhoCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3283 fhRhoCen->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
3286 // Background Subtracted Plots
3287 PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,20);
3288 fh020BSPt = new TH1D("fh020BSPt",PtString,fPtBins,fPtLow,fPtUp);
3289 fh020BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3290 fh020BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3293 PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",80,100);
3294 fh80100BSPt = new TH1D("fh80100BSPt",PtString,fPtBins,fPtLow,fPtUp);
3295 fh80100BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3296 fh80100BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3297 fh80100BSPt->Sumw2();
3299 PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,100);
3300 fhBSPt = new TH1D("fhBSPt",PtString,fPtBins,fPtLow,fPtUp);
3301 fhBSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3302 fhBSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3305 PtString = "Background Subtracted Jet Spectrum vs Centrality";
3306 fhBSPtCen = new TH2D("fhBSPtCen",PtString,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3307 fhBSPtCen->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3308 fhBSPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3309 fhBSPtCen->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3312 PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,20);
3313 fh020BSPtSignal = new TH1D("fh020BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
3314 fh020BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3315 fh020BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3316 fh020BSPtSignal->Sumw2();
3318 PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",80,100);
3319 fh80100BSPtSignal = new TH1D("fh80100BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
3320 fh80100BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3321 fh80100BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3322 fh80100BSPtSignal->Sumw2();
3324 PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,100);
3325 fhBSPtSignal = new TH1D("fhBSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
3326 fhBSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3327 fhBSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3328 fhBSPtSignal->Sumw2();
3330 PtString = "Background Subtracted Signal Jet Spectrum vs Centrality";
3331 fhBSPtCenSignal = new TH2D("fhBSPtCenSignal",PtString,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3332 fhBSPtCenSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
3333 fhBSPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3334 fhBSPtCenSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3335 fhBSPtCenSignal->Sumw2();
3337 // Delta Pt Plots with RC at least 2R away from Leading Signal
3338 DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
3339 fh020DeltaPt = new TH1D("fh020DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3340 fh020DeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3341 fh020DeltaPt->GetYaxis()->SetTitle("Probability Density");
3342 fh020DeltaPt->Sumw2();
3344 DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
3345 fh80100DeltaPt = new TH1D("fh80100DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3346 fh80100DeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3347 fh80100DeltaPt->GetYaxis()->SetTitle("Probability Density");
3348 fh80100DeltaPt->Sumw2();
3350 DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
3351 fhDeltaPt = new TH1D("fhDeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3352 fhDeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3353 fhDeltaPt->GetYaxis()->SetTitle("Probability Density");
3356 DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
3357 fhDeltaPtCen = new TH2D("fhDeltaPtCen",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3358 fhDeltaPtCen->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3359 fhDeltaPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3360 fhDeltaPtCen->GetZaxis()->SetTitle("Probability Density");
3361 fhDeltaPtCen->Sumw2();
3363 // Delta Pt Plots with no spatial restrictions on RC
3364 DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
3365 fh020DeltaPtSignal = new TH1D("fh020DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3366 fh020DeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3367 fh020DeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
3368 fh020DeltaPtSignal->Sumw2();
3370 DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
3371 fh80100DeltaPtSignal = new TH1D("fh80100DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3372 fh80100DeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3373 fh80100DeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
3374 fh80100DeltaPtSignal->Sumw2();
3376 DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
3377 fhDeltaPtSignal = new TH1D("fhDeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
3378 fhDeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3379 fhDeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
3380 fhDeltaPtSignal->Sumw2();
3382 DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
3383 fhDeltaPtCenSignal = new TH2D("fhDeltaPtCenSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3384 fhDeltaPtCenSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
3385 fhDeltaPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3386 fhDeltaPtCenSignal->GetZaxis()->SetTitle("Probability Density");
3387 fhDeltaPtCenSignal->Sumw2();
3389 // Background Fluctuations Pt Plots
3390 BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,20);
3391 fh020BckgFlucPt = new TH1D("fh020BckgFlucPt",PtString,fPtBins,fPtLow,fPtUp);
3392 fh020BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
3393 fh020BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3394 fh020BckgFlucPt->Sumw2();
3396 BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",80,100);
3397 fh80100BckgFlucPt = new TH1D("fh80100BckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
3398 fh80100BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
3399 fh80100BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3400 fh80100BckgFlucPt->Sumw2();
3402 BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,100);
3403 fhBckgFlucPt = new TH1D("fhBckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
3404 fhBckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
3405 fhBckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3406 fhBckgFlucPt->Sumw2();
3408 BckgFlucPtString = "Background Fluctuation p_{T} Spectrum vs Centrality";
3409 fhBckgFlucPtCen = new TH2D("fhBckgFlucPtCen",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
3410 fhBckgFlucPtCen->GetXaxis()->SetTitle("#p_{T} (GeV/c)");
3411 fhBckgFlucPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
3412 fhBckgFlucPtCen->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
3413 fhBckgFlucPtCen->Sumw2();
3415 // Background Density vs Centrality Profile
3416 RhoString = "Background Density vs Centrality";
3417 fpRho = new TProfile("fpRho",RhoString,fCentralityBins,fCentralityLow,fCentralityUp);
3418 fpRho->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
3419 fpRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
3421 // Background Density vs Leading Jet Profile
3422 fpLJetRho = new TProfile("fpLJetRho","#rho vs Leading Jet p_{T}",fLJetPtBins,fLJetPtLow,fLJetPtUp);
3423 fpLJetRho->GetXaxis()->SetTitle("Leading Jet p_{T}");
3424 fpLJetRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
3426 // Add Histos & Profiles to List
3427 fOutput->Add(fh020Rho);
3428 fOutput->Add(fh80100Rho);
3429 fOutput->Add(fhRho);
3430 fOutput->Add(fhRhoCen);
3431 fOutput->Add(fh020BSPt);
3432 fOutput->Add(fh80100BSPt);
3433 fOutput->Add(fhBSPt);
3434 fOutput->Add(fhBSPtCen);
3435 fOutput->Add(fh020BSPtSignal);
3436 fOutput->Add(fh80100BSPtSignal);
3437 fOutput->Add(fhBSPtSignal);
3438 fOutput->Add(fhBSPtCenSignal);
3439 fOutput->Add(fh020DeltaPt);
3440 fOutput->Add(fh80100DeltaPt);
3441 fOutput->Add(fhDeltaPt);
3442 fOutput->Add(fhDeltaPtCen);
3443 fOutput->Add(fh020DeltaPtSignal);
3444 fOutput->Add(fh80100DeltaPtSignal);
3445 fOutput->Add(fhDeltaPtSignal);
3446 fOutput->Add(fhDeltaPtCenSignal);
3447 fOutput->Add(fh020BckgFlucPt);
3448 fOutput->Add(fh80100BckgFlucPt);
3449 fOutput->Add(fhBckgFlucPt);
3450 fOutput->Add(fhBckgFlucPtCen);
3451 fOutput->Add(fpRho);
3452 fOutput->Add(fpLJetRho);
3455 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetName(const char *name)
3460 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetCentralityTag(const char *name)
3462 fCentralityTag = name;
3465 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetCentralityRange(Int_t bins, Double_t low, Double_t up)
3467 fCentralityBins=bins;
3472 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetPtRange(Int_t bins, Double_t low, Double_t up)
3479 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetRhoPtRange(Int_t bins, Double_t low, Double_t up)
3486 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetDeltaPtRange(Int_t bins, Double_t low, Double_t up)
3493 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetBackgroundFluctuationsPtRange(Int_t bins, Double_t low, Double_t up)
3495 fBckgFlucPtBins=bins;
3500 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetLeadingJetPtRange(Int_t bins, Double_t low, Double_t up)
3507 TList* AliAnalysisTaskFullpAJets::AlipAJetHistos::GetOutputHistos()
3512 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillRho(Double_t eventCentrality, Double_t rho)
3515 fhRhoCen->Fill(rho,eventCentrality);
3516 fpRho->Fill(eventCentrality,rho);
3518 if (eventCentrality<=20)
3520 fh020Rho->Fill(rho);
3522 else if (eventCentrality>=80)
3524 fh80100Rho->Fill(rho);
3528 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillBSJS(Double_t eventCentrality, Double_t rho, Double_t signalCut, TClonesArray *jetList, Int_t *indexJetList, Int_t nIndexJetList)
3531 Double_t tempPt=0.0;
3533 for (i=0;i<nIndexJetList;i++)
3535 AliEmcalJet *myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
3536 tempPt=myJet->Pt()-rho*myJet->Area();
3538 fhBSPt->Fill(tempPt);
3539 fhBSPtCen->Fill(tempPt,eventCentrality);
3540 if (eventCentrality<=20)
3542 fh020BSPt->Fill(tempPt);
3544 else if (eventCentrality>=80)
3546 fh80100BSPt->Fill(tempPt);
3549 if (myJet->Pt()>=signalCut)
3551 fhBSPtSignal->Fill(tempPt);
3552 fhBSPtCenSignal->Fill(tempPt,eventCentrality);
3553 if (eventCentrality<=20)
3555 fh020BSPtSignal->Fill(tempPt);
3557 else if (eventCentrality>=80)
3559 fh80100BSPtSignal->Fill(tempPt);
3566 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillDeltaPt(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
3569 Double_t tempPt=0.0;
3573 tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
3574 fhDeltaPt->Fill(tempPt);
3575 fhDeltaPtCen->Fill(tempPt,eventCentrality);
3576 if (eventCentrality<=20)
3578 fh020DeltaPt->Fill(tempPt);
3580 else if (eventCentrality>=80)
3582 fh80100DeltaPt->Fill(tempPt);
3588 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillDeltaPtSignal(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
3591 Double_t tempPt=0.0;
3595 tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
3596 fhDeltaPtSignal->Fill(tempPt);
3597 fhDeltaPtCenSignal->Fill(tempPt,eventCentrality);
3598 if (eventCentrality<=20)
3600 fh020DeltaPtSignal->Fill(tempPt);
3602 else if (eventCentrality>=80)
3604 fh80100DeltaPtSignal->Fill(tempPt);
3610 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillBackgroundFluctuations(Double_t eventCentrality, Double_t rho, Double_t jetRadius)
3612 Double_t tempPt=0.0;
3614 tempPt=rho*TMath::Power(jetRadius,2);
3615 fhBckgFlucPt->Fill(tempPt);
3616 fhBckgFlucPtCen->Fill(tempPt,eventCentrality);
3617 if (eventCentrality<=20)
3619 fh020BckgFlucPt->Fill(tempPt);
3621 else if (eventCentrality>=80)
3623 fh80100BckgFlucPt->Fill(tempPt);
3627 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillLeadingJetPtRho(Double_t jetPt, Double_t rho)
3629 fpLJetRho->Fill(jetPt,rho);