AliPhysics  vAN-20150630 (513c479)
 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 <map>
16 #include <vector>
17 
18 #include <TArrayD.h>
19 #include <TArrayI.h>
20 #include <TClonesArray.h>
21 #include <TMath.h>
22 #include <TObjArray.h>
23 #include <TString.h>
24 #include <TTree.h>
25 
26 #include "AliAnalysisUtils.h"
27 #include "AliAODEvent.h"
28 #include "AliAODMCParticle.h"
29 #include "AliAODTrack.h"
30 #include "AliCentrality.h"
31 #include "AliESDEvent.h"
32 #include "AliESDtrack.h"
33 #include "AliInputEventHandler.h"
34 #include "AliLog.h"
35 #include "AliMCEvent.h"
36 #include "AliPicoTrack.h"
37 #include "AliStack.h"
38 #include "AliVCluster.h"
39 #include "AliVTrack.h"
40 #include "AliVVertex.h"
42 
44 #include "AliReducedEmcalCluster.h"
46 #include "AliReducedHighPtEvent.h"
47 #include "AliReducedMCHeader.h"
49 #include "AliReducedPatchInfo.h"
52 
54 ClassImp(HighPtTracks::AliReducedTrackSelectionContainer)
55 ClassImp(HighPtTracks::AliReducedHighPtEventCreator)
57 
58 namespace HighPtTracks {
59 
63 AliReducedHighPtEventCreator::AliReducedHighPtEventCreator():
65  fOutputTree(NULL),
66  fOutputEvent(NULL),
67  fTrackSelections(NULL),
68  fEventSelectionBits(AliVEvent::kAny),
69  fSwapTriggerThresholds(kFALSE),
70  fMinClusterE(-1),
71  fMaxClusterE(1000),
72  fMinPt(-1),
73  fMaxPt(1000),
74  fMinEta(-1000),
75  fMaxEta(1000),
76  fApplyCentralitySelection(kFALSE),
77  fCentralityMethod("V0A"),
78  fTriggerSetup("4classes"),
79  fMinBiasSelection(AliVEvent::kINT7)
80 {
81  fSelectCentralityRange[0] = 0.;
82  fSelectCentralityRange[1] = 100.;
83 }
84 
89 AliReducedHighPtEventCreator::AliReducedHighPtEventCreator(const char* name):
90  AliAnalysisTaskEmcal(name, kTRUE),
91  fOutputTree(NULL),
92  fOutputEvent(NULL),
93  fTrackSelections(NULL),
94  fEventSelectionBits(AliVEvent::kAny),
95  fSwapTriggerThresholds(kFALSE),
96  fMinClusterE(-1),
97  fMaxClusterE(1000),
98  fMinPt(-1),
99  fMaxPt(1000),
100  fMinEta(-1000),
101  fMaxEta(1000),
102  fApplyCentralitySelection(kFALSE),
103  fCentralityMethod("V0A"),
104  fTriggerSetup("4classes"),
105  fMinBiasSelection(AliVEvent::kINT7)
106 {
107  DefineOutput(2, TTree::Class());
108 
109  fSelectCentralityRange[0] = 0.;
110  fSelectCentralityRange[1] = 100.;
111 
112  fTrackSelections = new TObjArray;
113  fTrackSelections->SetOwner(kTRUE);
114 }
115 
119 AliReducedHighPtEventCreator::~AliReducedHighPtEventCreator() {
120  if(fTrackSelections) delete fTrackSelections;
121 }
122 
126 void AliReducedHighPtEventCreator::UserCreateOutputObjects() {
128  OpenFile(2);
129 
130  fOutputTree = new TTree("ReducedEvent", "A simple reduced event");
131  fOutputEvent = new AliReducedHighPtEvent();
132  fOutputTree->Branch("ReducedEvent", "AliReducedHighPtEvent", fOutputEvent, 128000, 0);
133 
134  PostData(1, fOutput);
135  PostData(2, fOutputTree);
136 }
137 
148 Bool_t AliReducedHighPtEventCreator::Run() {
149  AliDebug(1, "Starting creation of the reconstructed event");
150  if(!fCaloClusters){
151  AliError("Cluster container missing");
152  return kFALSE;
153  }
154  if(!fTracks){
155  AliError("Track container missing");
156  return kFALSE;
157  }
158  if(!(this->fInputHandler->IsEventSelected() & fEventSelectionBits)){
159  AliDebug(1, "Trigger not selected");
160  return kFALSE;
161  } else {
162  AliDebug(1, "Trigger selected");
163  }
164  if(!SelectEvent(fInputEvent)) return kFALSE;
165  new(fOutputEvent) AliReducedHighPtEvent(kTRUE);
166 
167  // Write event specific information
168  fOutputEvent->SetVertexZ(fInputEvent->GetPrimaryVertex()->GetZ());
169  AliCentrality *centralityHandler = fInputEvent->GetCentrality();
170  if(centralityHandler) fOutputEvent->SetCentralityPercentile(centralityHandler->GetCentralityPercentile(fCentralityMethod.Data()));
171  ConvertTriggerPatches(fTriggerPatchInfo, fOutputEvent->GetPatchContainer());
172  TString triggerString(fInputEvent->GetFiredTriggerClasses());
173  if(!fTriggerSetup.CompareTo("4classes"))
174  fOutputEvent->SetDecisionFromTriggerString(triggerString.Contains("EG1"), triggerString.Contains("EG2"), triggerString.Contains("EJ1"), triggerString.Contains("EJ2"));
175  else
176  fOutputEvent->SetDecisionFromTriggerString(triggerString.Contains("EGA"), kFALSE, triggerString.Contains("EJE"), kFALSE);
177  fOutputEvent->SetMinBiasEvent(fInputHandler->IsEventSelected() & fMinBiasSelection);
178  fOutputEvent->SetRunNumber(fInputEvent->GetRunNumber());
179 
180  std::map<int, int> mcindexmap;
181  if(MCEvent()){
182  // Generate Monte-Carlo header
183  AliReducedMCHeader *mcheader = new AliReducedMCHeader();
184  if(fIsPythia){
185  mcheader->SetCrossSection(fXsection);
186  mcheader->SetNumberOfTrials(fNTrials);
187  mcheader->SetPtHard(fPtHard);
188  }
189  fOutputEvent->SetMonteCarloHeader(mcheader);
190 
191  // Convert Monte Carlo particles, fill index map
192  AliMCEvent *mcev = MCEvent();
193  int npart(0);
194  for(Int_t ipart = 0; ipart < mcev->GetNumberOfTracks(); ipart++){
195  AliVParticle *part = mcev->GetTrack(ipart);
196  Double_t pt(TMath::Abs(part->Pt())), eta(part->Eta());
197  if(pt < fMinPt || pt > fMaxPt) continue;
198  if(eta < fMinEta || eta > fMaxEta) continue;
199  if(!part->Charge()) continue;
200  if(part->IsA() == AliAODMCParticle::Class()){
201  AliAODMCParticle *aodpart = static_cast<AliAODMCParticle *>(part);
202  if(!aodpart->IsPhysicalPrimary()) continue;
203  } else {
204  // ESD - get information from stack
205  if(!mcev->Stack()->IsPhysicalPrimary(ipart)) continue;
206  }
207 
208  AliReducedGeneratedParticle *reducedgen = new AliReducedGeneratedParticle(npart, part->PdgCode(), part->Px(), part->Py(), part->Pz(), part->E());
209  fOutputEvent->AddReducedGeneratedParticle(reducedgen);
210  mcindexmap.insert(std::pair<int,int>(ipart, npart));
211  npart++;
212  }
213  }
214 
215  // Filter Clusters
216  std::map<int, int> clusterindexmap;
217  int ncluster(0);
218  Double_t vtxpos[3];fInputEvent->GetPrimaryVertex()->GetXYZ(vtxpos);
219  for(TIter cliter = TIter(fCaloClusters).Begin(); cliter != TIter::End(); ++cliter){
220  AliVCluster *incluster = static_cast<AliVCluster *>(*cliter);
221  if(!SelectCluster(incluster)) continue;
222  TLorentzVector clustervec;
223  incluster->GetMomentum(clustervec, vtxpos);
224  AliReducedEmcalCluster *redcluster = new AliReducedEmcalCluster(ncluster, incluster->E(), clustervec.Eta(), clustervec.Phi(), incluster->GetM02(), incluster->GetM20());
225  // Get leading 3 cell energies
226  TArrayD cellEnergies;
227  GetCellEnergies(incluster, cellEnergies);
228  TArrayI indices(cellEnergies.GetSize());
229  TMath::Sort(cellEnergies.GetSize(), cellEnergies.GetArray(), indices.GetArray(), kTRUE);
230  redcluster->SetLeadingCellEnergies(
231  cellEnergies[indices[0]],
232  cellEnergies.GetSize() > 1 ? cellEnergies[indices[1]] : 0,
233  cellEnergies.GetSize() > 2 ? cellEnergies[indices[2]] : 0
234  );
235  // Assing MC particles
236  if(MCEvent()){
237  for(Int_t ilab = 0; ilab < incluster->GetNLabels(); ilab++){
238  AliVParticle *assigned = MCEvent()->GetTrack(TMath::Abs(incluster->GetLabels()[ilab]));
239  if(!assigned) continue;
240  redcluster->AddTrueContributor(assigned->PdgCode(), assigned->Px(), assigned->Py(), assigned->Pz(), assigned->E());
241  }
242  }
243  fOutputEvent->AddReducedCluster(redcluster);
244  clusterindexmap.insert(std::pair<int,int>(incluster->GetID(), ncluster));
245  ncluster++;
246  }
247 
248  // Filter tracks
249  int nsel(0);
250  AliDebug(1, Form("Number of tracks in container: %d", fTracks->GetEntries()));
251  for(TIter trkiter = TIter(fTracks).Begin(); trkiter != TIter::End(); ++trkiter){
252  AliVTrack *trackRaw = static_cast<AliVTrack *>(*trkiter), *trackToCheck(NULL);
253  // handlinf for pico tracks as input
254  if(trackRaw->IsA() == AliPicoTrack::Class()){
255  AliPicoTrack *picotrack = static_cast<AliPicoTrack *>(trackRaw);
256  trackToCheck = picotrack->GetTrack();
257  } else {
258  trackToCheck = trackRaw;
259  }
260  FixTrackInputEvent(trackToCheck);
261  TArrayI cutIndices;
262  if(!SelectTrack(trackToCheck,cutIndices)) continue;
263  AliReducedReconstructedTrack *rectrack = new AliReducedReconstructedTrack;
264  double pvec[3];
265  trackToCheck->GetPxPyPz(pvec);
266  rectrack->SetMomentumVector(pvec[0], pvec[1], pvec[2]);
267  rectrack->SetCharge(static_cast<Char_t>(trackToCheck->Charge()));
268  for(Int_t icut = 0; icut < cutIndices.GetSize(); icut++) rectrack->SetTrackCuts(cutIndices[icut]);
269  if(MCEvent()){
270  // handle label
271  Int_t label = TMath::Abs(trackToCheck->GetLabel());
272  if(label < fMCEvent->GetNumberOfTracks()){
273  std::map<int,int>::iterator found = mcindexmap.find(label);
274  if(found != mcindexmap.end()){
275  rectrack->SetMatchedParticleIndex(found->second);
276  }
277  if(label >= 0) rectrack->SetGoodTrackLabel(kTRUE);
278  rectrack->SetTPCClusters(trackToCheck->GetTPCNcls());
279  rectrack->SetTPCCrossedRows(GetTPCCrossedRows(trackToCheck));
280  rectrack->SetTPCSharedClusters(trackToCheck->GetTPCSharedMapPtr()->CountBits());
281  rectrack->SetTPCFindableClusters(trackToCheck->GetTPCNclsF());
282  }
283  }
284  // handle clusters
285  Int_t clusterIndex = trackToCheck->GetEMCALcluster();
286  if(clusterIndex >= 0 && clusterIndex < fCaloClusters->GetEntries()){
287  std::map<int,int>::iterator found = clusterindexmap.find(clusterIndex);
288  if(found != clusterindexmap.end()){
289  rectrack->SetMatchedClusterIndex(found->second);
290  }
291  }
292  //std::cout << "Track selected" << std::endl;
293  fOutputEvent->AddReducedReconstructedParticle(rectrack);
294  nsel++;
295  }
296  AliDebug(1, Form("Number of selected tracks :%d ", nsel ));
297 
298  fOutputTree->Fill();
299  PostData(2, fOutputTree);
300  return kTRUE;
301 }
302 
308 void AliReducedHighPtEventCreator::AddVirtualTrackSelection(
310  fTrackSelections->Add(new AliReducedTrackSelectionContainer(cutindex, sel));
311 }
312 
318 Bool_t AliReducedHighPtEventCreator::SelectEvent(AliVEvent* event) const {
319  if(fApplyCentralitySelection && InputEvent()->GetCentrality()){
320  AliDebug(1, Form("Centrality selection applied in range %f - %f\n", fSelectCentralityRange[0], fSelectCentralityRange[1]));
321  Float_t centrality = InputEvent()->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
322  AliDebug(1, Form("Centrality before: %f, estimator %s\n", centrality, fCentralityMethod.Data()));
323  if(centrality < fSelectCentralityRange[0] || centrality > fSelectCentralityRange[1]) return kFALSE;
324  AliDebug(1, "Centrality selected");
325  }
326  const AliVVertex *primvtx = event->GetPrimaryVertex();
327  if(!primvtx || !primvtx->GetNContributors()) return kFALSE;
328  if(TMath::Abs(primvtx->GetZ()) > 10.) return kFALSE;
329  if(event->IsPileupFromSPD(3, 0.8, 3., 2., 5.)) return kFALSE;
330  AliAnalysisUtils eventSelUtil;
331  if(!eventSelUtil.IsVertexSelected2013pA(event)) return kFALSE;
332  return kTRUE;
333 }
334 
340 Bool_t AliReducedHighPtEventCreator::SelectCluster(const AliVCluster* clust) const {
341  if(!clust->IsEMCAL()) return kFALSE;
342  if(clust->E() < fMinClusterE || clust->E() > fMaxClusterE) return kFALSE;
343  return kTRUE;
344 }
345 
352 Int_t AliReducedHighPtEventCreator::SelectTrack(AliVTrack* track, TArrayI& cutindices) const {
353  std::vector<int> selected;
354 
355  double pt = TMath::Abs(track->Pt()), eta = track->Eta();
356  if(pt < fMinPt || pt > fMaxPt) return 0;
357  if(eta < fMinEta || eta > fMaxEta) return 0;
358 
359  for(TIter cutiter = TIter(fTrackSelections).Begin(); cutiter != TIter::End(); ++cutiter){
360  AliReducedTrackSelectionContainer *mycut = static_cast<AliReducedTrackSelectionContainer *>(*cutiter);
361  if(!mycut->GetTrackSelection()->IsTrackAccepted(track)) continue;
362  selected.push_back(mycut->GetIndex());
363  }
364  if(!selected.size()) return 0;
365  cutindices.Set(selected.size());
366  int counter = 0;
367  for(std::vector<int>::iterator inditer = selected.begin(); inditer != selected.end(); ++inditer){
368  cutindices[counter++] = *inditer;
369  }
370  return selected.size();
371 }
372 
378 Int_t AliReducedHighPtEventCreator::GetTPCCrossedRows(const AliVTrack* trk) const {
379  if(trk->IsA() == AliESDtrack::Class()){
380  return (static_cast<const AliESDtrack *>(trk))->GetTPCCrossedRows();
381  } else if(trk->IsA() == AliAODTrack::Class()){
382  return (static_cast<const AliAODTrack *>(trk))->GetTPCNCrossedRows();
383  }
384  return 0;
385 }
386 
392 void AliReducedHighPtEventCreator::GetCellEnergies(AliVCluster* emccluster, TArrayD& energies) const {
393  if(!fInputEvent->GetEMCALCells()) {
394  AliError("No EMCAL cells");
395  return;
396  }
397  AliDebug(2, Form("Number of cells: %d, array: %p", emccluster->GetNCells(), emccluster->GetCellsAbsId()));
398  energies.Set(emccluster->GetNCells());
399  for(int icell = 0; icell < emccluster->GetNCells(); icell++){
400  // printf("Cell ID: %d\n", emccluster->GetCellsAbsId()[icell]);
401  energies[icell] = fInputEvent->GetEMCALCells()->GetCellAmplitude(emccluster->GetCellsAbsId()[icell]);
402  }
403 }
404 
410 void AliReducedHighPtEventCreator::ConvertTriggerPatches(TClonesArray* patches,
411  AliReducedPatchContainer* cont) {
412  if(!patches){
413  AliError("Trigger patch container not found\n");
414  return;
415  }
416  for(TIter patchIter = TIter(patches).Begin(); patchIter != TIter::End(); ++patchIter){
417  AliEmcalTriggerPatchInfo *mypatch = static_cast<AliEmcalTriggerPatchInfo *>(*patchIter);
418  if(!mypatch->IsOfflineSimple() && mypatch->IsLevel0()) continue;
420  Bool_t isDefined(kFALSE);
421  if(mypatch->IsOfflineSimple()){
422  if(mypatch->IsGammaHighSimple()) { triggertype = AliReducedPatchContainer::kEMCGammaHigh; isDefined = kTRUE; }
423  if(mypatch->IsGammaLowSimple()) { triggertype = AliReducedPatchContainer::kEMCGammaLow; isDefined = kTRUE; }
424  if(mypatch->IsJetHighSimple()) { triggertype = AliReducedPatchContainer::kEMCJetHigh; isDefined = kTRUE; }
425  if(mypatch->IsJetLowSimple()) { triggertype = AliReducedPatchContainer::kEMCJetLow; isDefined = kTRUE; }
426  if(!isDefined){
427  AliDebug(2, "Unknown offline patch type");
428  continue;
429  }
430  AliDebug(2, Form("Adding offline patch of type %d", int(triggertype)));
431  cont->AddTriggerPatch(kTRUE, triggertype, mypatch->GetPatchE(), mypatch->GetADCAmp(), mypatch->GetEtaGeo(), mypatch->GetPhiGeo());
432  } else {
433  if(mypatch->IsGammaHigh()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCGammaLow : AliReducedPatchContainer::kEMCGammaHigh; isDefined=kTRUE; }
434  if(mypatch->IsGammaLow()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCGammaHigh : AliReducedPatchContainer::kEMCGammaLow; isDefined=kTRUE; }
435  if(mypatch->IsJetHigh()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCJetLow : AliReducedPatchContainer::kEMCJetHigh; isDefined=kTRUE; }
436  if(mypatch->IsJetLow()) { triggertype = fSwapTriggerThresholds ? AliReducedPatchContainer::kEMCJetHigh : AliReducedPatchContainer::kEMCJetLow; isDefined=kTRUE; }
437  if(!isDefined){
438  AliDebug(2, "Unknown online patch type");
439  continue;
440  }
441  AliDebug(2, Form("Adding online patch of type %d", int(triggertype)));
442  cont->AddTriggerPatch(kFALSE, triggertype, mypatch->GetPatchE(), mypatch->GetADCAmp(), mypatch->GetEtaGeo(), mypatch->GetPhiGeo());
443  }
444  }
445 }
446 
452 void AliReducedHighPtEventCreator::FixTrackInputEvent(AliVTrack* trk) {
453  if(!trk->GetEvent()){
454  if(trk->IsA() == AliESDtrack::Class())
455  (static_cast<AliESDtrack *>(trk))->SetESDEvent(static_cast<AliESDEvent *>(fInputEvent));
456  else if(trk->IsA() == AliAODTrack::Class())
457  (static_cast<AliAODTrack *>(trk))->SetAODEvent(static_cast<AliAODEvent *>(fInputEvent));
458  else if(trk->IsA() == AliPicoTrack::Class()){
459  AliPicoTrack *mytrk = static_cast<AliPicoTrack *>(trk);
460  if(!mytrk->GetEvent()){
461  if(mytrk->GetTrack()->IsA() == AliESDtrack::Class())
462  (static_cast<AliESDtrack *>(mytrk->GetTrack()))->SetESDEvent(static_cast<AliESDEvent *>(fInputEvent));
463  else
464  (static_cast<AliAODTrack *>(mytrk->GetTrack()))->SetAODEvent(static_cast<AliAODEvent *>(fInputEvent));
465  }
466  }
467  }
468 }
469 
473 AliReducedTrackSelectionContainer::AliReducedTrackSelectionContainer():
474  fIndex(-1),
475  fTrackSelection(NULL)
476 {
477 }
478 
484 AliReducedTrackSelectionContainer::AliReducedTrackSelectionContainer(
486  fIndex(index),
487  fTrackSelection(sel)
488 {
489 }
490 
494 AliReducedTrackSelectionContainer::~AliReducedTrackSelectionContainer() {
495  if(fTrackSelection) delete fTrackSelection;
496 }
497 
498 } /* namespace HighPtTracks */
499 
Reduced event structure for high- analysis.
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Declaration of class AliReducedHighPtEventCreator and AliReducedTrackSelectionContainer.
Declaration of class AlliReducedPatchInfo.
Declaration of class AliReducedPatchContainer, a container for reduced trigger patches.
Declaration of a reduced MC event header.
Declaration of class AliReducedGeneratedParticle.
AliVTrack * GetTrack() const
Definition: AliPicoTrack.h:59
Reduced information about reconstructed EMCAL clusters.
Declartion of class AliEMCalPtTaskVTrackSelection.