AliPhysics  ff1d528 (ff1d528)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalJetFlavourTagExample.cxx
Go to the documentation of this file.
1 /**************************************************************************
2 * Copyright(c) 1998-2007, 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 
16 /*
17 * This sample task demostrates the basics of tagging a jet. The PID response
18 * is retrived for both the TPC. A Hadronic tag is implemented for
19 * clusters matched to tracks in a jet with a E/P < 0.2
20 *
21 * Author: Andrew Castro (UTK) and Joel Mazer (UTK)
22 */
23 
25 
26 // general ROOT includes
27 #include <TCanvas.h>
28 #include <TChain.h>
29 #include <TClonesArray.h>
30 #include <TH1F.h>
31 #include <TH2F.h>
32 #include <TH3F.h>
33 #include <THnSparse.h>
34 #include <TList.h>
35 #include <TLorentzVector.h>
36 #include <TParameter.h>
37 #include <TParticle.h>
38 #include <TTree.h>
39 #include <TVector3.h>
40 #include <TObjArray.h>
41 
42 // AliROOT includes
43 #include "AliAODEvent.h"
44 #include "AliESDEvent.h"
45 #include "AliAnalysisManager.h"
46 #include "AliAnalysisTask.h"
47 #include "AliCentrality.h"
48 #include "AliEmcalJet.h"
49 #include "AliAODJet.h"
50 #include "AliVCluster.h"
51 #include "AliVTrack.h"
52 #include <AliVEvent.h>
53 #include <AliVParticle.h>
54 #include "AliRhoParameter.h"
55 #include "AliLog.h"
56 #include "AliJetContainer.h"
57 #include "AliParticleContainer.h"
58 #include "AliClusterContainer.h"
59 #include "AliEmcalParticle.h"
60 #include "AliESDCaloCluster.h"
61 #include <AliESDtrackCuts.h>
62 #include "AliPID.h"
63 
64 // event handler (and pico's) includes
65 #include <AliInputEventHandler.h>
66 #include <AliVEventHandler.h>
67 #include "AliESDInputHandler.h"
68 #include "AliPicoTrack.h"
69 #include "AliEventPoolManager.h"
70 #include "AliAODTrack.h"
71 #include "AliESDtrack.h"
72 
73 // PID includes
74 #include "AliPIDResponse.h"
75 #include "AliTPCPIDResponse.h"
76 #include "AliESDpid.h"
77 
78 // magnetic field includes
79 #include "TGeoGlobalMagField.h"
80 #include "AliMagF.h"
81 
82 using std::cout;
83 using std::endl;
84 
86 
87 //________________________________________________________________________
89  AliAnalysisTaskEmcalJet("heavyF",kFALSE),
90  event(0),
91  fCuts(0),
92  fPhimin(-10), fPhimax(10),
93  fEtamin(-0.9), fEtamax(0.9),
94  fAreacut(0.0),
95  fJetHIpt(50.0),
96  fTrackPtCut(2.0),
97  fTrackEta(0.9),
98  fesdTrackCuts(0),
99  fPIDResponse(0x0), fTPCResponse(),
100  fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fTracksJetCont(0), fCaloClustersJetCont(0),
101  fESD(0), fAOD(0),
102  fHistJetPhi(0),
103  fHistCorJetPt(0), fHistJetPt(0),
104  fHistHighJetPt(0),
105  fHistnSigElecPt(0),
106  fHistdEdXvsPt(0),
107  fHistnJetTrackvnJetClusters(0)
108 {
109  // Default constructor.
110 
111 
113 
114 }
115 
116 //________________________________________________________________________
118  AliAnalysisTaskEmcalJet(name,kTRUE),
119  event(0),
120  fCuts(0),
121  fPhimin(-10), fPhimax(10),
122  fEtamin(-0.9), fEtamax(0.9),
123  fAreacut(0.0),
124  fJetHIpt(50.0),
125  fTrackPtCut(2.0),
126  fTrackEta(0.9),
127  fesdTrackCuts(0),
128  fPIDResponse(0x0), fTPCResponse(),
129  fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fTracksJetCont(0), fCaloClustersJetCont(0),
130  fESD(0), fAOD(0),
131  fHistJetPhi(0),
132  fHistCorJetPt(0), fHistJetPt(0),
133  fHistHighJetPt(0),
134  fHistnSigElecPt(0),
135  fHistdEdXvsPt(0),
136  fHistnJetTrackvnJetClusters(0)
137 {
138 
140 
141  DefineInput(0,TChain::Class());
142  DefineOutput(1, TList::Class());
143 }
144 
145 //_______________________________________________________________________
147 {
148  // destructor
149  //
150  if (fOutput) {
151  delete fOutput;
152  fOutput = 0;
153  }
154 }
155 
156 //________________________________________________________________________
158 {
159  if (! fCreateHisto)
160  return;
162 
163  //fJetsCont = GetJetContainer(0);
164  if(fJetsCont) { //get particles and clusters connected to jets
167  }
168  else { //no jets, just analysis tracks and clusters
171 }
172 fTracksCont->SetClassName("AliVTrack");
173 fCaloClustersCont->SetClassName("AliVCluster");
174 
175  fHistJetPhi = new TH1F("NjetvsPhi", "NjetvsPhi", 288,-2*TMath::Pi(),2*TMath::Pi());
176  fHistJetPt = new TH1F("NjetvsJetPt", "NjetvsJetPt", 300, 0, 300);
177  fOutput->Add(fHistJetPhi);
178  fOutput->Add(fHistJetPt);
179 
180 
181  TString histname;
182 
183  fHistCorJetPt = new TH1F("CorrJetPt", "CorrJetPt", 300, -100, 200);
184  fHistnSigElecPt = new TH2F("nsigvsPt(TPC)","nSigmaElectronTPC_v_Pt",60,0,60,40,-10,10);
185  fHistdEdXvsPt = new TH2F("dEdXvsTrackPt","dEdXvsTrackPt",60,0,60,80,0,80);
186  fHistnJetTrackvnJetClusters = new TH2F("NumbJetTracksvJetClusters","NumbJetTracksvJetClusters",21,0,20,21,0,20);
187  fHistHighJetPt = new TH1F("HighestPtJetPerEvent","HighJetPt",80,0,80);
188 
189  TString name;
190  TString title;
191 
192  fOutput->Add(fHistCorJetPt);
193  fOutput->Add(fHistnSigElecPt);
194  fOutput->Add(fHistdEdXvsPt);
196  fOutput->Add(fHistHighJetPt);
197 
198 
199  // ****************************** PID *****************************************************
200  // set up PID handler
201  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
202  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
203  if(!inputHandler) {
204  AliFatal("Input handler needed");
205  }
206 
207  // PID response object
208  //fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
209  // inputHandler->CreatePIDResponse(fIsMC); // needed to create object, why though?
210  fPIDResponse = inputHandler->GetPIDResponse();
211  if (!fPIDResponse) {
212  AliError("PIDResponse object was not created");
213  }
214  // ***************************************************************************************
215 
216  PostData(1, fOutput);
217 
218 }
219 
220 //________________________________________________________
222 {
223  // Initialize the analysis
225 
226  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
227  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
228  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
229 
230 
231 } // end of ExecOnce
232 
233 //________________________________________________________________________
235 {
236  // check to see if we have any tracks
237  if (!fTracks) return kTRUE;
238  if (!fJets) return kTRUE;
239 
240  // what kind of event do we have: AOD or ESD?
241  Bool_t useAOD;
242  if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
243  else useAOD = kFALSE;
244 
245  // if we have ESD event, set up ESD object
246  if(!useAOD){
247  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
248  if (!fESD) {
249  AliError(Form("ERROR: fESD not available\n"));
250  return kTRUE;
251  }
252  }
253 
254  // if we have AOD event, set up AOD object
255  if(useAOD){
256  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
257  if(!fAOD) {
258  AliError(Form("ERROR: fAOD not available\n"));
259  return kTRUE;
260  }
261  }
262 
263  // get magnetic field info for DCA
264  Double_t MagF = fESD->GetMagneticField();
265  Double_t MagSign = 1.0;
266  if(MagF<0)MagSign = -1.0;
267  // set magnetic field
268  if (!TGeoGlobalMagField::Instance()->GetField()) {
269  AliMagF* field = new AliMagF("Maps","Maps", MagSign, MagSign, AliMagF::k5kG);
270  TGeoGlobalMagField::Instance()->SetField(field);
271  }
272 
273  // get vertex information
274  Double_t fvertex[3]={0,0,0};
275  InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
276  //Double_t zVtx=fvertex[2];
277 
278  // create pointer to list of input event
279  TList *list = InputEvent()->GetList();
280  if(!list) {
281  AliError(Form("ERROR: list not attached\n"));
282  return kTRUE;
283  }
284 
285  // background density
286  fRhoVal = fRho->GetVal();
287 
288  // initialize TClonesArray pointers to jets and tracks
289  TClonesArray *jets = 0;
290  //TClonesArray *tracks = 0;
291  //TClonesArray *clusters = 0;
292 
293  // get Jets object
294  jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
295  if(!jets){
296  AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
297  return kTRUE;
298  } // verify existence of jets
299 
300  // get number of jets and tracks
301  const Int_t Njets = jets->GetEntries();
302  if(Njets<1) return kTRUE;
303 
304  if (fTracksCont) {
305  fTracksCont->ResetCurrentID();
306  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
307  while(track) {
308  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
309  }
310  }
311  if (fCaloClustersCont) {
312  fCaloClustersCont->ResetCurrentID();
313  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
314  while(cluster) {
315  TLorentzVector nPart;
316  cluster->GetMomentum(nPart, fVertex);
318  }
319  }
320 
321  // Start Jet Analysis
322  // initialize jet parameters
323  Int_t ijethi=-1;
324  Double_t highestjetpt=0.0;
325 
326  // loop over jets in an event - to find highest jet pT and apply some cuts && JetQA Sparse
327  for (Int_t ijet = 0; ijet < Njets; ijet++){
328  // get our jets
329  AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
330  if (!jet) continue;
331 
332  // apply jet cuts
333  if(!AcceptMyJet(jet)) continue;
334 
335  if(highestjetpt<jet->Pt()){
336  ijethi=ijet;
337  highestjetpt=jet->Pt();
338  }
339  } // end of looping over jets
340 
341  fHistHighJetPt->Fill(ijethi);
342 
343  // **********************************************************************
344  // JET LOOP
345  // **********************************************************************
346 
347  // loop over jets in the event and make appropriate cuts
348  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
349  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
350  if (!jet) // see if we have a jet
351  continue;
352 
353  // apply jet cuts
354  if(!AcceptMyJet(jet)) continue;
355 
356  //AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
357  //jet->AddFlavourTag(tag);
358 
359 
360  Int_t JetClusters = jet->GetNumberOfClusters();
361  Int_t JetTracks = jet -> GetNumberOfTracks();
362  fHistnJetTrackvnJetClusters->Fill(JetClusters,JetTracks);
363  // Initializations and Calculations
364  Double_t jetPt = -500; // initialize corr jet pt LOCAL
365  Double_t jetarea = -500; // initialize jet area
366  jetPt = jet->Pt() - jetarea*fRhoVal; // semi-corrected pT of jet from GLOBAL rho value
367  fHistCorJetPt->Fill(jetPt);
368 
369  Bool_t bkgrnd1 = kFALSE;
370 
371  if(jet->Pt() > fJetHIpt) {
372  if(!fTracksCont || !fCaloClustersCont) continue;
373 
374  //******************************Cluster Matched To Closest Track
375  //**************************************************************
376  Int_t NumbTrackContainer = -999;
377  NumbTrackContainer = fTracksCont->GetNParticles();
378  for(int iTracks = 0; iTracks <= NumbTrackContainer; iTracks++){
379  AliVTrack *AcceptedTrack =static_cast<AliVTrack*>(fTracksCont->GetParticle(iTracks));
380  if(!AcceptedTrack){
381  AliError(Form("Couldn't get AliVTrack Container %d\n", iTracks));
382  continue;
383  }
384  if(!IsJetTrack(jet,iTracks,kFALSE))continue;
385  //Get matched cluster
386  Int_t emc1 = AcceptedTrack->GetEMCALcluster();
387 
388  Double_t acceptTrackP = AcceptedTrack->P();
389  Double_t acceptTrackPt = AcceptedTrack->Pt();
390  Double_t nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(AcceptedTrack,AliPID::kElectron);
391  fHistnSigElecPt->Fill(acceptTrackPt,nSigmaElectron_TPC);
392 
393  AliESDtrack *ESDacceptedTrack = static_cast<AliESDtrack*>(AcceptedTrack);
394  if(!ESDacceptedTrack){
395  AliError(Form("Couldn't get AliESDTrack %d\n", iTracks));
396  continue;
397  }
398  Double_t dEdx = AcceptedTrack->GetTPCsignal();
399  fHistdEdXvsPt->Fill(acceptTrackPt,dEdx);
400  if(fCaloClustersCont && emc1>=0) {
401  AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
402  if(!clusMatch){
403  AliError(Form("Couldn't get matched AliVCluster %d\n", emc1));
404  continue;
405  }
406 
407  Double_t mClusterE = clusMatch->E();
408  Double_t EovP_mc = -999;
409  EovP_mc = mClusterE/acceptTrackP;
410 
411  if(EovP_mc < 0.2) bkgrnd1 = kTRUE;
412  }
413 
414  } //loop over tracks for matching to closest cluster
415 
416 
417 
418 
419  } // highest pt jet cut
420 
421 
422  if(bkgrnd1 == kTRUE) {
424  jet->AddFlavourTag(tag);
425  }
426 
427  } // LOOP over JETS in event
428 
429 
430  return kTRUE;
431 
432 }
433 //________________________________________________________________________
435 {
436  cout<<"###########################"<<endl;
437  cout<<"#### Task Finished ####"<<endl;
438  cout<<"###########################"<<endl;
439  cout<<"###########################"<<endl;
440 } // end of terminate
441 
442 //________________________________________________________________________
444  //applies all jet cuts except pt
445  if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
446  if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
447  if (jet->Area()<fAreacut) return 0;
448  // prevents 0 area jets from sneaking by when area cut == 0
449  if (jet->Area()==0) return 0;
450  //exclude jets with extremely high pt tracks which are likely misreconstructed
451  if(jet->MaxTrackPt()>100) return 0;
452 
453  //passed all above cuts
454  return 1;
455 }
456 
457 
458 
459 
460 
void AddFlavourTag(Int_t tag)
Definition: AliEmcalJet.h:245
Double_t Area() const
Definition: AliEmcalJet.h:123
virtual AliVParticle * GetNextAcceptParticle()
double Double_t
Definition: External.C:58
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
Double_t Eta() const
Definition: AliEmcalJet.h:114
Double_t Phi() const
Definition: AliEmcalJet.h:110
Int_t GetNParticles() const
TList * list
AliClusterContainer * GetClusterContainer() const
TH1F * fHistJetPt
// (Njets) vs Corrected Jet Pt (local rho)
void ExecOnce()
Perform steps needed to initialize the analysis.
Bool_t IsJetTrack(AliEmcalJet *jet, Int_t itrack, Bool_t sorted=kFALSE) const
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
AliParticleContainer * GetParticleContainer() const
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:131
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
EFlavourTag
Bit definition for the flavor tagging.
Definition: AliEmcalJet.h:77
AliRhoParameter * fRho
! event rho
virtual AliVParticle * GetParticle(Int_t i=-1) const
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Double_t MaxTrackPt() const
Definition: AliEmcalJet.h:148
Generic background 1.
Definition: AliEmcalJet.h:82
ClassImp(AliAnalysisTaskEmcalJetFlavourTagExample) AliAnalysisTaskEmcalJetFlavourTagExample
AliVCluster * GetCluster(Int_t i) const
TClonesArray * fJets
! jets
Double_t Pt() const
Definition: AliEmcalJet.h:102
AliEmcalList * fOutput
!output list
TClonesArray * fTracks
!tracks
Double_t fVertex[3]
!event vertex
Bool_t fCreateHisto
whether or not create histograms
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
Double_t fRhoVal
! event rho value, same for local rho
AliVCluster * GetNextAcceptCluster()