AliPhysics  ec7afe5 (ec7afe5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskTracksInJet.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 <TChain.h>
16 #include <TDatabasePDG.h>
17 #include <TFile.h>
18 #include <THashList.h>
19 #include <THistManager.h>
20 #include <TKey.h>
21 #include <TList.h>
22 #include <TMath.h>
23 #include <TParticlePDG.h>
24 #include <TPDGCode.h>
25 #include <TProfile.h>
26 #include <TSystem.h>
27 #include <TTree.h>
28 
29 #include "AliAnalysisManager.h"
30 #include "AliAnalysisUtils.h"
31 #include "AliAODHeader.h"
32 #include "AliAODMCHeader.h"
33 #include "AliAODMCParticle.h"
34 #include "AliAODTrack.h"
35 #include "AliESDtrack.h"
36 #include "AliESDtrackCuts.h"
37 #include "AliExternalTrackParam.h"
38 #include "AliGenPythiaEventHeader.h"
39 #include "AliLog.h"
40 #include "AliMCEvent.h"
41 #include "AliMCParticle.h"
42 #include "AliVParticle.h"
43 #include "AliVTrack.h"
44 #include "AliVEvent.h"
45 
46 #include "fastjet/PseudoJet.hh"
47 #include "fastjet/ClusterSequence.hh"
48 #include "fastjet/JetDefinition.hh"
49 
51 
53 
54 namespace EMCalTriggerPtAnalysis {
55 
59 AliAnalysisTaskTracksInJet::AliAnalysisTaskTracksInJet() :
61  fJetStructure(),
62  fJetTree(NULL),
63  fAnalysisUtils(NULL),
64  fTrackCutsDefault(NULL),
65  fHybridCutsCat1(NULL),
66  fHybridCutsCat2(NULL),
67  fIsMC(kFALSE),
68  fFracPtHard(-1),
69  fHistosMC(NULL)
70 {
71 }
72 
73 AliAnalysisTaskTracksInJet::AliAnalysisTaskTracksInJet(const char *taskname) :
74  AliAnalysisTaskSE(taskname),
75  fJetStructure(),
76  fJetTree(NULL),
77  fAnalysisUtils(NULL),
78  fTrackCutsDefault(NULL),
79  fHybridCutsCat1(NULL),
80  fHybridCutsCat2(NULL),
81  fIsMC(kFALSE),
82  fFracPtHard(-1),
83  fHistosMC(NULL)
84 {
85  DefineOutput(1, TList::Class());
86  DefineOutput(2, TTree::Class());
87 }
88 
89 AliAnalysisTaskTracksInJet::~AliAnalysisTaskTracksInJet() {}
90 
91 void AliAnalysisTaskTracksInJet::UserCreateOutputObjects(){
92  fHistosMC = new THistManager("MCHistos");
93  fHistosMC->CreateTH1("hNtrials", "Number of trials", 1., 0.5, 1.5);
94  fHistosMC->CreateTH1("hNtrialsEvent", "Number of trials (from header, after event selection)", 1, 0.5, 1.5);
95  fHistosMC->CreateTH1("hPtHard", "Pt of the hard interaction", 1000, 0., 500);
96  fHistosMC->CreateTProfile("hCrossSection", "PYTHIA Cross section", 1, 0.5, 1.5);
97 
98  OpenFile(2);
99  fJetTree = new TTree("JetTreeRec", "Jet Tree on reconstructed data");
100  fJetTree->Branch("JetsRec", &fJetStructure, "pjet[3]/D:plead[3]/D:psub[3]/D:isData/I");
101 
102  fAnalysisUtils = new AliAnalysisUtils;
103 
104  fTrackCutsDefault = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(true, 1);
105  fTrackCutsDefault->SetName("Standard Track cuts");
106  fTrackCutsDefault->SetMinNCrossedRowsTPC(120);
107  fTrackCutsDefault->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
108 
109  // First class of hybrid track cuts
110  fHybridCutsCat1 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE,1);
111  fHybridCutsCat1->SetMaxDCAToVertexXY(2.4);
112  fHybridCutsCat1->SetMaxDCAToVertexZ(3.2);
113  fHybridCutsCat1->SetDCAToVertex2D(kTRUE);
114  fHybridCutsCat1->SetMaxChi2TPCConstrainedGlobal(36);
115  fHybridCutsCat1->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
116  fHybridCutsCat1->SetMaxFractionSharedTPCClusters(0.4);
117 
118  // Second class of hybrid track cuts
119  fHybridCutsCat2 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE,1);
120  fHybridCutsCat2->SetMaxDCAToVertexXY(2.4);
121  fHybridCutsCat2->SetMaxDCAToVertexZ(3.2);
122  fHybridCutsCat2->SetDCAToVertex2D(kTRUE);
123  fHybridCutsCat2->SetMaxChi2TPCConstrainedGlobal(36);
124  fHybridCutsCat2->SetMaxFractionSharedTPCClusters(0.4);
125  fHybridCutsCat2->SetRequireITSRefit(kFALSE);
126  fHybridCutsCat2->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
127 
128  PostData(1, fHistosMC->GetListOfHistograms());
129  PostData(2, fJetTree);
130 }
131 
136 Bool_t AliAnalysisTaskTracksInJet::UserNotify()
137 {
138  // Called when file changes.
139 
140  if(fIsMC){
141  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
142  if (!tree) {
143  AliError(Form("%s - UserNotify: No current tree!",GetName()));
144  return kFALSE;
145  }
146 
147  Float_t xsection = 0;
148  Float_t trials = 0;
149  Int_t pthardbin = 0;
150 
151  TFile *curfile = tree->GetCurrentFile();
152  if (!curfile) {
153  AliError(Form("%s - UserNotify: No current file!",GetName()));
154  return kFALSE;
155  }
156 
157  TChain *chain = dynamic_cast<TChain*>(tree);
158  if (chain) tree = chain->GetTree();
159 
160  Int_t nevents = tree->GetEntriesFast();
161 
162  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
163 
164  fHistosMC->FillTH1("hNtrials", 1., trials);
165  fHistosMC->FillProfile("hCrossSection", 1., xsection);
166  //fHistos->FillTH1(pthardbin, nevents);
167  }
168 
169  return kTRUE;
170 }
171 
176 void AliAnalysisTaskTracksInJet::UserExec(Option_t *){
177  // 1st Select event
178  if(!fAnalysisUtils->IsVertexSelected2013pA(InputEvent())) return;
179  if(fAnalysisUtils->IsPileUpEvent(InputEvent())) return;
180 
181  if(MCEvent()){
182  // Apply outlier cut
183  AliGenPythiaEventHeader *pyheader = GetPythiaHeader();
184  if(fFracPtHard > 0){
185  if(pyheader && IsOutlier(pyheader))
186  return;
187  }
188  fHistosMC->FillTH1("hNtrialsEvent", pyheader->Trials());
189  fHistosMC->FillTH1("hPtHard", pyheader->GetPtHard());
190  }
191 
192  Double_t pvecjet[3], pveclead[3], pvecsublead[3];
193  if(MCEvent()){
194  // Prepare MC jet finding
195  std::vector<fastjet::PseudoJet> mcpseudo;
196  for(int itrk = 0; itrk < fMCEvent->GetNumberOfTracks(); itrk++){
197  AliVParticle *mctrack = fMCEvent->GetTrack(itrk);
198  if(!mctrack->Charge()) continue;
199  if(TMath::Abs(mctrack->Pt()) < 0.15) continue;
200  if(TMath::Abs(mctrack->Eta()) > 0.8) continue;
201  if(!IsPhysicalPrimary(mctrack, MCEvent())) continue;
202 
203  fastjet::PseudoJet inputjet(mctrack->Px(), mctrack->Py(), mctrack->Pz(), mctrack->E());
204  inputjet.set_user_index(itrk);
205  mcpseudo.push_back(inputjet);
206  }
207 
208  // Run jet finding
209  fastjet::ClusterSequence jetfinder(mcpseudo, fastjet::JetDefinition(fastjet::antikt_algorithm, 0.4));
210  std::vector<fastjet::PseudoJet> genjets = sorted_by_pt(jetfinder.inclusive_jets());
211 
212  for(std::vector<fastjet::PseudoJet>::const_iterator jetit = genjets.begin(); jetit != genjets.end(); ++jetit){
213  if(TMath::Abs(jetit->eta()) > 0.5) continue; // expect jet to be within |eta| < +- 0.5
214  std::vector<fastjet::PseudoJet> constituentssorted = sorted_by_pt(jetit->constituents());
215  memset(pveclead, 0, sizeof(Double_t) * 3);
216  memset(pvecsublead, 0, sizeof(Double_t) * 3);
217  pvecjet[0] = jetit->px();
218  pvecjet[1] = jetit->py();
219  pvecjet[2] = jetit->pz();
220  int icounter = 0;
221  for(std::vector<fastjet::PseudoJet>::const_iterator cit = constituentssorted.begin(); cit != constituentssorted.end(); ++cit){
222  AliVParticle *basepart = fMCEvent->GetTrack(cit->user_index());
223  if(icounter == 0) {
224  basepart->PxPyPz(pveclead);
225  } else if(icounter == 1) {
226  basepart->PxPyPz(pvecsublead);
227  } else {
228  break;
229  }
230  icounter++;
231  }
232 
233  if(icounter >= 2){
234  // At least 2 particles, accept jet
235  fJetStructure.Reset();
236  fJetStructure.fIsData = 0;
237  memcpy(fJetStructure.fPvecJet, pvecjet, sizeof(Double_t) * 3);
238  memcpy(fJetStructure.fPvecLead, pveclead, sizeof(Double_t) * 3);
239  memcpy(fJetStructure.fPvecSubLead, pvecsublead, sizeof(Double_t) * 3);
240  fJetTree->Fill();
241  }
242  }
243  }
244 
245  // Select tracks
246  TVector3 pvectrack;
247  std::vector<fastjet::PseudoJet> datapseudo;
248  AliESDtrack *esdtrack(NULL);
249  AliAODTrack *aodtrack(NULL);
250  double mpion = TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass();
251  double mass(0);
252  for(int ipart = 0; ipart < fInputEvent->GetNumberOfTracks(); ipart++){
253  AliVParticle *recpart = fInputEvent->GetTrack(ipart);
254  if(MCEvent()){
255  AliVParticle *assocMC = fMCEvent->GetTrack(TMath::Abs(recpart->GetLabel()));
256  if(!assocMC) continue;
257  if(!IsPhysicalPrimary(assocMC, fMCEvent)) continue;
258  mass = TDatabasePDG::Instance()->GetParticle(assocMC->PdgCode())->Mass();
259  } else {
260  mass = mpion;
261  }
262 
263  // Select track
264  if(TMath::Abs(recpart->Eta()) > 0.8) continue;
265  if((esdtrack = dynamic_cast<AliESDtrack *>(recpart))){
266  AliESDtrack copytrack(*esdtrack);
267  if(!TrackSelectionESDHybrid(&copytrack)) continue;
268  pvectrack.SetXYZ(copytrack.Px(), copytrack.Py(), copytrack.Pz());
269  } else if((aodtrack = dynamic_cast<AliAODTrack *>(recpart))){
270  if(!TrackSelectionAODHybrid(aodtrack)) continue;
271  pvectrack.SetXYZ(recpart->Px(), recpart->Py(), recpart->Pz());
272  } else continue;
273 
274  // Create pseudojet under the assumption the particle is a pion
275  fastjet::PseudoJet inputparticle(pvectrack.Px(), pvectrack.Py(), pvectrack.Pz(), TMath::Sqrt(pvectrack.Perp()*pvectrack.Perp() + mass * mass));
276  inputparticle.set_user_index(ipart);
277  datapseudo.push_back(inputparticle);
278  }
279 
280  // Run jetfinder
281  fastjet::ClusterSequence jetfinder(datapseudo, fastjet::JetDefinition(fastjet::antikt_algorithm, 0.4));
282  std::vector<fastjet::PseudoJet> recjets = sorted_by_pt(jetfinder.inclusive_jets());
283 
284  // Loop over jets, select jets with at least 2 tracks, store momentum vector of track, leading jet and subleading jet
285  for(std::vector<fastjet::PseudoJet>::const_iterator jetit = recjets.begin(); jetit != recjets.end(); ++jetit){
286  if(TMath::Abs(jetit->eta()) > 0.5) continue; // expect jet to be within |eta| < +- 0.5
287  std::vector<fastjet::PseudoJet> constituentssorted = sorted_by_pt(jetit->constituents());
288  memset(pveclead, 0, sizeof(Double_t) * 3);
289  memset(pvecsublead, 0, sizeof(Double_t) * 3);
290  pvecjet[0] = jetit->px();
291  pvecjet[1] = jetit->py();
292  pvecjet[2] = jetit->pz();
293  int icounter = 0;
294  for(std::vector<fastjet::PseudoJet>::const_iterator cit = constituentssorted.begin(); cit != constituentssorted.end(); ++cit){
295  AliVParticle *basepart = fInputEvent->GetTrack(cit->user_index());
296  // Now apply stronger selection: default track cuts
297  bool isSelected = false;
298  if(basepart->IsA() == AliESDtrack::Class()){
299  isSelected = TrackSelectionESDDefault(static_cast<AliESDtrack *>(basepart));
300  } else {
301  isSelected = TrackSelectionAODDefault(static_cast<AliAODTrack *>(basepart));
302  }
303  if(!isSelected) continue;
304  if(icounter == 0) {
305  basepart->PxPyPz(pveclead);
306  } else if(icounter == 1) {
307  basepart->PxPyPz(pvecsublead);
308  } else {
309  break;
310  }
311  icounter++;
312  }
313 
314  if(icounter >= 2){
315  // At least 2 particles, accept jet
316  fJetStructure.Reset();
317  fJetStructure.fIsData = 1;
318  memcpy(fJetStructure.fPvecJet, pvecjet, sizeof(Double_t) * 3);
319  memcpy(fJetStructure.fPvecLead, pveclead, sizeof(Double_t) * 3);
320  memcpy(fJetStructure.fPvecSubLead, pvecsublead, sizeof(Double_t) * 3);
321  fJetTree->Fill();
322  }
323  }
324 
325  PostData(1, fHistosMC->GetListOfHistograms());
326  PostData(2, fJetTree);
327 }
328 
335 Bool_t AliAnalysisTaskTracksInJet::TrackSelectionESDHybrid(AliESDtrack* track) const {
336  if (fHybridCutsCat1->AcceptTrack(track)) {
337  track->SetBit(BIT(22),0);
338  track->SetBit(BIT(23),0);
339  return true;
340  } else if (fHybridCutsCat2->AcceptTrack(track)) {
341  if (!track->GetConstrainedParam())
342  return false;
343  UInt_t status = track->GetStatus();
344  if (((status&AliESDtrack::kITSrefit)==0))
345  return false;
346  const AliExternalTrackParam* constrainParam = track->GetConstrainedParam();
347  track->Set(constrainParam->GetX(),
348  constrainParam->GetAlpha(),
349  constrainParam->GetParameter(),
350  constrainParam->GetCovariance());
351  if ((status&AliESDtrack::kITSrefit)==0) {
352  track->SetBit(BIT(22),0); //type 2
353  track->SetBit(BIT(23),1);
354  } else {
355  track->SetBit(BIT(22),1); //type 1
356  track->SetBit(BIT(23),0);
357  }
358  return true;
359  }
360 
361  return false;
362 }
363 
369 Bool_t AliAnalysisTaskTracksInJet::TrackSelectionESDDefault(AliESDtrack* track) const {
370  return fTrackCutsDefault->AcceptTrack(track);
371 }
372 
378 Bool_t AliAnalysisTaskTracksInJet::TrackSelectionAODHybrid(AliAODTrack* track) const {
379  // @TODO: Change to hybrid track cuts (256,512)
380  if(!(track->TestFilterBit(256) || track->TestFilterBit(512))) return false;
381  return true;
382 }
383 
389 Bool_t AliAnalysisTaskTracksInJet::TrackSelectionAODDefault(AliAODTrack* track) const {
390  if(!track->TestFilterBit(AliAODTrack::kTrkGlobal)) return false;
391  if(track->GetTPCNCrossedRows() < 120) return false;
392  return true;
393 }
394 
407 Bool_t AliAnalysisTaskTracksInJet::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard) const {
408 
409  TString file(currFile);
410  fXsec = 0;
411  fTrials = 1;
412 
413  if (file.Contains(".zip#")) {
414  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
415  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
416  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
417  file.Replace(pos+1,pos2-pos1,"");
418  } else {
419  // not an archive take the basename....
420  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
421  }
422  AliDebug(1,Form("File name: %s",file.Data()));
423 
424  // Get the pt hard bin
425  TString strPthard(file);
426 
427  strPthard.Remove(strPthard.Last('/'));
428  strPthard.Remove(strPthard.Last('/'));
429  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
430  strPthard.Remove(0,strPthard.Last('/')+1);
431  if (strPthard.IsDec())
432  pthard = strPthard.Atoi();
433  else
434  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
435 
436  // 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
437  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
438 
439  if (!fxsec) {
440  // next trial fetch the histgram file
441  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
442  if (!fxsec) {
443  // not a severe condition but inciate that we have no information
444  return kFALSE;
445  } else {
446  // find the tlist we want to be independtent of the name so use the Tkey
447  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
448  if (!key) {
449  fxsec->Close();
450  return kFALSE;
451  }
452  TList *list = dynamic_cast<TList*>(key->ReadObj());
453  if (!list) {
454  fxsec->Close();
455  return kFALSE;
456  }
457  fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
458  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
459  fxsec->Close();
460  }
461  } else { // no tree pyxsec.root
462  TTree *xtree = (TTree*)fxsec->Get("Xsection");
463  if (!xtree) {
464  fxsec->Close();
465  return kFALSE;
466  }
467  UInt_t ntrials = 0;
468  Double_t xsection = 0;
469  xtree->SetBranchAddress("xsection",&xsection);
470  xtree->SetBranchAddress("ntrials",&ntrials);
471  xtree->GetEntry(0);
472  fTrials = ntrials;
473  fXsec = xsection;
474  fxsec->Close();
475  }
476  return kTRUE;
477 }
478 
487 Bool_t AliAnalysisTaskTracksInJet::IsPhysicalPrimary(const AliVParticle* const part, AliMCEvent* const mcevent) const {
488  Bool_t physprim = false;
489  const AliAODMCParticle *aodmc = dynamic_cast<const AliAODMCParticle *>(part);
490  if(aodmc){
491  physprim = aodmc->IsPhysicalPrimary();
492  } else {
493  physprim = mcevent->IsPhysicalPrimary(part->GetLabel());
494  }
495  return physprim;
496 }
497 
502 AliGenPythiaEventHeader *AliAnalysisTaskTracksInJet::GetPythiaHeader() const {
503  AliGenPythiaEventHeader *pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
504  if (!pythiaHeader) {
505  // Check if AOD
506  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
507  if (aodMCH) {
508  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
509  pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
510  if (pythiaHeader) break;
511  }
512  }
513  }
514  return pythiaHeader;
515 }
516 
522 Bool_t AliAnalysisTaskTracksInJet::IsOutlier(AliGenPythiaEventHeader * const header) const {
523  Bool_t hasOutlier = kFALSE;
524  Float_t pbuf[4];
525  TLorentzVector jetvec;
526  for(int ijet = 0; ijet < header->NTriggerJets(); ijet++){
527  memset(pbuf, 0, sizeof(Float_t) * 4);
528  header->TriggerJet(ijet, pbuf);
529  jetvec.SetPxPyPzE(pbuf[0], pbuf[1], pbuf[2], pbuf[3]);
530  if(TMath::Abs(jetvec.Pt()) >= fFracPtHard * header->GetPtHard()){
531  hasOutlier = true;
532  break;
533  }
534  }
535  return hasOutlier;
536 }
537 
538 } /* namespace EMCalTriggerPtAnalysis */
double Double_t
Definition: External.C:58
Double_t mass
TSystem * gSystem
TList * list
ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskTracksInJet) namespace EMCalTriggerPtAnalysis
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Stores p-vector of jet, leading track and subleading track.
TFile * file
Int_t nevents[nsamples]
Container class for histograms.
Definition: THistManager.h:99
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65