1 // $Id: AliJetEmbeddingFromAODTask.cxx $
3 // Jet embedding from AOD task.
5 // Author: S.Aiola, C.Loizides
7 #include "AliJetEmbeddingFromAODTask.h"
9 // C++ standard library
15 #include <TClonesArray.h>
16 #include <TObjArray.h>
17 #include <TObjString.h>
21 #include <TStreamerInfo.h>
24 #include <TLorentzVector.h>
27 #include "AliVEvent.h"
28 #include "AliAODTrack.h"
29 #include "AliESDtrack.h"
30 #include "AliPicoTrack.h"
31 #include "AliVTrack.h"
32 #include "AliVCluster.h"
33 #include "AliVCaloCells.h"
34 #include "AliAODMCParticle.h"
35 #include "AliCentrality.h"
36 #include "AliVHeader.h"
37 #include "AliVVertex.h"
38 #include "AliAODHeader.h"
39 #include "AliFJWrapper.h"
41 #include "AliInputEventHandler.h"
43 ClassImp(AliJetEmbeddingFromAODTask)
45 //________________________________________________________________________
46 AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
47 AliJetModelBaseTask("AliJetEmbeddingFromAODTask"),
49 fRandomAccess(kFALSE),
56 fAODMCParticlesName(),
59 fTriggerMask(AliVEvent::kAny),
67 fJetConstituentMinPt(0),
71 fJetParticleLevel(kTRUE),
72 fIncludeNoITS(kFALSE),
73 fCutMaxFractionSharedTPCClusters(0.4),
74 fUseNegativeLabels(kTRUE),
96 fHistRejectedEvents(0)
98 // Default constructor.
99 SetSuffix("AODEmbedding");
100 fAODfilterBits[0] = -1;
101 fAODfilterBits[1] = -1;
110 //________________________________________________________________________
111 AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa) :
112 AliJetModelBaseTask(name, drawqa),
114 fRandomAccess(kFALSE),
115 fAODTreeName("aodTree"),
116 fAODHeaderName("header"),
117 fAODVertexName("vertices"),
118 fAODTrackName("tracks"),
120 fAODCellsName("emcalCells"),
121 fAODMCParticlesName(AliAODMCParticle::StdBranchName()),
124 fTriggerMask(AliVEvent::kAny),
132 fJetConstituentMinPt(0),
136 fJetParticleLevel(kTRUE),
137 fIncludeNoITS(kFALSE),
138 fCutMaxFractionSharedTPCClusters(0.4),
139 fUseNegativeLabels(kTRUE),
144 fEsdTreeMode(kFALSE),
146 fCurrentAODFileID(0),
148 fPicoTrackVersion(0),
156 fCurrentAODEntry(-1),
157 fHistFileMatching(0),
158 fHistAODFileError(0),
161 fHistRejectedEvents(0)
163 // Standard constructor.
164 SetSuffix("AODEmbedding");
165 fAODfilterBits[0] = -1;
166 fAODfilterBits[1] = -1;
175 //________________________________________________________________________
176 AliJetEmbeddingFromAODTask::~AliJetEmbeddingFromAODTask()
180 if (fCurrentAODFile) {
181 fCurrentAODFile->Close();
182 delete fCurrentAODFile;
186 //________________________________________________________________________
187 void AliJetEmbeddingFromAODTask::UserCreateOutputObjects()
192 AliJetModelBaseTask::UserCreateOutputObjects();
194 fHistFileMatching = new TH2C("fHistFileMatching", "fHistFileMatching", fTotalFiles, -0.5, fTotalFiles - 0.5, fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
195 fHistFileMatching->GetXaxis()->SetTitle("File no. (PYTHIA)");
196 fHistFileMatching->GetYaxis()->SetTitle("File no. (Embedded AOD)");
197 fHistFileMatching->GetZaxis()->SetTitle("counts");
198 fOutput->Add(fHistFileMatching);
200 fHistAODFileError = new TH1C("fHistAODFileError", "fHistAODFileError", fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
201 fHistAODFileError->GetXaxis()->SetTitle("File no. (Embedded AOD)");
202 fHistAODFileError->GetYaxis()->SetTitle("counts");
203 fOutput->Add(fHistAODFileError);
205 fHistNotEmbedded = new TH1C("fHistNotEmbedded", "fHistNotEmbedded", fTotalFiles, -0.5, fTotalFiles - 0.5);
206 fHistNotEmbedded->GetXaxis()->SetTitle("File no. (PYTHIA)");
207 fHistNotEmbedded->GetYaxis()->SetTitle("counts");
208 fOutput->Add(fHistNotEmbedded);
210 fHistEmbeddingQA = new TH1F("fHistEmbeddingQA", "fHistEmbeddingQA", 2, 0, 2);
211 fHistEmbeddingQA->GetXaxis()->SetTitle("Event state");
212 fHistEmbeddingQA->GetYaxis()->SetTitle("counts");
213 fHistEmbeddingQA->GetXaxis()->SetBinLabel(1, "OK");
214 fHistEmbeddingQA->GetXaxis()->SetBinLabel(2, "Not embedded");
215 fOutput->Add(fHistEmbeddingQA);
217 fHistRejectedEvents = new TH1F("fHistRejectedEvents", "fHistRejectedEvents", 500, -0.5, 499.5);
218 fHistRejectedEvents->GetXaxis()->SetTitle("# of rejected events");
219 fHistRejectedEvents->GetYaxis()->SetTitle("counts");
220 fOutput->Add(fHistRejectedEvents);
222 PostData(1, fOutput);
225 //________________________________________________________________________
226 Bool_t AliJetEmbeddingFromAODTask::UserNotify()
228 TString path(fInputHandler->GetTree()->GetCurrentFile()->GetPath());
229 if (path.EndsWith("/"))
230 path.Remove(path.Length()-1);
231 path.Remove(path.Last('/'));
232 path.Remove(0, path.Last('/')+1);
234 AliWarning(Form("Could not extract file number from path %s", path.Data()));
237 fCurrentFileID = path.Atoi();
238 if (!fRandomAccess) {
239 fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles-1;
240 AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
245 //________________________________________________________________________
246 Bool_t AliJetEmbeddingFromAODTask::ExecOnce()
248 if (fAODTreeName.Contains("aod"))
249 fEsdTreeMode = kFALSE;
251 fEsdTreeMode = kTRUE;
253 return AliJetModelBaseTask::ExecOnce();
256 //________________________________________________________________________
257 Bool_t AliJetEmbeddingFromAODTask::OpenNextFile()
259 if (fCurrentAODFile) {
260 fCurrentAODFile->Close();
261 delete fCurrentAODFile;
267 while ((!fCurrentAODFile || fCurrentAODFile->IsZombie()) && i < fAttempts) {
268 if (i > 0 && fHistAODFileError) {
269 fHistAODFileError->Fill(fCurrentAODFileID);
272 fCurrentAODFile = GetNextFile();
276 if (!fCurrentAODFile || fCurrentAODFile->IsZombie()) {
277 AliError("Could not open AOD file to embed!");
281 const TList *clist = fCurrentAODFile->GetStreamerInfoCache();
283 TStreamerInfo *cinfo = static_cast<TStreamerInfo*>(clist->FindObject("AliPicoTrack"));
285 fPicoTrackVersion = cinfo->GetClassVersion();
287 fPicoTrackVersion = 0;
290 fCurrentAODTree = static_cast<TTree*>(fCurrentAODFile->Get(fAODTreeName));
291 if (!fCurrentAODTree) {
292 AliError(Form("Could not get tree %s from file %s", fAODTreeName.Data(), fCurrentAODFile->GetName()));
296 if (!fAODHeaderName.IsNull())
297 fCurrentAODTree->SetBranchAddress(fAODHeaderName, &fAODHeader);
299 if (!fAODVertexName.IsNull())
300 fCurrentAODTree->SetBranchAddress(fAODVertexName, &fAODVertex);
302 if (!fAODTrackName.IsNull())
303 fCurrentAODTree->SetBranchAddress(fAODTrackName, &fAODTracks);
305 if (!fAODClusName.IsNull())
306 fCurrentAODTree->SetBranchAddress(fAODClusName, &fAODClusters);
308 if (!fAODCellsName.IsNull())
309 fCurrentAODTree->SetBranchAddress(fAODCellsName, &fAODCaloCells);
311 if (!fAODMCParticlesName.IsNull())
312 fCurrentAODTree->SetBranchAddress(fAODMCParticlesName, &fAODMCParticles);
315 fCurrentAODEntry = TMath::Nint(gRandom->Rndm()*fCurrentAODTree->GetEntries());
317 fCurrentAODEntry = -1;
319 AliDebug(3,Form("Will start embedding from entry %d", fCurrentAODEntry+1));
321 if (fHistFileMatching)
322 fHistFileMatching->Fill(fCurrentFileID, fCurrentAODFileID-1);
327 //________________________________________________________________________
328 TFile* AliJetEmbeddingFromAODTask::GetNextFile()
331 fCurrentAODFileID = TMath::Nint(gRandom->Rndm()*fFileList->GetEntriesFast());
335 if (fCurrentAODFileID >= fFileList->GetEntriesFast()) {
336 AliError("No more file in the list!");
340 TObjString *objFileName = static_cast<TObjString*>(fFileList->At(fCurrentAODFileID));
341 TString fileName(objFileName->GetString());
343 if (fileName.BeginsWith("alien://") && !gGrid) {
344 AliInfo("Trying to connect to AliEn ...");
345 TGrid::Connect("alien://");
348 TString baseFileName(fileName);
349 if (baseFileName.Contains(".zip#")) {
350 Ssiz_t pos = baseFileName.Last('#');
351 baseFileName.Remove(pos);
354 if (gSystem->AccessPathName(baseFileName)) {
355 AliError(Form("File %s does not exist!", baseFileName.Data()));
359 AliDebug(3,Form("Trying to open file %s...", fileName.Data()));
360 TFile *file = TFile::Open(fileName);
362 if (!file || file->IsZombie()) {
363 AliError(Form("Unable to open file: %s!", fileName.Data()));
370 //________________________________________________________________________
371 Bool_t AliJetEmbeddingFromAODTask::GetNextEntry()
376 if (!fCurrentAODFile || !fCurrentAODTree || fCurrentAODEntry+1 >= fCurrentAODTree->GetEntries()) {
377 if (!OpenNextFile()) {
378 AliError("Could not open the next file!");
383 if (!fCurrentAODTree) {
384 AliError("Could not get the tree!");
389 fCurrentAODTree->GetEntry(fCurrentAODEntry);
392 if (attempts == 1000)
393 AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
395 } while (!IsAODEventSelected());
397 fHistRejectedEvents->Fill(attempts);
399 if (!fCurrentAODTree)
405 //________________________________________________________________________
406 Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
408 // AOD event selection.
410 if (!fEsdTreeMode && fAODHeader) {
411 AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);
414 if (fTriggerMask != AliVEvent::kAny) {
415 UInt_t offlineTrigger = aodHeader->GetOfflineTrigger();
417 if ((offlineTrigger & fTriggerMask) == 0) {
418 AliDebug(2,Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.",
419 offlineTrigger, fTriggerMask));
424 // Centrality selection
425 if (fMinCentrality >= 0) {
426 AliCentrality *cent = aodHeader->GetCentralityP();
427 Float_t centVal = cent->GetCentralityPercentile("V0M");
428 if (centVal < fMinCentrality || centVal >= fMaxCentrality) {
429 AliDebug(2,Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f",
430 centVal, fMinCentrality, fMaxCentrality));
438 Double_t vert[3]={0};
439 ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
440 if (TMath::Abs(vert[2]) > fZVertexCut) {
441 AliDebug(2,Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f",
442 vert[2], fZVertexCut));
445 Double_t dist = TMath::Sqrt((vert[0]-fVertex[0])*(vert[0]-fVertex[0])+(vert[1]-fVertex[1])*(vert[1]-fVertex[1])+(vert[2]-fVertex[2])*(vert[2]-fVertex[2]));
446 if (dist > fMaxVertexDist) {
447 AliDebug(2,Form("Event rejected because the distance between the current and embedded verteces is > %f. "
448 "Current event vertex (%f, %f, %f), embedded event vertex (%f, %f, %f). Distance = %f",
449 fMaxVertexDist, fVertex[0], fVertex[1], fVertex[2], vert[0], vert[1], vert[2], dist));
459 if (fJetParticleLevel) {
461 jet = GetLeadingJet(fAODMCParticles);
463 AliWarning("Particle level jets selected, but not MC particles found. The jet event selection will be skipped.");
468 if (fAODTracks || fAODClusters)
469 jet = GetLeadingJet(fAODTracks, fAODClusters);
471 AliWarning("Detector level jets selected, but not tracks or clusters found. The jet event selection will be skipped.");
476 AliDebug(1, Form("Leading jet pt = %f", jet.Pt()));
477 if (jet.Pt() < fJetMinPt)
484 //________________________________________________________________________
485 void AliJetEmbeddingFromAODTask::Run()
487 if (!GetNextEntry()) {
488 if (fHistNotEmbedded)
489 fHistNotEmbedded->Fill(fCurrentFileID);
490 if (fHistEmbeddingQA)
491 fHistEmbeddingQA->Fill("Not embedded", 1);
492 AliError("Unable to get the AOD event to embed. Nothing will be embedded.");
496 if (fHistEmbeddingQA)
497 fHistEmbeddingQA->Fill("OK", 1);
499 if (fOutMCParticles) {
501 if (fCopyArray && fMCParticles)
504 if (fAODMCParticles) {
505 AliDebug(3, Form("%d MC particles will be processed for embedding.", fAODMCParticles->GetEntriesFast()));
506 for (Int_t i = 0; i < fAODMCParticles->GetEntriesFast(); i++) {
507 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fAODMCParticles->At(i));
509 AliError(Form("Could not find MC particle %d in branch %s of tree %s!", i, fAODMCParticlesName.Data(), fAODTreeName.Data()));
513 AliDebug(3, Form("Processing MC particle %d with pT = %f, eta = %f, phi = %f", i, part->Pt(), part->Eta(), part->Phi()));
515 if (!part->IsPhysicalPrimary())
518 if (part->Pt() < fPtMin || part->Pt() > fPtMax ||
519 part->Eta() < fEtaMin || part->Eta() > fEtaMax ||
520 part->Phi() < fPhiMin || part->Phi() > fPhiMax)
523 AddMCParticle(part, i);
524 AliDebug(3, "Embedded!");
531 if (fCopyArray && fTracks)
534 AliDebug(3, Form("Start embedding with %d tracks.", fOutTracks->GetEntriesFast()));
537 AliDebug(3, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast()));
538 for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) {
539 AliVTrack *track = static_cast<AliVTrack*>(fAODTracks->At(i));
541 AliError(Form("Could not find track %d in branch %s of tree %s!", i, fAODTrackName.Data(), fAODTreeName.Data()));
545 AliDebug(3, Form("Processing track %d with pT = %f, eta = %f, phi = %f, label = %d", i, track->Pt(), track->Eta(), track->Phi(), track->GetLabel()));
548 Bool_t isEmc = kFALSE;
550 AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
551 if (aodtrack->TestFilterBit(fAODfilterBits[0])) {
554 else if (aodtrack->TestFilterBit(fAODfilterBits[1])) {
555 if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
559 AliDebug(3, "Track not embedded because ITS refit failed.");
567 else { /*not a good track*/
568 AliDebug(3, "Track not embedded because not an hybrid track.");
571 if (fCutMaxFractionSharedTPCClusters > 0) {
572 Double_t frac = Double_t(aodtrack->GetTPCnclsS()) / Double_t(aodtrack->GetTPCncls());
573 if (frac > fCutMaxFractionSharedTPCClusters)
576 if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 &&
577 aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() &&
578 aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad())
581 else if (fPicoTrackVersion > 0) { /*not AOD mode, let's see if it is PicoTrack*/
582 AliPicoTrack *ptrack = dynamic_cast<AliPicoTrack*>(track);
584 if (fPicoTrackVersion >= 3)
585 type = ptrack->GetTrackType();
587 type = ptrack->GetLabel();
588 isEmc = ptrack->IsEMCAL();
592 if (track->Pt() < fPtMin || track->Pt() > fPtMax ||
593 track->Eta() < fEtaMin || track->Eta() > fEtaMax ||
594 track->Phi() < fPhiMin || track->Phi() > fPhiMax) {
595 AliDebug(3, "Track not embedded because out of limits.");
599 if (fTrackEfficiency < 1) {
600 Double_t r = gRandom->Rndm();
601 if (fTrackEfficiency < r) {
602 AliDebug(3, "Track not embedded because of artificial inefiiciency.");
609 if (fUseNegativeLabels)
610 label = track->GetLabel();
612 label = TMath::Abs(track->GetLabel());
615 AliDebug(1,Form("%s: Track %d with label==0", GetName(), i));
618 AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal(), isEmc, label, track->Charge());
619 AliDebug(3, "Track embedded!");
626 if (fCopyArray && fClusters)
630 for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) {
631 AliVCluster *clus = static_cast<AliVCluster*>(fAODClusters->At(i));
633 AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data()));
637 if (!clus->IsEMCAL())
641 Double_t vert[3] = {0,0,0};
642 clus->GetMomentum(vect,vert);
644 if (vect.Pt() < fPtMin || vect.Pt() > fPtMax ||
645 vect.Eta() < fEtaMin || vect.Eta() > fEtaMax ||
646 vect.Phi() < fPhiMin || vect.Phi() > fPhiMax)
651 label = clus->GetLabel();
653 AliDebug(1,Form("%s: Clus %d with label<=0", GetName(), i));
662 Double_t totalEnergy = 0;
663 Int_t totalCells = 0;
666 totalCells += fCaloCells->GetNumberOfCells();
668 totalCells += fAODCaloCells->GetNumberOfCells();
670 SetNumberOfOutCells(totalCells);
676 for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
680 Short_t cellNum = -1;
683 fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
687 AliDebug(1,Form("%s: Cell %d with label<=0", GetName(), i));
693 AliDebug(2,Form("Adding cell with amplitude %f, absolute ID %d, time %f, mc label %d", amp, cellNum, time, mclabel));
694 AddCell(amp, cellNum, time, mclabel);
698 AliDebug(2,Form("Added cells = %d (energy = %f), total cells = %d", fAddedCells, totalEnergy, totalCells));
702 //________________________________________________________________________
703 TLorentzVector AliJetEmbeddingFromAODTask::GetLeadingJet(TClonesArray *tracks, TClonesArray *clusters)
706 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
709 jalgo = fastjet::antikt_algorithm;
713 AliFJWrapper fjw(name, name);
714 fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
715 fjw.SetGhostArea(1); // set a very large ghost area to speed up jet finding
716 fjw.SetR(fJetRadius);
717 fjw.SetAlgorithm(jalgo);
722 const Int_t Ntracks = tracks->GetEntries();
723 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
724 AliVParticle *t = static_cast<AliVParticle*>(tracks->At(iTracks));
728 AliAODMCParticle *aodmcpart = dynamic_cast<AliAODMCParticle*>(t);
729 if (aodmcpart && !aodmcpart->IsPhysicalPrimary())
732 if ((fJetType == 1 && t->Charge() == 0) ||
733 (fJetType == 2 && t->Charge() != 0))
736 if (t->Pt() < fJetConstituentMinPt)
739 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
743 if (clusters && fJetType != 1) {
744 Double_t vert[3]={0};
746 ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
748 const Int_t Nclus = clusters->GetEntries();
749 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
750 AliVCluster *c = static_cast<AliVCluster*>(clusters->At(iClus));
758 c->GetMomentum(nP, vert);
760 if (nP.Pt() < fJetConstituentMinPt)
763 fjw.AddInputVector(nP.Px(), nP.Py(), nP.Pz(), nP.P(), -iClus - 100);
770 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
771 AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
775 Int_t njets = jets_incl.size();
778 //std::vector<fastjet::PseudoJet> jets_incl_sorted = fastjet::sorted_by_pt(jets_incl);
779 for (Int_t i = 0; i < njets; i++) {
780 if (jet.Pt() >= jets_incl[i].perp())
782 if (jets_incl[i].eta() < fJetMinEta || jets_incl[i].eta() > fJetMaxEta || jets_incl[i].phi() < fJetMinPhi || jets_incl[i].phi() > fJetMaxPhi)
784 jet.SetPxPyPzE(jets_incl[i].px(), jets_incl[i].py(), jets_incl[i].pz(), jets_incl[i].E());