first supposedly working version
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliEsdJetTask.cxx
1 // $Id$
2
3 #include "AliEsdJetTask.h"
4 #include <TClonesArray.h>
5 #include <TParticle.h>
6 #include "AliESDJet.h"
7 #include "AliAnalysisManager.h"
8 #include "AliESDtrack.h"
9 #include "AliFJWrapper.h"
10 #include "AliESDCaloCluster.h"
11
12 ClassImp(AliEsdJetTask)
13
14 #ifdef __HAVE_FJINTERFACE__
15 //________________________________________________________________________
16 AliEsdJetTask::AliEsdJetTask(const char *name) : 
17   AliAnalysisTaskSE("AliEsdJetTask"),
18   fTracksName("Tracks"),
19   fCaloName("CaloClusters"),
20   fJetsName("Jets"),
21   fAlgo(1),
22   fRadius(0.4),
23   fType(0),
24   fHadCorr(0),
25   fJets(0)
26 {
27   // Standard constructor.
28
29   if (!name)
30     return;
31
32   SetName(name);
33   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
34 }
35
36 //________________________________________________________________________
37 AliEsdJetTask::~AliEsdJetTask()
38 {
39   // Destructor
40 }
41
42 //________________________________________________________________________
43 void AliEsdJetTask::UserCreateOutputObjects()
44 {
45   // Create user objects.
46
47   fJets = new TClonesArray("AliESDJet");
48   fJets->SetName(fJetsName);
49 }
50
51 //________________________________________________________________________
52 void AliEsdJetTask::UserExec(Option_t *) 
53 {
54   // Main loop, called for each event.
55   // Add jets to event if not yet there
56
57   if (!(InputEvent()->FindListObject(fJetsName)))
58     InputEvent()->AddObject(fJets);
59
60   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
61   TClonesArray *tracks = 0;
62   TClonesArray *clus   = 0;
63   TList *l = InputEvent()->GetList();
64   if ((fType==0)||(fType==1)) {
65     if (fTracksName == "Tracks")
66       am->LoadBranch("Tracks");
67     tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
68     if (!tracks) {
69       AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
70       return;
71     }
72   }
73   if ((fType==0)||(fType==2)) {
74     if (fCaloName == "CaloClusters")
75       am->LoadBranch("CaloClusters");
76     clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
77     if (!clus) {
78       AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
79       return;
80     }
81   }
82       
83   FindJets(tracks, clus, fAlgo, fRadius);
84 }
85
86 //________________________________________________________________________
87 void AliEsdJetTask::Terminate(Option_t *) 
88 {
89   // Called once at the end of the analysis.
90
91 }
92
93 //________________________________________________________________________
94 void AliEsdJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius)
95 {
96   // Find jets.
97
98   TString name("kt");
99   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
100   if (algo>=1) {
101     name  = "antikt";
102     jalgo = fastjet::antikt_algorithm;
103   }
104
105   AliFJWrapper fjw(name, name);
106   fjw.SetR(radius);
107   fjw.SetAlgorithm(jalgo);
108   fjw.SetMaxRap(0.9);
109   fjw.Clear();
110
111   if (tracks) {
112     const Int_t Ntracks = tracks->GetEntries();
113     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
114       AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
115       if (!t)
116         continue;
117       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
118     }
119   }
120   if (clus) {
121     Double_t vertex[3] = {0, 0, 0};
122     InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
123     const Int_t Nclus = clus->GetEntries();
124     for (Int_t iClus = 0, iN = 0; iClus < Nclus; ++iClus) {
125       AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
126       if (!c->IsEMCAL())
127         continue;
128       TLorentzVector nPart;
129       c->GetMomentum(nPart, vertex);
130       Double_t energy = nPart.P();
131       if (fHadCorr>0) {
132         Int_t imin = static_cast<Int_t>(c->GetEmcCpvDistance());
133         if (imin>=0) {
134           Double_t dPhiMin = c->GetTrackDx();
135           Double_t dEtaMin = c->GetTrackDz();
136           if (dPhiMin<0.05 && dEtaMin<0.025) { // pp cuts!!!
137             AliVTrack *t = dynamic_cast<AliESDtrack*>(tracks->At(imin));
138             if (t) {
139               energy -= fHadCorr*t->P();
140               if (energy<0)
141                 continue;
142             }
143           }
144         }
145       }
146       fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), energy, -iN-1);
147     }
148   }
149
150   // run jet finder
151   fjw.Run();
152
153   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
154   for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
155     if (jets_incl[ij].perp()<1) 
156       continue;
157 #if 0
158     AliAODJet *jet = new ((*fJets)[jetCount]) 
159       AliAODJet(jets_incl[ij].px(), jets_incl[ij].py(), jets_incl[ij].pz(), jets_incl[ij].E());
160     jet->SetEffArea(fjw.GetJetArea(ij),0);
161     jet->GetRefTracks()->Clear();
162     vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
163     Double_t neutralE = 0;
164     for(UInt_t ic=0; ic<constituents.size(); ++ic) {
165       Int_t uid = constituents[ic].user_index();
166       if (uid>=0)
167         jet->AddTrack(tracks->At(uid));
168       else {
169         TLorentzVector *nP = static_cast<TLorentzVector*>(fNeutrals->At(-(uid+1)));
170         neutralE += nP->E();
171         jet->AddTrack(nP);
172       }
173     }
174     jet->SetNEF(neutralE/jet->E());
175 #endif
176     //jet->Print("");
177     jetCount++;
178   }
179 }
180 #else
181 //________________________________________________________________________
182 AliEsdJetTask::AliEsdJetTask(const char *name) : 
183   AliAnalysisTaskSE("AliEsdJetTask"),
184   fPrimTrCuts(0),
185   fPrimTracksName("Tracks"),
186   fJetsName("Jets"),
187   fCaloName("CaloClusters"),
188   fAlgo(1),
189   fRadius(0.4),
190   fType(0),
191   fHadCorr(0),
192   fJets(0),
193   fNeutrals(0)
194 {
195   // Standard constructor.
196   AliFatal("Compiled without FASTJET package. Task cannot be used!");
197 }
198
199 //________________________________________________________________________
200 AliEsdJetTask::~AliEsdJetTask()
201 {
202   // Destructor
203 }
204
205 //________________________________________________________________________
206 void AliEsdJetTask::UserCreateOutputObjects()
207 {
208   // Create user objects.
209 }
210
211 //________________________________________________________________________
212 void AliEsdJetTask::UserExec(Option_t *) 
213 {
214 }
215
216 //________________________________________________________________________
217 void AliEsdJetTask::Terminate(Option_t *) 
218 {
219   // Called once at the end of the analysis.
220 }
221
222 //________________________________________________________________________
223 void AliEsdJetTask::FindJets(TObjArray *, TObjArray *, Int_t, Double_t)
224 {
225   // Find jets.
226 }
227 #endif