AliPhysics  b76e98e (b76e98e)
AliReducedHighPtEventCreator.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2015, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 #include <cfloat>
16 #include <map>
17 #include <vector>
18 
19 #include <TArrayD.h>
20 #include <TArrayI.h>
21 #include <TClonesArray.h>
22 #include <TMath.h>
23 #include <TObjArray.h>
24 #include <TRandom.h>
25 #include <TString.h>
26 #include <TTree.h>
27 
28 #include "AliAnalysisUtils.h"
29 #include "AliAODEvent.h"
30 #include "AliAODMCParticle.h"
31 #include "AliAODTrack.h"
32 #include "AliCentrality.h"
33 #include "AliESDEvent.h"
34 #include "AliESDtrack.h"
35 #include "AliEmcalTrackSelection.h"
36 #include "AliInputEventHandler.h"
37 #include "AliLog.h"
38 #include "AliMCEvent.h"
39 #include "AliPicoTrack.h"
40 #include "AliStack.h"
41 #include "AliVCluster.h"
42 #include "AliVTrack.h"
43 #include "AliVVertex.h"
44 #include "AliEMCALTriggerPatchInfo.h"
45 
46 #include "AliReducedEmcalCluster.h"
48 #include "AliReducedHighPtEvent.h"
49 #include "AliReducedMCHeader.h"
51 #include "AliReducedPatchInfo.h"
54 
56 ClassImp(HighPtTracks::AliReducedTrackSelectionContainer)
57 ClassImp(HighPtTracks::AliReducedHighPtEventCreator)
59 
60 namespace HighPtTracks {
61 
65 AliReducedHighPtEventCreator::AliReducedHighPtEventCreator():
67  fOutputTree(NULL),
68  fOutputEvent(NULL),
69  fTrackSelections(NULL),
70  fEventSelectionBits(AliVEvent::kAny),
71  fSwapTriggerThresholds(kFALSE),
72  fMinClusterE(-1),
73  fMaxClusterE(1000),
74  fMinPt(-1),
75  fMaxPt(1000),
76  fMinEta(-1000),
77  fMaxEta(1000),
78  fKeepFractionEvents(1.),
79  fApplyCentralitySelection(kFALSE),
80  fCentralityMethod("V0A"),
81  fTriggerSetup("4classes"),
82  fMinBiasSelection(AliVEvent::kINT7)
83 {
84  fSelectCentralityRange[0] = 0.;
85  fSelectCentralityRange[1] = 100.;
86 }
87 
92 AliReducedHighPtEventCreator::AliReducedHighPtEventCreator(const char* name):
93  AliAnalysisTaskEmcal(name, kTRUE),
94  fOutputTree(NULL),
95  fOutputEvent(NULL),
96  fTrackSelections(NULL),
97  fEventSelectionBits(AliVEvent::kAny),
98  fSwapTriggerThresholds(kFALSE),
99  fMinClusterE(-1),
100  fMaxClusterE(1000),
101  fMinPt(-1),
102  fMaxPt(1000),
103  fMinEta(-1000),
104  fMaxEta(1000),
105  fKeepFractionEvents(1.),
106  fApplyCentralitySelection(kFALSE),
107  fCentralityMethod("V0A"),
108  fTriggerSetup("4classes"),
109  fMinBiasSelection(AliVEvent::kINT7)
110 {
111  DefineOutput(2, TTree::Class());
112 
113  fSelectCentralityRange[0] = 0.;
114  fSelectCentralityRange[1] = 100.;
115 
116  fTrackSelections = new TObjArray;
117  fTrackSelections->SetOwner(kTRUE);
118 }
119 
123 AliReducedHighPtEventCreator::~AliReducedHighPtEventCreator() {
124  if(fTrackSelections) delete fTrackSelections;
125 }
126 
130 void AliReducedHighPtEventCreator::UserCreateOutputObjects() {
132  OpenFile(2);
133 
134  fOutputTree = new TTree("ReducedEvent", "A simple reduced event");
135  fOutputEvent = new AliReducedHighPtEvent();
136  fOutputTree->Branch("ReducedEvent", "AliReducedHighPtEvent", fOutputEvent, 128000, 0);
137 
138  PostData(1, fOutput);
139  PostData(2, fOutputTree);
140 }
141 
152 Bool_t AliReducedHighPtEventCreator::Run() {
153  AliDebug(1, "Starting creation of the reconstructed event");
154  if(!fCaloClusters){
155  AliError("Cluster container missing");
156  return kFALSE;
157  }
158  if(!fTracks){
159  AliError("Track container missing");
160  return kFALSE;
161  }
162  if(!(this->fInputHandler->IsEventSelected() & fEventSelectionBits)){
163  AliDebug(1, "Trigger not selected");
164  return kFALSE;
165  } else {
166  AliDebug(1, "Trigger selected");
167  }
168  if(!SelectEvent(fInputEvent)) return kFALSE;
169  new(fOutputEvent) AliReducedHighPtEvent(kTRUE);
170 
171  // Write event specific information
172  fOutputEvent->SetVertexZ(fInputEvent->GetPrimaryVertex()->GetZ());
173  AliCentrality *centralityHandler = fInputEvent->GetCentrality();
174  if(centralityHandler) fOutputEvent->SetCentralityPercentile(centralityHandler->GetCentralityPercentile(fCentralityMethod.Data()));
175  ConvertTriggerPatches(fTriggerPatchInfo, fOutputEvent->GetPatchContainer());
176  TString triggerString(fInputEvent->GetFiredTriggerClasses());
177  if(!fTriggerSetup.CompareTo("4classes"))
178  fOutputEvent->SetDecisionFromTriggerString(triggerString.Contains("EG1"), triggerString.Contains("EG2"), triggerString.Contains("EJ1"), triggerString.Contains("EJ2"));
179  else
180  fOutputEvent->SetDecisionFromTriggerString(triggerString.Contains("EGA"), kFALSE, triggerString.Contains("EJE"), kFALSE);
181  fOutputEvent->SetMinBiasEvent(fInputHandler->IsEventSelected() & fMinBiasSelection);
182  fOutputEvent->SetRunNumber(fInputEvent->GetRunNumber());
183 
184  std::map<int, int> mcindexmap;
185  if(MCEvent()){
186  // Generate Monte-Carlo header
187  AliReducedMCHeader *mcheader = new AliReducedMCHeader();
188  if(fIsPythia){
189  mcheader->SetCrossSection(fXsection);
190  mcheader->SetNumberOfTrials(fNTrials);
191  mcheader->SetPtHard(fPtHard);
192  }
193  fOutputEvent->SetMonteCarloHeader(mcheader);
194 
195  // Convert Monte Carlo particles, fill index map
196  AliMCEvent *mcev = MCEvent();
197  int npart(0);
198  for(Int_t ipart = 0; ipart < mcev->GetNumberOfTracks(); ipart++){
199  AliVParticle *part = mcev->GetTrack(ipart);
200  Double_t pt(TMath::Abs(part->Pt())), eta(part->Eta());
201  if(pt < fMinPt || pt > fMaxPt) continue;
202  if(eta < fMinEta || eta > fMaxEta) continue;
203  if(!part->Charge()) continue;
204  if(part->IsA() == AliAODMCParticle::Class()){
205  AliAODMCParticle *aodpart = static_cast<AliAODMCParticle *>(part);
206  if(!aodpart->IsPhysicalPrimary()) continue;
207  } else {
208  // ESD - get information from stack
209  if(!mcev->Stack()->IsPhysicalPrimary(ipart)) continue;
210  }
211 
212  AliReducedGeneratedParticle *reducedgen = new AliReducedGeneratedParticle(npart, part->PdgCode(), part->Px(), part->Py(), part->Pz(), part->E());
213  fOutputEvent->AddReducedGeneratedParticle(reducedgen);
214  mcindexmap.insert(std::pair<int,int>(ipart, npart));
215  npart++;
216  }
217  }
218 
219  // Filter Clusters
220  std::map<int, int> clusterindexmap;
221  int ncluster(0);
222  Double_t vtxpos[3];fInputEvent->GetPrimaryVertex()->GetXYZ(vtxpos);
223  for(TIter cliter = TIter(fCaloClusters).Begin(); cliter != TIter::End(); ++cliter){
224  AliVCluster *incluster = static_cast<AliVCluster *>(*cliter);
225  if(!SelectCluster(incluster)) continue;
226  TLorentzVector clustervec;
227  incluster->GetMomentum(clustervec, vtxpos);
228  AliReducedEmcalCluster *redcluster = new AliReducedEmcalCluster(ncluster, incluster->E(), clustervec.Eta(), clustervec.Phi(), incluster->GetM02(), incluster->GetM20());
229  // Get leading 3 cell energies
230  TArrayD cellEnergies;
231  GetCellEnergies(incluster, cellEnergies);
232  TArrayI indices(cellEnergies.GetSize());
233  TMath::Sort(cellEnergies.GetSize(), cellEnergies.GetArray(), indices.GetArray(), kTRUE);
234  redcluster->SetLeadingCellEnergies(
235  cellEnergies[indices[0]],
236  cellEnergies.GetSize() > 1 ? cellEnergies[indices[1]] : 0,
237  cellEnergies.GetSize() > 2 ? cellEnergies[indices[2]] : 0
238  );
239  // Assing MC particles
240  if(MCEvent()){
241  for(Int_t ilab = 0; ilab < incluster->GetNLabels(); ilab++){
242  AliVParticle *assigned = MCEvent()->GetTrack(TMath::Abs(incluster->GetLabels()[ilab]));
243  if(!assigned) continue;
244  redcluster->AddTrueContributor(assigned->PdgCode(), assigned->Px(), assigned->Py(), assigned->Pz(), assigned->E());
245  }
246  }
247  fOutputEvent->AddReducedCluster(redcluster);
248  clusterindexmap.insert(std::pair<int,int>(incluster->GetID(), ncluster));
249  ncluster++;
250  }
251 
252  // Filter tracks
253  int nsel(0);
254  AliDebug(1, Form("Number of tracks in container: %d", fTracks->GetEntries()));
255  for(TIter trkiter = TIter(fTracks).Begin(); trkiter != TIter::End(); ++trkiter){
256  AliVTrack *trackRaw = static_cast<AliVTrack *>(*trkiter), *trackToCheck(NULL);
257  // handlinf for pico tracks as input
258  if(trackRaw->IsA() == AliPicoTrack::Class()){
259  AliPicoTrack *picotrack = static_cast<AliPicoTrack *>(trackRaw);
260  trackToCheck = picotrack->GetTrack();
261  } else {
262  trackToCheck = trackRaw;
263  }
264  FixTrackInputEvent(trackToCheck);
265  TArrayI cutIndices;
266  if(!SelectTrack(trackToCheck,cutIndices)) continue;
267  AliReducedReconstructedTrack *rectrack = new AliReducedReconstructedTrack;
268  double pvec[3];
269  trackToCheck->GetPxPyPz(pvec);
270  rectrack->SetMomentumVector(pvec[0], pvec[1], pvec[2]);
271  rectrack->SetCharge(static_cast<Char_t>(trackToCheck->Charge()));
272  for(Int_t icut = 0; icut < cutIndices.GetSize(); icut++) rectrack->SetTrackCuts(cutIndices[icut]);
273  if(MCEvent()){
274  // handle label
275  Int_t label = TMath::Abs(trackToCheck->GetLabel());
276  if(label < fMCEvent->GetNumberOfTracks()){
277  std::map<int,int>::iterator found = mcindexmap.find(label);
278  if(found != mcindexmap.end()){
279  rectrack->SetMatchedParticleIndex(found->second);
280  }
281  if(label >= 0) rectrack->SetGoodTrackLabel(kTRUE);
282  rectrack->SetTPCClusters(trackToCheck->GetTPCNcls());
283  rectrack->SetTPCCrossedRows(GetTPCCrossedRows(trackToCheck));
284  rectrack->SetTPCSharedClusters(trackToCheck->GetTPCSharedMapPtr()->CountBits());
285  rectrack->SetTPCFindableClusters(trackToCheck->GetTPCNclsF());
286  }
287  }
288  // handle clusters
289  Int_t clusterIndex = trackToCheck->GetEMCALcluster();
290  if(clusterIndex >= 0 && clusterIndex < fCaloClusters->GetEntries()){
291  std::map<int,int>::iterator found = clusterindexmap.find(clusterIndex);
292  if(found != clusterindexmap.end()){
293  rectrack->SetMatchedClusterIndex(found->second);
294  }
295  }
296  //std::cout << "Track selected" << std::endl;
297  fOutputEvent->AddReducedReconstructedParticle(rectrack);
298  nsel++;
299  }
300  AliDebug(1, Form("Number of selected tracks :%d ", nsel ));
301 
302  fOutputTree->Fill();
303  PostData(2, fOutputTree);
304  return kTRUE;
305 }
306 
312 void AliReducedHighPtEventCreator::AddVirtualTrackSelection(AliEmcalTrackSelection* sel, int cutindex) {
313  fTrackSelections->Add(new AliReducedTrackSelectionContainer(cutindex, sel));
314 }
315 
322 Bool_t AliReducedHighPtEventCreator::SelectEvent(AliVEvent* event) const {
323  if(fApplyCentralitySelection && InputEvent()->GetCentrality()){
324  AliDebug(1, Form("Centrality selection applied in range %f - %f\n", fSelectCentralityRange[0], fSelectCentralityRange[1]));
325  Float_t centrality = InputEvent()->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
326  AliDebug(1, Form("Centrality before: %f, estimator %s\n", centrality, fCentralityMethod.Data()));
327  if(centrality < fSelectCentralityRange[0] || centrality > fSelectCentralityRange[1]) return kFALSE;
328  AliDebug(1, "Centrality selected");
329  }
330  const AliVVertex *primvtx = event->GetPrimaryVertex();
331  if(!primvtx || !primvtx->GetNContributors()) return kFALSE;
332  if(TMath::Abs(primvtx->GetZ()) > 10.) return kFALSE;
333  if(event->IsPileupFromSPD(3, 0.8, 3., 2., 5.)) return kFALSE;
334  AliAnalysisUtils eventSelUtil;
335  if(!eventSelUtil.IsVertexSelected2013pA(event)) return kFALSE;
336  if(fKeepFractionEvents < (1. - FLT_EPSILON)){
337  // Apply downscaling
338  Float_t choice = gRandom->Uniform(0., 1.);
339  if(choice > fKeepFractionEvents) return kFALSE; // event rejected by downscaling
340  }
341  return kTRUE;
342 }
343 
349 Bool_t AliReducedHighPtEventCreator::SelectCluster(const AliVCluster* clust) const {
350  if(!clust->IsEMCAL()) return kFALSE;
351  if(clust->E() < fMinClusterE || clust->E() > fMaxClusterE) return kFALSE;
352  return kTRUE;
353 }
354 
361 Int_t AliReducedHighPtEventCreator::SelectTrack(AliVTrack* track, TArrayI& cutindices) const {
362  std::vector<int> selected;
363 
364  double pt = TMath::Abs(track->Pt()), eta = track->Eta();
365  if(pt < fMinPt || pt > fMaxPt) return 0;
366  if(eta < fMinEta || eta > fMaxEta) return 0;
367 
368  for(TIter cutiter = TIter(fTrackSelections).Begin(); cutiter != TIter::End(); ++cutiter){
369  AliReducedTrackSelectionContainer *mycut = static_cast<AliReducedTrackSelectionContainer *>(*cutiter);
370  if(!mycut->GetTrackSelection()->IsTrackAccepted(track)) continue;
371  selected.push_back(mycut->GetIndex());
372  }
373  if(!selected.size()) return 0;
374  cutindices.Set(selected.size());
375  int counter = 0;
376  for(std::vector<int>::iterator inditer = selected.begin(); inditer != selected.end(); ++inditer){
377  cutindices[counter++] = *inditer;
378  }
379  return selected.size();
380 }
381 
387 Int_t AliReducedHighPtEventCreator::GetTPCCrossedRows(const AliVTrack* trk) const {
388  if(trk->IsA() == AliESDtrack::Class()){
389  return (static_cast<const AliESDtrack *>(trk))->GetTPCCrossedRows();
390  } else if(trk->IsA() == AliAODTrack::Class()){
391  return (static_cast<const AliAODTrack *>(trk))->GetTPCNCrossedRows();
392  }
393  return 0;
394 }
395 
401 void AliReducedHighPtEventCreator::GetCellEnergies(AliVCluster* emccluster, TArrayD& energies) const {
402  if(!fInputEvent->GetEMCALCells()) {
403  AliError("No EMCAL cells");
404  return;
405  }
406  AliDebug(2, Form("Number of cells: %d, array: %p", emccluster->GetNCells(), emccluster->GetCellsAbsId()));
407  energies.Set(emccluster->GetNCells());
408  for(int icell = 0; icell < emccluster->GetNCells(); icell++){
409  // printf("Cell ID: %d\n", emccluster->GetCellsAbsId()[icell]);
410  energies[icell] = fInputEvent->GetEMCALCells()->GetCellAmplitude(emccluster->GetCellsAbsId()[icell]);
411  }
412 }
413 
419 void AliReducedHighPtEventCreator::ConvertTriggerPatches(TClonesArray* patches,
420  AliReducedPatchContainer* cont) {
421  if(!patches){
422  AliError("Trigger patch container not found\n");
423  return;
424  }
425  for(TIter patchIter = TIter(patches).Begin(); patchIter != TIter::End(); ++patchIter){
426  AliEMCALTriggerPatchInfo *mypatch = static_cast<AliEMCALTriggerPatchInfo *>(*patchIter);
427  if(!mypatch->IsOfflineSimple() && mypatch->IsLevel0()) continue;
429  Bool_t isDefined(kFALSE);
430  if(mypatch->IsOfflineSimple()){
431  if(mypatch->IsGammaHighSimple()) { triggertype = AliReducedPatchContainer::kEMCGammaHigh; isDefined = kTRUE; }
432  if(mypatch->IsGammaLowSimple()) { triggertype = AliReducedPatchContainer::kEMCGammaLow; isDefined = kTRUE; }
433  if(mypatch->IsJetHighSimple()) { triggertype = AliReducedPatchContainer::kEMCJetHigh; isDefined = kTRUE; }
434  if(mypatch->IsJetLowSimple()) { triggertype = AliReducedPatchContainer::kEMCJetLow; isDefined = kTRUE; }
435  if(!isDefined){
436  AliDebug(2, "Unknown offline patch type");
437  continue;
438  }
439  AliDebug(2, Form("Adding offline patch of type %d", int(triggertype)));
440  cont->AddTriggerPatch(kTRUE, triggertype, mypatch->GetPatchE(), mypatch->GetADCAmp(), mypatch->GetEtaGeo(), mypatch->GetPhiGeo());
441  } else {
442  if(mypatch->IsGammaHigh()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCGammaLow : AliReducedPatchContainer::kEMCGammaHigh; isDefined=kTRUE; }
443  if(mypatch->IsGammaLow()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCGammaHigh : AliReducedPatchContainer::kEMCGammaLow; isDefined=kTRUE; }
444  if(mypatch->IsJetHigh()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCJetLow : AliReducedPatchContainer::kEMCJetHigh; isDefined=kTRUE; }
445  if(mypatch->IsJetLow()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCJetHigh : AliReducedPatchContainer::kEMCJetLow; isDefined=kTRUE; }
446  if(!isDefined){
447  AliDebug(2, "Unknown online patch type");
448  continue;
449  }
450  AliDebug(2, Form("Adding online patch of type %d", int(triggertype)));
451  cont->AddTriggerPatch(kFALSE, triggertype, mypatch->GetPatchE(), mypatch->GetADCAmp(), mypatch->GetEtaGeo(), mypatch->GetPhiGeo());
452  }
453  }
454 }
455 
461 void AliReducedHighPtEventCreator::FixTrackInputEvent(AliVTrack* trk) {
462  if(!trk->GetEvent()){
463  if(trk->IsA() == AliESDtrack::Class())
464  (static_cast<AliESDtrack *>(trk))->SetESDEvent(static_cast<AliESDEvent *>(fInputEvent));
465  else if(trk->IsA() == AliAODTrack::Class())
466  (static_cast<AliAODTrack *>(trk))->SetAODEvent(static_cast<AliAODEvent *>(fInputEvent));
467  else if(trk->IsA() == AliPicoTrack::Class()){
468  AliPicoTrack *mytrk = static_cast<AliPicoTrack *>(trk);
469  if(!mytrk->GetEvent()){
470  if(mytrk->GetTrack()->IsA() == AliESDtrack::Class())
471  (static_cast<AliESDtrack *>(mytrk->GetTrack()))->SetESDEvent(static_cast<AliESDEvent *>(fInputEvent));
472  else
473  (static_cast<AliAODTrack *>(mytrk->GetTrack()))->SetAODEvent(static_cast<AliAODEvent *>(fInputEvent));
474  }
475  }
476  }
477 }
478 
482 AliReducedTrackSelectionContainer::AliReducedTrackSelectionContainer():
483  fIndex(-1),
484  fTrackSelection(NULL)
485 {
486 }
487 
493 AliReducedTrackSelectionContainer::AliReducedTrackSelectionContainer(Int_t index, AliEmcalTrackSelection* sel):
494  fIndex(index),
495  fTrackSelection(sel)
496 {
497 }
498 
502 AliReducedTrackSelectionContainer::~AliReducedTrackSelectionContainer() {
503  if(fTrackSelection) delete fTrackSelection;
504 }
505 
506 } /* namespace HighPtTracks */
507 
Interface for virtual track selection.
Reduced event structure for high- analysis.
double Double_t
Definition: External.C:58
Base task in the EMCAL framework.
Declaration of class AliReducedHighPtEventCreator and AliReducedTrackSelectionContainer.
Declaration of class AlliReducedPatchInfo.
centrality
Declaration of class AliReducedPatchContainer, a container for reduced trigger patches.
TRandom * gRandom
int Int_t
Definition: External.C:63
Declaration of a reduced MC event header.
float Float_t
Definition: External.C:68
Namespace for classes creating trees of events with jets.
Declaration of class AliReducedGeneratedParticle.
AliVTrack * GetTrack() const
Definition: AliPicoTrack.h:60
Reduced information about reconstructed EMCAL clusters.
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65