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