Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / HMPID / AliAnalysisTaskJetsHMPID.cxx
CommitLineData
bfcb155c 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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// AliAnalysysTaskJetsHMPID - Class performing PID analysis in jets with the HMPID
18// A set of histograms is created.
19//=================================================================================
20//
21// By means of AnalysisTrainHMPID.C macro it is possible to use this class
22// to perform the analysis on local data, on data on alien using local machine
23// and on CAF.
24
25#include <TClonesArray.h>
26#include <TRandom2.h>
27#include "AliAnalysisTaskJetsHMPID.h"
28#include "AliAnalysisManager.h"
29#include "AliAODHandler.h"
30#include "AliAODInputHandler.h"
31#include "AliAODEvent.h"
32#include "AliAODJet.h"
33#include "AliAODJetEventBackground.h"
34
35ClassImp(AliAnalysisTaskJetsHMPID)
36
37//__________________________________________________________________________
38AliAnalysisTaskJetsHMPID::AliAnalysisTaskJetsHMPID() :
39 fJetBranch("jets"),
40 fBkgBranch(""),
41 fJetPtCut(10.),
42 fAOD(0x0),
43 fHistList(0x0),
44 fThetaChJet(0x0),
45 fThetaChBkg(0x0),
46 fThetaChRndCone(0x0),
47 fEvSelDPhi(0x0),
48 fJetsPt(0x0),
49 fRndConePt(0x0),
50 fAwayJetPt(0x0),
51 fJetsEtaPhi(0x0),
52 fTrksEtaPhiJet(0x0),
53 fTrksEtaPhiBkg(0x0),
54 fTree(0x0),
55 fTrackPt(0),
56 fJetPt(0),
57 fPionBkg(1.1),
58 fKaonBkg(1.1),
59 fProtBkg(1.1),
60 fPionJet(1.1),
61 fKaonJet(1.1),
62 fProtJet(1.1)
63{
64 //
65 //Default ctor
66 //
67}
68
69//__________________________________________________________________________
70AliAnalysisTaskJetsHMPID::AliAnalysisTaskJetsHMPID(const Char_t* name) :
71 AliAnalysisTaskSE(name),
72 fJetBranch("jets"),
73 fBkgBranch(""),
74 fJetPtCut(10.),
75 fAOD(0x0),
76 fHistList(0x0),
77 fThetaChJet(0x0),
78 fThetaChBkg(0x0),
79 fThetaChRndCone(0x0),
80 fEvSelDPhi(0x0),
81 fJetsPt(0x0),
82 fRndConePt(0x0),
83 fAwayJetPt(0x0),
84 fJetsEtaPhi(0x0),
85 fTrksEtaPhiJet(0x0),
86 fTrksEtaPhiBkg(0x0),
87 fTree(0x0),
88 fTrackPt(0),
89 fJetPt(0),
90 fPionBkg(1.1),
91 fKaonBkg(1.1),
92 fProtBkg(1.1),
93 fPionJet(1.1),
94 fKaonJet(1.1),
95 fProtJet(1.1)
96{
97 //
98 // Constructor. Initialization of Inputs and Outputs
99 //
100
101 DefineOutput(1,TList::Class());
102 DefineOutput(2,TTree::Class());
103}
104
105//___________________________________________________________________________
106AliAnalysisTaskJetsHMPID& AliAnalysisTaskJetsHMPID::operator=(const AliAnalysisTaskJetsHMPID& c)
107{
108 //
109 // Assignment operator
110 //
111 if (this!=&c) {
112 AliAnalysisTaskSE::operator=(c);
113 fJetBranch = c.fJetBranch;
114 fBkgBranch = c.fBkgBranch;
115 fJetPtCut = c.fJetPtCut;
116 fAOD = c.fAOD;
117 fHistList = c.fHistList;
118 fThetaChJet = c.fThetaChJet;
119 fThetaChBkg = c.fThetaChBkg;
120 fThetaChRndCone = c.fThetaChRndCone;
121 fEvSelDPhi = c.fEvSelDPhi;
122 fJetsPt = c.fJetsPt;
123 fRndConePt = c.fRndConePt;
124 fAwayJetPt = c.fAwayJetPt;
125 fJetsEtaPhi = c.fJetsEtaPhi;
126 fTrksEtaPhiJet = c.fTrksEtaPhiJet;
127 fTrksEtaPhiBkg = c.fTrksEtaPhiBkg;
128 fTree = c.fTree;
129 fTrackPt = c.fTrackPt;
130 fJetPt = c.fJetPt;
131 fPionBkg = c.fPionBkg;
132 fKaonBkg = c.fKaonBkg;
133 fProtBkg = c.fProtBkg;
134 fPionJet = c.fPionJet;
135 fKaonJet = c.fKaonJet;
136 fProtJet = c.fProtJet;
137 }
138 return *this;
139}
140
141//___________________________________________________________________________
142AliAnalysisTaskJetsHMPID::AliAnalysisTaskJetsHMPID(const AliAnalysisTaskJetsHMPID& c) :
143 AliAnalysisTaskSE(c),
144 fJetBranch(c.fJetBranch),
145 fBkgBranch(c.fBkgBranch),
146 fJetPtCut(c.fJetPtCut),
147 fAOD(c.fAOD),
148 fHistList(c.fHistList),
149 fThetaChJet(c.fThetaChJet),
150 fThetaChBkg(c.fThetaChBkg),
151 fThetaChRndCone(c.fThetaChRndCone),
152 fEvSelDPhi(c.fEvSelDPhi),
153 fJetsPt(c.fJetsPt),
154 fRndConePt(c.fRndConePt),
155 fAwayJetPt(c.fAwayJetPt),
156 fJetsEtaPhi(c.fJetsEtaPhi),
157 fTrksEtaPhiJet(c.fTrksEtaPhiJet),
158 fTrksEtaPhiBkg(c.fTrksEtaPhiBkg),
159 fTree(c.fTree),
160 fTrackPt(c.fTrackPt),
161 fJetPt(c.fJetPt),
162 fPionBkg(c.fPionBkg),
163 fKaonBkg(c.fKaonBkg),
164 fProtBkg(c.fProtBkg),
165 fPionJet(c.fPionJet),
166 fKaonJet(c.fKaonJet),
167 fProtJet(c.fProtJet)
168{
169 //
170 // Copy Constructor
171 //
172}
173
174//___________________________________________________________________________
175AliAnalysisTaskJetsHMPID::~AliAnalysisTaskJetsHMPID()
176{
177 //
178 //destructor
179 //
180 Info("~AliAnalysisTaskJetsHMPID","Calling Destructor");
181 if (fHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHistList;
182}
183
184//__________________________________________________________________________
185void AliAnalysisTaskJetsHMPID::ConnectInputData(Option_t *)
186{
187 if (fDebug) AliInfo("ConnectInputData() ");
188
189 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
190
191 if(handler && handler->InheritsFrom("AliAODInputHandler"))
192 { // input AOD
193 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
194 if (fDebug) AliInfo("Tracks and Jets from AliAODInputHandler");
195 }
196 else
197 { //output AOD
198 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
199 if (handler && handler->InheritsFrom("AliAODHandler"))
200 {
201 fAOD = ((AliAODHandler*)handler)->GetAOD();
202 if (fDebug > 1) AliInfo("Tracks and Jets from AliAODHandler");
203 }
204 else
205 { // no AOD
206 AliWarning("No AOD event in the input or in the output");
207 }
208 }
209
210}
211
212//___________________________________________________________________________
213void AliAnalysisTaskJetsHMPID::UserCreateOutputObjects()
214{
215 if (!fHistList) fHistList = new TList();
216 fHistList->SetOwner();
217
218 fThetaChJet = new TH2F("ThetaChJet","Theta Cherenkov distribution in the jets",1000,0.,10.,400,0.,0.8);
219 fHistList->Add(fThetaChJet);
220
221 fThetaChBkg = new TH2F("ThetaChBkg","Theta Cherenkov distribution out of the jets",1000,0.,10.,400,0.,0.8);
222 fHistList->Add(fThetaChBkg);
223
224 fThetaChRndCone = new TH2F("ThetaChRndCone","Theta Cherenkov distribution in random cones",1000,0.,10.,400,0.,0.8);
225 fHistList->Add(fThetaChRndCone);
226
227 fEvSelDPhi = new TH1F("PhiJets","Jets phi in different selections",600,0.,6.);
228 fHistList->Add(fEvSelDPhi);
229
230 fJetsPt = new TH1F("JetsPt","Jets spectrum in the HMPID",200,0.,200.);
231 fHistList->Add(fJetsPt);
232
233 fRndConePt = new TH1F("RndConePt","Pt of random cones in the HMPID",200,0.,200.);
234 fHistList->Add(fRndConePt);
235
236 fAwayJetPt = new TH1F("AwayJetPt","Pt of jets opposite to the HMPID",200,0.,200.);
237 fHistList->Add(fAwayJetPt);
238
239 fJetsEtaPhi = new TH2F("JetsEtaPhi","Eta-Phi of jets in the HMPID",180,-0.9,0.9,200,-0.25*TMath::Pi(),0.75*TMath::Pi());
240 fHistList->Add(fJetsEtaPhi);
241
242 fTrksEtaPhiJet = new TH2F("TrksEtaPhiJet","Eta-Phi of jet tracks in the HMPID",180,-0.9,0.9,200,-0.25*TMath::Pi(),0.75*TMath::Pi());
243 fHistList->Add(fTrksEtaPhiJet);
244
245 fTrksEtaPhiBkg = new TH2F("TrksEtaPhiBkg","Eta-Phi of bkg tracks in the HMPID",180,-0.9,0.9,200,-0.25*TMath::Pi(),0.75*TMath::Pi());
246 fHistList->Add(fTrksEtaPhiBkg);
247
248 fTree = new TTree("Tree","Tree with data");
249 fTree->Branch("Pt",&fTrackPt);
250 fTree->Branch("PtJet",&fJetPt);
251 fTree->Branch("PionBkg",&fPionBkg);
252 fTree->Branch("KaonBkg",&fKaonBkg);
253 fTree->Branch("ProtBkg",&fProtBkg);
254 fTree->Branch("PionJet",&fPionJet);
255 fTree->Branch("KaonJet",&fKaonJet);
256 fTree->Branch("ProtJet",&fProtJet);
257
258 PostData(1,fHistList);
259 PostData(2,fTree);
260}
261
262//___________________________________________________________________________
263void AliAnalysisTaskJetsHMPID::UserExec(Option_t *)
264{
265 if (!fAOD){
266 PostData(1,fHistList);
267 PostData(2,fTree);
268 return;
269 }
270
271 const Float_t pi = TMath::Pi(), jetR = 0.4;
272 const Float_t hmpPhiMid = 30*TMath::DegToRad();
273 const Float_t hmpPhiWid = 25*TMath::DegToRad();
274
275 TClonesArray* jets = (TClonesArray*) fAOD->FindListObject(fJetBranch.Data());
276 if (fDebug) printf("Jet Branch: %s\n",fJetBranch.Data());
277 if (fDebug) printf("Bkg Branch: %s\n",fBkgBranch.Data());
278
279 if (!jets) AliError("Jet branch not found");
280
281 Int_t nj = jets->GetEntriesFast();
282 Float_t *jetPt = new Float_t[nj];
283 Int_t *iPerm = new Int_t[nj];
284// Correct jets Pt
285 static AliAODJetEventBackground* externalBackground = 0;
286 if(fBkgBranch.Length()) externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBkgBranch.Data()));
287 Float_t pTbkg = 0;
288 AliAODJet *jet = 0;
289 for(Int_t iJet=0; iJet<nj; iJet++){
290 jet = (AliAODJet*)jets->At(iJet);
291 jetPt[iJet] = jet->Pt();
292 pTbkg = (externalBackground) ? externalBackground->GetBackground(2)*jet->EffectiveAreaCharged() : 0;
293 jetPt[iJet] -= pTbkg;
294 }
295 TMath::Sort(nj,jetPt,iPerm,kTRUE);
296
297// Event characterisation:
298// iEvSel = 0: no jets in the HMPID, consider it as background
299// iEvSel = 1: jet in the HMPID, use it for PID in jets
300// iEvSel = 2: HMPID is in the away side of a leading jet, but no jet found around
301// the detector; dismiss the event as it might alter the analysis
302
303 Int_t iEvSel = -1, iJetSel = -1;
304 Float_t phiTrig = 0;
305 if (nj && jetPt[iPerm[0]] > fJetPtCut){
306 phiTrig = ((AliAODJet*)jets->At(iPerm[0]))->Phi();
307 if (phiTrig < pi+hmpPhiMid) phiTrig-=hmpPhiMid;
308 else if (phiTrig > pi+hmpPhiMid) phiTrig-=(2*pi+hmpPhiMid);
309 if (TMath::Abs(phiTrig) < hmpPhiWid+jetR) iEvSel = 1, iJetSel = 0;
310 else if (TMath::Abs(phiTrig) > pi/3 && TMath::Abs(phiTrig) < 2*pi/3) iEvSel = 0;
311 else if (TMath::Abs(phiTrig) > 2*pi/3) iEvSel = 2, iJetSel = 0;
312 }
313 else iEvSel = 0;
314
315 for (Int_t iJet=1; iJet<jets->GetEntriesFast(); iJet++){
316 if (jetPt[iPerm[iJet]] < fJetPtCut || iEvSel == 1) break;
317 jet = (AliAODJet*)jets->At(iPerm[iJet]);
318 phiTrig = jet->Phi();
319 if (phiTrig < pi+hmpPhiMid) phiTrig-=hmpPhiMid;
320 else if (phiTrig > pi+hmpPhiMid) phiTrig-=(2*pi+hmpPhiMid);
321 if (TMath::Abs(phiTrig) < hmpPhiWid+jetR) iEvSel = 1, iJetSel = iJet;
322 }
323
324 TRandom2 rndm(0);
325 Float_t etaRnd = rndm.Uniform(-0.8,0.8);
326 Float_t phiRnd = rndm.Uniform(hmpPhiMid-hmpPhiWid-jetR,hmpPhiMid+hmpPhiWid+jetR);
327 AliAODJet *rndJet = new AliAODJet();
328 rndJet->SetPtEtaPhiM(10.,etaRnd,phiRnd,0.);
329 rndJet->SetEffArea(pi*jetR*jetR,pi*jetR*jetR);
330 Float_t ptSum=0., phi=0;
331 TRefArray *jetTracks = 0;
332 Double_t probs[5]={0, 0, 0, 0, 0};
333 switch(iEvSel){
334 case 0:
335 for (Int_t iTr = 0; iTr < fAOD->GetNumberOfTracks(); iTr++){
336 AliAODTrack *tr=(AliAODTrack*)fAOD->GetTrack(iTr);
337 if (rndJet->DeltaR(tr) < jetR) ptSum+=tr->Pt();
338 AliAODPid *pid = tr->GetDetPid();
339 if (!pid) continue;
340 if (pid->GetHMPIDsignal() > 0.){
341 phi = tr->Phi()>pi ? tr->Phi()-2*pi : tr->Phi();
342 fThetaChBkg->Fill(tr->Pt(),pid->GetHMPIDsignal());
343 if (rndJet->DeltaR(tr) < jetR) fThetaChRndCone->Fill(tr->Pt(),pid->GetHMPIDsignal());
344 fTrksEtaPhiBkg->Fill(tr->Eta(),phi);
345 pid->GetHMPIDprobs(probs);
346 Short_t charge = tr->Charge();
347 fTrackPt = tr->Pt();
348 fJetPt = 0;
349 fPionBkg = charge*(probs[0] + probs[1] + probs[2]);
350 fKaonBkg = charge*probs[3];
351 fProtBkg = charge*probs[4];
352 fPionJet = charge*1.2;
353 fKaonJet = charge*1.2;
354 fProtJet = charge*1.2;
355 fTree->Fill();
356 }
357 }
358 pTbkg = (externalBackground) ? externalBackground->GetBackground(2)*jet->EffectiveAreaCharged() : 0;
359 fRndConePt->Fill(ptSum-pTbkg);
360 rndJet->Delete();
361 break;
362 case 1:
363 jet = (AliAODJet*)jets->At(iPerm[iJetSel]);
364 fJetsPt->Fill(jetPt[iPerm[iJetSel]]);
365 phi = jet->Phi()>pi ? jet->Phi()-2*pi : jet->Phi();
366 fJetsEtaPhi->Fill(jet->Eta(),phi);
367 jetTracks = (TRefArray*)jet->GetRefTracks();
368 for (Int_t iTr=0; iTr<jetTracks->GetEntriesFast(); iTr++){
369 AliAODTrack *tr = (AliAODTrack*)jetTracks->At(iTr);
370 AliAODPid *pid = tr->GetDetPid();
371 if (!pid) continue;
372 if (pid->GetHMPIDsignal() > 0.){
373 phi = tr->Phi()>pi ? tr->Phi()-2*pi : tr->Phi();
374 fThetaChJet->Fill(tr->Pt(),pid->GetHMPIDsignal());
375 fTrksEtaPhiJet->Fill(tr->Eta(),phi);
376 pid->GetHMPIDprobs(probs);
377 Short_t charge = tr->Charge();
378 fTrackPt = tr->Pt();
379 fJetPt = jetPt[iPerm[iJetSel]];
380 fPionBkg = charge*1.2;
381 fKaonBkg = charge*1.2;
382 fProtBkg = charge*1.2;
383 fPionJet = charge*(probs[0] + probs[1] + probs[2]);
384 fKaonJet = charge*probs[3];
385 fProtJet = charge*probs[4];
386 fTree->Fill();
387 }
388 }
389 break;
390 case 2:
391 fAwayJetPt->Fill(jetPt[iPerm[iJetSel]]);
392 for (Int_t iTr = 0; iTr < fAOD->GetNumberOfTracks(); iTr++){
393 AliAODTrack *tr=(AliAODTrack*)fAOD->GetTrack(iTr);
394 AliAODPid *pid = tr->GetDetPid();
395 if (!pid) continue;
396 if (pid->GetHMPIDsignal() > 0.){
397 fThetaChBkg->Fill(tr->Pt(),pid->GetHMPIDsignal());
398 pid->GetHMPIDprobs(probs);
399 Short_t charge = tr->Charge();
400 fTrackPt = tr->Pt();
401 fJetPt = jetPt[iPerm[iJetSel]];
402 fPionBkg = charge*(2 + probs[0] + probs[1] + probs[2]);
403 fKaonBkg = charge*(2 + probs[3]);
404 fProtBkg = charge*(2 + probs[4]);
405 fPionJet = charge*1.2;
406 fKaonJet = charge*1.2;
407 fProtJet = charge*1.2;
408 fTree->Fill();
409 }
410 }
411 break;
412 }
413
414 if (iEvSel > -1){
415 if (iJetSel > -1){
416 jet = (AliAODJet*)jets->At(iPerm[iJetSel]);
417 fEvSelDPhi->Fill(jet->Phi()/pi+2*iEvSel);
418 } else {
419 for (Int_t iJet=0; iJet<jets->GetEntriesFast(); iJet++){
420 if (jetPt[iPerm[iJet]] < fJetPtCut) break;
421 jet = (AliAODJet*)jets->At(iPerm[iJet]);
422 fEvSelDPhi->Fill(jet->Phi()/pi+2*iEvSel);
423 }
424 }
425 }
426
427 delete [] jetPt;
428 delete [] iPerm;
429
430 PostData(1,fHistList);
431 PostData(2,fTree);
432 return;
433}
434
435//___________________________________________________________________________
436void AliAnalysisTaskJetsHMPID::Terminate(Option_t*)
437{
438 AliAnalysisTaskSE::Terminate();
439
440}