1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets
4 // *******************************************
7 /**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
13 * Permission to use, copy, modify and distribute this software and its *
14 * documentation strictly for non-commercial purposes is hereby granted *
15 * without fee, provided that the above copyright notice appears in all *
16 * copies and that both the copyright notice and this permission notice *
17 * appear in the supporting documentation. The authors make no claims *
18 * about the suitability of this software for any purpose. It is *
19 * provided "as is" without express or implied warranty. *
20 **************************************************************************/
26 #include <TInterpreter.h>
35 #include <TLorentzVector.h>
36 #include <TClonesArray.h>
37 #include "TDatabasePDG.h"
39 #include "AliAnalysisTaskJetSpectrum2.h"
40 #include "AliAnalysisManager.h"
41 #include "AliJetFinder.h"
42 #include "AliJetHeader.h"
43 #include "AliJetReader.h"
44 #include "AliJetReaderHeader.h"
45 #include "AliUA1JetHeaderV1.h"
46 #include "AliESDEvent.h"
47 #include "AliAODEvent.h"
48 #include "AliAODHandler.h"
49 #include "AliAODTrack.h"
50 #include "AliAODJet.h"
51 #include "AliAODMCParticle.h"
52 #include "AliMCEventHandler.h"
53 #include "AliMCEvent.h"
55 #include "AliGenPythiaEventHeader.h"
56 #include "AliJetKineReaderHeader.h"
57 #include "AliGenCocktailEventHeader.h"
58 #include "AliInputEventHandler.h"
61 #include "AliAnalysisHelperJetTasks.h"
63 ClassImp(AliAnalysisTaskJetSpectrum2)
65 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
72 fUseAODJetInput(kFALSE),
73 fUseAODTrackInput(kFALSE),
74 fUseAODMCInput(kFALSE),
75 fUseGlobalSelection(kFALSE),
76 fUseExternalWeightOnly(kFALSE),
77 fLimitGenJetEta(kFALSE),
80 fTrackTypeRec(kTrackUndef),
81 fTrackTypeGen(kTrackUndef),
93 fh1SumPtTrackRec(0x0),
94 fh1SumPtTrackAreaRec(0x0),
96 fh1PtJetsLeadingRecIn(0x0),
97 fh1PtTracksRecIn(0x0),
98 fh1PtTracksLeadingRecIn(0x0),
99 fh1PtTracksGenIn(0x0),
101 fh2NRecTracksPt(0x0),
102 fh2JetsLeadingPhiEta(0x0),
103 fh2JetsLeadingPhiPt(0x0),
104 fh2TracksLeadingPhiEta(0x0),
105 fh2TracksLeadingPhiPt(0x0),
106 fh2TracksLeadingJetPhiPt(0x0),
109 for(int i = 0;i < kMaxStep*2;++i){
110 fhnJetContainer[i] = 0;
112 for(int i = 0;i < kMaxJets;++i){
113 fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
114 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
119 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
120 AliAnalysisTaskSE(name),
127 fUseAODJetInput(kFALSE),
128 fUseAODTrackInput(kFALSE),
129 fUseAODMCInput(kFALSE),
130 fUseGlobalSelection(kFALSE),
131 fUseExternalWeightOnly(kFALSE),
132 fLimitGenJetEta(kFALSE),
135 fTrackTypeRec(kTrackUndef),
136 fTrackTypeGen(kTrackUndef),
144 fh1PtHardTrials(0x0),
148 fh1SumPtTrackRec(0x0),
149 fh1SumPtTrackAreaRec(0x0),
151 fh1PtJetsLeadingRecIn(0x0),
152 fh1PtTracksRecIn(0x0),
153 fh1PtTracksLeadingRecIn(0x0),
154 fh1PtTracksGenIn(0x0),
156 fh2NRecTracksPt(0x0),
157 fh2JetsLeadingPhiEta(0x0),
158 fh2JetsLeadingPhiPt(0x0),
159 fh2TracksLeadingPhiEta(0x0),
160 fh2TracksLeadingPhiPt(0x0),
161 fh2TracksLeadingJetPhiPt(0x0),
165 for(int i = 0;i < kMaxStep*2;++i){
166 fhnJetContainer[i] = 0;
168 for(int i = 0;i < kMaxJets;++i){
169 fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
170 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
172 DefineOutput(1, TList::Class());
177 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
180 // Implemented Notify() to read the cross sections
181 // and number of trials from pyxsec.root
184 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
185 Float_t xsection = 0;
190 TFile *curfile = tree->GetCurrentFile();
192 Error("Notify","No current file");
195 if(!fh1Xsec||!fh1Trials){
196 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
199 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
200 fh1Xsec->Fill("<#sigma>",xsection);
201 // construct a poor man average trials
202 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
203 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
208 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
212 // Create the output container
222 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
225 if(!fHistList)fHistList = new TList();
226 fHistList->SetOwner(kTRUE);
227 Bool_t oldStatus = TH1::AddDirectoryStatus();
228 TH1::AddDirectory(kFALSE);
233 const Int_t nBinPt = 100;
234 Double_t binLimitsPt[nBinPt+1];
235 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
237 binLimitsPt[iPt] = 0.0;
240 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
244 const Int_t nBinPhi = 90;
245 Double_t binLimitsPhi[nBinPhi+1];
246 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
248 binLimitsPhi[iPhi] = -1.*TMath::Pi();
251 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
257 const Int_t nBinEta = 40;
258 Double_t binLimitsEta[nBinEta+1];
259 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
261 binLimitsEta[iEta] = -2.0;
264 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
268 const Int_t nBinFrag = 25;
271 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
272 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
274 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
275 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
277 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
278 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
279 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
281 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
282 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
284 fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
285 fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
286 fh1SumPtTrackAreaRec = new TH1F("fh1SumPtTrackAreaRec","Sum Rec track P_T #eta <0.9 / 1.8 * 2 * 0.4*0.4;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
288 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
289 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
290 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
291 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
292 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
294 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
295 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
298 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
299 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
300 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
301 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
302 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
303 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
304 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
305 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
307 for(int ij = 0;ij < kMaxJets;++ij){
308 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
309 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
311 fh2PhiPt[ij] = new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
312 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
314 fh2PhiEta[ij] = new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
315 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
318 fh2FragRec[ij] = new TH2F(Form("fh2FragRec_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
319 nBinFrag,0.,1.,nBinPt,binLimitsPt);
320 fh2FragLnRec[ij] = new TH2F(Form("fh2FragLnRec_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
321 nBinFrag,0.,10.,nBinPt,binLimitsPt);
323 fh2FragGen[ij] = new TH2F(Form("fh2FragGen_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
324 nBinFrag,0.,1.,nBinPt,binLimitsPt);
325 fh2FragLnGen[ij] = new TH2F(Form("fh2FragLnGen_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
326 nBinFrag,0.,10.,nBinPt,binLimitsPt);
330 const Int_t saveLevel = 3; // large save level more histos
332 fHistList->Add(fh1Xsec);
333 fHistList->Add(fh1Trials);
334 fHistList->Add(fh1PtHard);
335 fHistList->Add(fh1PtHardNoW);
336 fHistList->Add(fh1PtHardTrials);
337 if(fBranchGen.Length()>0){
338 fHistList->Add(fh1NGenJets);
339 fHistList->Add(fh1PtTracksGenIn);
341 fHistList->Add(fh1PtJetsRecIn);
342 fHistList->Add(fh1PtJetsLeadingRecIn);
343 fHistList->Add(fh1PtTracksRecIn);
344 fHistList->Add(fh1PtTracksLeadingRecIn);
345 fHistList->Add(fh1NRecJets);
346 fHistList->Add(fh1PtTrackRec);
347 fHistList->Add(fh1SumPtTrackRec);
348 fHistList->Add(fh1SumPtTrackAreaRec);
349 fHistList->Add(fh2NRecJetsPt);
350 fHistList->Add(fh2NRecTracksPt);
351 fHistList->Add(fh2JetsLeadingPhiEta );
352 fHistList->Add(fh2JetsLeadingPhiPt );
353 fHistList->Add(fh2TracksLeadingPhiEta);
354 fHistList->Add(fh2TracksLeadingPhiPt);
355 for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
356 for(int ij = 0;ij<kMaxJets;++ij){
357 fHistList->Add( fh1PtRecIn[ij]);
358 if(fBranchGen.Length()>0){
359 fHistList->Add( fh1PtGenIn[ij]);
360 fHistList->Add( fh2FragGen[ij]);
361 fHistList->Add( fh2FragLnGen[ij]);
363 fHistList->Add( fh2PhiPt[ij]);
364 fHistList->Add( fh2PhiEta[ij]);
365 fHistList->Add( fh2FragRec[ij]);
366 fHistList->Add( fh2FragLnRec[ij]);
368 fHistList->Add(fhnCorrelation);
371 // =========== Switch on Sumw2 for all histos ===========
372 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
373 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
378 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
381 TH1::AddDirectory(oldStatus);
384 void AliAnalysisTaskJetSpectrum2::Init()
390 Printf(">>> AnalysisTaskJetSpectrum2::Init() debug level %d\n",fDebug);
391 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
395 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
398 if(!AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
399 // no selection by the service task, we continue
400 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
401 PostData(1, fHistList);
406 // Execute analysis for current event
408 AliESDEvent *fESD = 0;
410 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
412 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
418 // assume that the AOD is in the general output...
421 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
425 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
432 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
435 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
438 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
442 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
443 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
445 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
449 // ==== General variables needed
452 // We use statice array, not to fragment the memory
453 AliAODJet genJets[kMaxJets];
455 AliAODJet recJets[kMaxJets];
457 ///////////////////////////
462 Double_t nTrials = 1; // Trials for MC trigger
464 if(fUseExternalWeightOnly){
465 eventW = fExternalWeight;
468 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
469 // if(fDebug>0)aodH->SetFillAOD(kFALSE);
470 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
471 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
472 // this is the part we only use when we have MC information
473 AliMCEvent* mcEvent = MCEvent();
474 // AliStack *pStack = 0;
476 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
479 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
482 nTrials = pythiaGenHeader->Trials();
483 ptHard = pythiaGenHeader->GetPtHard();
484 int iProcessType = pythiaGenHeader->ProcessType();
486 // 12 f+barf -> f+barf
491 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
492 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
494 // fetch the pythia generated jets only to be used here
495 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
496 AliAODJet pythiaGenJets[kMaxJets];
497 for(int ip = 0;ip < nPythiaGenJets;++ip){
498 if(iCount>=kMaxJets)continue;
500 pythiaGenHeader->TriggerJet(ip,p);
501 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
505 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
506 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
510 if(fBranchGen.Length()==0){
511 // if we have MC particles and we do not read from the aod branch
512 // use the pythia jets
513 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
518 if(fBranchGen.Length()==0)nGenJets = iCount;
519 }// (fAnalysisType&kMCESD)==kMCESD)
522 // we simply fetch the tracks/mc particles as a list of AliVParticles
530 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
531 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
532 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
533 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
536 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
537 fh1PtHard->Fill(ptHard,eventW);
538 fh1PtHardNoW->Fill(ptHard,1);
539 fh1PtHardTrials->Fill(ptHard,nTrials);
541 // If we set a second branch for the input jets fetch this
542 if(fBranchGen.Length()>0){
543 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
546 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
547 if(iCount>=kMaxJets)continue;
548 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
552 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
553 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
556 genJets[iCount] = *tmp;
562 if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
563 if(fDebug>2)fAOD->Print();
567 fh1NGenJets->Fill(nGenJets);
568 // We do not want to exceed the maximum jet number
569 nGenJets = TMath::Min(nGenJets,kMaxJets);
571 // Fetch the reconstructed jets...
573 nRecJets = aodRecJets->GetEntries();
575 nRecJets = aodRecJets->GetEntries();
576 fh1NRecJets->Fill(nRecJets);
578 // Do something with the all rec jets
579 Int_t nRecOver = nRecJets;
581 // check that the jets are sorted
582 Float_t ptOld = 999999;
583 for(int ir = 0;ir < nRecJets;ir++){
584 AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
585 Float_t tmpPt = tmp->Pt();
587 Printf("%s:%d Jets Not Sorted!! %d:%.3E %d%.3E",(char*)__FILE__,__LINE__,ir,tmpPt,ir-1,ptOld);
594 TIterator *recIter = aodRecJets->MakeIterator();
595 AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
596 Float_t pt = tmpRec->Pt();
598 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
599 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
600 while(pt<ptCut&&tmpRec){
602 tmpRec = (AliAODJet*)(recIter->Next());
607 if(nRecOver<=0)break;
608 fh2NRecJetsPt->Fill(ptCut,nRecOver);
612 AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
613 Float_t phi = leading->Phi();
614 if(phi<0)phi+=TMath::Pi()*2.;
615 Float_t eta = leading->Eta();
617 while((tmpRec = (AliAODJet*)(recIter->Next()))){
618 Float_t tmpPt = tmpRec->Pt();
619 fh1PtJetsRecIn->Fill(tmpPt);
621 fh1PtJetsLeadingRecIn->Fill(tmpPt);
625 Float_t tmpPhi = tmpRec->Phi();
627 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
628 Float_t dPhi = phi - tmpRec->Phi();
629 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
630 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
631 Float_t dEta = eta - tmpRec->Eta();
632 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
633 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
638 Int_t nTrackOver = recParticles.GetSize();
639 // do the same for tracks and jets
641 TIterator *recIter = recParticles.MakeIterator();
642 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
643 Float_t pt = tmpRec->Pt();
644 // Printf("Leading track p_t %3.3E",pt);
645 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
646 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
647 while(pt<ptCut&&tmpRec){
649 tmpRec = (AliVParticle*)(recIter->Next());
654 if(nTrackOver<=0)break;
655 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
659 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
660 Float_t phi = leading->Phi();
661 if(phi<0)phi+=TMath::Pi()*2.;
662 Float_t eta = leading->Eta();
664 while((tmpRec = (AliVParticle*)(recIter->Next()))){
665 Float_t tmpPt = tmpRec->Pt();
666 fh1PtTracksRecIn->Fill(tmpPt);
668 fh1PtTracksLeadingRecIn->Fill(tmpPt);
672 Float_t tmpPhi = tmpRec->Phi();
674 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
675 Float_t dPhi = phi - tmpRec->Phi();
676 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
677 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
678 Float_t dEta = eta - tmpRec->Eta();
679 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
680 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
685 if(genParticles.GetSize()){
686 TIterator *genIter = genParticles.MakeIterator();
687 AliVParticle *tmpGen = 0;
688 while((tmpGen = (AliVParticle*)(genIter->Next()))){
689 if(TMath::Abs(tmpGen->Eta())<0.9){
690 Float_t tmpPt = tmpGen->Pt();
691 fh1PtTracksGenIn->Fill(tmpPt);
697 fh1NRecJets->Fill(nRecJets);
698 nRecJets = TMath::Min(nRecJets,kMaxJets);
700 for(int ir = 0;ir < nRecJets;++ir){
701 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
707 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
709 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
710 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
713 for(int i = 0;i<kMaxJets;++i){
714 iGenIndex[i] = iRecIndex[i] = -1;
717 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
718 iGenIndex,iRecIndex,fDebug);
719 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
722 for(int i = 0;i<kMaxJets;++i){
723 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
724 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
731 Double_t container[6];
733 // loop over generated jets
735 // radius; tmp, get from the jet header itself
736 // at some point, how todeal woht FastJet on MC?
737 Float_t radiusGen =0.4;
738 Float_t radiusRec =0.4;
740 for(int ig = 0;ig < nGenJets;++ig){
741 Double_t ptGen = genJets[ig].Pt();
742 Double_t phiGen = genJets[ig].Phi();
743 if(phiGen<0)phiGen+=TMath::Pi()*2.;
744 Double_t etaGen = genJets[ig].Eta();
746 container[3] = ptGen;
747 container[4] = etaGen;
748 container[5] = phiGen;
749 fhnJetContainer[kStep0]->Fill(&container[3],eventW);
750 Int_t ir = iRecIndex[ig];
752 if(TMath::Abs(etaGen)<fRecEtaWindow){
753 fhnJetContainer[kStep1]->Fill(&container[3],eventW);
754 fh1PtGenIn[ig]->Fill(ptGen,eventW);
755 // fill the fragmentation function
756 for(int it = 0;it<genParticles.GetEntries();++it){
757 AliVParticle *part = (AliVParticle*)genParticles.At(it);
758 if(genJets[ig].DeltaR(part)<radiusGen){
759 Float_t z = part->Pt()/ptGen;
760 Float_t lnz = -1.*TMath::Log(z);
761 fh2FragGen[ig]->Fill(z,ptGen,eventW);
762 fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
765 if(ir>=0&&ir<nRecJets){
766 fhnJetContainer[kStep3]->Fill(&container[3],eventW);
770 if(ir>=0&&ir<nRecJets){
771 fhnJetContainer[kStep2]->Fill(&container[3],eventW);
772 Double_t etaRec = recJets[ir].Eta();
773 if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
775 }// loop over generated jets
779 for(int it = 0;it<recParticles.GetEntries();++it){
780 AliVParticle *part = (AliVParticle*)recParticles.At(it);
781 // fill sum pt and P_t of all paritles
782 if(TMath::Abs(part->Eta())<0.9){
783 Float_t pt = part->Pt();
784 fh1PtTrackRec->Fill(pt,eventW);
788 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
789 fh1SumPtTrackRec->Fill(sumPt,eventW);
792 // loop over reconstructed jets
793 for(int ir = 0;ir < nRecJets;++ir){
794 Double_t etaRec = recJets[ir].Eta();
795 Double_t ptRec = recJets[ir].Pt();
796 Double_t phiRec = recJets[ir].Phi();
797 if(phiRec<0)phiRec+=TMath::Pi()*2.;
798 container[0] = ptRec;
799 container[1] = etaRec;
800 container[2] = phiRec;
802 if(ptRec>10.&&fDebug>0){
803 // need to cast to int, otherwise the printf overwrites
804 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
805 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
806 // aodH->SetFillAOD(kTRUE);
807 fAOD->GetHeader()->Print();
808 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
809 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
810 AliAODTrack *tr = fAOD->GetTrack(it);
811 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
815 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
823 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
824 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
825 if(TMath::Abs(etaRec)<fRecEtaWindow){
826 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
827 fh1PtRecIn[ir]->Fill(ptRec,eventW);
828 // fill the fragmentation function
829 for(int it = 0;it<recParticles.GetEntries();++it){
830 AliVParticle *part = (AliVParticle*)recParticles.At(it);
831 Float_t eta = part->Eta();
832 if(TMath::Abs(eta)<0.9){
833 Float_t phi = part->Phi();
834 if(phi<0)phi+=TMath::Pi()*2.;
835 Float_t dPhi = phi - phiRec;
836 Float_t dEta = eta - etaRec;
837 if(dPhi<(-1.*TMath::Pi()))phiRec+=TMath::Pi()*2.;
838 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
839 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
841 if(recJets[ir].DeltaR(part)<radiusRec){
842 Float_t z = part->Pt()/ptRec;
843 Float_t lnz = -1.*TMath::Log(z);
844 fh2FragRec[ir]->Fill(z,ptRec,eventW);
845 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
852 Int_t ig = iGenIndex[ir];
853 if(ig>=0 && ig<nGenJets){
854 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
855 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
856 Double_t ptGen = genJets[ig].Pt();
857 Double_t phiGen = genJets[ig].Phi();
858 if(phiGen<0)phiGen+=TMath::Pi()*2.;
859 Double_t etaGen = genJets[ig].Eta();
861 container[3] = ptGen;
862 container[4] = etaGen;
863 container[5] = phiGen;
866 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
869 if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
870 if(TMath::Abs(etaRec)<fRecEtaWindow){
871 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
872 fhnCorrelation->Fill(container);
873 }// if etarec in window
876 ////////////////////////////////////////////////////
880 }// loop over reconstructed jets
883 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
884 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
885 PostData(1, fHistList);
889 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
891 // Create the particle container for the correction framework manager and
894 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
895 const Double_t kPtmin = 5.0, kPtmax = 105.; // we do not want to have empty bins at the beginning...
896 const Double_t kEtamin = -3.0, kEtamax = 3.0;
897 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
899 // can we neglect migration in eta and phi?
900 // phi should be no problem since we cover full phi and are phi symmetric
901 // eta migration is more difficult due to needed acceptance correction
902 // in limited eta range
905 //arrays for the number of bins in each dimension
907 iBin[0] = 100; //bins in pt
908 iBin[1] = 1; //bins in eta
909 iBin[2] = 1; // bins in phi
911 //arrays for lower bounds :
912 Double_t* binEdges[kNvar];
913 for(Int_t ivar = 0; ivar < kNvar; ivar++)
914 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
916 //values for bin lower bounds
917 // for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
918 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
919 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
920 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
923 for(int i = 0;i < kMaxStep*2;++i){
924 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
925 for (int k=0; k<kNvar; k++) {
926 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
929 //create correlation matrix for unfolding
930 Int_t thnDim[2*kNvar];
931 for (int k=0; k<kNvar; k++) {
932 //first half : reconstructed
935 thnDim[k+kNvar] = iBin[k];
938 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
939 for (int k=0; k<kNvar; k++) {
940 fhnCorrelation->SetBinEdges(k,binEdges[k]);
941 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
943 fhnCorrelation->Sumw2();
945 // Add a histogram for Fake jets
946 // thnDim[3] = AliPID::kSPECIES;
947 // fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
948 // for(Int_t idim = 0; idim < kNvar; idim++)
949 // fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
950 for(Int_t ivar = 0; ivar < kNvar; ivar++)
951 delete [] binEdges[ivar];
955 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
957 // Terminate analysis
959 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
963 Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
965 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
969 AliAODEvent *aod = 0;
970 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
971 else aod = AODEvent();
975 for(int it = 0;it < aod->GetNumberOfTracks();++it){
976 AliAODTrack *tr = aod->GetTrack(it);
977 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
978 if(TMath::Abs(tr->Eta())>0.9)continue;
981 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
982 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
985 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
987 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
997 else if (type == kTrackKineAll||type == kTrackKineCharged){
998 AliMCEvent* mcEvent = MCEvent();
999 if(!mcEvent)return iCount;
1000 // we want to have alivpartilces so use get track
1001 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1002 if(!mcEvent->IsPhysicalPrimary(it))continue;
1003 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1004 if(type == kTrackKineAll){
1008 else if(type == kTrackKineCharged){
1009 if(part->Particle()->GetPDG()->Charge()==0)continue;
1015 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1016 AliAODEvent *aod = 0;
1017 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1018 else aod = AODEvent();
1019 if(!aod)return iCount;
1020 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1021 if(!tca)return iCount;
1022 for(int it = 0;it < tca->GetEntriesFast();++it){
1023 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1024 if(!part->IsPhysicalPrimary())continue;
1025 if(type == kTrackAODMCAll){
1029 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1030 if(part->Charge()==0)continue;
1031 if(kTrackAODMCCharged){
1035 if(TMath::Abs(part->Eta())>0.9)continue;