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