AliPhysics  6cf2591 (6cf2591)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliEmcalMCTreeWriter.cxx
Go to the documentation of this file.
1 /*
2  * AliEmcalMCTreeWriter.cxx
3  *
4  * Created on: 03.10.2014
5  * Author: markusfasel
6  */
7 
8 #include <iostream>
9 #include <string>
10 #include <TArrayI.h>
11 #include <TClonesArray.h>
12 #include <TLorentzVector.h>
13 #include <TPDGCode.h>
14 #include <TMath.h>
15 #include <TTree.h>
16 #include <TVector3.h>
17 
18 #include "AliAODMCParticle.h"
19 #include "AliMCEvent.h"
20 #include "AliVParticle.h"
21 #include "AliVTrack.h"
22 #include "AliVCluster.h"
23 
24 #include "AliEmcalMCTreeWriter.h"
25 
28  fOutputTree(NULL),
29  fOutputInfo()
30 {
31  /*
32  * Dummy constructor
33  */
34 }
35 
37  AliAnalysisTaskEmcal(name,true),
38  fOutputTree(NULL),
39  fOutputInfo()
40 {
41  /*
42  * Constructor
43  */
44  DefineOutput(2,TTree::Class());
45 }
46 
48  /*
49  * Destructor
50  */
51  if(fOutputTree) delete fOutputTree;
52 }
53 
55  /*
56  * Create output tree, with two branches, one for the tracks matched and one for the clusters
57  */
59 
60  // set the listnames
61  std::string clusterlist(""), tracklist("");
62  if(fIsEsd){
63  tracklist = "ESDFilterTracks";
64  //clusterlist = "EmcCaloClusters";
65  clusterlist = "CaloClustersCorr";
66  } else {
67  tracklist = "AODFilterTracks";
68  //clusterlist = "EmcCaloClusters";
69  clusterlist = "CaloClustersCorr";
70  }
71  this->SetTracksName(tracklist.c_str());
72  this->SetClusName(clusterlist.c_str());
73 
74  // Build the tree
75  OpenFile(2);
76  fOutputTree = new TTree("EMCalTree", "A tree with emcal information");
77  fOutputTree->Branch("pdg", &fOutputInfo.pdg, "pdg/I");
78  fOutputTree->Branch("energy", &fOutputInfo.E, "energy/D");
79  fOutputTree->Branch("energyGen", &fOutputInfo.Egen, "energyGen/D");
80  fOutputTree->Branch("isUnique", &fOutputInfo.isUnique, "isUnique/B");
81  fOutputTree->Branch("isPhysicalPrimary", &fOutputInfo.isPhysicalPrimary, "isPhysicalPrimary/B");
82  fOutputTree->Branch("pgx", &fOutputInfo.pgen[0], "pgx/D");
83  fOutputTree->Branch("pgy", &fOutputInfo.pgen[1], "pgy/D");
84  fOutputTree->Branch("pgz", &fOutputInfo.pgen[2], "pgz/D");
85  fOutputTree->Branch("prx", &fOutputInfo.prec[0], "prx/D");
86  fOutputTree->Branch("pry", &fOutputInfo.prec[1], "pry/D");
87  fOutputTree->Branch("prz", &fOutputInfo.prec[2], "prz/D");
88  fOutputTree->Branch("erx", &fOutputInfo.erec[0], "erx/D");
89  fOutputTree->Branch("ery", &fOutputInfo.erec[1], "ery/D");
90  fOutputTree->Branch("erz", &fOutputInfo.erec[2], "erz/D");
91  fOutputTree->Branch("m02", &fOutputInfo.showershape[0], "m02/D");
92  fOutputTree->Branch("m20", &fOutputInfo.showershape[1], "m20/D");
93  fOutputTree->Branch("cells", "TVectorD", &fOutputInfo.fCellEnergies);
94  fOutputTree->Branch("indices", "TVectorD", &fOutputInfo.fCellIndices);
95  PostData(1, fOutput);
96  PostData(2, fOutputTree);
97 }
98 
100  /*
101  * Build the tree
102  */
103 
104  const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
105  if(!vertex) return kFALSE;
106 
107  Double_t vpos[3];
108  vertex->GetXYZ(vpos);
109  TIter trackIter(fTracks);
110  AliVTrack *rectrack(NULL);
111  AliVParticle *particle(NULL);
112  AliAODMCParticle *aodmc(NULL);
113  AliVCluster *assoccluster(NULL);
114  for(unsigned int ipart = 0; ipart < static_cast<unsigned int>(fMCEvent->GetNumberOfTracks()); ipart++){
115  particle = fMCEvent->GetTrack(static_cast<int>(ipart));
116  if(!AcceptParticle(particle)) continue;
117  if(TMath::Abs(particle->Eta()) > 2) continue; // Reject particles far outside the alice acceptance
118 
119  fOutputInfo.Reset();
120  fOutputInfo.pdg = particle->PdgCode();
121  if((aodmc = dynamic_cast<AliAODMCParticle *>(particle)))
122  fOutputInfo.isPhysicalPrimary = aodmc->IsPhysicalPrimary();
123  else
124  fOutputInfo.isPhysicalPrimary = fMCEvent->IsPhysicalPrimary(static_cast<int>(ipart));
125  particle->PxPyPz(fOutputInfo.pgen);
126  fOutputInfo.Egen = particle->E();
127 
128  // First try to find a matching track
129  std::vector<AliVTrack *> tracks;
130  FindTracks(ipart, tracks);
131  if(tracks.size() > 1) fOutputInfo.isUnique = false;
132  for(std::vector<AliVTrack *>::iterator trackIter = tracks.begin(); trackIter != tracks.end(); trackIter++){
133  rectrack = *trackIter;
134  rectrack->GetPxPyPz(fOutputInfo.prec);
135  memset(fOutputInfo.erec, 0, sizeof(Double_t) * 3);
136  fOutputInfo.E = 0;
137  assoccluster = dynamic_cast<AliVCluster *>(fCaloClusters->At(rectrack->GetEMCALcluster()));
138  if(assoccluster){
139  fOutputInfo.E = GetClusterEnergy(assoccluster, vpos, fOutputInfo.erec);
140  GetShowerShape(assoccluster, fOutputInfo.showershape);
142  }
143  fOutputTree->Fill();
144  }
145 
146  if(!tracks.size()){
147  // No track found , neutral particle? Look for EMCal clusters
148  std::vector<AliVCluster *> clusters;
149  FindClusters(ipart, clusters);
150  if(clusters.size() > 1) fOutputInfo.isUnique = false;
151  for(std::vector<AliVCluster *>::iterator clustIter = clusters.begin(); clustIter != clusters.end(); clustIter++){
152  assoccluster = *clustIter;
153  fOutputInfo.E = GetClusterEnergy(assoccluster, vpos, fOutputInfo.erec);
154  GetShowerShape(assoccluster, fOutputInfo.showershape);
156  fOutputTree->Fill();
157  }
158  }
159  }
160 
161  PostData(2, fOutputTree);
162  return kTRUE;
163 }
164 
165 double AliEmcalMCTreeWriter::GetClusterEnergy(AliVCluster* clust, Double_t *vpos, Double_t* evec) const {
166  /*
167  * Get energy from a cluster. Fill also as directed quantity
168  */
169  TLorentzVector clustervector;
170  clust->GetMomentum(clustervector, vpos);
171  TVector3 vec = clustervector.Vect().Unit();
172  vec *= clustervector.E();
173  if(evec)
174  vec.GetXYZ(evec);
175  return clustervector.E();
176 }
177 
178 void AliEmcalMCTreeWriter::FindTracks(unsigned int label, std::vector<AliVTrack*>& tracks) {
179  /*
180  * Reverse order finding: Find all tracks matched to a particle
181  */
182  tracks.clear();
183  TIter trackIter(fTracks);
184  AliVTrack *track;
185  while((track = dynamic_cast<AliVTrack *>(trackIter()))){
186  if(static_cast<unsigned int>(TMath::Abs(track->GetLabel())) == label){
187  // we have found a matching track
188  tracks.push_back(track);
189  }
190  }
191 }
192 
193 void AliEmcalMCTreeWriter::FindClusters(unsigned int label, std::vector<AliVCluster*>& clusters) {
194  /*
195  * Reverse order: Find all clusters matched to a particle
196  */
197  clusters.clear();
198  TIter clusterIter(fCaloClusters);
199  AliVCluster *cluster;
200  while((cluster = dynamic_cast<AliVCluster *>(clusterIter()))){
201  TArrayI labels;
202  labels.Set(cluster->GetNLabels(), cluster->GetLabels());
203  for(int *labIter = labels.GetArray(); labIter < labels.GetArray() + labels.GetSize(); labIter++){
204  if(static_cast<unsigned int>(TMath::Abs(*labIter)) == label){
205  // we have found a matching track
206  clusters.push_back(cluster);
207  break;
208  }
209  }
210  }
211 }
212 
213 void AliEmcalMCTreeWriter::GetShowerShape(AliVCluster* clust, Double_t* vector) const {
214  vector[0] = clust->GetM02();
215  vector[1] = clust->GetM20();
216 }
217 
218 void AliEmcalMCTreeWriter::GetCellEnergies(AliVCluster* clust, TVectorD& cells, TVectorF &indices) const {
219  cells.ResizeTo(clust->GetNCells());
220  indices.ResizeTo(clust->GetNCells());
221  UShort_t *cellIndices = clust->GetCellsAbsId();
222  for(int icell = 0; icell < clust->GetNCells(); icell++){
223  indices[icell] = cellIndices[icell];
224  cells[icell] = fInputEvent->GetEMCALCells()->GetCellAmplitude(cellIndices[icell]);
225  }
226 }
227 
228 bool AliEmcalMCTreeWriter::AcceptParticle(const AliVParticle* const particle) const {
229  /*
230  * accept particle in case it is a photon, pi0, charged pion, kaon, proton or electron
231  */
232  const int kAcceptPdg[7] = {kGamma, kPi0, kPiPlus, kKPlus, kProton, kElectron};
233  bool result = false;
234  for(const int *pdgiter = kAcceptPdg; pdgiter < kAcceptPdg + 7; pdgiter++){
235  if(TMath::Abs(particle->PdgCode()) == *pdgiter){
236  result = true;
237  break;
238  }
239  }
240  return result;
241 }
double GetClusterEnergy(AliVCluster *clust, Double_t *vpos, Double_t *evec) const
double Double_t
Definition: External.C:58
bool AcceptParticle(const AliVParticle *const particle) const
Base task in the EMCAL framework.
void FindTracks(unsigned int label, std::vector< AliVTrack * > &tracks)
TClonesArray * fCaloClusters
!clusters
TrackInfo fOutputInfo
Output tree with tracks.
void GetShowerShape(AliVCluster *clust, Double_t *vector) const
AliEmcalList * fOutput
!output list
TClonesArray * fTracks
!tracks
Bool_t fIsEsd
!whether it's an ESD analysis
unsigned short UShort_t
Definition: External.C:28
void GetCellEnergies(AliVCluster *clust, TVectorD &cells, TVectorF &indices) const
bool Bool_t
Definition: External.C:53
virtual void UserCreateOutputObjects()
void FindClusters(unsigned int label, std::vector< AliVCluster * > &clusters)
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65