]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalDiJetBase.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
from marta
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalDiJetBase.cxx
CommitLineData
6e8b6371 1//
2// Dijet base analysis task.
3//
4// Author: M.Verweij
5
6#include <TClonesArray.h>
7#include <TH1F.h>
8#include <TH2F.h>
9#include <TH3F.h>
10#include <THnSparse.h>
11#include <TList.h>
12#include <TLorentzVector.h>
13#include <TProfile.h>
14#include <TChain.h>
15#include <TSystem.h>
16#include <TFile.h>
17#include <TKey.h>
18
19#include "AliVCluster.h"
20#include "AliVTrack.h"
21#include "AliEmcalJet.h"
22#include "AliRhoParameter.h"
23#include "AliLog.h"
24#include "AliAnalysisUtils.h"
25#include "AliEmcalParticle.h"
26#include "AliMCEvent.h"
27#include "AliGenPythiaEventHeader.h"
28#include "AliAODMCHeader.h"
29#include "AliMCEvent.h"
30#include "AliAnalysisManager.h"
31#include "AliJetContainer.h"
32
33#include "AliAnalysisTaskEmcalDiJetBase.h"
34
35ClassImp(AliAnalysisTaskEmcalDiJetBase)
36
37//________________________________________________________________________
38AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase() :
39 AliAnalysisTaskEmcalJetDev("AliAnalysisTaskEmcalDiJetBase", kTRUE),
40 fDebug(kFALSE),
c66364a6 41 fJetCorrelationType(kCorrelateAll),
6e8b6371 42 fJetFullChargedMatchingType(kFraction),
43 fTriggerClass(""),
44 fContainerCharged(1),
45 fContainerFull(0),
46 fContainerChargedMC(3),
47 fContainerFullMC(2),
48 fRhoType(0),
49 fRhoChVal(0),
50 fRhoFullVal(0),
51 fDoChargedCharged(kTRUE),
52 fDoFullCharged(kTRUE),
bb84b374 53 fDoFullFull(kFALSE),
6e8b6371 54 fUseAnaUtils(kTRUE),
55 fAnalysisUtils(0),
56 fPtMinTriggerJet(10.),
57 fMinFractionShared(0.5),
58 fMatchingDone(kFALSE),
59 faFullFracIndex(0),
60 faFullFracIndexMC(0),
61 fIsPythiaPtHard(0),
62 fPythiaHeader(0),
63 fPtHard(0),
64 fPtHardBin(0),
65 fNTrials(0),
66 fhNEvents(0),
67 fHistTrials(0),
68 fHistTrialsSelEvents(0),
69 fHistXsection(0),
70 fHistEvents(0)
71{
72 // Default constructor.
73
74
75 SetMakeGeneralHistograms(kTRUE);
76
77}
78
79//________________________________________________________________________
80AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase(const char *name) :
81 AliAnalysisTaskEmcalJetDev(name, kTRUE),
82 fDebug(kFALSE),
c66364a6 83 fJetCorrelationType(kCorrelateAll),
6e8b6371 84 fJetFullChargedMatchingType(kFraction),
85 fTriggerClass(""),
86 fContainerCharged(1),
87 fContainerFull(0),
88 fContainerChargedMC(3),
89 fContainerFullMC(2),
90 fRhoType(0),
91 fRhoChVal(0),
92 fRhoFullVal(0),
93 fDoChargedCharged(kTRUE),
94 fDoFullCharged(kTRUE),
bb84b374 95 fDoFullFull(kFALSE),
6e8b6371 96 fUseAnaUtils(kTRUE),
97 fAnalysisUtils(0),
98 fPtMinTriggerJet(10.),
99 fMinFractionShared(0.5),
100 fMatchingDone(kFALSE),
101 faFullFracIndex(0),
102 faFullFracIndexMC(0),
103 fIsPythiaPtHard(0),
104 fPythiaHeader(0),
105 fPtHard(0),
106 fPtHardBin(0),
107 fNTrials(0),
108 fhNEvents(0),
109 fHistTrials(0),
110 fHistTrialsSelEvents(0),
111 fHistXsection(0),
112 fHistEvents(0)
113{
114 // Standard constructor.
115
116 SetMakeGeneralHistograms(kTRUE);
117
118}
119
120//________________________________________________________________________
121AliAnalysisTaskEmcalDiJetBase::~AliAnalysisTaskEmcalDiJetBase()
122{
123 // Destructor.
124}
125
126//----------------------------------------------------------------------------------------------
127void AliAnalysisTaskEmcalDiJetBase::InitOnce() {
128 //
129 // only initialize once
130 //
131
132 // Initialize analysis util class for vertex selection
133 if(fUseAnaUtils) {
134 fAnalysisUtils = new AliAnalysisUtils();
135 fAnalysisUtils->SetMinVtxContr(2);
136 fAnalysisUtils->SetMaxVtxZ(10.);
137 }
138}
139
140//________________________________________________________________________
141Bool_t AliAnalysisTaskEmcalDiJetBase::SelectEvent() {
142 //
143 // Decide if event should be selected for analysis
144 //
145
146 fhNEvents->Fill(0.5);
147
148 if(!fAnalysisUtils && fUseAnaUtils)
149 InitOnce();
150
151
152 if(fAnalysisUtils) {
153 if(!fAnalysisUtils->IsVertexSelected2013pA(InputEvent()))
154 return kFALSE;
155
156 fhNEvents->Fill(2.5);
157
158 if(fAnalysisUtils->IsPileUpEvent(InputEvent()))
159 return kFALSE;
160 }
161 else{
162 if(fUseAnaUtils)
163 AliError("fAnalysisUtils not initialized. Call AliAnalysisTaskEmcalJetTriggerQA::InitOnce()");
164 }
165
166 fhNEvents->Fill(3.5);
167
168 if(!fTriggerClass.IsNull()) {
169 //Check if requested trigger was fired
170 TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
171 if(!firedTrigClass.Contains(fTriggerClass))
172 return kFALSE;
173 }
174
175 fhNEvents->Fill(1.5);
176
177 return kTRUE;
178
179}
180
181//________________________________________________________________________
182void AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects()
183{
184 // Create user output.
185
186 InitOnce();
187
188 AliAnalysisTaskEmcalJetDev::UserCreateOutputObjects();
189
190 Bool_t oldStatus = TH1::AddDirectoryStatus();
191 TH1::AddDirectory(kFALSE);
192
193 fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
194 fOutput->Add(fhNEvents);
195
196
197 //Pythia info
198
199 fHistTrials = new TH1F("fHistTrials", "fHistTrials", 12, 0, 12);
200 fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
201 fHistTrials->GetYaxis()->SetTitle("trials");
202 fOutput->Add(fHistTrials);
203
204 fHistTrialsSelEvents = new TH1F("fHistTrialsSelEvents", "fHistTrialsSelEvents", 12, 0, 12);
205 fHistTrialsSelEvents->GetXaxis()->SetTitle("p_{T} hard bin");
206 fHistTrialsSelEvents->GetYaxis()->SetTitle("trials");
207 fOutput->Add(fHistTrialsSelEvents);
208
209 fHistXsection = new TProfile("fHistXsection", "fHistXsection", 12, 0, 12);
210 fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
211 fHistXsection->GetYaxis()->SetTitle("xsection");
212 fOutput->Add(fHistXsection);
213
214 fHistEvents = new TH1F("fHistEvents", "fHistEvents", 12, 0, 12);
215 fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
216 fHistEvents->GetYaxis()->SetTitle("total events");
217 fOutput->Add(fHistEvents);
218
219
220 TH1::AddDirectory(oldStatus);
221
222 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
223}
224
225//________________________________________________________________________
226Bool_t AliAnalysisTaskEmcalDiJetBase::IsSameJet(const Int_t ijt, const Int_t ija, const Int_t type, const Bool_t isMC) {
227 //check if two jets are the same one
228
229 Bool_t bSame = kFALSE;
230
231 if(type<2 && ijt==ija)
232 bSame = kTRUE;
233 if(type==2) {
234
235 if(fJetFullChargedMatchingType==kFraction) {
236 if(isMC && faFullFracIndexMC.At(ijt)==ija)
237 bSame = kTRUE;
238 else if(!isMC && faFullFracIndex.At(ijt)==ija)
239 bSame = kTRUE;
240 }
241 else if(fJetFullChargedMatchingType==kGeo) {
242 Int_t contFull = fContainerFull;
243 Int_t contChar = fContainerCharged;
244 if(isMC) {
245 contFull = fContainerFullMC;
246 contChar = fContainerChargedMC;
247 }
248 AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, contFull));
249 AliEmcalJet *chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, contChar));
250 AliEmcalJet *matchJet = fullJet->ClosestJet();
251 if(chJet==matchJet)
252 bSame = kTRUE;
253 }
254 else if(fJetFullChargedMatchingType==kNoMatching) {
255 return kFALSE;
256 }
257 else{
258 AliWarning(Form("%s: matching type unknown", GetName()));
259 }
260
261
262 }
263
264 return bSame;
265 }
266
267
268//________________________________________________________________________
269Double_t AliAnalysisTaskEmcalDiJetBase::GetJetPt(const AliEmcalJet *jet, const Int_t type) {
270
271 if(!jet) return -99;
272
273 if(type==0)
274 return jet->Pt() - fRhoFullVal * jet->Area();
275 else if(type==1)
276 return jet->Pt() - fRhoChVal * jet->Area();
277 else
278 return jet->Pt();
279
280}
281
282//________________________________________________________________________
283Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const AliVParticle *trk, const AliEmcalJet *jet) const
284{
285 // Get Z of constituent trk
286
287 return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
288}
289
290//________________________________________________________________________
291Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
292{
293 //
294 // Get the z of a constituent inside of a jet
295 //
296
297 Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
298
299 if(pJetSq>0.)
300 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
301 else {
302 AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
303 return 0;
304 }
305}
306
307//________________________________________________________________________
308Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const {
309 //
310 // Helper function to calculate the distance between two jets
311 //
312
313 Double_t dPhi = jet1->Phi() - jet2->Phi();
314 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
315 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
316 Double_t dEta = jet1->Eta() - jet2->Eta();
317 Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
318 return dR;
319}
320
321//________________________________________________________________________
322Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2) {
323 //
324 // Calculate azimuthal angle between the axises of the jets
325 //
326
327 return GetDeltaPhi(jet1->Phi(),jet2->Phi());
328
329}
330
331//________________________________________________________________________
332Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaPhi(Double_t phi1,Double_t phi2) {
333 //
334 // Calculate azimuthal angle between the axises of the jets
335 //
336
337 Double_t dPhi = phi1-phi2;
338
339 if(dPhi > 2*TMath::Pi()) dPhi -= TMath::TwoPi();
340 if(dPhi <-2*TMath::Pi()) dPhi += TMath::TwoPi();
341 if(dPhi <-0.5*TMath::Pi()) dPhi += TMath::TwoPi();
342 if(dPhi > 1.5*TMath::Pi()) dPhi -= TMath::TwoPi();
343
344 return dPhi;
345}
346
347//________________________________________________________________________
348Double_t AliAnalysisTaskEmcalDiJetBase::GetFractionSharedPt(const AliEmcalJet *jetFull, const AliEmcalJet *jetCharged) const
349{
350 //
351 // Get fraction of shared pT between matched full and charged jet
352 // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
353 //
354
355 Double_t fraction = 0.;
356 Double_t jetPtCh = jetCharged->Pt();
357
358 if(jetPtCh>0) {
359
360 if(fDebug>10) AliInfo(Form("%s: nConstituents: %d, ch: %d chne: %d ne: %d",GetName(),jetFull->GetNumberOfConstituents(),jetCharged->GetNumberOfTracks(),jetFull->GetNumberOfTracks(),jetFull->GetNumberOfClusters()));
361
362 Double_t sumPt = 0.;
363 AliVParticle *vpf = 0x0;
364 Int_t iFound = 0;
365 for(Int_t icc=0; icc<jetCharged->GetNumberOfTracks(); icc++) {
366 Int_t idx = (Int_t)jetCharged->TrackAt(icc);
367 iFound = 0;
368 for(Int_t icf=0; icf<jetFull->GetNumberOfTracks(); icf++) {
369 if(idx == jetFull->TrackAt(icf) && iFound==0 ) {
370 iFound=1;
371 vpf = static_cast<AliVParticle*>(jetFull->TrackAt(icf, fTracks));
372 if(!vpf) continue;
373 if(vpf->Charge()!=0)
374 sumPt += vpf->Pt();
375 continue;
376 }
377 }
378 }
379
380 fraction = sumPt/jetPtCh;
381 }
382
383 if(fDebug>10) AliInfo(Form("%s: charged shared fraction: %.2f",GetName(),fraction));
384
385 return fraction;
386
387}
388
389//________________________________________________________________________
390void AliAnalysisTaskEmcalDiJetBase::MatchJetsGeo(Int_t cFull, Int_t cCharged,
391 Int_t iDebug, Float_t maxDist, Int_t type) {
392
393 //
394 // Match the full jets to the corresponding charged jets
395 // Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
396 // type: 0 = use acceptance cuts of container 1 = allow 0.1 one more for cCharged(MC) in eta 2 = allow 0.1 more in eta and phi for cCharged(MC)
397
398 const int kMode = 3;
399
400 const Int_t nChJets = GetNJets(cCharged);
401 const Int_t nFullJets = GetNJets(cFull);
402
403 if(nFullJets==0 || nChJets==0) return;
404
405 TArrayI faChNeMatchIndex;
406 faChNeMatchIndex.Set(nChJets+1);
407
408 TArrayI faChMatchIndex;
409 faChMatchIndex.Set(nFullJets+1);
410
411 static TArrayS iFlag(nChJets*nFullJets);
412 if(iFlag.GetSize()<(nChJets*nFullJets)){
413 iFlag.Set(nChJets*nFullJets+1);
414 }
415 iFlag.Reset(0);
416
417 AliJetContainer *contCh = GetJetContainer(cCharged);
418
419 // find the closest distance to the full jet
420 for(int ifu = 0;ifu<nFullJets;ifu++){
421
422 AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, cFull));
423 if(!fullJet) continue;
424
425 Float_t dist = maxDist;
426
427 for(int ich = 0;ich<nChJets;ich++){
428 AliEmcalJet *chJet = 0x0;
429 if(type==0)
430 chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
431 else {
432 chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
433 if(!chJet) continue;;
434 if(type>0) {
435 if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) || chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
436 continue;
437 if(type==2) {
438 if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) || chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
439 continue;
440 }
441 }
442 }
443 if(!chJet)
444 continue;
445
446 Double_t frac = GetFractionSharedPt(fullJet,chJet);
447 Double_t dR = GetDeltaR(fullJet,chJet);
448 if(dR<dist){
449 faChMatchIndex[ifu]=ich;
450 dist = dR;
451 }
452 if(iDebug>10) Printf("Distance (%d)--(%d) %3.3f frac:%.2f",ifu,ich,dR,frac);
453 }
454 if(faChMatchIndex[ifu]>=0) iFlag[ifu*nChJets+faChMatchIndex[ifu]]+=1;//ich closest to ifu
455 if(iDebug>10) Printf("Full Distance (%d)--(%d) %3.3f flag[%d] = %d",ifu,faChMatchIndex[ifu],dist,ifu*nChJets+faChMatchIndex[ifu],iFlag[ifu*nChJets+faChMatchIndex[ifu]]);
456
457 // reset...
458 faChMatchIndex[ifu]=-1;
459
460
461 }//full jet loop
462
463
464 // other way around
465 for(int ich = 0;ich<nChJets;ich++){
466 AliEmcalJet *chJet = 0x0;
467 if(type==0)
468 chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
469 else {
470 chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
471 if(!chJet) continue;;
472 if(type>0) {
473 if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) || chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
474 continue;
475 if(type==2) {
476 if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) || chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
477 continue;
478 }
479 }
480 }
481 if(!chJet)
482 continue;
483
484 Float_t dist = maxDist;
485 for(int ifu = 0;ifu<nFullJets;ifu++){
486 AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, cFull));
487 if(!fullJet)
488 continue;
489
490 Double_t dR = GetDeltaR(fullJet,chJet);
491 if(dR<dist){
492 faChNeMatchIndex[ich]=ifu;
493 dist = dR;
494 }
495 }
496 if(faChNeMatchIndex[ich]>=0) iFlag[faChNeMatchIndex[ich]*nChJets+ich]+=2;//ifu closest to ich
497 if(iDebug>10) Printf("Other way Distance (%d)--(%d) %3.3f flag[%d] = %d",faChNeMatchIndex[ich],ich,dist,faChNeMatchIndex[ich]*nChJets+ich,iFlag[faChNeMatchIndex[ich]*nChJets+ich]);
498
499 // reset...
500 faChNeMatchIndex[ich]=-1;
501
502 }
503
504
505 // check for "true" correlations
506 for(int ifu = 0;ifu<nFullJets;ifu++){
507 for(int ich = 0;ich<nChJets;ich++){
508 if(iDebug>10) AliInfo(Form("%s: Flag[%d][%d] %d ",GetName(),ifu,ich,iFlag[ifu*nChJets+ich]));
509
510 if(kMode==3){
511 // we have a uniqe correlation
512 if(iFlag[ifu*nChJets+ich]==3){
513
514 AliEmcalJet *chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
515 AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetJetFromArray(ifu, cFull));
516 Double_t dR = GetDeltaR(fullJet,chJet);
517
518 if(iDebug>10) Printf("closest jets %d %d dR = %f",ich,ifu,dR);
519
520 chJet->SetClosestJet(fullJet,dR);
521 fullJet->SetClosestJet(chJet,dR);
522
523 }
524 }
525 }
526 }
527
528 fMatchingDone = kTRUE;
529
530}
531
532//________________________________________________________________________
533void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndex() {
534
535 // take each full jet and set the index of the charged jet with the largest shared charged fraction
536
537 const Int_t nJetsCh = GetNJets(fContainerCharged);
538 const Int_t nJetsFull = GetNJets(fContainerFull);
539 faFullFracIndex.Set(nJetsFull+1);
540 faFullFracIndex.Reset(-1);
541
542 AliJetContainer *cont = GetJetContainer(fContainerFull);
543 Float_t radius = cont->GetJetRadius();
544
545 for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
546 Double_t frac = 0.;
547 Double_t dist = 10.;
548 AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, fContainerFull));
549 if(!jetFull) {
550 faFullFracIndex.SetAt(-1,ifu);
551 continue;
552 }
553
554 for(Int_t ich = 0; ich<nJetsCh; ich++) {
555 AliEmcalJet *jetCh = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, fContainerCharged));
556 if(!jetCh)
557 continue;
558 Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
559 dist = GetDeltaR(jetFull,jetCh);
560 if(tmpFrac>frac && dist<radius) {
561 frac = tmpFrac;
562 faFullFracIndex.SetAt(ich,ifu);
563 }
564 }
565
566 }
567
568}
569
570//________________________________________________________________________
571void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndexMC() {
572
573 // take each full jet and set the index of the charged jet with the largest shared charged fraction
574
575 const Int_t nJetsCh = GetNJets(fContainerChargedMC);
576 const Int_t nJetsFull = GetNJets(fContainerFullMC);
577 faFullFracIndexMC.Set(nJetsFull);
578 faFullFracIndexMC.Reset(-1);
579
580 AliJetContainer *cont = GetJetContainer(fContainerFullMC);
581 Float_t radius = cont->GetJetRadius();
582
583 for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
584 Double_t frac = 0.;
585 Double_t dist = 10.;
586 AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, fContainerFullMC));
587 if(!jetFull) {
588 faFullFracIndexMC.SetAt(-1,ifu);
589 continue;
590 }
591
592 for(Int_t ich = 0; ich<nJetsCh; ich++) {
593 AliEmcalJet *jetCh = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, fContainerChargedMC));
594 if(!jetCh)
595 continue;
596 Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
597 dist = GetDeltaR(jetFull,jetCh);
598 if(tmpFrac>frac && dist<radius) {
599 frac = tmpFrac;
600 faFullFracIndexMC.SetAt(ich,ifu);
601 }
602 }
603
604 }
605
606}
607
6ab30d5f 608//_______________________________________________________________________
609AliEmcalJet* AliAnalysisTaskEmcalDiJetBase::GetLeadingJetOppositeHemisphere(const Int_t type, const Int_t typea, const AliEmcalJet *jetTrig) {
610
611 // Get leading jet in opposite hemisphere from trigger jet
612 // type = correlation type
613 // typea = container of associated jets
614
615 Int_t nJetsAssoc = GetNJets(typea);
616 Double_t ptLead = -999;
617 Int_t iJetLead = -1;
618 for(Int_t ija=0; ija<nJetsAssoc; ija++) {
619
620 AliEmcalJet *jetAssoc = NULL;
621 if(type==0) {
622 jetAssoc = static_cast<AliEmcalJet*>(GetJetFromArray(ija, typea));
623 if(TMath::Abs(jetAssoc->Eta())>0.5)
624 jetAssoc = NULL;
625 }
626 else
627 jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));
628
629 if(!jetAssoc)
630 continue;
631
632 Double_t dPhi = GetDeltaPhi(jetTrig,jetAssoc);
633 Double_t phiMin = 0.5*TMath::Pi();
634 Double_t phiMax = 1.5*TMath::Pi();
3c709670 635 if(dPhi<phiMin || dPhi>phiMax)
6ab30d5f 636 continue;
637
638 Double_t jetAssocPt = GetJetPt(jetAssoc,typea);
639
640 if(jetAssocPt>ptLead) {
641 ptLead = jetAssocPt;
642 iJetLead = ija;
643 }
644
645 }
646
647 AliEmcalJet *jetAssocLead = static_cast<AliEmcalJet*>(GetJetFromArray(iJetLead, typea));
648
649 return jetAssocLead;
650
651}
652
6e8b6371 653//________________________________________________________________________
654Bool_t AliAnalysisTaskEmcalDiJetBase::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
655{
656 //
657 // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
658 // Get the pt hard bin from the file path
659 // This is to called in Notify and should provide the path to the AOD/ESD file
660 // (Partially copied from AliAnalysisHelperJetTasks)
661
662 TString file(currFile);
663 xsec = 0;
664 trials = 1;
665
666 if(file.Contains(".zip#")){
667 Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
668 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
669 Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
670 file.Replace(pos+1,pos2-pos1,"");
671 }
672 else {
673 // not an archive take the basename....
674 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
675 }
676 //Printf("%s",file.Data());
677
678 // Get the pt hard bin
679 TString strPthard(file);
680
681 strPthard.Remove(strPthard.Last('/'));
682 strPthard.Remove(strPthard.Last('/'));
683 if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
684 strPthard.Remove(0,strPthard.Last('/')+1);
685 if (strPthard.IsDec())
686 pthard = strPthard.Atoi();
687 else
688 AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
689
690 TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
691 if(!fxsec){
692 // next trial fetch the histgram file
693 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
694 if(!fxsec){
695 // not a severe condition but inciate that we have no information
696 return kFALSE;
697 }
698 else{
699 // find the tlist we want to be independtent of the name so use the Tkey
700 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
701 if(!key){
702 fxsec->Close();
703 return kFALSE;
704 }
705 TList *list = dynamic_cast<TList*>(key->ReadObj());
706 if(!list){
707 fxsec->Close();
708 return kFALSE;
709 }
710 xsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
711 trials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
712 fxsec->Close();
713 }
714 } // no tree pyxsec.root
715 else {
716 TTree *xtree = (TTree*)fxsec->Get("Xsection");
717 if(!xtree){
718 fxsec->Close();
719 return kFALSE;
720 }
721 UInt_t ntrials = 0;
722 Double_t xsection = 0;
723 xtree->SetBranchAddress("xsection",&xsection);
724 xtree->SetBranchAddress("ntrials",&ntrials);
725 xtree->GetEntry(0);
726 trials = ntrials;
727 xsec = xsection;
728 fxsec->Close();
729 }
730 return kTRUE;
731}
732
733//________________________________________________________________________
734Bool_t AliAnalysisTaskEmcalDiJetBase::UserNotify()
735{
736 // Get pythia info
737
738 if (!fIsPythiaPtHard) {
739
740 return kTRUE;
741 }
742
743 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
744 if (!tree) {
745 AliError(Form("%s - UserNotify: No current tree!",GetName()));
746 return kFALSE;
747 }
748
749 Float_t xsection = 0;
750 Float_t trials = 0;
751 Int_t pthard = 0;
752
753 TFile *curfile = tree->GetCurrentFile();
754 if (!curfile) {
755 AliError(Form("%s - UserNotify: No current file!",GetName()));
756 return kFALSE;
757 }
758
759 TChain *chain = dynamic_cast<TChain*>(tree);
760 if (chain)
761 tree = chain->GetTree();
762
763 Int_t nevents = tree->GetEntriesFast();
764
765 PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthard);
766
767 fPtHardBin = pthard;
768 // fXsec = xsection;
769 // fTrials = trials;
770
771 fHistTrials->Fill(pthard, trials);
772 fHistXsection->Fill(pthard, xsection);
773 fHistEvents->Fill(pthard, nevents);
774
775 return kTRUE;
776}
777
778//________________________________________________________________________
779Bool_t AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects() {
780 //
781 // get charged jets
782 //
783
784 if (!AliAnalysisTaskEmcalJetDev::RetrieveEventObjects())
785 return kFALSE;
786
787 if(fRhoType==0) {
788 fRhoFullVal = 0.;
789 fRhoChVal = 0.;
790 }
791 if(fRhoType==1) {
792 fRhoFullVal = GetRhoVal(fContainerFull);
793 fRhoChVal = GetRhoVal(fContainerCharged);
794 }
795
796 if (MCEvent()) {
797 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
798 if (!fPythiaHeader) {
799 // Check if AOD
800 AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
801
802 if (aodMCH) {
803 for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
804 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
805 if (fPythiaHeader) break;
806 }
807 }
808 }
809 }
810
811 if (fPythiaHeader) {
812 fPtHard = fPythiaHeader->GetPtHard();
813
814 fNTrials = fPythiaHeader->Trials();
815
816 }
817
818 return kTRUE;
819
820}
821
822//_______________________________________________________________________
823void AliAnalysisTaskEmcalDiJetBase::Terminate(Option_t *)
824{
825 // Called once at the end of the analysis.
826}