2 // Dijet base analysis task.
6 #include <TClonesArray.h>
10 #include <THnSparse.h>
12 #include <TLorentzVector.h>
19 #include "AliVCluster.h"
20 #include "AliVTrack.h"
21 #include "AliEmcalJet.h"
22 #include "AliRhoParameter.h"
24 #include "AliAnalysisUtils.h"
25 #include "AliEmcalParticle.h"
26 #include "AliMCEvent.h"
27 #include "AliGenPythiaEventHeader.h"
28 #include "AliAODMCHeader.h"
29 #include "AliMCEvent.h"
30 #include "AliAnalysisManager.h"
31 #include "AliJetContainer.h"
33 #include "AliAnalysisTaskEmcalDiJetBase.h"
35 ClassImp(AliAnalysisTaskEmcalDiJetBase)
37 //________________________________________________________________________
38 AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase() :
39 AliAnalysisTaskEmcalJetDev("AliAnalysisTaskEmcalDiJetBase", kTRUE),
41 fJetCorrelationType(kCorrelateAll),
42 fJetFullChargedMatchingType(kFraction),
46 fContainerChargedMC(3),
51 fDoChargedCharged(kTRUE),
52 fDoFullCharged(kTRUE),
56 fPtMinTriggerJet(10.),
57 fMinFractionShared(0.5),
58 fMatchingDone(kFALSE),
68 fHistTrialsSelEvents(0),
72 // Default constructor.
75 SetMakeGeneralHistograms(kTRUE);
79 //________________________________________________________________________
80 AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase(const char *name) :
81 AliAnalysisTaskEmcalJetDev(name, kTRUE),
83 fJetCorrelationType(kCorrelateAll),
84 fJetFullChargedMatchingType(kFraction),
88 fContainerChargedMC(3),
93 fDoChargedCharged(kTRUE),
94 fDoFullCharged(kTRUE),
98 fPtMinTriggerJet(10.),
99 fMinFractionShared(0.5),
100 fMatchingDone(kFALSE),
102 faFullFracIndexMC(0),
110 fHistTrialsSelEvents(0),
114 // Standard constructor.
116 SetMakeGeneralHistograms(kTRUE);
120 //________________________________________________________________________
121 AliAnalysisTaskEmcalDiJetBase::~AliAnalysisTaskEmcalDiJetBase()
126 //----------------------------------------------------------------------------------------------
127 void AliAnalysisTaskEmcalDiJetBase::InitOnce() {
129 // only initialize once
132 // Initialize analysis util class for vertex selection
134 fAnalysisUtils = new AliAnalysisUtils();
135 fAnalysisUtils->SetMinVtxContr(2);
136 fAnalysisUtils->SetMaxVtxZ(10.);
140 //________________________________________________________________________
141 Bool_t AliAnalysisTaskEmcalDiJetBase::SelectEvent() {
143 // Decide if event should be selected for analysis
146 fhNEvents->Fill(0.5);
148 if(!fAnalysisUtils && fUseAnaUtils)
153 if(!fAnalysisUtils->IsVertexSelected2013pA(InputEvent()))
156 fhNEvents->Fill(2.5);
158 if(fAnalysisUtils->IsPileUpEvent(InputEvent()))
163 AliError("fAnalysisUtils not initialized. Call AliAnalysisTaskEmcalJetTriggerQA::InitOnce()");
166 fhNEvents->Fill(3.5);
168 if(!fTriggerClass.IsNull()) {
169 //Check if requested trigger was fired
170 TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
171 if(!firedTrigClass.Contains(fTriggerClass))
175 fhNEvents->Fill(1.5);
181 //________________________________________________________________________
182 void AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects()
184 // Create user output.
188 AliAnalysisTaskEmcalJetDev::UserCreateOutputObjects();
190 Bool_t oldStatus = TH1::AddDirectoryStatus();
191 TH1::AddDirectory(kFALSE);
193 fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
194 fOutput->Add(fhNEvents);
199 fHistTrials = new TH1F("fHistTrials", "fHistTrials", 12, 0, 12);
200 fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
201 fHistTrials->GetYaxis()->SetTitle("trials");
202 fOutput->Add(fHistTrials);
204 fHistTrialsSelEvents = new TH1F("fHistTrialsSelEvents", "fHistTrialsSelEvents", 12, 0, 12);
205 fHistTrialsSelEvents->GetXaxis()->SetTitle("p_{T} hard bin");
206 fHistTrialsSelEvents->GetYaxis()->SetTitle("trials");
207 fOutput->Add(fHistTrialsSelEvents);
209 fHistXsection = new TProfile("fHistXsection", "fHistXsection", 12, 0, 12);
210 fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
211 fHistXsection->GetYaxis()->SetTitle("xsection");
212 fOutput->Add(fHistXsection);
214 fHistEvents = new TH1F("fHistEvents", "fHistEvents", 12, 0, 12);
215 fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
216 fHistEvents->GetYaxis()->SetTitle("total events");
217 fOutput->Add(fHistEvents);
220 TH1::AddDirectory(oldStatus);
222 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
225 //________________________________________________________________________
226 Bool_t AliAnalysisTaskEmcalDiJetBase::IsSameJet(const Int_t ijt, const Int_t ija, const Int_t type, const Bool_t isMC) {
227 //check if two jets are the same one
229 Bool_t bSame = kFALSE;
231 if(type<2 && ijt==ija)
235 if(fJetFullChargedMatchingType==kFraction) {
236 if(isMC && faFullFracIndexMC.At(ijt)==ija)
238 else if(!isMC && faFullFracIndex.At(ijt)==ija)
241 else if(fJetFullChargedMatchingType==kGeo) {
242 Int_t contFull = fContainerFull;
243 Int_t contChar = fContainerCharged;
245 contFull = fContainerFullMC;
246 contChar = fContainerChargedMC;
248 AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, contFull));
249 AliEmcalJet *chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, contChar));
250 AliEmcalJet *matchJet = fullJet->ClosestJet();
254 else if(fJetFullChargedMatchingType==kNoMatching) {
258 AliWarning(Form("%s: matching type unknown", GetName()));
268 //________________________________________________________________________
269 Double_t AliAnalysisTaskEmcalDiJetBase::GetJetPt(const AliEmcalJet *jet, const Int_t type) {
274 return jet->Pt() - fRhoFullVal * jet->Area();
276 return jet->Pt() - fRhoChVal * jet->Area();
282 //________________________________________________________________________
283 Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const AliVParticle *trk, const AliEmcalJet *jet) const
285 // Get Z of constituent trk
287 return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
290 //________________________________________________________________________
291 Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
294 // Get the z of a constituent inside of a jet
297 Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
300 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
302 AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
307 //________________________________________________________________________
308 Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const {
310 // Helper function to calculate the distance between two jets
313 Double_t dPhi = jet1->Phi() - jet2->Phi();
314 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
315 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
316 Double_t dEta = jet1->Eta() - jet2->Eta();
317 Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
321 //________________________________________________________________________
322 Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2) {
324 // Calculate azimuthal angle between the axises of the jets
327 return GetDeltaPhi(jet1->Phi(),jet2->Phi());
331 //________________________________________________________________________
332 Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaPhi(Double_t phi1,Double_t phi2) {
334 // Calculate azimuthal angle between the axises of the jets
337 Double_t dPhi = phi1-phi2;
339 if(dPhi > 2*TMath::Pi()) dPhi -= TMath::TwoPi();
340 if(dPhi <-2*TMath::Pi()) dPhi += TMath::TwoPi();
341 if(dPhi <-0.5*TMath::Pi()) dPhi += TMath::TwoPi();
342 if(dPhi > 1.5*TMath::Pi()) dPhi -= TMath::TwoPi();
347 //________________________________________________________________________
348 Double_t AliAnalysisTaskEmcalDiJetBase::GetFractionSharedPt(const AliEmcalJet *jetFull, const AliEmcalJet *jetCharged) const
351 // Get fraction of shared pT between matched full and charged jet
352 // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
355 Double_t fraction = 0.;
356 Double_t jetPtCh = jetCharged->Pt();
360 if(fDebug>10) AliInfo(Form("%s: nConstituents: %d, ch: %d chne: %d ne: %d",GetName(),jetFull->GetNumberOfConstituents(),jetCharged->GetNumberOfTracks(),jetFull->GetNumberOfTracks(),jetFull->GetNumberOfClusters()));
363 AliVParticle *vpf = 0x0;
365 for(Int_t icc=0; icc<jetCharged->GetNumberOfTracks(); icc++) {
366 Int_t idx = (Int_t)jetCharged->TrackAt(icc);
368 for(Int_t icf=0; icf<jetFull->GetNumberOfTracks(); icf++) {
369 if(idx == jetFull->TrackAt(icf) && iFound==0 ) {
371 vpf = static_cast<AliVParticle*>(jetFull->TrackAt(icf, fTracks));
380 fraction = sumPt/jetPtCh;
383 if(fDebug>10) AliInfo(Form("%s: charged shared fraction: %.2f",GetName(),fraction));
389 //________________________________________________________________________
390 void AliAnalysisTaskEmcalDiJetBase::MatchJetsGeo(Int_t cFull, Int_t cCharged,
391 Int_t iDebug, Float_t maxDist, Int_t type) {
394 // Match the full jets to the corresponding charged jets
395 // Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
396 // type: 0 = use acceptance cuts of container 1 = allow 0.1 one more for cCharged(MC) in eta 2 = allow 0.1 more in eta and phi for cCharged(MC)
400 const Int_t nChJets = GetNJets(cCharged);
401 const Int_t nFullJets = GetNJets(cFull);
403 if(nFullJets==0 || nChJets==0) return;
405 TArrayI faChNeMatchIndex;
406 faChNeMatchIndex.Set(nChJets+1);
408 TArrayI faChMatchIndex;
409 faChMatchIndex.Set(nFullJets+1);
411 static TArrayS iFlag(nChJets*nFullJets);
412 if(iFlag.GetSize()<(nChJets*nFullJets)){
413 iFlag.Set(nChJets*nFullJets+1);
417 AliJetContainer *contCh = GetJetContainer(cCharged);
419 // find the closest distance to the full jet
420 for(int ifu = 0;ifu<nFullJets;ifu++){
422 AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, cFull));
423 if(!fullJet) continue;
425 Float_t dist = maxDist;
427 for(int ich = 0;ich<nChJets;ich++){
428 AliEmcalJet *chJet = 0x0;
430 chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
432 chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
433 if(!chJet) continue;;
435 if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) || chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
438 if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) || chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
446 Double_t frac = GetFractionSharedPt(fullJet,chJet);
447 Double_t dR = GetDeltaR(fullJet,chJet);
449 faChMatchIndex[ifu]=ich;
452 if(iDebug>10) Printf("Distance (%d)--(%d) %3.3f frac:%.2f",ifu,ich,dR,frac);
454 if(faChMatchIndex[ifu]>=0) iFlag[ifu*nChJets+faChMatchIndex[ifu]]+=1;//ich closest to ifu
455 if(iDebug>10) Printf("Full Distance (%d)--(%d) %3.3f flag[%d] = %d",ifu,faChMatchIndex[ifu],dist,ifu*nChJets+faChMatchIndex[ifu],iFlag[ifu*nChJets+faChMatchIndex[ifu]]);
458 faChMatchIndex[ifu]=-1;
465 for(int ich = 0;ich<nChJets;ich++){
466 AliEmcalJet *chJet = 0x0;
468 chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
470 chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
471 if(!chJet) continue;;
473 if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) || chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
476 if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) || chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
484 Float_t dist = maxDist;
485 for(int ifu = 0;ifu<nFullJets;ifu++){
486 AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, cFull));
490 Double_t dR = GetDeltaR(fullJet,chJet);
492 faChNeMatchIndex[ich]=ifu;
496 if(faChNeMatchIndex[ich]>=0) iFlag[faChNeMatchIndex[ich]*nChJets+ich]+=2;//ifu closest to ich
497 if(iDebug>10) Printf("Other way Distance (%d)--(%d) %3.3f flag[%d] = %d",faChNeMatchIndex[ich],ich,dist,faChNeMatchIndex[ich]*nChJets+ich,iFlag[faChNeMatchIndex[ich]*nChJets+ich]);
500 faChNeMatchIndex[ich]=-1;
505 // check for "true" correlations
506 for(int ifu = 0;ifu<nFullJets;ifu++){
507 for(int ich = 0;ich<nChJets;ich++){
508 if(iDebug>10) AliInfo(Form("%s: Flag[%d][%d] %d ",GetName(),ifu,ich,iFlag[ifu*nChJets+ich]));
511 // we have a uniqe correlation
512 if(iFlag[ifu*nChJets+ich]==3){
514 AliEmcalJet *chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
515 AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetJetFromArray(ifu, cFull));
516 Double_t dR = GetDeltaR(fullJet,chJet);
518 if(iDebug>10) Printf("closest jets %d %d dR = %f",ich,ifu,dR);
520 chJet->SetClosestJet(fullJet,dR);
521 fullJet->SetClosestJet(chJet,dR);
528 fMatchingDone = kTRUE;
532 //________________________________________________________________________
533 void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndex() {
535 // take each full jet and set the index of the charged jet with the largest shared charged fraction
537 const Int_t nJetsCh = GetNJets(fContainerCharged);
538 const Int_t nJetsFull = GetNJets(fContainerFull);
539 faFullFracIndex.Set(nJetsFull+1);
540 faFullFracIndex.Reset(-1);
542 AliJetContainer *cont = GetJetContainer(fContainerFull);
543 Float_t radius = cont->GetJetRadius();
545 for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
548 AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, fContainerFull));
550 faFullFracIndex.SetAt(-1,ifu);
554 for(Int_t ich = 0; ich<nJetsCh; ich++) {
555 AliEmcalJet *jetCh = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, fContainerCharged));
558 Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
559 dist = GetDeltaR(jetFull,jetCh);
560 if(tmpFrac>frac && dist<radius) {
562 faFullFracIndex.SetAt(ich,ifu);
570 //________________________________________________________________________
571 void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndexMC() {
573 // take each full jet and set the index of the charged jet with the largest shared charged fraction
575 const Int_t nJetsCh = GetNJets(fContainerChargedMC);
576 const Int_t nJetsFull = GetNJets(fContainerFullMC);
577 faFullFracIndexMC.Set(nJetsFull);
578 faFullFracIndexMC.Reset(-1);
580 AliJetContainer *cont = GetJetContainer(fContainerFullMC);
581 Float_t radius = cont->GetJetRadius();
583 for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
586 AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, fContainerFullMC));
588 faFullFracIndexMC.SetAt(-1,ifu);
592 for(Int_t ich = 0; ich<nJetsCh; ich++) {
593 AliEmcalJet *jetCh = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, fContainerChargedMC));
596 Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
597 dist = GetDeltaR(jetFull,jetCh);
598 if(tmpFrac>frac && dist<radius) {
600 faFullFracIndexMC.SetAt(ich,ifu);
608 //________________________________________________________________________
609 Bool_t AliAnalysisTaskEmcalDiJetBase::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
612 // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
613 // Get the pt hard bin from the file path
614 // This is to called in Notify and should provide the path to the AOD/ESD file
615 // (Partially copied from AliAnalysisHelperJetTasks)
617 TString file(currFile);
621 if(file.Contains(".zip#")){
622 Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
623 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
624 Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
625 file.Replace(pos+1,pos2-pos1,"");
628 // not an archive take the basename....
629 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
631 //Printf("%s",file.Data());
633 // Get the pt hard bin
634 TString strPthard(file);
636 strPthard.Remove(strPthard.Last('/'));
637 strPthard.Remove(strPthard.Last('/'));
638 if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
639 strPthard.Remove(0,strPthard.Last('/')+1);
640 if (strPthard.IsDec())
641 pthard = strPthard.Atoi();
643 AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
645 TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
647 // next trial fetch the histgram file
648 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
650 // not a severe condition but inciate that we have no information
654 // find the tlist we want to be independtent of the name so use the Tkey
655 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
660 TList *list = dynamic_cast<TList*>(key->ReadObj());
665 xsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
666 trials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
669 } // no tree pyxsec.root
671 TTree *xtree = (TTree*)fxsec->Get("Xsection");
677 Double_t xsection = 0;
678 xtree->SetBranchAddress("xsection",&xsection);
679 xtree->SetBranchAddress("ntrials",&ntrials);
688 //________________________________________________________________________
689 Bool_t AliAnalysisTaskEmcalDiJetBase::UserNotify()
693 if (!fIsPythiaPtHard) {
698 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
700 AliError(Form("%s - UserNotify: No current tree!",GetName()));
704 Float_t xsection = 0;
708 TFile *curfile = tree->GetCurrentFile();
710 AliError(Form("%s - UserNotify: No current file!",GetName()));
714 TChain *chain = dynamic_cast<TChain*>(tree);
716 tree = chain->GetTree();
718 Int_t nevents = tree->GetEntriesFast();
720 PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthard);
726 fHistTrials->Fill(pthard, trials);
727 fHistXsection->Fill(pthard, xsection);
728 fHistEvents->Fill(pthard, nevents);
733 //________________________________________________________________________
734 Bool_t AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects() {
739 if (!AliAnalysisTaskEmcalJetDev::RetrieveEventObjects())
747 fRhoFullVal = GetRhoVal(fContainerFull);
748 fRhoChVal = GetRhoVal(fContainerCharged);
752 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
753 if (!fPythiaHeader) {
755 AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
758 for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
759 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
760 if (fPythiaHeader) break;
767 fPtHard = fPythiaHeader->GetPtHard();
769 fNTrials = fPythiaHeader->Trials();
777 //_______________________________________________________________________
778 void AliAnalysisTaskEmcalDiJetBase::Terminate(Option_t *)
780 // Called once at the end of the analysis.