AliPhysics  64f4410 (64f4410)
AliHighPtReconstructionEfficiency.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 <iostream>
16 #include <sstream>
17 
18 #include <TFile.h>
19 #include <TKey.h>
20 #include <TList.h>
21 #include <TMath.h>
22 #include <TProfile.h>
23 #include <TNtuple.h>
24 #include <TString.h>
25 #include <TSystem.h>
26 #include <TTree.h>
27 
28 #include "AliAnalysisManager.h"
29 #include "AliAODMCParticle.h"
30 #include "AliAODTrack.h"
31 #include "AliESDtrackCuts.h"
32 #include "AliESDtrack.h"
33 #include "AliLog.h"
34 #include "AliMCEvent.h"
35 #include "AliMCParticle.h"
36 #include "AliTrackReference.h"
37 #include "AliVEvent.h"
38 #include "AliVParticle.h"
39 #include "AliVTrack.h"
40 
41 #include <fastjet/PseudoJet.hh>
42 #include <fastjet/JetDefinition.hh>
43 #include <fastjet/ClusterSequence.hh>
44 
46 #include "AliReducedJetEvent.h"
47 #include "AliReducedJetInfo.h"
48 #include "AliReducedJetParticle.h"
49 #include "AliReducedMatchedTrack.h"
50 
52 
56 
57 namespace HighPtTracks {
58 
62 AliHighPtReconstructionEfficiency::AliHighPtReconstructionEfficiency() :
64  fParticleMap(NULL),
65  fJetTree(NULL),
66  fJetEvent(NULL),
67  fMaxEtaJets(0.5),
68  fMaxEtaParticles(0.8),
69  fMinPtParticles(0.15),
70  fMaxDR(0.5),
71  fCrossSection(0.),
72  fNtrials(0),
73  fPtHardBin(0),
74  fTaskDebugMode(false)
75 {
76  memset(fTrackCuts, 0, sizeof(AliESDtrackCuts *) * 2);
77 }
78 
86  AliAnalysisTaskSE(name),
87  fParticleMap(NULL),
88  fJetTree(NULL),
89  fJetEvent(NULL),
90  fMaxEtaJets(0.5),
91  fMaxEtaParticles(0.8),
92  fMinPtParticles(0.15),
93  fMaxDR(0.5),
94  fCrossSection(0.),
95  fNtrials(0),
96  fPtHardBin(0),
97  fTaskDebugMode(false)
98 {
99  DefineOutput(1, TNtuple::Class());
100  memset(fTrackCuts, 0, sizeof(AliESDtrackCuts *) * 2);
101 }
102 
107  for(int icut = 0; icut < 2; icut++)
108  if(fTrackCuts[icut]) delete fTrackCuts[icut];
109  if(fParticleMap) delete fParticleMap;
110 }
111 
117  OpenFile(1);
118 
119  fJetTree = new TTree("JetEvent","A tree with jet information");
121  fJetTree->Branch("JetEvent", "AliReducedJetEvent", fJetEvent, 128000,0);
122 
123  PostData(1, fJetTree);
124 }
125 
132  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
133  if (!tree) {
134  AliError(Form("%s - UserNotify: No current tree!",GetName()));
135  return kFALSE;
136  }
137 
138  TFile *curfile = tree->GetCurrentFile();
139  if (!curfile) {
140  AliError(Form("%s - UserNotify: No current file!",GetName()));
141  return kFALSE;
142  }
143 
144  PythiaInfoFromFile(curfile->GetName(), fCrossSection, fNtrials, fPtHardBin);
145 
146  return kTRUE;
147 }
148 
162  TList listforjets;
163  SelectParticlesForJetfinding(listforjets);
164 
165  // Find Jets
166  AliVParticle *part(NULL);
167  TIter partIter(&listforjets);
168  std::vector<fastjet::PseudoJet> inputjets;
169  while((part = dynamic_cast<AliVParticle *>(partIter()))){
170  fastjet::PseudoJet inputpart(part->Px(), part->Py(), part->Pz(), part->E());
171  inputpart.set_user_index(part->GetLabel());
172  inputjets.push_back(inputpart);
173  }
174  fastjet::ClusterSequence jetfinder(inputjets, fastjet::JetDefinition(fastjet::antikt_algorithm, 0.4));
175  std::vector<fastjet::PseudoJet> recjets = fastjet::sorted_by_pt(jetfinder.inclusive_jets());
177 
179 
180  std::vector<AliReconstructedParticlePair> selectedParticles = SelectParticles();
181 
182  if(fTaskDebugMode){
183  std::stringstream debugmessage;
184  debugmessage << "Number of particles: gen[" << selectedParticles.size() << "], rec[" << fParticleMap->GetNumberOfParticles() << "]";
185  AliDebug(1, debugmessage.str().c_str());
186 
187  int nparticlesGen = 0, nparticlesRec = 0;
188  for(std::vector<AliReconstructedParticlePair>::iterator pairiter = selectedParticles.begin(); pairiter != selectedParticles.end(); ++pairiter){
189  if(pairiter->GetMCTrueParticle()) nparticlesGen++;
190  if(pairiter->GetRecTracks().GetNumberOfParticles()) nparticlesRec++;
191  }
192  AliDebug(1, Form("Among %d selected particles we find %d reconstructed particles.", nparticlesGen, nparticlesRec));
193  }
194 
195  AliReducedJetInfo *recjet(NULL);
196  for(std::vector<fastjet::PseudoJet>::iterator jetiter = recjets.begin(); jetiter != recjets.end(); ++jetiter){
197  if(TMath::Abs(jetiter->eta()) > fMaxEtaJets) continue;
198  recjet = new AliReducedJetInfo(jetiter->px(), jetiter->py(), jetiter->pz(), jetiter->E());
199  ProcessJet(recjet, selectedParticles);
200  ConvertConstituents(recjet, *jetiter);
202  }
203  fJetTree->Fill();
204  PostData(1, fJetTree);
205 }
206 
215 bool AliHighPtReconstructionEfficiency::IsSelected(const AliVTrack* const track, CutType_t cuttype) const {
216  /*
217  * Select reconstructed particle by applying track cuts
218  */
219  if(!fTrackCuts[cuttype]) return kTRUE;
220  const AliESDtrack *inputtrack = dynamic_cast<const AliESDtrack *>(track);
221  if(inputtrack){
222  return fTrackCuts[cuttype]->AcceptTrack(inputtrack);
223  } else {
224  // AOD track: Use copy
225  AliESDtrack copytrack(track);
226  return fTrackCuts[cuttype]->AcceptTrack(&copytrack);
227  }
228 }
229 
237 bool AliHighPtReconstructionEfficiency::IsTrueSelected(const AliVParticle* const track) const {
238  if(TMath::Abs(track->Eta()) > fMaxEtaParticles) return false;
239  if(TMath::Abs(track->Pt()) < fMinPtParticles) return false;
240  if(!track->Charge()) return false;
241  return true;
242 }
243 
252  AliVParticle *part = 0;
253  particles.Clear();
254  for(int ipart = 0; ipart < fMCEvent->GetNumberOfTracks(); ipart++){
255  part = fMCEvent->GetTrack(ipart);
256  if(!IsPhysicalPrimary(part)) continue;
257  if(TMath::Abs(part->Eta()) > fMaxEtaParticles) continue;
258  if(TMath::Abs(part->Pt()) < fMinPtParticles) continue;
259  particles.Add(part);
260  }
261 }
262 
273 std::vector<AliReconstructedParticlePair> AliHighPtReconstructionEfficiency::SelectParticles() const {
274  AliVParticle *part(NULL);
275  AliParticleList tracks;
276  std::vector<AliReconstructedParticlePair> result;
277  for(int ipart = 0; ipart <fMCEvent->GetNumberOfTracks(); ipart++){
278  part = fMCEvent->GetTrack(ipart);
279  if(!IsPhysicalPrimary(part)) continue;
280  if(!IsTrueSelected(part)) continue;
281  const AliParticleList *tmptracks = FindReconstructedParticleFast(ipart);
282  if(tmptracks) tracks = *tmptracks;
283  if(fTaskDebugMode){
284  std::stringstream debugmessage;
285  debugmessage << "Accepted particle with label " << ipart;
286  if(tracks.GetNumberOfParticles())
287  debugmessage << " - found reconstructed track";
288  else
289  debugmessage << " - did not find reconstructed track";
290  AliDebug(1, debugmessage.str().c_str());
291  }
292 
294  foundpair.SetMCTrueParticle(part);
295  foundpair.SetRecParticles(tracks);
296  result.push_back(foundpair);
297  }
298  return result;
299 }
300 
308 double AliHighPtReconstructionEfficiency::GetDR(const fastjet::PseudoJet& recjet, const AliVParticle* const inputtrack) const {
309  return recjet.delta_R(fastjet::PseudoJet(inputtrack->Px(), inputtrack->Py(), inputtrack->Pz(), inputtrack->E()));
310 }
311 
322  const std::vector<AliReconstructedParticlePair>& particles) const
323 {
324  double pvec[4];
325  recjet->GetPxPyPxE(pvec[0], pvec[1], pvec[2], pvec[3]);
326  fastjet::PseudoJet frecjet(pvec[0], pvec[1], pvec[2], pvec[3]); // For distance to the main jet axis
327  AliVParticle *mcpart(NULL);
328  const AliVTrack *rectrack(NULL);
329  for(std::vector<AliReconstructedParticlePair>::const_iterator partiter = particles.begin(); partiter != particles.end(); ++partiter){
330  mcpart = partiter->GetMCTrueParticle();
331  const AliParticleList &matchedTracks = partiter->GetRecTracks();
332  double dr = GetDR(frecjet, mcpart);
333  if(dr < fMaxDR){
334  // Create new Particle and add it to the jet reconstructed jet
335  AliReducedJetParticle * part = new AliReducedJetParticle(mcpart->Px(), mcpart->Py(), mcpart->Pz(), mcpart->E(), mcpart->PdgCode());
337  part->SetDistanceToMainJetAxis(dr);
338  for(int itrk = 0; itrk < matchedTracks.GetNumberOfParticles(); itrk++){
339  rectrack = matchedTracks.GetParticle(itrk);
340  AliReducedMatchedTrack *matchedTrack = new AliReducedMatchedTrack(rectrack->Px(), rectrack->Py(), rectrack->Pz());
341  matchedTrack->SetNumberOfClustersTPC(rectrack->GetTPCNcls());
342  if(rectrack->GetLabel() >= 0) matchedTrack->SetGoodTrackLabel(true);
343  // FilterTrackCuts
346  part->AddMatchedTrack(matchedTrack);
347  }
348  recjet->AddParticleInCone(part);
349  }
350  }
351 }
352 
357  if(fParticleMap) delete fParticleMap;
359  for(int ipart = 0; ipart < fInputEvent->GetNumberOfTracks(); ipart++){
360  AliVTrack *mytrack = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(ipart));
361  // Only apply a minimum selection: TPC refit and ITS refit
362  // check particle cuts on tree level
363  if(!((mytrack->GetStatus() & AliVTrack::kITSrefit) && (mytrack->GetStatus() & AliVTrack::kTPCrefit))) continue;
364  //if(!IsSelected(mytrack)) continue;
365  fParticleMap->AddParticle(mytrack);
366  }
367 }
368 
379  return fParticleMap->GetParticles(label);
380 }
381 
390 bool AliHighPtReconstructionEfficiency::IsPhysicalPrimary(const AliVParticle* const part) const {
391  const AliAODMCParticle *aodpart = dynamic_cast<const AliAODMCParticle *>(part);
392  if(aodpart) return aodpart->IsPhysicalPrimary(); // AOD MC particle
393  else if (dynamic_cast<const AliMCParticle *>(part)) return fMCEvent->IsPhysicalPrimary(part->GetLabel()); // ESD MC particle
394  return false;
395 }
396 
404 void AliHighPtReconstructionEfficiency::ConvertConstituents(AliReducedJetInfo* const recjet, const fastjet::PseudoJet& inputjet) {
405  std::vector<fastjet::PseudoJet> constituents = inputjet.constituents();
406  for(std::vector<fastjet::PseudoJet>::const_iterator citer = constituents.begin(); citer != constituents.end(); ++citer){
407  AliVParticle *mcpart = fMCEvent->GetTrack(citer->user_index());
408  AliReducedJetConstituent *jetconst = new AliReducedJetConstituent(citer->px(), citer->py(), citer->px(), citer->E(), mcpart->PdgCode());
409  recjet->AddConstituent(jetconst);
410  }
411 }
412 
424 bool AliHighPtReconstructionEfficiency::PythiaInfoFromFile(const char* currFile, double &fXsec, double &fTrials, int &pthard) const {
425  TString file(currFile);
426  fXsec = 0;
427  fTrials = 1;
428 
429  if (file.Contains(".zip#")) {
430  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
431  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
432  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
433  file.Replace(pos+1,pos2-pos1,"");
434  } else {
435  // not an archive take the basename....
436  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
437  }
438 
439  // Get the pt hard bin
440  TString strPthard(file);
441 
442  strPthard.Remove(strPthard.Last('/'));
443  strPthard.Remove(strPthard.Last('/'));
444  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
445  strPthard.Remove(0,strPthard.Last('/')+1);
446  if (strPthard.IsDec())
447  pthard = strPthard.Atoi();
448 
449  // problem that we cannot really test the existance of a file in a archive so we have to live with open error message from root
450  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
451 
452  if (!fxsec) {
453  // next trial fetch the histgram file
454  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
455  if (!fxsec) {
456  // not a severe condition but inciate that we have no information
457  return kFALSE;
458  } else {
459  // find the tlist we want to be independtent of the name so use the Tkey
460  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
461  if (!key) {
462  fxsec->Close();
463  return kFALSE;
464  }
465  TList *list = dynamic_cast<TList*>(key->ReadObj());
466  if (!list) {
467  fxsec->Close();
468  return kFALSE;
469  }
470  fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
471  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
472  fxsec->Close();
473  }
474  } else { // no tree pyxsec.root
475  TTree *xtree = (TTree*)fxsec->Get("Xsection");
476  if (!xtree) {
477  fxsec->Close();
478  return kFALSE;
479  }
480  UInt_t ntrials = 0;
481  Double_t xsection = 0;
482  xtree->SetBranchAddress("xsection",&xsection);
483  xtree->SetBranchAddress("ntrials",&ntrials);
484  xtree->GetEntry(0);
485  fTrials = ntrials;
486  fXsec = xsection;
487  fxsec->Close();
488  }
489  return kTRUE;
490  }
491 
498 unsigned short AliHighPtReconstructionEfficiency::GetNumberOfTPCTrackReferences(AliVParticle* const trk) const {
499  unsigned short nref = 0;
500  if(trk->IsA() == AliMCParticle::Class()){
501  AliMCParticle *part = dynamic_cast<AliMCParticle *>(trk);
502  AliTrackReference *tref(NULL);
503  for(int iref = 0; iref < part->GetNumberOfTrackReferences(); ++iref){
504  tref = part->GetTrackReference(iref);
505  if(tref->DetectorId() == AliTrackReference::kTPC) nref++;
506  }
507  }
508  return nref;
509 }
510 
511 }
Definition of class AliReducedJetEvent, an event structure containing reduced jet information...
Minimal stucture for jet constituents associated to a jet by the jet clustering algorithm.
double Double_t
Definition: External.C:58
double fMaxDR
maximum distance of a particle to the main jet axis
void SetNumberOfTPCtrackReferences(unsigned short nref)
Definintion of class AliReducedJetParticle, a structure for a reduced information set of particles as...
void AddParticle(AliVTrack *track)
TSystem * gSystem
AliParticleList * GetParticles(int label) const
void ConvertConstituents(AliReducedJetInfo *const recjet, const fastjet::PseudoJet &inputjet)
bool IsSelected(const AliVTrack *const track, CutType_t type) const
void SetSurvivedTrackCuts(ETrackCutsType_t cuts)
void AddMatchedTrack(AliReducedMatchedTrack *trk)
Event structure containing reduced jet information.
Container of reconstructed particles.
bool PythiaInfoFromFile(const char *currFile, double &fXsec, double &fTrials, int &pthard) const
void AddParticleInCone(AliReducedJetParticle *part)
Analysis task producing filtered trees with reconstructed jets at generator level.
void AddConstituent(AliReducedJetConstituent *con)
bool IsTrueSelected(const AliVParticle *const track) const
unsigned int UInt_t
Definition: External.C:33
unsigned short GetNumberOfTPCTrackReferences(AliVParticle *const trk) const
void ProcessJet(AliReducedJetInfo *const jet, const std::vector< AliReconstructedParticlePair > &particles) const
Namespace for classes creating trees of events with jets.
void SetMCTrueParticle(AliVParticle *const part)
Reduced information about a reconstructed jet.
Pair of a Monte-Carlo true particle and the associated reconstructed information. ...
void SetNumberOfClustersTPC(unsigned char ncls)
AliParticleMap * fParticleMap
Map of reconstructed particles associate to a Monte-Carlo label.
void GetPxPyPxE(double &px, double &py, double &pz, double &e)
double GetDR(const fastjet::PseudoJet &recjet, const AliVParticle *const inputtrack) const
Reduced information set of particles associated with a jet.
const AliParticleList * FindReconstructedParticleFast(int label) const
Find reconstructed particles for a given label.
Definition of class AliReducedTrack, a structure with reduced track information at reconstruction lev...
TFile * file
TList with histograms for a given trigger.
Definition of class AliReducedJetConstituent, a minimal stucture for jet constituents associated to a...
void AddReconstructedJet(AliReducedJetInfo *jet)
const char Option_t
Definition: External.C:48
Class with reduced track information at reconstruction level.
bool IsPhysicalPrimary(const AliVParticle *const part) const
Checks if the particle is a physical primary particle.
AliVTrack * GetParticle(int itrack) const
std::vector< AliReconstructedParticlePair > SelectParticles() const
Definition of class AliReducedJetInfo, a structure for reduced information about a reconstructed jet...
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
Definition of the analysis task producing filtered trees with reconstructed jets at generator level...
Map of reconstructed particles which share the same Monte-Carlo label.