6588f6b984a60d5909570c897fc32ec6c84ce488
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetSpectra.cxx
1 // $Id$
2 //
3 // Jet spectrum task.
4 //
5 // Author: R.Reed, M.Connors
6
7 #include "AliAnalysisTaskEmcalJetSpectra.h"
8
9 #include <TCanvas.h>
10 #include <TChain.h>
11 #include <TClonesArray.h>
12 #include <TH1F.h>
13 #include <TH2F.h>
14 #include <TList.h>
15 #include <TLorentzVector.h>
16 #include <TParameter.h>
17 #include <TParticle.h>
18 #include <TTree.h>
19 #include <TVector3.h>
20
21 #include "AliAODEvent.h"
22 #include "AliAnalysisManager.h"
23 #include "AliAnalysisTask.h"
24 #include "AliCentrality.h"
25 #include "AliESDEvent.h"
26 #include "AliESDInputHandler.h"
27 #include "AliEmcalJet.h"
28 #include "AliVCluster.h"
29
30 ClassImp(AliAnalysisTaskEmcalJetSpectra)
31
32 //________________________________________________________________________
33 AliAnalysisTaskEmcalJetSpectra::AliAnalysisTaskEmcalJetSpectra() : 
34   AliAnalysisTaskSE(), 
35   fTracksName("tracks"),
36   fJetsName("jets"),
37   fClustersName("clusters"),
38   fRhos1Name("fArrRhos"),
39   fRhos2Name("fArrRhos"),
40   fRhos3Name("fArrRhos"),
41   fPhimin(-10), 
42   fPhimax(10),
43   fEtamin(-0.9), 
44   fEtamax(0.9),
45   fAreacut(0.0),
46   fESD(0), 
47   fOutputList(0), 
48   fHistCentrality(0),  
49   fHistDeltaRho12vsCent(0), 
50   fHistDeltaRho13vsCent(0), 
51   fHistDeltaRho23vsCent(0), 
52   fHistDeltaJetPt12vsCent(0), 
53   fHistDeltaJetPt13vsCent(0), 
54   fHistDeltaJetPt23vsCent(0), 
55   fHistRho1vsCent(0), 
56   fHistRho2vsCent(0), 
57   fHistRho3vsCent(0)
58 {
59   // Default constructor.
60
61   for (Int_t i = 0;i<6;++i) {
62     fHistRawJetPt[i]       = 0;
63     fHistAreavsRawPt[i]    = 0;
64     for (Int_t j = 0;j<4;++j) {
65       fHistNEFvsPt[i][j]   = 0;
66        fHistZvsPt[i][j]     = 0;
67        fHistZchvsPt[i][j]   = 0;
68        fHistZemvsPt[i][j]   = 0;
69        fHistJetPt[i][j]     = 0;
70        fHistNconsvsPt[i][j] = 0;
71        fHistJetPt5[i][j]    = 0;
72        fHistJetPt6[i][j]    = 0;
73        fHistJetPt7[i][j]    = 0;
74        fHistJetPt8[i][j]    = 0;
75     }
76   }
77 }
78
79 //________________________________________________________________________
80 AliAnalysisTaskEmcalJetSpectra::AliAnalysisTaskEmcalJetSpectra(const char *name) :
81   AliAnalysisTaskSE(name), 
82   fTracksName("tracks"),
83   fJetsName("jets"),
84   fClustersName("clusters"),
85   fRhos1Name("fArrRhos"),
86   fRhos2Name("fArrRhos"),
87   fRhos3Name("fArrRhos"),
88   fPhimin(-10), 
89   fPhimax(10),
90   fEtamin(-0.9), 
91   fEtamax(0.9),
92   fAreacut(0.0),
93   fESD(0), 
94   fOutputList(0), 
95   fHistCentrality(0),  
96   fHistDeltaRho12vsCent(0), 
97   fHistDeltaRho13vsCent(0), 
98   fHistDeltaRho23vsCent(0), 
99   fHistDeltaJetPt12vsCent(0), 
100   fHistDeltaJetPt13vsCent(0), 
101   fHistDeltaJetPt23vsCent(0), 
102   fHistRho1vsCent(0), 
103   fHistRho2vsCent(0), 
104   fHistRho3vsCent(0)
105 {
106   // Constructor
107
108   for (Int_t i = 0;i<6;++i) {
109     fHistRawJetPt[i]       = 0;
110     fHistAreavsRawPt[i]    = 0;
111     for (Int_t j = 0;j<4;++j) {
112       fHistNEFvsPt[i][j]   = 0;
113       fHistZvsPt[i][j]     = 0;
114       fHistZchvsPt[i][j]   = 0;
115       fHistZemvsPt[i][j]   = 0;
116       fHistJetPt[i][j]     = 0;
117       fHistNconsvsPt[i][j] = 0;
118       fHistJetPt5[i][j]    = 0;
119       fHistJetPt6[i][j]    = 0;
120       fHistJetPt7[i][j]    = 0;
121       fHistJetPt8[i][j]    = 0;
122     }
123   }
124   
125   DefineInput(0, TChain::Class());
126   DefineOutput(1, TList::Class());
127 }
128
129 //________________________________________________________________________
130 void AliAnalysisTaskEmcalJetSpectra::UserCreateOutputObjects()
131 {
132   // Called once at the beginning of the analysis.
133
134   AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
135   if (!handler) {
136     AliError("Input handler not available!");
137     return;
138   }
139
140   OpenFile(1);
141   fOutputList = new TList();
142   fOutputList->SetOwner();
143
144   fHistCentrality            = new TH1F("Centrality",             "Centrality",            101, -1,  100);
145   fHistDeltaRho12vsCent      = new TH2F("DeltaRho12vsCent",       "DeltaRho12vsCent",      100, 0.0, 100.0, 500,-250,250);
146   fHistDeltaRho13vsCent      = new TH2F("DeltaRho13vsCent",       "DeltaRho13vsCent",      100, 0.0, 100.0, 500,-250,250);
147   fHistDeltaRho23vsCent      = new TH2F("DeltaRho23vsCent",       "DeltaRho23vsCent",      100, 0.0, 100.0, 500,-250,250);
148   fHistDeltaJetPt12vsCent    = new TH2F("DeltaJetPt12vsCent",     "DeltaJetPt12vsCent",    100, 0.0, 100.0, 500,-250,250);
149   fHistDeltaJetPt13vsCent    = new TH2F("DeltaJetPt13vsCent",     "DeltaJetPt13vsCent",    100, 0.0, 100.0, 500,-250,250);
150   fHistDeltaJetPt23vsCent    = new TH2F("DeltaJetPt23vsCent",     "DeltaJetPt23vsCent",    100, 0.0, 100.0, 500,-250,250);
151   fHistRho1vsCent            = new TH2F("Rho1vsCent",             "Rho1vsCent",            100, 0.0, 100.0, 500, 0, 500);
152   fHistRho2vsCent            = new TH2F("Rho2vsCent",             "Rho2vsCent",            100, 0.0, 100.0, 500, 0, 500);
153   fHistRho3vsCent            = new TH2F("Rho3vsCent",             "Rho3vsCent",            100, 0.0, 100.0, 500, 0, 500);
154
155   for (Int_t i = 0;i<6;++i) {
156     TString name00(Form("fHistRawJetPt_%i",i));
157     fHistRawJetPt[i] = new TH1F(name00,name00,250,0,500);
158     fOutputList->Add(fHistRawJetPt[i]);
159     TString name01(Form("fHistAreavsRawPt_%i",i));
160     fHistAreavsRawPt[i] = new TH2F(name01,name01,250,0,500,100,0,1);
161     fOutputList->Add(fHistAreavsRawPt[i]);
162     for (Int_t j = 0;j<4;j++) {
163       TString name0(Form("fHistNEFvsPt_%i_%i",i,j));
164       fHistNEFvsPt[i][j] = new TH2F(name0,name0,250,-250,250,200,0,2);
165       fOutputList->Add(fHistNEFvsPt[i][j]);
166       TString name1(Form("fHistZvsPt_%i_%i",i,j));
167       fHistZvsPt[i][j] = new TH2F(name1,name1,250,-250,250,200,0,2);
168       fOutputList->Add(fHistZvsPt[i][j]);
169       TString name2(Form("fHistZchvsPt_%i_%i",i,j));
170       fHistZchvsPt[i][j] = new TH2F(name2,name2,250,-250,250,200,0,2);
171       fOutputList->Add(fHistZchvsPt[i][j]);
172       TString name3(Form("fHistZemvsPt_%i_%i",i,j));
173       fHistZemvsPt[i][j] = new TH2F(name3,name3,250,-250,250,200,0,2);
174       fOutputList->Add(fHistZemvsPt[i][j]);
175       TString name5(Form("fHistJetPt_%i_%i",i,j));
176       fHistJetPt[i][j] = new TH1F(name5,name5,250,-250,250);
177       fOutputList->Add(fHistJetPt[i][j]);
178       TString name6(Form("fHistNconsvsPt_%i_%i",i,j));
179       fHistNconsvsPt[i][j] = new TH2F(name6,name6,250,-250,250,500,0,500);
180       fOutputList->Add(fHistNconsvsPt[i][j]);
181       TString name7(Form("fHistJetPt5_%i_%i",i,j));
182       fHistJetPt5[i][j] = new TH1F(name7,name7,250,-250,250);
183       fOutputList->Add(fHistJetPt5[i][j]);
184       TString name8(Form("fHistJetPt6_%i_%i",i,j));
185       fHistJetPt6[i][j] = new TH1F(name8,name8,250,-250,250);
186       fOutputList->Add(fHistJetPt6[i][j]);
187       TString name9(Form("fHistJetPt7_%i_%i",i,j));
188       fHistJetPt7[i][j] = new TH1F(name9,name9,250,-250,250);
189       fOutputList->Add(fHistJetPt7[i][j]);
190       TString name10(Form("fHistJetPt8_%i_%i",i,j));
191       fHistJetPt8[i][j] = new TH1F(name10,name10,250,-250,250);
192       fOutputList->Add(fHistJetPt8[i][j]);
193     }
194   }
195   fOutputList->Add(fHistCentrality);
196   fOutputList->Add(fHistDeltaRho12vsCent);
197   fOutputList->Add(fHistDeltaRho13vsCent);
198   fOutputList->Add(fHistDeltaRho23vsCent);
199   fOutputList->Add(fHistDeltaJetPt12vsCent);
200   fOutputList->Add(fHistDeltaJetPt13vsCent);
201   fOutputList->Add(fHistDeltaJetPt23vsCent);
202   fOutputList->Add(fHistRho1vsCent);
203   fOutputList->Add(fHistRho2vsCent);
204   fOutputList->Add(fHistRho3vsCent);
205   
206   PostData(1, fOutputList);
207 }
208
209 //________________________________________________________________________
210 Int_t AliAnalysisTaskEmcalJetSpectra::GetCentBin(Double_t cent) const 
211 {
212   // Get centrality bin.
213
214   Int_t centbin = -1;
215   if (cent>=0 && cent<10)
216     centbin = 0;
217   else if (cent>=10 && cent<20)
218     centbin = 1;
219   else if (cent>=20 && cent<30)
220     centbin = 2;
221   else if (cent>=30 && cent<40)
222     centbin = 3;
223   else if (cent>=40 && cent<50)
224     centbin = 4;
225   else if (cent>=50 && cent<90)
226     centbin = 5;
227   return centbin;
228 }
229
230 //________________________________________________________________________
231 void AliAnalysisTaskEmcalJetSpectra::UserExec(Option_t *) 
232 {
233   // Main loop, called for each event.
234
235   // esd or aod mode
236   Bool_t esdMode = kTRUE;
237   if (dynamic_cast<AliAODEvent*>(InputEvent()))
238     esdMode = kFALSE;
239
240   if (esdMode) {
241     // optimization in case autobranch loading is off
242     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
243     if (fTracksName == "Tracks")
244       am->LoadBranch("Tracks");
245   }
246
247   // get centrality 
248   Double_t fCent = -1; 
249   TList *l = InputEvent()->GetList();
250   AliCentrality *centrality = InputEvent()->GetCentrality() ;
251   if (centrality)
252     fCent = centrality->GetCentralityPercentile("V0M");
253   else
254     fCent=99; // probably pp data
255   if (fCent<0) {
256     AliError(Form("Centrality negative: %f", fCent));
257     return;
258   }
259
260   Int_t centbin = GetCentBin(fCent);
261
262   //don't want to analyze events 90-100
263   if (centbin < 0)
264     return;
265
266   TClonesArray *jets = 0;
267   TClonesArray *tracks = 0;
268
269   tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
270   if (!tracks) {
271     AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
272     return;
273   }
274     
275   jets = dynamic_cast<TClonesArray*>(l->FindObject(fJetsName));
276   if (!jets) {
277     AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
278     return;
279   }
280
281   Double_t rho1 = -1;
282   TParameter<Double_t> *Rho1Param = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhos1Name));
283   if (Rho1Param)
284     rho1 = Rho1Param->GetVal();
285
286   Double_t rho2 = -1;
287   TParameter<Double_t> *Rho2Param = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhos2Name));
288   if (Rho2Param)
289     rho2 = Rho2Param->GetVal();
290
291   Double_t rho3 = -1;
292   TParameter<Double_t> *Rho3Param = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhos3Name));
293   if (Rho3Param)
294     rho3 = Rho3Param->GetVal();
295
296   if (rho1>0)
297     fHistRho1vsCent->Fill(fCent,rho1);
298   if (rho2>0)
299     fHistRho2vsCent->Fill(fCent,rho2);
300   if (rho3>0)
301     fHistRho3vsCent->Fill(fCent,rho3);
302
303   if (( rho1>0 ) && ( rho2>0 )) 
304     fHistDeltaRho12vsCent->Fill(fCent,rho1-rho2);
305   if (( rho1>0 ) && ( rho3>0 )) 
306     fHistDeltaRho13vsCent->Fill(fCent,rho1-rho3);
307   if (( rho2>0 ) && ( rho3>0 )) 
308     fHistDeltaRho23vsCent->Fill(fCent,rho2-rho3);
309   fHistCentrality->Fill(fCent);
310
311   Double_t lJetPt   = -500;
312   Double_t lJetPtun = -500;
313   const Int_t Njets = jets->GetEntries();
314   Int_t NjetAcc = 0;
315  
316   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
317     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(iJets));
318     if (!jet)
319       continue; 
320     if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax))
321       continue;
322     if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax))
323       continue;
324     fHistAreavsRawPt[centbin]->Fill(jet->Pt(),jet->Area());
325     if (jet->Area()<fAreacut)
326       continue;
327     //prevents 0 area jets from sneaking by when area cut == 0
328     if (jet->Area()==0)
329       continue;
330     fHistRawJetPt[centbin]->Fill(jet->Pt());
331     
332     NjetAcc++;
333     Double_t Z; 
334     Double_t jetPt1 = -500;
335     Double_t jetPt2 = -500;
336     Double_t jetPt3 = -500;
337     if (rho1 > 0)
338       jetPt1 = jet->Pt()-jet->Area()*rho1;
339     if (rho2 > 0)
340       jetPt2 = jet->Pt()-jet->Area()*rho2;
341     if (rho3 > 0)
342       jetPt3 = jet->Pt()-jet->Area()*rho3;
343     if (( rho1>0 ) && ( rho2>0 )) 
344       fHistDeltaJetPt12vsCent->Fill(fCent,jetPt1-jetPt2);
345     if (( rho1>0 ) && ( rho3>0 )) 
346       fHistDeltaJetPt13vsCent->Fill(fCent,jetPt1-jetPt3);
347     if (( rho2>0 ) && ( rho3>0 )) 
348       fHistDeltaJetPt23vsCent->Fill(fCent,jetPt2-jetPt3);
349
350     if (lJetPt < jetPt2) {
351       lJetPt= jetPt2;
352       lJetPtun = jet->Pt();
353     }
354
355     if (jet->MaxTrackPt()>jet->MaxClusterPt()) {
356       Z = jet->MaxTrackPt();
357     }
358     else{
359       Z = jet->MaxClusterPt();
360     }
361
362     fHistNEFvsPt[centbin][0]->Fill(jet->Pt(),jet->NEF());
363     fHistZvsPt[centbin][0]->Fill(jet->Pt(),Z/jet->Pt());
364     fHistZchvsPt[centbin][0]->Fill(jet->Pt(),jet->MaxTrackPt()/jet->Pt());
365     fHistZemvsPt[centbin][0]->Fill(jet->Pt(),jet->MaxClusterPt()/jet->Pt());
366     fHistJetPt[centbin][0]->Fill(jet->Pt());
367     fHistNconsvsPt[centbin][0]->Fill(jet->Pt(),jet->N());
368     if ((jet->MaxTrackPt()>5) || (jet->MaxClusterPt()>5))
369       fHistJetPt5[centbin][0]->Fill(jet->Pt());
370     if ((jet->MaxTrackPt()>6) || (jet->MaxClusterPt()>6))
371       fHistJetPt6[centbin][0]->Fill(jet->Pt());
372     if ((jet->MaxTrackPt()>7) || (jet->MaxClusterPt()>7))
373       fHistJetPt7[centbin][0]->Fill(jet->Pt());
374     if ((jet->MaxTrackPt()>8) || (jet->MaxClusterPt()>8))
375       fHistJetPt8[centbin][0]->Fill(jet->Pt());
376
377     fHistNEFvsPt[centbin][1]->Fill(jetPt1,jet->NEF());
378     fHistZvsPt[centbin][1]->Fill(jetPt1,Z/jetPt1);
379     fHistZchvsPt[centbin][1]->Fill(jetPt1,jet->MaxTrackPt()/jetPt1);
380     fHistZemvsPt[centbin][1]->Fill(jetPt1,jet->MaxClusterPt()/jetPt1);
381     fHistJetPt[centbin][1]->Fill(jetPt1);
382     fHistNconsvsPt[centbin][1]->Fill(jetPt1,jet->N());
383     if ((jet->MaxTrackPt()>5) || (jet->MaxClusterPt()>5))
384       fHistJetPt5[centbin][1]->Fill(jetPt1);
385     if ((jet->MaxTrackPt()>6) || (jet->MaxClusterPt()>6))
386       fHistJetPt6[centbin][1]->Fill(jetPt1);
387     if ((jet->MaxTrackPt()>7) || (jet->MaxClusterPt()>7))
388       fHistJetPt7[centbin][1]->Fill(jetPt1);
389     if ((jet->MaxTrackPt()>8) || (jet->MaxClusterPt()>8))
390       fHistJetPt8[centbin][1]->Fill(jetPt1);
391
392     fHistNEFvsPt[centbin][2]->Fill(jetPt2,jet->NEF());
393     fHistZvsPt[centbin][2]->Fill(jetPt2,Z/jetPt2);
394     fHistZchvsPt[centbin][2]->Fill(jetPt2,jet->MaxTrackPt()/jetPt2);
395     fHistZemvsPt[centbin][2]->Fill(jetPt2,jet->MaxClusterPt()/jetPt2);
396     fHistJetPt[centbin][2]->Fill(jetPt2);
397     fHistNconsvsPt[centbin][2]->Fill(jetPt2,jet->N());
398     if ((jet->MaxTrackPt()>5) || (jet->MaxClusterPt()>5))
399       fHistJetPt5[centbin][2]->Fill(jetPt2);
400     if ((jet->MaxTrackPt()>6) || (jet->MaxClusterPt()>6))
401       fHistJetPt6[centbin][2]->Fill(jetPt2);
402     if ((jet->MaxTrackPt()>7) || (jet->MaxClusterPt()>7))
403       fHistJetPt7[centbin][2]->Fill(jetPt2);
404     if ((jet->MaxTrackPt()>8) || (jet->MaxClusterPt()>8))
405       fHistJetPt8[centbin][2]->Fill(jetPt2);
406
407     fHistNEFvsPt[centbin][3]->Fill(jetPt3,jet->NEF());
408     fHistZvsPt[centbin][3]->Fill(jetPt3,Z/jetPt3);
409     fHistZchvsPt[centbin][3]->Fill(jetPt3,jet->MaxTrackPt()/jetPt3);
410     fHistZemvsPt[centbin][3]->Fill(jetPt3,jet->MaxClusterPt()/jetPt3);
411     fHistJetPt[centbin][3]->Fill(jetPt3);
412     fHistNconsvsPt[centbin][3]->Fill(jetPt3,jet->N());
413     if ((jet->MaxTrackPt()>5) || (jet->MaxClusterPt()>5))
414       fHistJetPt5[centbin][3]->Fill(jetPt3);
415     if ((jet->MaxTrackPt()>6) || (jet->MaxClusterPt()>6))
416       fHistJetPt6[centbin][3]->Fill(jetPt3);
417     if ((jet->MaxTrackPt()>7) || (jet->MaxClusterPt()>7))
418       fHistJetPt7[centbin][3]->Fill(jetPt3);
419     if ((jet->MaxTrackPt()>8) || (jet->MaxClusterPt()>8))
420       fHistJetPt8[centbin][3]->Fill(jetPt3);
421   }
422
423   PostData(1, fOutputList);
424 }      
425
426 //________________________________________________________________________
427 void AliAnalysisTaskEmcalJetSpectra::Terminate(Option_t *) 
428 {
429   // Called once at the end of the analysis.
430 }
431
432