]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskJetCore.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New task for analysis of jet core (L. Cunqueiro)
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetCore.cxx
CommitLineData
75bf77e3 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 "THnSparse.h"
8#include "TCanvas.h"
9
10#include "AliLog.h"
11
12#include "AliAnalysisTask.h"
13#include "AliAnalysisManager.h"
14
15#include "AliVEvent.h"
16#include "AliESDEvent.h"
17#include "AliESDInputHandler.h"
18#include "AliCentrality.h"
19#include "AliAnalysisHelperJetTasks.h"
20#include "AliInputEventHandler.h"
21#include "AliAODJetEventBackground.h"
22#include "AliAnalysisTaskFastEmbedding.h"
23
24#include "AliAODEvent.h"
25#include "AliAODJet.h"
26
27#include "AliAnalysisTaskJetCore.h"
28
29ClassImp(AliAnalysisTaskJetCore)
30
31AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
32AliAnalysisTaskSE(),
33fESD(0x0),
34fAOD(0x0),
35fBackgroundBranch(""),
36fIsPbPb(kTRUE),
37fOfflineTrgMask(AliVEvent::kAny),
38fMinContribVtx(1),
39fVtxZMin(-8.),
40fVtxZMax(8.),
41fEvtClassMin(0),
42fEvtClassMax(4),
43fCentMin(0.),
44fCentMax(100.),
45fNInputTracksMin(0),
46fNInputTracksMax(-1),
47fJetEtaMin(-.5),
48fJetEtaMax(.5),
49fJetPtMin(20.),
50fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
51fJetPtFractionMin(0.5),
52fNMatchJets(4),
53fMatchMaxDist(0.8),
54fKeepJets(kFALSE),
55fRadioFrac(0.2),
56fMinDist(0.1),
57fkNbranches(2),
58fkEvtClasses(12),
59fOutputList(0x0),
60fbEvent(kTRUE),
61fHistEvtSelection(0x0),
62fHistJetSelection(0x0),
63fh2JetSelection(0x0),
64fh2JetCoreMethod1C10(0x0),
65fh2JetCoreMethod2C10(0x0),
66fh2JetCoreMethod3C10(0x0),
67fh2JetCoreMethod1C30(0x0),
68fh2JetCoreMethod2C30(0x0),
69fh2JetCoreMethod3C30(0x0),
70fh2JetCoreMethod1C60(0x0),
71fh2JetCoreMethod2C60(0x0),
72fh2JetCoreMethod3C60(0x0)
73
74
75{
76 // default Constructor
77
78 fJetBranchName[0] = "";
79 fJetBranchName[1] = "";
80
81 fListJets[0] = new TList;
82 fListJets[1] = new TList;
83}
84
85AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
86AliAnalysisTaskSE(name),
87fESD(0x0),
88fAOD(0x0),
89fBackgroundBranch(""),
90fIsPbPb(kTRUE),
91fOfflineTrgMask(AliVEvent::kAny),
92fMinContribVtx(1),
93fVtxZMin(-8.),
94fVtxZMax(8.),
95fEvtClassMin(0),
96fEvtClassMax(4),
97fCentMin(0.),
98fCentMax(100.),
99fNInputTracksMin(0),
100fNInputTracksMax(-1),
101fJetEtaMin(-.5),
102fJetEtaMax(.5),
103fJetPtMin(20.),
104fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
105fJetPtFractionMin(0.5),
106fNMatchJets(4),
107fMatchMaxDist(0.8),
108fKeepJets(kFALSE),
109fRadioFrac(0.2),
110fMinDist(0.1),
111fkNbranches(2),
112fkEvtClasses(12),
113fOutputList(0x0),
114fbEvent(kTRUE),
115fHistEvtSelection(0x0),
116fHistJetSelection(0x0),
117fh2JetSelection(0x0),
118fh2JetCoreMethod1C10(0x0),
119fh2JetCoreMethod2C10(0x0),
120fh2JetCoreMethod3C10(0x0),
121fh2JetCoreMethod1C30(0x0),
122fh2JetCoreMethod2C30(0x0),
123fh2JetCoreMethod3C30(0x0),
124fh2JetCoreMethod1C60(0x0),
125fh2JetCoreMethod2C60(0x0),
126fh2JetCoreMethod3C60(0x0)
127
128
129 {
130 // Constructor
131
132 fJetBranchName[0] = "";
133 fJetBranchName[1] = "";
134
135 fListJets[0] = new TList;
136 fListJets[1] = new TList;
137
138 DefineOutput(1, TList::Class());
139}
140
141AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
142{
143 delete fListJets[0];
144 delete fListJets[1];
145}
146
147void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
148{
149 fJetBranchName[0] = branch1;
150 fJetBranchName[1] = branch2;
151}
152
153void AliAnalysisTaskJetCore::Init()
154{
155
156 // check for jet branches
157 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
158 AliError("Jet branch name not set.");
159 }
160
161}
162
163void AliAnalysisTaskJetCore::UserCreateOutputObjects()
164{
165 // Create histograms
166 // Called once
167 OpenFile(1);
168 if(!fOutputList) fOutputList = new TList;
169 fOutputList->SetOwner(kTRUE);
170
171 Bool_t oldStatus = TH1::AddDirectoryStatus();
172 TH1::AddDirectory(kFALSE);
173
174
175 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
176 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
177 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
178 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
179 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
180 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
181 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
182
183 fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
184 fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
185 fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
186 fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
187 fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
188 fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
189 fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
190 fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
191 fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
192
193 fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
194 fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
195 fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
196 fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
197 fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
198 fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
199 fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
200 fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
201 fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
202
203
204 // UInt_t entries = 0; // bit coded, see GetDimParams() below
205 // UInt_t opt = 0; // bit coded, default (0) or high resolution (1)
206
207 fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
208 fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
209 fh2JetCoreMethod3C10 = new TH2F("JetCoreMethod3C10","",150, 0., 150.,100, 0., 1.5);
210 fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
211 fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
212 fh2JetCoreMethod3C30 = new TH2F("JetCoreMethod3C30","",150, 0., 150.,100, 0., 1.5);
213 fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
214 fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);
215 fh2JetCoreMethod3C60 = new TH2F("JetCoreMethod3C60","",150, 0., 150.,100, 0., 1.5);
216
217
218 fOutputList->Add(fHistEvtSelection);
219 fOutputList->Add(fHistJetSelection);
220 fOutputList->Add(fh2JetSelection);
221
222
223 fOutputList->Add(fh2JetCoreMethod1C10);
224 fOutputList->Add(fh2JetCoreMethod2C10);
225 fOutputList->Add(fh2JetCoreMethod3C10);
226 fOutputList->Add(fh2JetCoreMethod1C30);
227 fOutputList->Add(fh2JetCoreMethod2C30);
228 fOutputList->Add(fh2JetCoreMethod3C30);
229 fOutputList->Add(fh2JetCoreMethod1C60);
230 fOutputList->Add(fh2JetCoreMethod2C60);
231 fOutputList->Add(fh2JetCoreMethod3C60);
232
233
234
235 // =========== Switch on Sumw2 for all histos ===========
236 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
237 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
238 if (h1){
239 h1->Sumw2();
240 continue;
241 }
242
243 }
244 TH1::AddDirectory(oldStatus);
245
246 PostData(1, fOutputList);
247}
248
249void AliAnalysisTaskJetCore::UserExec(Option_t *)
250{
251
252
253 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
254 AliError("Jet branch name not set.");
255 return;
256 }
257
258 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
259 if (!fESD) {
260 AliError("ESD not available");
261 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
262 } else {
263 fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
264 }
265 if (!fAOD) {
266 AliError("AOD not available");
267 return;
268 }
269
270 // -- event selection --
271 fHistEvtSelection->Fill(1); // number of events before event selection
272
273 // physics selection
274 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
275 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
276 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
277 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
278 fHistEvtSelection->Fill(2);
279 PostData(1, fOutputList);
280 return;
281 }
282
283 // vertex selection
284 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
285 Int_t nTracksPrim = primVtx->GetNContributors();
286 if ((nTracksPrim < fMinContribVtx) ||
287 (primVtx->GetZ() < fVtxZMin) ||
288 (primVtx->GetZ() > fVtxZMax) ){
289 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
290 fHistEvtSelection->Fill(3);
291 PostData(1, fOutputList);
292 return;
293 }
294
295 // event class selection (from jet helper task)
296 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
297 if(fDebug) Printf("Event class %d", eventClass);
298 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
299 fHistEvtSelection->Fill(4);
300 PostData(1, fOutputList);
301 return;
302 }
303
304 // centrality selection
305 AliCentrality *cent = 0x0;
306 Float_t centValue = 0.;
307 if(fESD) cent = fESD->GetCentrality();
308 if(cent) centValue = cent->GetCentralityPercentile("V0M");
309 if(fDebug) printf("centrality: %f\n", centValue);
310 if (centValue < fCentMin || centValue > fCentMax){
311 fHistEvtSelection->Fill(4);
312 PostData(1, fOutputList);
313 return;
314 }
315
316
317 // multiplicity due to input tracks
318 Int_t nInputTracks = GetNInputTracks();
319
320 if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
321 fHistEvtSelection->Fill(5);
322 PostData(1, fOutputList);
323 return;
324 }
325
326
327 fHistEvtSelection->Fill(0); // accepted events
328 // -- end event selection --
329
330
331 // get background
332 AliAODJetEventBackground* externalBackground = 0;
333 if(!externalBackground&&fBackgroundBranch.Length()){
334 externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
335 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
336 }
337 Float_t rho = 0;
338 if(externalBackground)rho = externalBackground->GetBackground(0);
339
340
341 // fetch jets
342 TClonesArray *aodJets[2];
343 aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); //
344 aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); //
345
346
347
348 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
349 fListJets[iJetType]->Clear();
350 if (!aodJets[iJetType]) continue;
351
352 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
353
354 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
355 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
356 if (jet) fListJets[iJetType]->Add(jet);
357 }
358 fListJets[iJetType]->Sort();
359 }
360
361 Double_t etabig=0;
362 Double_t ptbig=0;
363 Double_t areabig=0;
364 Double_t phibig=0.;
365 Double_t etasmall=0;
366 Double_t ptsmall=0;
367 Double_t areasmall=0;
368 Double_t distr=0.;
369 Double_t phismall=0.;
370 for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
371 AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
372 etabig = jetbig->Eta();
373 phibig = jetbig->Phi();
374 ptbig = jetbig->Pt();
375 if(ptbig==0) continue;
376 areabig = jetbig->EffectiveAreaCharged();
377 if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
378 Double_t dismin=100.;
379 Double_t ptmax=-10.;
380 Int_t index1=-1;
381 Int_t index2=-1;
382 Double_t fracin=0.;
383 //compute sum of the pt of the tracks in a concentric cone
384 TRefArray *genTrackList = jetbig->GetRefTracks();
385 Int_t nTracksGenJet = genTrackList->GetEntriesFast();
386 AliAODTrack* genTrack;
387 for(Int_t ir=0; ir<nTracksGenJet; ++ir){
388 genTrack = (AliAODTrack*)(genTrackList->At(ir));
389 Float_t etr=genTrack->Eta();
390 Float_t phir=genTrack->Phi();
391 distr=(etr-etabig)*(etr-etabig)+(phir-phibig)*(phir-phibig);
392 distr=TMath::Sqrt(distr);
393 if(distr<=fRadioFrac){ fracin=fracin+genTrack->Pt();}}
394 if(centValue<10) fh2JetCoreMethod3C10->Fill(ptbig-rho*areabig,fracin/ptbig);
395 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod3C30->Fill(ptbig-rho*areabig,fracin/ptbig);
396 if(centValue>60) fh2JetCoreMethod3C60->Fill(ptbig-rho*areabig,fracin/ptbig);
397 //cout<<"method3 "<<fracin/ptbig<<" "<<ptbig-rho*areabig<<endl;
398 /////////////////////////////////////////////////////////////
399
400 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
401 AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
402 etasmall = jetsmall->Eta();
403 phismall = jetsmall->Phi();
404 ptsmall = jetsmall->Pt();
405 areasmall = jetsmall->EffectiveAreaCharged();
406
407 if((etasmall<fJetEtaMin)||(etasmall>fJetEtaMax)) continue;
408
409 Float_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
410 tmpDeltaR=TMath::Sqrt(tmpDeltaR);
411
412 if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
413 index2=j;}
414
415 if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
416 index1=j;}}
417
418 //method1:most concentric jet=core
419 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
420
421 if(centValue<10) fh2JetCoreMethod1C10->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig);
422 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig);
423 if(centValue>60) fh2JetCoreMethod1C60->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig);
424
425 //cout<<"method1 "<<jetmethod1->Pt()/ptbig<<" "<<ptbig<<endl;
426 }
427 //method2:hardest contained jet=core
428 if(index2!=-1){
429 AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
430 if(centValue<10) fh2JetCoreMethod2C10->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig);
431 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig);
432 if(centValue>60) fh2JetCoreMethod2C60->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig);
433
434 //cout<<"method2 "<<jetmethod2->Pt()/ptbig<<" "<<ptbig<<endl;
435 }
436
437
438
439
440 }
441
442
443
444
445 PostData(1, fOutputList);
446}
447
448void AliAnalysisTaskJetCore::Terminate(const Option_t *)
449{
450 // Draw result to the screen
451 // Called once at the end of the query
452
453 if (!GetOutputData(1))
454 return;
455}
456
457Int_t AliAnalysisTaskJetCore::GetNInputTracks()
458{
459
460 Int_t nInputTracks = 0;
461
462 TString jbname(fJetBranchName[1]);
463 //needs complete event, use jets without background subtraction
464 for(Int_t i=1; i<=3; ++i){
465 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
466 }
467 // use only HI event
468 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
469 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
470
471 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
472 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
473 if(!tmpAODjets){
474 Printf("Jet branch %s not found", jbname.Data());
475 Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
476 return -1;
477 }
478
479 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
480 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
481 if(!jet) continue;
482 TRefArray *trackList = jet->GetRefTracks();
483 Int_t nTracks = trackList->GetEntriesFast();
484 nInputTracks += nTracks;
485 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
486 }
487 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
488
489 return nInputTracks;
490}
491
492
493
494
495