Changes to fast embedding and jetresponse: QA plot for delta AOD, Jet eta phi selecti...
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetResponse.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TMath.h"
4 #include "TH1F.h"
5 #include "TH2F.h"
6 #include "TH3F.h"
7 #include "TCanvas.h"
8
9 #include "AliLog.h"
10
11 #include "AliAnalysisTask.h"
12 #include "AliAnalysisManager.h"
13
14 #include "AliVEvent.h"
15 #include "AliESDEvent.h"
16 #include "AliESDInputHandler.h"
17 #include "AliCentrality.h"
18 #include "AliAnalysisHelperJetTasks.h"
19 #include "AliInputEventHandler.h"
20
21 #include "AliAODEvent.h"
22 #include "AliAODJet.h"
23
24 #include "AliAnalysisTaskJetResponse.h"
25
26 ClassImp(AliAnalysisTaskJetResponse)
27
28 AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse() :
29   AliAnalysisTaskSE(),
30   fESD(0x0),
31   fAOD(0x0),
32   fOfflineTrgMask(AliVEvent::kAny),
33   fMinContribVtx(1),
34   fVtxZMin(-8.),
35   fVtxZMax(8.),
36   fEvtClassMin(1),
37   fEvtClassMax(4),
38   fCentMin(0.),
39   fCentMax(100.),
40   fJetEtaMin(-.5),
41   fJetEtaMax(.5),
42   fJetPtMin(20.),
43   fJetPtFractionMin(0.5),
44   fNMatchJets(4),
45   //fJetDeltaEta(.4),
46   //fJetDeltaPhi(.4),
47   fkNbranches(2),
48   fkEvtClasses(5),
49   fOutputList(0x0),
50   fHistEvtSelection(0x0),
51   fHistEvtClass(0x0),
52   fHistCentrality(0x0),
53   fHistPtJet(new TH1F*[fkNbranches]),
54   fHistEtaPhiJet(new TH2F*[fkNbranches]),
55   fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
56   fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
57   fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
58   fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
59   fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
60   fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
61   fHistPtFraction(new TH2F*[fkEvtClasses]),
62   fHistPtPtExtra(0x0),
63   fHistPtResponse(new TH2F*[fkEvtClasses]),
64   fHistPtSmearing(new TH2F*[fkEvtClasses]),
65   fHistDeltaR(new TH2F*[fkEvtClasses]),
66   fHistArea(new TH2F*[fkEvtClasses]),
67   fHistDeltaArea(new TH2F*[fkEvtClasses])
68 {
69   // default Constructor
70
71   fJetBranchName[0] = "";
72   fJetBranchName[1] = "";
73
74   fListJets[0] = new TList;
75   fListJets[1] = new TList;
76 }
77
78 AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse(const char *name) :
79   AliAnalysisTaskSE(name),
80   fESD(0x0),
81   fAOD(0x0),
82   fOfflineTrgMask(AliVEvent::kAny),
83   fMinContribVtx(1),
84   fVtxZMin(-8.),
85   fVtxZMax(8.),
86   fEvtClassMin(1),
87   fEvtClassMax(4),
88   fCentMin(0.),
89   fCentMax(100.),
90   fJetEtaMin(-.5),
91   fJetEtaMax(.5),
92   fJetPtMin(20.),
93   fJetPtFractionMin(0.5),
94   fNMatchJets(4),
95   //fJetDeltaEta(.4),
96   //fJetDeltaPhi(.4),
97   fkNbranches(2),
98   fkEvtClasses(5),
99   fOutputList(0x0),
100   fHistEvtSelection(0x0),
101   fHistEvtClass(0x0),
102   fHistCentrality(0x0),
103   fHistPtJet(new TH1F*[fkNbranches]),
104   fHistEtaPhiJet(new TH2F*[fkNbranches]),
105   fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
106   fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
107   fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
108   fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
109   fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
110   fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
111   fHistPtFraction(new TH2F*[fkEvtClasses]),
112   fHistPtPtExtra(0x0),
113   fHistPtResponse(new TH2F*[fkEvtClasses]),
114   fHistPtSmearing(new TH2F*[fkEvtClasses]),
115   fHistDeltaR(new TH2F*[fkEvtClasses]),
116   fHistArea(new TH2F*[fkEvtClasses]),
117   fHistDeltaArea(new TH2F*[fkEvtClasses])
118 {
119   // Constructor
120
121   fJetBranchName[0] = "";
122   fJetBranchName[1] = "";
123
124   fListJets[0] = new TList;
125   fListJets[1] = new TList;
126
127   DefineOutput(1, TList::Class());
128 }
129
130 AliAnalysisTaskJetResponse::~AliAnalysisTaskJetResponse()
131 {
132   delete fListJets[0];
133   delete fListJets[1];
134 }
135
136 void AliAnalysisTaskJetResponse::SetBranchNames(const TString &branch1, const TString &branch2)
137 {
138   fJetBranchName[0] = branch1;
139   fJetBranchName[1] = branch2;
140 }
141
142 void AliAnalysisTaskJetResponse::Init()
143 {
144
145    // check for jet branches
146    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
147       AliError("Jet branch name not set.");
148    }
149
150 }
151
152 void AliAnalysisTaskJetResponse::UserCreateOutputObjects()
153 {
154   // Create histograms
155   // Called once
156   OpenFile(1);
157   if(!fOutputList) fOutputList = new TList;
158   fOutputList->SetOwner(kTRUE);
159
160   Bool_t oldStatus = TH1::AddDirectoryStatus();
161   TH1::AddDirectory(kFALSE);
162
163
164   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 5, -0.5, 4.5);
165   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
166   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
167   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
168   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
169   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
170   
171   fHistEvtClass   = new TH1I("fHistEvtClass", "event classes", fkEvtClasses, -0.5, fkEvtClasses-0.5);
172   fHistCentrality = new TH1F("fHistCentrality", "event centrality", 100, 0., 100.);
173
174   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
175     fHistPtJet[iBranch] = new TH1F(Form("fHistPtJet%i", iBranch),
176                                           Form("p_{T} of jets from branch %i;p_{T} (GeV/c);counts", iBranch),
177                                           250, 0., 250.);
178     fHistPtJet[iBranch]->SetMarkerStyle(kFullCircle);
179
180     fHistEtaPhiJet[iBranch]    = new TH2F(Form("fHistEtaPhiJet%i", iBranch),
181                                                          Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
182                                                          48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
183         fHistEtaPhiJetCut[iBranch] = new TH2F(Form("fHistEtaPhiJetCut%i", iBranch),
184                                                          Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
185                                                          48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
186
187   }
188   
189
190   fHistPtPtExtra = new TH2F("fHistPtPtExtra", "p_{T} response;p_{T} (GeV/c);p_{T} (GeV/c)",
191                             250, 0., 250., 250, 0., 250.);
192
193   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
194     fHistDeltaEtaDeltaPhiJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJet%i",iEvtClass),
195                                                                       "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
196                                                                       101, -1.01, 1.01, 101, -1.01, 1.01);
197     fHistDeltaEtaDeltaPhiJetCut[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJetCut%i",iEvtClass),
198                                                                       "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
199                                                                       101, -1.01, 1.01, 101, -1.01, 1.01);                                                                                                
200     fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass] = new TH2F("fHistDeltaEtaDeltaPhiJetNOMatching",
201                                                                                  "#Delta#eta - #Delta#phi of jets which do not match;#Delta#eta;#Delta#phi",
202                                                                                  100, -2., 2., 100, TMath::Pi(), TMath::Pi());
203     fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
204                                                               250, 0., 250., 250, 0., 250.);
205     fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;#Deltap_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
206                                                               200, -50., 150., 250, 0., 250.);
207     fHistDeltaR[iEvtClass]     = new TH2F(Form("hist_DeltaR%i",iEvtClass), "#DeltaR of matched jets;#Deltap_{T};#DeltaR", 200, -50.,150., 60, 0.,.6);
208     fHistArea[iEvtClass]       = new TH2F(Form("hist_Area%i",iEvtClass), "jet area;#Deltap_{T};jet area", 200, -50.,150., 100, 0.,1.0);
209     fHistDeltaArea[iEvtClass]  = new TH2F(Form("hist_DeltaArea%i",iEvtClass), "#DeltaArea of matched jets;#Deltap_{T};#Deltaarea", 200, -50., 150., 60, 0., .3); 
210         
211     fHistDeltaEtaEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaJet%i", iEvtClass),
212                                                 Form("#eta - #Delta#eta of matched jets from event class %i;#eta;#Delta#eta", iEvtClass),
213                                                 60, -.6, .6, 100, -.5, .5);
214     fHistDeltaPtEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaJet%i", iEvtClass),
215                                                 Form("#eta - #Deltap_{T} of matched jets from event class %i;#eta;#Deltap_{T}", iEvtClass),
216                                                 60, -.6, .6, 200, -50., 150);
217     fHistPtFraction[iEvtClass] = new TH2F(Form("fHistPtFraction%i", iEvtClass),
218                                               Form("fraction from embedded jet in reconstructed jet;fraction (event class %i);p_{T}^{emb}", iEvtClass),
219                                                                                   100, 0., 1., 250, 0., 250.);
220   }
221
222   
223   fOutputList->Add(fHistEvtSelection);
224   fOutputList->Add(fHistEvtClass);
225   fOutputList->Add(fHistCentrality);
226
227   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
228     fOutputList->Add(fHistPtJet[iBranch]);
229     fOutputList->Add(fHistEtaPhiJet[iBranch]);
230         fOutputList->Add(fHistEtaPhiJetCut[iBranch]);
231   }
232   
233   fOutputList->Add(fHistPtPtExtra);
234   
235   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
236     fOutputList->Add(fHistDeltaEtaDeltaPhiJet[iEvtClass]);
237     fOutputList->Add(fHistDeltaEtaDeltaPhiJetCut[iEvtClass]);
238     fOutputList->Add(fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass]);
239         fOutputList->Add(fHistDeltaEtaEtaJet[iEvtClass]);
240     fOutputList->Add(fHistDeltaPtEtaJet[iEvtClass]);
241         fOutputList->Add(fHistPtFraction[iEvtClass]);
242         
243     fOutputList->Add(fHistPtResponse[iEvtClass]);
244     fOutputList->Add(fHistPtSmearing[iEvtClass]);
245     fOutputList->Add(fHistDeltaR[iEvtClass]);
246     fOutputList->Add(fHistArea[iEvtClass]);
247     fOutputList->Add(fHistDeltaArea[iEvtClass]);
248   }
249
250   // =========== Switch on Sumw2 for all histos ===========
251   for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
252     TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
253     if (h1){
254       h1->Sumw2();
255       continue;
256     }
257   }
258   TH1::AddDirectory(oldStatus);
259
260   PostData(1, fOutputList);
261 }
262
263 void AliAnalysisTaskJetResponse::UserExec(Option_t *)
264 {
265   // load events, apply event cuts, then compare leading jets
266
267   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
268       AliError("Jet branch name not set.");
269       return;
270    }
271
272   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
273   if (!fESD) {
274     AliError("ESD not available");
275     return;
276   }
277   fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
278   if (!fAOD) {
279     AliError("AOD not available");
280     return;
281   }
282
283   // -- event selection --
284   fHistEvtSelection->Fill(1); // number of events before event selection
285
286   // physics selection
287   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
288     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
289   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
290     if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
291     fHistEvtSelection->Fill(2);
292     PostData(1, fOutputList);
293     return;
294   }
295
296   // vertex selection
297   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
298   Int_t nTracksPrim = primVtx->GetNContributors();
299   if ((nTracksPrim < fMinContribVtx) ||
300       (primVtx->GetZ() < fVtxZMin) ||
301       (primVtx->GetZ() > fVtxZMax) ){
302     if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
303     fHistEvtSelection->Fill(3);
304     PostData(1, fOutputList);
305     return;
306   }
307
308   // event class selection (from jet helper task)
309   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
310   if(fDebug) Printf("Event class %d", eventClass);
311   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
312     fHistEvtSelection->Fill(4);
313     PostData(1, fOutputList);
314     return;
315   }
316
317   // centrality selection
318   AliCentrality *cent = fESD->GetCentrality();
319   Float_t centValue = cent->GetCentralityPercentile("TRK");
320   if(fDebug) printf("centrality: %f\n", centValue);
321   if (centValue < fCentMin || centValue > fCentMax){
322     fHistEvtSelection->Fill(4);
323     PostData(1, fOutputList);
324     return;
325   }
326
327   fHistEvtSelection->Fill(0); // accepted events
328   fHistEvtClass->Fill(eventClass);
329   fHistCentrality->Fill(centValue);
330   // -- end event selection --
331
332   // fetch jets
333   TClonesArray *aodJets[2];
334   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
335   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
336
337   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
338     fListJets[iJetType]->Clear();
339     if (!aodJets[iJetType]) continue;
340         
341     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
342       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
343       if (jet) fListJets[iJetType]->Add(jet);
344     }
345     fListJets[iJetType]->Sort();
346   }
347   
348   // jet matching 
349   static TArrayI aMatchIndex(fListJets[0]->GetEntries());
350   static TArrayF aPtFraction(fListJets[0]->GetEntries());
351   if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
352   if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
353   
354   // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
355   AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
356                                             fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
357                                             aMatchIndex, aPtFraction, fDebug);
358                                                                                         
359   // loop over matched jets
360   Int_t ir = -1; // index of matched reconstruced jet
361   Float_t fraction = 0.;
362   AliAODJet *jet[2]  = { 0x0, 0x0 };
363   Float_t jetEta[2]  = { 0., 0. };
364   Float_t jetPhi[2]  = { 0., 0. };
365   Float_t jetPt[2]   = { 0., 0. };
366   Float_t jetArea[2] = { 0., 0. };
367    
368   for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
369      ir = aMatchIndex[ig];
370          if(ir<0) continue;
371          fraction = aPtFraction[ig];
372          
373          // fetch jets
374          jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
375          jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
376          if(!jet[0] || !jet[1]) continue;
377          
378          for(Int_t i=0; i<fkNbranches; ++i){
379             jetEta[i]  = jet[i]->Eta();
380                 jetPhi[i]  = jet[i]->Phi();
381                 jetPt[i]   = jet[i]->Pt();
382                 jetArea[i] = jet[i]->EffectiveAreaCharged();
383          }
384          if(eventClass > -1 && eventClass < fkEvtClasses){
385         fHistPtFraction[eventClass] -> Fill(fraction, jetPt[0]);
386          }
387          
388          if(fraction<fJetPtFractionMin) continue;
389          
390          // calculate parameters of associated jets
391          Float_t deltaPt    = jetPt[1]-jetPt[0];
392          Float_t deltaEta   = jetEta[1]-jetEta[0];
393          Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
394          Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
395          Float_t deltaArea  = jetArea[1]-jetArea[0];
396          
397          // fill histograms before acceptance cut
398          for(Int_t i=0; i<fkNbranches; ++i){
399              fHistEtaPhiJet[i] -> Fill(jetEta[i], jetPhi[i]);
400          }
401          if(eventClass > -1 && eventClass < fkEvtClasses){
402              fHistDeltaEtaDeltaPhiJet[eventClass] -> Fill(deltaEta, deltaPhi);
403          }
404          
405          // jet acceptance + minimum pT check
406          if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
407             jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
408                 
409                 if(fDebug){
410                 Printf("Jet not in eta acceptance.");
411                         Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
412                         Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
413                 }
414                 continue;
415      }
416          if(jetPt[1] < fJetPtMin){
417             if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]);
418             continue;
419          }
420          
421          
422          
423          // fill histograms
424          for(Int_t i=0; i<fkNbranches; ++i){
425              fHistPtJet[i]        -> Fill(jetPt[i]);
426              fHistEtaPhiJetCut[i] -> Fill(jetEta[i], jetPhi[i]);
427          }
428          
429          fHistPtPtExtra->Fill(jetPt[0], jetPt[1]);
430          
431          if(eventClass > -1 && eventClass < fkEvtClasses){
432                  fHistDeltaEtaDeltaPhiJetCut[eventClass] -> Fill(deltaEta, deltaPhi);
433                  
434                  fHistPtResponse[eventClass]             -> Fill(jetPt[1], jetPt[0]);
435                  fHistPtSmearing[eventClass]             -> Fill(deltaPt,  jetPt[0]);
436                  
437                  fHistDeltaPtEtaJet[eventClass]          -> Fill(jetEta[0], deltaPt);
438                  fHistDeltaEtaEtaJet[eventClass]         -> Fill(jetEta[0], deltaEta);
439                  
440                  fHistDeltaR[eventClass]                 -> Fill(deltaPt, deltaR);
441                  fHistArea[eventClass]                   -> Fill(deltaPt, jetArea[1]);
442                  fHistDeltaArea[eventClass]              -> Fill(deltaPt, deltaArea);
443                  
444          }
445   }
446
447   PostData(1, fOutputList);
448 }
449
450 void AliAnalysisTaskJetResponse::Terminate(const Option_t *)
451 {
452   // Draw result to the screen
453   // Called once at the end of the query
454
455   if (!GetOutputData(1))
456     return;
457 }
458