proof of principle
[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     //cout << jets_incl[ij].perp() << " " << jets_incl[ij].eta() << " " << jets_incl[ij].phi() << " " <<  jets_incl[ij].m() << endl;
158     AliESDJet *jet = new ((*fJets)[jetCount]) 
159       AliESDJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
160     jet->SetArea(fjw.GetJetArea(ij));
161 #if 0 /*tobedone*/
162     jet->GetRefTracks()->Clear();
163     vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
164     Double_t neutralE = 0;
165     for(UInt_t ic=0; ic<constituents.size(); ++ic) {
166       Int_t uid = constituents[ic].user_index();
167       if (uid>=0)
168         jet->AddTrack(tracks->At(uid));
169       else {
170         TLorentzVector *nP = static_cast<TLorentzVector*>(fNeutrals->At(-(uid+1)));
171         neutralE += nP->E();
172         jet->AddTrack(nP);
173       }
174     }
175     jet->SetNEF(neutralE/jet->E());
176 #endif
177     //jet->Print("");
178     jetCount++;
179   }
180 }
181 #else
182 //________________________________________________________________________
183 AliEsdJetTask::AliEsdJetTask(const char *name) : 
184   AliAnalysisTaskSE("AliEsdJetTask"),
185   fPrimTrCuts(0),
186   fPrimTracksName("Tracks"),
187   fJetsName("Jets"),
188   fCaloName("CaloClusters"),
189   fAlgo(1),
190   fRadius(0.4),
191   fType(0),
192   fHadCorr(0),
193   fJets(0),
194   fNeutrals(0)
195 {
196   // Standard constructor.
197   AliFatal("Compiled without FASTJET package. Task cannot be used!");
198 }
199
200 //________________________________________________________________________
201 AliEsdJetTask::~AliEsdJetTask()
202 {
203   // Destructor
204 }
205
206 //________________________________________________________________________
207 void AliEsdJetTask::UserCreateOutputObjects()
208 {
209   // Create user objects.
210 }
211
212 //________________________________________________________________________
213 void AliEsdJetTask::UserExec(Option_t *) 
214 {
215 }
216
217 //________________________________________________________________________
218 void AliEsdJetTask::Terminate(Option_t *) 
219 {
220   // Called once at the end of the analysis.
221 }
222
223 //________________________________________________________________________
224 void AliEsdJetTask::FindJets(TObjArray *, TObjArray *, Int_t, Double_t)
225 {
226   // Find jets.
227 }
228 #endif