1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
21 * Author: Andrew Castro (UTK) and Joel Mazer (UTK)
24 #include "AliAnalysisTaskEmcalJetFlavourTagExample.h"
26 // general ROOT includes
29 #include <TClonesArray.h>
33 #include <THnSparse.h>
35 #include <TLorentzVector.h>
36 #include <TParameter.h>
37 #include <TParticle.h>
40 #include <TObjArray.h>
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"
56 #include "AliJetContainer.h"
57 #include "AliParticleContainer.h"
58 #include "AliClusterContainer.h"
59 #include "AliEmcalParticle.h"
60 #include "AliESDCaloCluster.h"
61 #include <AliESDtrackCuts.h>
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"
74 #include "AliPIDResponse.h"
75 #include "AliTPCPIDResponse.h"
76 #include "AliESDpid.h"
78 // magnetic field includes
79 #include "TGeoGlobalMagField.h"
85 ClassImp(AliAnalysisTaskEmcalJetFlavourTagExample)
87 //________________________________________________________________________
88 AliAnalysisTaskEmcalJetFlavourTagExample::AliAnalysisTaskEmcalJetFlavourTagExample() :
89 AliAnalysisTaskEmcalJet("heavyF",kFALSE),
92 fPhimin(-10), fPhimax(10),
93 fEtamin(-0.9), fEtamax(0.9),
99 fPIDResponse(0x0), fTPCResponse(),
100 fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fTracksJetCont(0), fCaloClustersJetCont(0),
103 fHistCorJetPt(0), fHistJetPt(0),
107 fHistnJetTrackvnJetClusters(0)
109 // Default constructor.
112 SetMakeGeneralHistograms(kTRUE);
116 //________________________________________________________________________
117 AliAnalysisTaskEmcalJetFlavourTagExample::AliAnalysisTaskEmcalJetFlavourTagExample(const char *name) :
118 AliAnalysisTaskEmcalJet(name,kTRUE),
121 fPhimin(-10), fPhimax(10),
122 fEtamin(-0.9), fEtamax(0.9),
128 fPIDResponse(0x0), fTPCResponse(),
129 fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fTracksJetCont(0), fCaloClustersJetCont(0),
132 fHistCorJetPt(0), fHistJetPt(0),
136 fHistnJetTrackvnJetClusters(0)
139 SetMakeGeneralHistograms(kTRUE);
141 DefineInput(0,TChain::Class());
142 DefineOutput(1, TList::Class());
145 //_______________________________________________________________________
146 AliAnalysisTaskEmcalJetFlavourTagExample::~AliAnalysisTaskEmcalJetFlavourTagExample()
156 //________________________________________________________________________
157 void AliAnalysisTaskEmcalJetFlavourTagExample::UserCreateOutputObjects()
161 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
163 //fJetsCont = GetJetContainer(0);
164 if(fJetsCont) { //get particles and clusters connected to jets
165 fTracksJetCont = fJetsCont->GetParticleContainer();
166 fCaloClustersJetCont = fJetsCont->GetClusterContainer();
168 else { //no jets, just analysis tracks and clusters
169 fTracksCont = GetParticleContainer(0);
170 fCaloClustersCont = GetClusterContainer(0);
172 fTracksCont->SetClassName("AliVTrack");
173 fCaloClustersCont->SetClassName("AliVCluster");
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);
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);
192 fOutput->Add(fHistCorJetPt);
193 fOutput->Add(fHistnSigElecPt);
194 fOutput->Add(fHistdEdXvsPt);
195 fOutput->Add(fHistnJetTrackvnJetClusters);
196 fOutput->Add(fHistHighJetPt);
199 // ****************************** PID *****************************************************
200 // set up PID handler
201 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
202 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
204 AliFatal("Input handler needed");
207 // PID response object
208 //fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
209 // inputHandler->CreatePIDResponse(fIsMC); // needed to create object, why though?
210 fPIDResponse = inputHandler->GetPIDResponse();
212 AliError("PIDResponse object was not created");
214 // ***************************************************************************************
216 PostData(1, fOutput);
220 //________________________________________________________
221 void AliAnalysisTaskEmcalJetFlavourTagExample::ExecOnce()
223 // Initialize the analysis
224 AliAnalysisTaskEmcalJet::ExecOnce();
226 if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
227 if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
228 if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
233 //________________________________________________________________________
234 Bool_t AliAnalysisTaskEmcalJetFlavourTagExample::Run()
236 // check to see if we have any tracks
237 if (!fTracks) return kTRUE;
238 if (!fJets) return kTRUE;
240 // what kind of event do we have: AOD or ESD?
242 if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
243 else useAOD = kFALSE;
245 // if we have ESD event, set up ESD object
247 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
249 AliError(Form("ERROR: fESD not available\n"));
254 // if we have AOD event, set up AOD object
256 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
258 AliError(Form("ERROR: fAOD not available\n"));
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);
273 // get vertex information
274 Double_t fvertex[3]={0,0,0};
275 InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
276 //Double_t zVtx=fvertex[2];
278 // create pointer to list of input event
279 TList *list = InputEvent()->GetList();
281 AliError(Form("ERROR: list not attached\n"));
285 // background density
286 fRhoVal = fRho->GetVal();
288 // initialize TClonesArray pointers to jets and tracks
289 TClonesArray *jets = 0;
290 //TClonesArray *tracks = 0;
291 //TClonesArray *clusters = 0;
294 jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
296 AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
298 } // verify existence of jets
300 // get number of jets and tracks
301 const Int_t Njets = jets->GetEntries();
302 if(Njets<1) return kTRUE;
305 AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
307 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
310 if (fCaloClustersCont) {
311 AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
313 TLorentzVector nPart;
314 cluster->GetMomentum(nPart, fVertex);
315 cluster = fCaloClustersCont->GetNextAcceptCluster();
319 // Start Jet Analysis
320 // initialize jet parameters
322 Double_t highestjetpt=0.0;
324 // loop over jets in an event - to find highest jet pT and apply some cuts && JetQA Sparse
325 for (Int_t ijet = 0; ijet < Njets; ijet++){
327 AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
331 if(!AcceptMyJet(jet)) continue;
333 if(highestjetpt<jet->Pt()){
335 highestjetpt=jet->Pt();
337 } // end of looping over jets
339 fHistHighJetPt->Fill(ijethi);
341 // **********************************************************************
343 // **********************************************************************
345 // loop over jets in the event and make appropriate cuts
346 for (Int_t iJets = 0; iJets < Njets; ++iJets) {
347 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
348 if (!jet) // see if we have a jet
352 if(!AcceptMyJet(jet)) continue;
354 //AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
355 //jet->AddFlavourTag(tag);
358 Int_t JetClusters = jet->GetNumberOfClusters();
359 Int_t JetTracks = jet -> GetNumberOfTracks();
360 fHistnJetTrackvnJetClusters->Fill(JetClusters,JetTracks);
361 // Initializations and Calculations
362 Double_t jetPt = -500; // initialize corr jet pt LOCAL
363 Double_t jetarea = -500; // initialize jet area
364 jetPt = jet->Pt() - jetarea*fRhoVal; // semi-corrected pT of jet from GLOBAL rho value
365 fHistCorJetPt->Fill(jetPt);
367 Bool_t bkgrnd1 = kFALSE;
369 if(jet->Pt() > fJetHIpt) {
370 if(!fTracksCont || !fCaloClustersCont) continue;
372 //******************************Cluster Matched To Closest Track
373 //**************************************************************
374 Int_t NumbTrackContainer = -999;
375 NumbTrackContainer = fTracksCont->GetNParticles();
376 for(int iTracks = 0; iTracks <= NumbTrackContainer; iTracks++){
377 AliVTrack *AcceptedTrack =static_cast<AliVTrack*>(fTracksCont->GetParticle(iTracks));
379 AliError(Form("Couldn't get AliVTrack Container %d\n", iTracks));
382 if(!IsJetTrack(jet,iTracks,kFALSE))continue;
383 //Get matched cluster
384 Int_t emc1 = AcceptedTrack->GetEMCALcluster();
386 Double_t acceptTrackP = AcceptedTrack->P();
387 Double_t acceptTrackPt = AcceptedTrack->Pt();
388 Double_t nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(AcceptedTrack,AliPID::kElectron);
389 fHistnSigElecPt->Fill(acceptTrackPt,nSigmaElectron_TPC);
391 AliESDtrack *ESDacceptedTrack = static_cast<AliESDtrack*>(AcceptedTrack);
392 if(!ESDacceptedTrack){
393 AliError(Form("Couldn't get AliESDTrack %d\n", iTracks));
396 Double_t dEdx = AcceptedTrack->GetTPCsignal();
397 fHistdEdXvsPt->Fill(acceptTrackPt,dEdx);
398 if(fCaloClustersCont && emc1>=0) {
399 AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
401 AliError(Form("Couldn't get matched AliVCluster %d\n", emc1));
405 Double_t mClusterE = clusMatch->E();
406 Double_t EovP_mc = -999;
407 EovP_mc = mClusterE/acceptTrackP;
409 if(EovP_mc < 0.2) bkgrnd1 = kTRUE;
412 } //loop over tracks for matching to closest cluster
417 } // highest pt jet cut
420 if(bkgrnd1 == kTRUE) {
421 AliEmcalJet::EFlavourTag tag=AliEmcalJet::kBckgrd1;
422 jet->AddFlavourTag(tag);
425 } // LOOP over JETS in event
431 //________________________________________________________________________
432 void AliAnalysisTaskEmcalJetFlavourTagExample::Terminate(Option_t *)
434 cout<<"###########################"<<endl;
435 cout<<"#### Task Finished ####"<<endl;
436 cout<<"###########################"<<endl;
437 cout<<"###########################"<<endl;
438 } // end of terminate
440 //________________________________________________________________________
441 Int_t AliAnalysisTaskEmcalJetFlavourTagExample::AcceptMyJet(AliEmcalJet *jet) {
442 //applies all jet cuts except pt
443 if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
444 if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
445 if (jet->Area()<fAreacut) return 0;
446 // prevents 0 area jets from sneaking by when area cut == 0
447 if (jet->Area()==0) return 0;
448 //exclude jets with extremely high pt tracks which are likely misreconstructed
449 if(jet->MaxTrackPt()>100) return 0;
451 //passed all above cuts