fix AliHLTGlobalEsdConverterComponent
[u/mrichter/AliRoot.git] / HMPID / AliAnalysisTaskJetsHMPID.cxx
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
35 ClassImp(AliAnalysisTaskJetsHMPID)
36
37 //__________________________________________________________________________
38 AliAnalysisTaskJetsHMPID::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 //__________________________________________________________________________
70 AliAnalysisTaskJetsHMPID::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 //___________________________________________________________________________
106 AliAnalysisTaskJetsHMPID& 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 //___________________________________________________________________________
142 AliAnalysisTaskJetsHMPID::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 //___________________________________________________________________________
175 AliAnalysisTaskJetsHMPID::~AliAnalysisTaskJetsHMPID()
176 {
177   //
178   //destructor
179   //
180   Info("~AliAnalysisTaskJetsHMPID","Calling Destructor");
181   if (fHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHistList;
182 }
183
184 //__________________________________________________________________________
185 void 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 //___________________________________________________________________________
213 void 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 //___________________________________________________________________________
263 void 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 //___________________________________________________________________________
436 void AliAnalysisTaskJetsHMPID::Terminate(Option_t*)
437 {
438   AliAnalysisTaskSE::Terminate();
439
440 }