1 // **************************************************************************************
3 // * Task for Jet properties and jet shape analysis in PWG4 Jet Task Force Train for pp *
5 // **************************************************************************************
8 /**************************************************************************
9 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
11 * Author: The ALICE Off-line Project. *
12 * Contributors are mentioned in the code where appropriate. *
14 * Permission to use, copy, modify and distribute this software and its *
15 * documentation strictly for non-commercial purposes is hereby granted *
16 * without fee, provided that the above copyright notice appears in all *
17 * copies and that both the copyright notice and this permission notice *
18 * appear in the supporting documentation. The authors make no claims *
19 * about the suitability of this software for any purpose. It is *
20 * provided "as is" without express or implied warranty. *
21 **************************************************************************/
30 #include "THnSparse.h"
36 #include "AliAODInputHandler.h"
37 #include "AliAODHandler.h"
38 #include "AliESDEvent.h"
39 #include "AliAODMCParticle.h"
40 #include "AliAODJet.h"
41 #include "AliGenPythiaEventHeader.h"
42 #include "AliGenHijingEventHeader.h"
43 #include "AliInputEventHandler.h"
45 #include "AliAnalysisHelperJetTasks.h"
46 #include "AliAnalysisManager.h"
47 #include "AliAnalysisTaskSE.h"
48 #include "AliVParticle.h"
49 #include "AliAODTrack.h"
50 #include "AliVEvent.h"
53 #include "AliAnalysisTaskJetProperties.h"
55 ClassImp(AliAnalysisTaskJetProperties)
57 //_________________________________________________________________________________//
59 AliAnalysisTaskJetProperties::AliAnalysisTaskJetProperties()
67 ,fTrackType(kTrackAOD)
69 ,fUseAODInputJets(kTRUE)
71 ,fUsePhysicsSelection(kTRUE)
85 ,fh1VertexNContributors(0)
105 ,fh2NtracksLeadingJet(0)
106 ,fProNtracksLeadingJet(0)
118 for(Int_t ii=0; ii<5; ii++){
119 fProDelRNchSum[ii] = NULL;
120 fProDelRPtSum[ii] = NULL;
121 fProDiffJetShapeA[ii] = NULL;
122 fProIntJetShapeA[ii] = NULL;
124 // default constructor
126 //_________________________________________________________________________________//
128 AliAnalysisTaskJetProperties::AliAnalysisTaskJetProperties(const char *name)
129 : AliAnalysisTaskSE(name)
136 ,fTrackType(kTrackAOD)
138 ,fUseAODInputJets(kTRUE)
140 ,fUsePhysicsSelection(kTRUE)
154 ,fh1VertexNContributors(0)
174 ,fh2NtracksLeadingJet(0)
175 ,fProNtracksLeadingJet(0)
187 for(Int_t ii=0; ii<5; ii++){
188 fProDelRNchSum[ii] = NULL;
189 fProDelRPtSum[ii] = NULL;
190 fProDiffJetShapeA[ii] = NULL;
191 fProIntJetShapeA[ii] = NULL;
194 DefineOutput(1,TList::Class());
196 //_________________________________________________________________________________//
198 AliAnalysisTaskJetProperties::~AliAnalysisTaskJetProperties()
201 if(fJetList) delete fJetList;
202 if(fTrackList) delete fTrackList;
204 //_________________________________________________________________________________//
206 Bool_t AliAnalysisTaskJetProperties::Notify()
209 // Implemented Notify() to read the cross sections
210 // and number of trials from pyxsec.root
211 // (taken from AliAnalysisTaskJetSpectrum2)
213 if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::Notify()");
215 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
216 Float_t xsection = 0;
221 TFile *curfile = tree->GetCurrentFile();
223 Error("Notify","No current file");
226 if(!fh1Xsec||!fh1Trials){
227 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
230 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
231 fh1Xsec->Fill("<#sigma>",xsection);
232 // construct a poor man average trials
233 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
234 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
238 //_________________________________________________________________________________//
240 void AliAnalysisTaskJetProperties::UserCreateOutputObjects()
242 //(Here, event selection part is taken from AliAnalysisTaskFramentationFunction)
243 // create output objects
244 if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::UserCreateOutputObjects()");
246 // create list of tracks and jets
247 fJetList = new TList();
248 fJetList->SetOwner(kFALSE);
249 fTrackList = new TList();
250 fTrackList->SetOwner(kFALSE);
252 // Create histograms / output container
254 fCommonHistList = new TList();
255 fCommonHistList->SetOwner();
257 Bool_t oldStatus = TH1::AddDirectoryStatus();
258 TH1::AddDirectory(kFALSE);
259 // Histograms/TProfiles
260 fh1EvtSelection = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
261 fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
262 fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
263 fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
264 fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
265 fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
266 fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
269 fh1VertexNContributors = new TH1F("fh1VertexNContributors", "Vertex N contributors", 2500,-.5, 2499.5);
270 fh1VertexZ = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
271 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
272 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
273 fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
274 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
275 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
276 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
279 Int_t kNbinsPtSlice=20; Float_t xMinPtSlice=0.0; Float_t xMaxPtSlice=100.0;
280 Int_t kNbinsEta=40; Float_t xMinEta=-2.0; Float_t xMaxEta=2.0;
281 Int_t kNbinsPhi=90; Float_t xMinPhi=0.0; Float_t xMaxPhi=TMath::TwoPi();
282 Int_t kNbinsPt=125; Float_t xMinPt=0.0; Float_t xMaxPt=250.0;
283 Int_t kNbinsNtracks=50; Float_t xMinNtracks=0.0; Float_t xMaxNtracks=50.0;
284 Int_t kNbinsFF=80; Float_t xMinFF=0.0; Float_t xMaxFF=2.0;
285 Int_t kNbinsDelR1D=50; Float_t xMinDelR1D=0.0; Float_t xMaxDelR1D=1.0;
288 fh2EtaJet = new TH2F("EtaJet","EtaJet",
289 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
290 kNbinsEta, xMinEta, xMaxEta);
291 fh2PhiJet = new TH2F("PhiJet","PhiJet",
292 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
293 kNbinsPhi, xMinPhi, xMaxPhi);
294 fh2PtJet = new TH2F("PtJet","PtJet",
295 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
296 kNbinsPt, xMinPt, xMaxPt);
297 fh1PtJet = new TH1F("PtJet1D","PtJet1D",
298 kNbinsPt, xMinPt, xMaxPt);
299 fh2NtracksJet = new TH2F("NtracksJet","NtracksJet",
300 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
301 kNbinsNtracks, xMinNtracks, xMaxNtracks);
302 fProNtracksJet = new TProfile("AvgNoOfTracksJet","AvgNoOfTracksJet",
303 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
304 xMinNtracks, xMaxNtracks);
305 fh2EtaTrack = new TH2F("EtaTrack","EtaTrack",
306 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
307 kNbinsEta, xMinEta, xMaxEta);
308 fh2PhiTrack = new TH2F("PhiTrack","PhiTrack",
309 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
310 kNbinsPhi, xMinPhi, xMaxPhi);
311 fh2PtTrack = new TH2F("PtTrack","PtTrack",
312 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
313 kNbinsPt, xMinPt, xMaxPt);
314 fh2FF = new TH2F("FF","FF",
315 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
316 kNbinsFF, xMinFF, xMaxFF);
317 fh2DelEta = new TH2F("DelEta","DelEta",
318 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
319 kNbinsEta, xMinEta, xMaxEta);
320 fh2DelPhi = new TH2F("DelPhi","DelPhi",
321 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
322 kNbinsPhi, xMinPhi, xMaxPhi);
323 fh2DelR = new TH2F("DelR","DelR",
324 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
325 kNbinsDelR1D, xMinDelR1D, xMaxDelR1D);
329 Int_t kNbinsDelR=100; Float_t xMinDelR=0.0; Float_t xMaxDelR=1.0;
330 Int_t kNbinsPtSliceJS=100; Float_t xMinPtSliceJS=0.0; Float_t xMaxPtSliceJS=250.0;
332 fh1PtLeadingJet = new TH1F("PtLeadingJet","PtLeadingJet",
333 kNbinsPt, xMinPt, xMaxPt);
334 fh2NtracksLeadingJet = new TH2F("NtracksLeadingJet","NtracksLeadingJet",
335 kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
336 kNbinsNtracks, xMinNtracks, xMaxNtracks);
337 fProNtracksLeadingJet = new TProfile("AvgNoOfTracksLeadingJet","AvgNoOfTracksLeadingJet",
338 kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
339 xMinNtracks, xMaxNtracks);
340 fh2DelR80pcNch = new TH2F("delR80pcNch","delR80pcNch",
341 kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
342 kNbinsDelR, xMinDelR, xMaxDelR);
343 fProDelR80pcNch = new TProfile("AvgdelR80pcNch","AvgdelR80pcNch",
344 kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
346 fh2DelR80pcPt = new TH2F("delR80pcPt","delR80pcPt",
347 kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
348 kNbinsDelR, xMinDelR, xMaxDelR);
349 fProDelR80pcPt = new TProfile("AvgdelR80pcPt","AvgdelR80pcPt",
350 kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
352 fh2AreaCh = new TH2F("Area","Area",
353 kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
354 72, 0.0, TMath::Pi());
355 fProAreaCh = new TProfile("AvgArea","AvgArea",
356 kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
358 fh3PtDelRNchSum = new TH3F("PtDelRNchSum","PtDelRNchSum",
359 kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
360 kNbinsDelR1D, xMinDelR1D, xMaxDelR1D,
361 kNbinsNtracks, xMinNtracks, xMaxNtracks);
362 fh3PtDelRPtSum = new TH3F("PtDelRPtSum","PtDelRPtSum",
363 kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
364 kNbinsDelR1D, xMinDelR1D, xMaxDelR1D,
365 kNbinsPt, xMinPt, xMaxPt);
366 fProDiffJetShape = new TProfile("DiffJetShape","DiffJetShape",
367 10,0.0,1.0,0.0,250.0);
368 fProIntJetShape = new TProfile("IntJetShape","IntJetShape",
369 10,0.0,1.0,0.0,250.0);
372 for(Int_t ii=0; ii<5; ii++){
373 if(ii==0)title = "_JetPt20to30";
374 if(ii==1)title = "_JetPt30to40";
375 if(ii==2)title = "_JetPt40to60";
376 if(ii==3)title = "_JetPt60to80";
377 if(ii==4)title = "_JetPt80to100";
378 fProDelRNchSum[ii] = new TProfile(Form("AvgNchSumDelR%s",title.Data()),Form("AvgNchSumDelR%s",title.Data()),
379 kNbinsDelR1D, xMinDelR1D, xMaxDelR1D,
380 xMinNtracks, xMaxNtracks);
382 fProDelRPtSum[ii] = new TProfile(Form("AvgPtSumDelR%s",title.Data()),Form("AvgPtSumDelR%s",title.Data()),
383 kNbinsDelR1D, xMinDelR1D, xMaxDelR1D,
385 fProDelRNchSum[ii] ->GetXaxis()->SetTitle("R");
386 fProDelRNchSum[ii] ->GetYaxis()->SetTitle("<NchSum>");
387 fProDelRPtSum[ii] ->GetXaxis()->SetTitle("R");
388 fProDelRPtSum[ii] ->GetYaxis()->SetTitle("<PtSum>");
390 fProDiffJetShapeA[ii] = new TProfile(Form("DiffJetShape%s",title.Data()),Form("DiffJetShape%s",title.Data()),
391 10,0.0,1.0,0.0,250.0);
392 fProIntJetShapeA[ii] = new TProfile(Form("IntJetShape%s",title.Data()),Form("IntJetShape%s",title.Data()),
393 10,0.0,1.0,0.0,250.0);
395 fProDiffJetShapeA[ii]->GetXaxis()->SetTitle("R");
396 fProDiffJetShapeA[ii]->GetYaxis()->SetTitle("Diff jet shape");
397 fProIntJetShapeA[ii]->GetXaxis()->SetTitle("R");
398 fProIntJetShapeA[ii]->GetYaxis()->SetTitle("Integrated jet shape");
400 fCommonHistList->Add(fProDelRNchSum[ii]);
401 fCommonHistList->Add(fProDelRPtSum[ii]);
402 fCommonHistList->Add(fProDiffJetShapeA[ii]);
403 fCommonHistList->Add(fProIntJetShapeA[ii]);
407 fh2EtaJet ->GetXaxis()->SetTitle("JetPt"); fh2EtaJet ->GetYaxis()->SetTitle("JetEta");
408 fh2PhiJet ->GetXaxis()->SetTitle("JetPt"); fh2PhiJet ->GetYaxis()->SetTitle("JetPhi");
409 fh2PtJet ->GetXaxis()->SetTitle("JetPt"); fh2PtJet ->GetYaxis()->SetTitle("JetPt");
410 fh1PtJet ->GetXaxis()->SetTitle("JetPt"); fh1PtJet ->GetYaxis()->SetTitle("#jets");
411 fh2NtracksJet ->GetXaxis()->SetTitle("JetPt"); fh2NtracksJet ->GetYaxis()->SetTitle("#tracks");
412 fProNtracksJet->GetXaxis()->SetTitle("JetPt"); fProNtracksJet->GetYaxis()->SetTitle("AgvNtracks");
413 fh2EtaTrack ->GetXaxis()->SetTitle("JetPt"); fh2EtaTrack ->GetYaxis()->SetTitle("TrackEta");
414 fh2PhiTrack ->GetXaxis()->SetTitle("JetPt"); fh2PhiTrack ->GetYaxis()->SetTitle("TrackPhi");
415 fh2PtTrack ->GetXaxis()->SetTitle("JetPt"); fh2PtTrack ->GetYaxis()->SetTitle("TrackPt");
416 fh2FF ->GetXaxis()->SetTitle("JetPt"); fh2FF ->GetYaxis()->SetTitle("FF");
417 fh2DelEta ->GetXaxis()->SetTitle("JetPt"); fh2DelEta ->GetYaxis()->SetTitle("DelEta");
418 fh2DelPhi ->GetXaxis()->SetTitle("JetPt"); fh2DelPhi ->GetYaxis()->SetTitle("DelPhi");
419 fh2DelR ->GetXaxis()->SetTitle("JetPt"); fh2DelR ->GetYaxis()->SetTitle("DelR");
421 fh1PtLeadingJet ->GetXaxis()->SetTitle("JetPt"); fh1PtLeadingJet ->GetYaxis()->SetTitle("#leading jets");
422 fh2NtracksLeadingJet ->GetXaxis()->SetTitle("JetPt"); fh2NtracksLeadingJet ->GetYaxis()->SetTitle("#tracks leading jet");
423 fProNtracksLeadingJet ->GetXaxis()->SetTitle("JetPt"); fProNtracksLeadingJet ->GetYaxis()->SetTitle("AvgNtracks leading jet");
424 fh2DelR80pcNch ->GetXaxis()->SetTitle("JetPt"); fh2DelR80pcNch ->GetYaxis()->SetTitle("R containing 80% of tracks");
425 fProDelR80pcNch ->GetXaxis()->SetTitle("JetPt"); fProDelR80pcNch ->GetYaxis()->SetTitle("<R> containing 80% of tracks");
426 fh2DelR80pcPt ->GetXaxis()->SetTitle("JetPt"); fh2DelR80pcPt ->GetYaxis()->SetTitle("R containing 80% of pT");
427 fProDelR80pcPt ->GetXaxis()->SetTitle("JetPt"); fProDelR80pcPt ->GetYaxis()->SetTitle("<R> containing 80% of pT");
428 fh2AreaCh ->GetXaxis()->SetTitle("JetPt"); fh2AreaCh ->GetYaxis()->SetTitle("Jet area");
429 fProAreaCh ->GetXaxis()->SetTitle("JetPt"); fProAreaCh ->GetYaxis()->SetTitle("<jet area>");
430 fh3PtDelRNchSum ->GetXaxis()->SetTitle("JetPt"); fh3PtDelRNchSum ->GetYaxis()->SetTitle("R"); fh3PtDelRNchSum->GetZaxis()->SetTitle("NchSum");
431 fh3PtDelRPtSum ->GetXaxis()->SetTitle("JetPt"); fh3PtDelRPtSum ->GetYaxis()->SetTitle("R"); fh3PtDelRPtSum ->GetZaxis()->SetTitle("PtSum");
432 fProDiffJetShape->GetXaxis()->SetTitle("R");
433 fProDiffJetShape->GetYaxis()->SetTitle("Diff jet shape");
434 fProIntJetShape->GetXaxis()->SetTitle("R");
435 fProIntJetShape->GetYaxis()->SetTitle("Integrated jet shape");
437 fCommonHistList->Add(fh1EvtSelection);
438 fCommonHistList->Add(fh1VertexNContributors);
439 fCommonHistList->Add(fh1VertexZ);
440 fCommonHistList->Add(fh1Xsec);
441 fCommonHistList->Add(fh1Trials);
442 fCommonHistList->Add(fh1PtHard);
443 fCommonHistList->Add(fh1PtHardTrials);
444 fCommonHistList->Add(fh2EtaJet);
445 fCommonHistList->Add(fh2PhiJet);
446 fCommonHistList->Add(fh2PtJet);
447 fCommonHistList->Add(fh1PtJet);
448 fCommonHistList->Add(fh2NtracksJet);
449 fCommonHistList->Add(fProNtracksJet);
450 fCommonHistList->Add(fh2EtaTrack);
451 fCommonHistList->Add(fh2PhiTrack);
452 fCommonHistList->Add(fh2PtTrack);
453 fCommonHistList->Add(fh2FF);
454 fCommonHistList->Add(fh2DelEta);
455 fCommonHistList->Add(fh2DelPhi);
456 fCommonHistList->Add(fh2DelR);
458 fCommonHistList->Add(fh1PtLeadingJet);
459 fCommonHistList->Add(fh2NtracksLeadingJet);
460 fCommonHistList->Add(fProNtracksLeadingJet);
461 fCommonHistList->Add(fh2DelR80pcNch);
462 fCommonHistList->Add(fProDelR80pcNch);
463 fCommonHistList->Add(fh2DelR80pcPt);
464 fCommonHistList->Add(fProDelR80pcPt);
465 fCommonHistList->Add(fh2AreaCh);
466 fCommonHistList->Add(fProAreaCh);
467 fCommonHistList->Add(fh3PtDelRNchSum);
468 fCommonHistList->Add(fh3PtDelRPtSum);
469 fCommonHistList->Add(fProDiffJetShape);
470 fCommonHistList->Add(fProIntJetShape);
471 // =========== Switch on Sumw2 for all histos ===========
472 for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
473 TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
476 TProfile *hPro = dynamic_cast<TProfile*>(fCommonHistList->At(i));
477 if(hPro) hPro->Sumw2();
481 TH1::AddDirectory(oldStatus);
482 PostData(1, fCommonHistList);
484 //_________________________________________________________________________________//
486 void AliAnalysisTaskJetProperties::Init()
489 if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::Init()");
492 //_________________________________________________________________________________//
494 void AliAnalysisTaskJetProperties::UserExec(Option_t *)
497 // Called for each event
498 if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::UserExec()");
500 if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
503 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
504 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
505 if(!(inputHandler->IsEventSelected() & AliVEvent::kMB)){
506 if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
507 fh1EvtSelection->Fill(1.);
508 if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
509 PostData(1, fCommonHistList);
514 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
516 if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
519 fMCEvent = MCEvent();
521 if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
524 // get AOD event from input/ouput
525 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
526 if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
527 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
528 if(fUseAODInputJets) fAODJets = fAOD;
529 if (fDebug > 1) Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
532 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
533 if( handler && handler->InheritsFrom("AliAODHandler") ) {
534 fAOD = ((AliAODHandler*)handler)->GetAOD();
536 if (fDebug > 1) Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
540 if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
541 TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
542 if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
543 fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
544 if (fDebug > 1) Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
548 if(fNonStdFile.Length()!=0){
549 // case we have an AOD extension - fetch the jets from the extended output
551 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
552 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
554 if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
559 Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
563 Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
568 // *** vertex cut ***
569 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
570 Int_t nTracksPrim = primVtx->GetNContributors();
571 fh1VertexNContributors->Fill(nTracksPrim);
573 if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
575 if(nTracksPrim<fNContributors){
576 if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
577 fh1EvtSelection->Fill(3.);
578 PostData(1, fCommonHistList);
582 fh1VertexZ->Fill(primVtx->GetZ());
584 if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
585 if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
586 fh1EvtSelection->Fill(4.);
587 PostData(1, fCommonHistList);
591 TString primVtxName(primVtx->GetName());
593 if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
594 if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
595 fh1EvtSelection->Fill(5.);
596 PostData(1, fCommonHistList);
600 if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
601 fh1EvtSelection->Fill(0.);
602 //___ get MC information __________________________________________________________________
603 Double_t ptHard = 0.;
604 Double_t nTrials = 1; // trials for MC trigger weight for real data
607 AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
609 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
610 AliGenHijingEventHeader* hijingGenHeader = 0x0;
613 if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
614 nTrials = pythiaGenHeader->Trials();
615 ptHard = pythiaGenHeader->GetPtHard();
617 fh1PtHard->Fill(ptHard);
618 fh1PtHardTrials->Fill(ptHard,nTrials);
621 } else { // no pythia, hijing?
623 if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
625 hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
626 if(!hijingGenHeader){
627 Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
629 if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
633 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
637 //___ fetch jets __________________________________________________________________________
638 Int_t nJ = GetListOfJets(fJetList);
640 if(nJ>=0) nJets = fJetList->GetEntries();
642 Printf("%s:%d Selected jets: %d %d",(char*)__FILE__,__LINE__,nJ,nJets);
643 if(nJ != nJets) Printf("%s:%d Mismatch Selected Jets: %d %d",(char*)__FILE__,__LINE__,nJ,nJets);
645 FillJetProperties(fJetList);
646 FillJetShape(fJetList);
649 PostData(1, fCommonHistList);
651 //_________________________________________________________________________________//
653 void AliAnalysisTaskJetProperties::Terminate(Option_t *)
657 if(fDebug > 1) printf("AliAnalysisTaskJetProperties::Terminate() \n");
659 //_________________________________________________________________________________//
661 Int_t AliAnalysisTaskJetProperties::GetListOfJets(TList *list)
663 //this functionality is motivated by AliAnalysisTaskFragmentationFunction
664 if(fDebug > 1) printf("AliAnalysisTaskJetProperties::GetListOfJets() \n");
665 // fill list of jets selected according to type
667 if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
671 if(fBranchJets.Length()==0){
672 Printf("%s:%d no jet branch specified", (char*)__FILE__,__LINE__);
673 if(fDebug>1)fAOD->Print();
676 TClonesArray *aodJets = 0;
677 if(fBranchJets.Length()) aodJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchJets.Data()));
678 if(!aodJets) aodJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchJets.Data()));
679 if(fAODExtension&&!aodJets) aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchJets.Data()));
682 if(fBranchJets.Length())Printf("%s:%d no jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchJets.Data());
683 if(fDebug>1)fAOD->Print();
688 for(Int_t ijet=0; ijet<aodJets->GetEntries(); ijet++){
689 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodJets->At(ijet));
691 if( tmp->Pt() < fJetPtCut ) continue;
692 if( tmp->Eta() < fJetEtaMin || tmp->Eta() > fJetEtaMax)continue;
693 if(fJetRejectType==kReject1Track && tmp->GetRefTracks()->GetEntriesFast()==1)continue;//rejecting 1track jet if...
700 //_________________________________________________________________________________//
702 Int_t AliAnalysisTaskJetProperties::GetListOfJetTracks(TList* list, const AliAODJet* jet)
704 //this functionality is motivated by AliAnalysisTaskFragmentationFunction
705 if(fDebug > 1) printf("AliAnalysisTaskJetProperties::GetListOfJetTracks() \n");
707 // list of jet tracks from trackrefs
708 Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
710 for (Int_t itrack=0; itrack<nTracks; itrack++) {
711 if(fTrackType==kTrackUndef){
712 if(fDebug>2)Printf("%s:%d unknown track type %d in AOD", (char*)__FILE__,__LINE__,kTrackUndef);
715 else if(fTrackType == kTrackAOD){
716 AliAODTrack* trackAod = dynamic_cast<AliAODTrack*>(jet->GetRefTracks()->At(itrack));
718 AliError("expected ref track not found ");
725 else if(fTrackType==kTrackAODMC){
726 AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(jet->GetRefTracks()->At(itrack));
728 AliError("expected ref trackmc not found ");
735 else if (fTrackType==kTrackKine){
736 AliVParticle* trackkine = dynamic_cast<AliVParticle*>(jet->GetRefTracks()->At(itrack));
738 AliError("expected ref trackkine not found ");
741 list->Add(trackkine);
748 //_________________________________________________________________________________//
750 void AliAnalysisTaskJetProperties::FillJetProperties(TList *jetlist){
751 //filling up the histograms jets and tracks inside jet
752 if(fDebug > 1) printf("AliAnalysisTaskJetProperties::FillJetProperties() \n");
754 for(Int_t iJet=0; iJet < jetlist->GetEntries(); iJet++){
755 Float_t JetPt;Float_t JetPhi;Float_t JetEta;Float_t JetE;
756 AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(iJet));
762 fh2EtaJet ->Fill(JetPt,JetEta);
763 fh2PhiJet ->Fill(JetPt,JetPhi);
764 fh2PtJet ->Fill(JetPt,JetPt);
765 fh1PtJet ->Fill(JetPt);
767 Int_t nJT = GetListOfJetTracks(fTrackList,jet);
768 Int_t nJetTracks = 0;
769 if(nJT>=0) nJetTracks = fTrackList->GetEntries();
771 Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
772 if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
775 fh2NtracksJet ->Fill(JetPt,fTrackList->GetEntries());
776 fProNtracksJet ->Fill(JetPt,fTrackList->GetEntries());
778 for (Int_t j =0; j< fTrackList->GetEntries(); j++){
779 if(fTrackType==kTrackUndef)continue;
780 Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
781 Float_t FF=-99.0; Float_t DelEta=-99.0; Float_t DelPhi=-99.0; Float_t DelR=-99.0; Float_t AreaJ=-99.0;
782 if(fTrackType==kTrackAOD){
783 AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackList->At(j));
784 if(!trackaod)continue;
785 TrackEta = trackaod->Eta();
786 TrackPhi = trackaod->Phi();
787 TrackPt = trackaod->Pt();
789 else if(fTrackType==kTrackAODMC){
790 AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackList->At(j));
791 if(!trackmc)continue;
792 TrackEta = trackmc->Eta();
793 TrackPhi = trackmc->Phi();
794 TrackPt = trackmc->Pt();
796 else if(fTrackType==kTrackKine){
797 AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackList->At(j));
798 if(!trackkine)continue;
799 TrackEta = trackkine->Eta();
800 TrackPhi = trackkine->Phi();
801 TrackPt = trackkine->Pt();
803 if(JetPt)FF = TrackPt/JetPt;
804 DelEta = TMath::Abs(JetEta - TrackEta);
805 DelPhi = TMath::Abs(JetPhi - TrackPhi);
806 if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
807 DelR = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
808 AreaJ = TMath::Pi()*DelR*DelR;
809 fh2EtaTrack ->Fill(JetPt,TrackEta);
810 fh2PhiTrack ->Fill(JetPt,TrackPhi);
811 fh2PtTrack ->Fill(JetPt,TrackPt);
812 fh2FF ->Fill(JetPt,FF);
813 fh2DelEta ->Fill(JetPt,DelEta);
814 fh2DelPhi ->Fill(JetPt,DelPhi);
815 fh2DelR ->Fill(JetPt,DelR);
820 //_________________________________________________________________________________//
822 void AliAnalysisTaskJetProperties::FillJetShape(TList *jetlist){
823 //filling up the histograms
824 if(fDebug > 1) printf("AliAnalysisTaskJetProperties::FillJetShape() \n");
825 Float_t JetEta; Float_t JetPhi; Float_t JetPt;
826 AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(0));//Leading jet
831 fh1PtLeadingJet->Fill(JetPt);
832 Float_t NchSumA[50] = {0.};
833 Float_t PtSumA[50] = {0.};
834 Float_t delRPtSum80pc = 0;
835 Float_t delRNtrkSum80pc = 0;
836 Float_t PtSumDiffShape[10] = {0.0};
837 Float_t PtSumIntShape[10] = {0.0};
840 Int_t nJT = GetListOfJetTracks(fTrackList,jet);
841 Int_t nJetTracks = 0;
842 if(nJT>=0) nJetTracks = fTrackList->GetEntries();
844 Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
845 if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
847 fh2NtracksLeadingJet->Fill(JetPt,nJetTracks);
848 fProNtracksLeadingJet->Fill(JetPt,nJetTracks);
849 Int_t *index = new Int_t [nJetTracks];//dynamic array containing index
850 Float_t *delRA = new Float_t [nJetTracks];//dynamic array containing delR
851 Float_t *delEtaA = new Float_t [nJetTracks];//dynamic array containing delEta
852 Float_t *delPhiA = new Float_t [nJetTracks];//dynamic array containing delPhi
853 Float_t *trackPtA = new Float_t [nJetTracks];//dynamic array containing pt-track
854 Float_t *trackEtaA = new Float_t [nJetTracks];//dynamic array containing eta-track
855 Float_t *trackPhiA = new Float_t [nJetTracks];//dynamic array containing phi-track
856 for(Int_t ii=0; ii<nJetTracks; ii++){
866 for (Int_t j =0; j< nJetTracks; j++){
867 if(fTrackType==kTrackUndef)continue;
869 Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
870 Float_t DelEta=-99.0; Float_t DelPhi=-99.0; Float_t DelR=-99.0; Float_t AreaJ=-99.0;
872 if(fTrackType==kTrackAOD){
873 AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackList->At(j));
874 if(!trackaod)continue;
875 TrackEta = trackaod->Eta();
876 TrackPhi = trackaod->Phi();
877 TrackPt = trackaod->Pt();
879 else if(fTrackType==kTrackAODMC){
880 AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackList->At(j));
881 if(!trackmc)continue;
882 TrackEta = trackmc->Eta();
883 TrackPhi = trackmc->Phi();
884 TrackPt = trackmc->Pt();
886 else if(fTrackType==kTrackKine){
887 AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackList->At(j));
888 if(!trackkine)continue;
889 TrackEta = trackkine->Eta();
890 TrackPhi = trackkine->Phi();
891 TrackPt = trackkine->Pt();
894 DelEta = TMath::Abs(JetEta - TrackEta);
895 DelPhi = TMath::Abs(JetPhi - TrackPhi);
896 if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
897 DelR = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
898 AreaJ = TMath::Pi()*DelR*DelR;
900 fh2AreaCh ->Fill(JetPt,AreaJ);
901 fProAreaCh->Fill(JetPt,AreaJ);
905 trackPtA[j] = TrackPt;
906 trackEtaA[j] = TrackEta;
907 trackPhiA[j] = TrackPhi;
909 //calculating diff and integrated jet shapes
910 Float_t kDeltaR = 0.1;
911 Float_t RMin = kDeltaR/2.0;
912 Float_t RMax = kDeltaR/2.0;
914 for(Int_t ii1=0; ii1<kNbinsR;ii1++){
915 if((DelR > (tmpR-RMin)) && (DelR <=(tmpR+RMax)))PtSumDiffShape[ii1]+= TrackPt;
916 if(DelR>0.0 && DelR <=(tmpR+RMax))PtSumIntShape[ii1]+= TrackPt;
920 for(Int_t ibin=1; ibin<=50; ibin++){
921 Float_t xlow = 0.02*(ibin-1);
922 Float_t xup = 0.02*ibin;
923 if( xlow <= DelR && DelR < xup){
925 PtSumA[ibin-1]+= TrackPt;
931 //---------------------
932 Float_t tmp1R = 0.05;
933 for(Int_t jj1=0; jj1<kNbinsR;jj1++){
934 if(JetPt>20 && JetPt<=100){
935 fProDiffJetShape->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
936 fProIntJetShape ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
938 Float_t jetPtMin0=20.0; Float_t jetPtMax0=30.0;
939 for(Int_t k=0; k<5; k++){
940 if(k==0){jetPtMin0=20.0;jetPtMax0=30.0;}
941 if(k==1){jetPtMin0=30.0;jetPtMax0=40.0;}
942 if(k==2){jetPtMin0=40.0;jetPtMax0=60.0;}
943 if(k==3){jetPtMin0=60.0;jetPtMax0=80.0;}
944 if(k==4){jetPtMin0=80.0;jetPtMax0=100.0;}
945 if(JetPt>jetPtMin0 && JetPt<=jetPtMax0){
946 fProDiffJetShapeA[k]->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
947 fProIntJetShapeA[k] ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
952 //----------------------//
955 Bool_t iflagPtSum = kFALSE;
956 Bool_t iflagNtrkSum = kFALSE;
957 TMath::Sort(nJetTracks,delRA,index,0);
958 for(Int_t ii=0; ii<nJetTracks; ii++){
960 PtSum += trackPtA[index[ii]];
962 cout << index[ii] << "\t" <<
966 trackPt[ii]<< "\t" <<
967 trackEta[ii]<< "\t" <<
968 trackPhi[ii]<< "\t DelR " <<
969 delR[index[ii]] << endl;
972 if((Float_t)NtrkSum/(Float_t)nJetTracks > 0.79){
973 delRNtrkSum80pc = delRA[index[ii]];
974 iflagNtrkSum = kTRUE;
978 if(PtSum/JetPt >= 0.8000){
979 delRPtSum80pc = delRA[index[ii]];
991 fh2DelR80pcNch ->Fill(JetPt,delRNtrkSum80pc);
992 fProDelR80pcNch->Fill(JetPt,delRNtrkSum80pc);
993 fh2DelR80pcPt ->Fill(JetPt,delRPtSum80pc);
994 fProDelR80pcPt ->Fill(JetPt,delRPtSum80pc);
995 for(Int_t ibin=0; ibin<50; ibin++){
996 Float_t iR = 0.02*ibin + 0.01;
997 fh3PtDelRNchSum->Fill(JetPt,iR,NchSumA[ibin]);
998 fh3PtDelRPtSum ->Fill(JetPt,iR,PtSumA[ibin]);
999 Float_t jetPtMin=20.0; Float_t jetPtMax=30.0;
1000 for(Int_t k=0; k<5; k++){
1001 if(k==0){jetPtMin=20.0;jetPtMax=30.0;}
1002 if(k==1){jetPtMin=30.0;jetPtMax=40.0;}
1003 if(k==2){jetPtMin=40.0;jetPtMax=60.0;}
1004 if(k==3){jetPtMin=60.0;jetPtMax=80.0;}
1005 if(k==4){jetPtMin=80.0;jetPtMax=100.0;}
1006 if(JetPt>jetPtMin && JetPt<jetPtMax){
1007 fProDelRNchSum[k]->Fill(iR,NchSumA[ibin]);
1008 fProDelRPtSum[k]->Fill(iR,PtSumA[ibin]);
1014 //_________________________________________________________________________________//