AliPhysics  vAN-20150827 (3e81cbb)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
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 "AliInputEventHandler.h"
36 #include "AliLog.h"
37 #include "AliMCEvent.h"
38 #include "AliPicoTrack.h"
39 #include "AliStack.h"
40 #include "AliVCluster.h"
41 #include "AliVTrack.h"
42 #include "AliVVertex.h"
44 
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(
314  fTrackSelections->Add(new AliReducedTrackSelectionContainer(cutindex, sel));
315 }
316 
323 Bool_t AliReducedHighPtEventCreator::SelectEvent(AliVEvent* event) const {
324  if(fApplyCentralitySelection && InputEvent()->GetCentrality()){
325  AliDebug(1, Form("Centrality selection applied in range %f - %f\n", fSelectCentralityRange[0], fSelectCentralityRange[1]));
326  Float_t centrality = InputEvent()->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
327  AliDebug(1, Form("Centrality before: %f, estimator %s\n", centrality, fCentralityMethod.Data()));
328  if(centrality < fSelectCentralityRange[0] || centrality > fSelectCentralityRange[1]) return kFALSE;
329  AliDebug(1, "Centrality selected");
330  }
331  const AliVVertex *primvtx = event->GetPrimaryVertex();
332  if(!primvtx || !primvtx->GetNContributors()) return kFALSE;
333  if(TMath::Abs(primvtx->GetZ()) > 10.) return kFALSE;
334  if(event->IsPileupFromSPD(3, 0.8, 3., 2., 5.)) return kFALSE;
335  AliAnalysisUtils eventSelUtil;
336  if(!eventSelUtil.IsVertexSelected2013pA(event)) return kFALSE;
337  if(fKeepFractionEvents < (1. - FLT_EPSILON)){
338  // Apply downscaling
339  Float_t choice = gRandom->Uniform(0., 1.);
340  if(choice > fKeepFractionEvents) return kFALSE; // event rejected by downscaling
341  }
342  return kTRUE;
343 }
344 
350 Bool_t AliReducedHighPtEventCreator::SelectCluster(const AliVCluster* clust) const {
351  if(!clust->IsEMCAL()) return kFALSE;
352  if(clust->E() < fMinClusterE || clust->E() > fMaxClusterE) return kFALSE;
353  return kTRUE;
354 }
355 
362 Int_t AliReducedHighPtEventCreator::SelectTrack(AliVTrack* track, TArrayI& cutindices) const {
363  std::vector<int> selected;
364 
365  double pt = TMath::Abs(track->Pt()), eta = track->Eta();
366  if(pt < fMinPt || pt > fMaxPt) return 0;
367  if(eta < fMinEta || eta > fMaxEta) return 0;
368 
369  for(TIter cutiter = TIter(fTrackSelections).Begin(); cutiter != TIter::End(); ++cutiter){
370  AliReducedTrackSelectionContainer *mycut = static_cast<AliReducedTrackSelectionContainer *>(*cutiter);
371  if(!mycut->GetTrackSelection()->IsTrackAccepted(track)) continue;
372  selected.push_back(mycut->GetIndex());
373  }
374  if(!selected.size()) return 0;
375  cutindices.Set(selected.size());
376  int counter = 0;
377  for(std::vector<int>::iterator inditer = selected.begin(); inditer != selected.end(); ++inditer){
378  cutindices[counter++] = *inditer;
379  }
380  return selected.size();
381 }
382 
388 Int_t AliReducedHighPtEventCreator::GetTPCCrossedRows(const AliVTrack* trk) const {
389  if(trk->IsA() == AliESDtrack::Class()){
390  return (static_cast<const AliESDtrack *>(trk))->GetTPCCrossedRows();
391  } else if(trk->IsA() == AliAODTrack::Class()){
392  return (static_cast<const AliAODTrack *>(trk))->GetTPCNCrossedRows();
393  }
394  return 0;
395 }
396 
402 void AliReducedHighPtEventCreator::GetCellEnergies(AliVCluster* emccluster, TArrayD& energies) const {
403  if(!fInputEvent->GetEMCALCells()) {
404  AliError("No EMCAL cells");
405  return;
406  }
407  AliDebug(2, Form("Number of cells: %d, array: %p", emccluster->GetNCells(), emccluster->GetCellsAbsId()));
408  energies.Set(emccluster->GetNCells());
409  for(int icell = 0; icell < emccluster->GetNCells(); icell++){
410  // printf("Cell ID: %d\n", emccluster->GetCellsAbsId()[icell]);
411  energies[icell] = fInputEvent->GetEMCALCells()->GetCellAmplitude(emccluster->GetCellsAbsId()[icell]);
412  }
413 }
414 
420 void AliReducedHighPtEventCreator::ConvertTriggerPatches(TClonesArray* patches,
421  AliReducedPatchContainer* cont) {
422  if(!patches){
423  AliError("Trigger patch container not found\n");
424  return;
425  }
426  for(TIter patchIter = TIter(patches).Begin(); patchIter != TIter::End(); ++patchIter){
427  AliEmcalTriggerPatchInfo *mypatch = static_cast<AliEmcalTriggerPatchInfo *>(*patchIter);
428  if(!mypatch->IsOfflineSimple() && mypatch->IsLevel0()) continue;
430  Bool_t isDefined(kFALSE);
431  if(mypatch->IsOfflineSimple()){
432  if(mypatch->IsGammaHighSimple()) { triggertype = AliReducedPatchContainer::kEMCGammaHigh; isDefined = kTRUE; }
433  if(mypatch->IsGammaLowSimple()) { triggertype = AliReducedPatchContainer::kEMCGammaLow; isDefined = kTRUE; }
434  if(mypatch->IsJetHighSimple()) { triggertype = AliReducedPatchContainer::kEMCJetHigh; isDefined = kTRUE; }
435  if(mypatch->IsJetLowSimple()) { triggertype = AliReducedPatchContainer::kEMCJetLow; isDefined = kTRUE; }
436  if(!isDefined){
437  AliDebug(2, "Unknown offline patch type");
438  continue;
439  }
440  AliDebug(2, Form("Adding offline patch of type %d", int(triggertype)));
441  cont->AddTriggerPatch(kTRUE, triggertype, mypatch->GetPatchE(), mypatch->GetADCAmp(), mypatch->GetEtaGeo(), mypatch->GetPhiGeo());
442  } else {
443  if(mypatch->IsGammaHigh()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCGammaLow : AliReducedPatchContainer::kEMCGammaHigh; isDefined=kTRUE; }
444  if(mypatch->IsGammaLow()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCGammaHigh : AliReducedPatchContainer::kEMCGammaLow; isDefined=kTRUE; }
445  if(mypatch->IsJetHigh()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCJetLow : AliReducedPatchContainer::kEMCJetHigh; isDefined=kTRUE; }
446  if(mypatch->IsJetLow()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCJetHigh : AliReducedPatchContainer::kEMCJetLow; isDefined=kTRUE; }
447  if(!isDefined){
448  AliDebug(2, "Unknown online patch type");
449  continue;
450  }
451  AliDebug(2, Form("Adding online patch of type %d", int(triggertype)));
452  cont->AddTriggerPatch(kFALSE, triggertype, mypatch->GetPatchE(), mypatch->GetADCAmp(), mypatch->GetEtaGeo(), mypatch->GetPhiGeo());
453  }
454  }
455 }
456 
462 void AliReducedHighPtEventCreator::FixTrackInputEvent(AliVTrack* trk) {
463  if(!trk->GetEvent()){
464  if(trk->IsA() == AliESDtrack::Class())
465  (static_cast<AliESDtrack *>(trk))->SetESDEvent(static_cast<AliESDEvent *>(fInputEvent));
466  else if(trk->IsA() == AliAODTrack::Class())
467  (static_cast<AliAODTrack *>(trk))->SetAODEvent(static_cast<AliAODEvent *>(fInputEvent));
468  else if(trk->IsA() == AliPicoTrack::Class()){
469  AliPicoTrack *mytrk = static_cast<AliPicoTrack *>(trk);
470  if(!mytrk->GetEvent()){
471  if(mytrk->GetTrack()->IsA() == AliESDtrack::Class())
472  (static_cast<AliESDtrack *>(mytrk->GetTrack()))->SetESDEvent(static_cast<AliESDEvent *>(fInputEvent));
473  else
474  (static_cast<AliAODTrack *>(mytrk->GetTrack()))->SetAODEvent(static_cast<AliAODEvent *>(fInputEvent));
475  }
476  }
477  }
478 }
479 
483 AliReducedTrackSelectionContainer::AliReducedTrackSelectionContainer():
484  fIndex(-1),
485  fTrackSelection(NULL)
486 {
487 }
488 
494 AliReducedTrackSelectionContainer::AliReducedTrackSelectionContainer(
496  fIndex(index),
497  fTrackSelection(sel)
498 {
499 }
500 
504 AliReducedTrackSelectionContainer::~AliReducedTrackSelectionContainer() {
505  if(fTrackSelection) delete fTrackSelection;
506 }
507 
508 } /* namespace HighPtTracks */
509 
Reduced event structure for high- analysis.
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Declaration of class AliReducedHighPtEventCreator and AliReducedTrackSelectionContainer.
Declaration of class AlliReducedPatchInfo.
centrality
Declaration of class AliReducedPatchContainer, a container for reduced trigger patches.
TRandom * gRandom
Declaration of a reduced MC event header.
Main data structure storing all relevant information of EMCAL/DCAL trigger patches.
Declaration of class AliReducedGeneratedParticle.
AliVTrack * GetTrack() const
Definition: AliPicoTrack.h:59
Reduced information about reconstructed EMCAL clusters.
Class to make array of trigger patch objects in AOD/ESD events.
Declartion of class AliEMCalPtTaskVTrackSelection.