2 // ******************************************
3 // This task computes several jet observables like
4 // the fraction of energy in inner and outer coronnas,
5 // jet-track correlations,triggered jet shapes and
6 // correlation strength distribution of particles inside jets.
7 // Author: lcunquei@cern.ch
8 // *******************************************
11 /**************************************************************************
12 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
14 * Author: The ALICE Off-line Project. *
15 * Contributors are mentioned in the code where appropriate. *
17 * Permission to use, copy, modify and distribute this software and its *
18 * documentation strictly for non-commercial purposes is hereby granted *
19 * without fee, provided that the above copyright notice appears in all *
20 * copies and that both the copyright notice and this permission notice *
21 * appear in the supporting documentation. The authors make no claims *
22 * about the suitability of this software for any purpose. It is *
23 * provided "as is" without express or implied warranty. *
24 **************************************************************************/
33 #include "THnSparse.h"
38 #include "AliAnalysisTask.h"
39 #include "AliAnalysisManager.h"
41 #include "AliVEvent.h"
42 #include "AliESDEvent.h"
43 #include "AliESDInputHandler.h"
44 #include "AliCentrality.h"
45 #include "AliAnalysisHelperJetTasks.h"
46 #include "AliInputEventHandler.h"
47 #include "AliAODJetEventBackground.h"
48 #include "AliAnalysisTaskFastEmbedding.h"
49 #include "AliAODEvent.h"
50 #include "AliAODHandler.h"
51 #include "AliAODJet.h"
53 #include "AliAnalysisTaskJetCore.h"
55 ClassImp(AliAnalysisTaskJetCore)
57 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
63 fBackgroundBranch(""),
66 fOfflineTrgMask(AliVEvent::kAny),
79 fAngStructCloseTracks(0),
94 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
95 fJetPtFractionMin(0.5),
103 fHistEvtSelection(0x0),
106 fh2JetCoreMethod1C10(0x0),
107 fh2JetCoreMethod2C10(0x0),
108 fh2JetCoreMethod1C20(0x0),
109 fh2JetCoreMethod2C20(0x0),
110 fh2JetCoreMethod1C30(0x0),
111 fh2JetCoreMethod2C30(0x0),
112 fh2JetCoreMethod1C60(0x0),
113 fh2JetCoreMethod2C60(0x0),
114 fh3JetTrackC3060(0x0),
116 fh3JetTrackC4080(0x0),
117 fh2AngStructpt1C10(0x0),
118 fh2AngStructpt2C10(0x0),
119 fh2AngStructpt3C10(0x0),
120 fh2AngStructpt4C10(0x0),
121 fh2AngStructpt1C20(0x0),
122 fh2AngStructpt2C20(0x0),
123 fh2AngStructpt3C20(0x0),
124 fh2AngStructpt4C20(0x0),
125 fh2AngStructpt1C30(0x0),
126 fh2AngStructpt2C30(0x0),
127 fh2AngStructpt3C30(0x0),
128 fh2AngStructpt4C30(0x0),
129 fh2AngStructpt1C60(0x0),
130 fh2AngStructpt2C60(0x0),
131 fh2AngStructpt3C60(0x0),
132 fh2AngStructpt4C60(0x0),
135 fh2JetDensityA4(0x0),
137 fh3spectriggeredC4080(0x0),
138 fh3spectriggeredC20(0x0),
139 fh3spectriggeredC3060(0x0)
143 // default Constructor
147 for(Int_t i=0; i<10; i++) {
148 for(Int_t j=0; j<6; j++) {
157 fJetBranchName[0] = "";
158 fJetBranchName[1] = "";
160 fListJets[0] = new TList;
161 fListJets[1] = new TList;
164 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
165 AliAnalysisTaskSE(name),
170 fBackgroundBranch(""),
173 fOfflineTrgMask(AliVEvent::kAny),
185 fNInputTracksMax(-1),
186 fAngStructCloseTracks(0),
201 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
202 fJetPtFractionMin(0.5),
210 fHistEvtSelection(0x0),
213 fh2JetCoreMethod1C10(0x0),
214 fh2JetCoreMethod2C10(0x0),
215 fh2JetCoreMethod1C20(0x0),
216 fh2JetCoreMethod2C20(0x0),
217 fh2JetCoreMethod1C30(0x0),
218 fh2JetCoreMethod2C30(0x0),
219 fh2JetCoreMethod1C60(0x0),
220 fh2JetCoreMethod2C60(0x0),
221 fh3JetTrackC3060(0x0),
223 fh3JetTrackC4080(0x0),
224 fh2AngStructpt1C10(0x0),
225 fh2AngStructpt2C10(0x0),
226 fh2AngStructpt3C10(0x0),
227 fh2AngStructpt4C10(0x0),
228 fh2AngStructpt1C20(0x0),
229 fh2AngStructpt2C20(0x0),
230 fh2AngStructpt3C20(0x0),
231 fh2AngStructpt4C20(0x0),
232 fh2AngStructpt1C30(0x0),
233 fh2AngStructpt2C30(0x0),
234 fh2AngStructpt3C30(0x0),
235 fh2AngStructpt4C30(0x0),
236 fh2AngStructpt1C60(0x0),
237 fh2AngStructpt2C60(0x0),
238 fh2AngStructpt3C60(0x0),
239 fh2AngStructpt4C60(0x0),
242 fh2JetDensityA4(0x0),
244 fh3spectriggeredC4080(0x0),
245 fh3spectriggeredC20(0x0),
246 fh3spectriggeredC3060(0x0)
252 for(Int_t i=0; i<10; i++) {
253 for(Int_t j=0; j<6; j++) {
260 fJetBranchName[0] = "";
261 fJetBranchName[1] = "";
263 fListJets[0] = new TList;
264 fListJets[1] = new TList;
266 DefineOutput(1, TList::Class());
269 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
275 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
277 fJetBranchName[0] = branch1;
278 fJetBranchName[1] = branch2;
281 void AliAnalysisTaskJetCore::Init()
284 // check for jet branches
285 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
286 AliError("Jet branch name not set.");
291 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
296 if(!fOutputList) fOutputList = new TList;
297 fOutputList->SetOwner(kTRUE);
299 Bool_t oldStatus = TH1::AddDirectoryStatus();
300 TH1::AddDirectory(kFALSE);
303 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
304 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
305 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
306 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
307 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
308 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
309 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
311 UInt_t entries = 0; // bit coded, see GetDimParams() below
312 entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7;
313 fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
315 //change binning in pTtrack
316 Double_t *xPt3 = new Double_t[10];
318 for(int i = 1; i<=9;i++){
319 if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.4; // 1 - 5
320 else if(xPt3[i-1]<11)xPt3[i] = xPt3[i-1] + 3; // 5 - 12
321 else xPt3[i] = xPt3[i-1] + 150.; // 18
323 fhnDeltaR->SetBinEdges(2,xPt3);
326 //change binning in HTI
327 Double_t *xPt4 = new Double_t[14];
329 for(int i = 1; i<=13;i++){
330 if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 1; // 1 - 10
331 else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 5; // 10 - 12
332 else xPt4[i] = xPt4[i-1] + 30.; // 13
334 fhnDeltaR->SetBinEdges(6,xPt4);
344 UInt_t cifras = 0; // bit coded, see GetDimParams() below
345 cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7;
346 fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);}
350 fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
351 fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
352 fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
353 fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
354 fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
355 fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
356 fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
357 fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
359 fh3JetTrackC3060=new TH3F("JetTrackC3060","",50,0,50,150,0.,150.,34,0,3.4);
360 fh3JetTrackC20=new TH3F("JetTrackC20","",50,0,50,150,0.,150.,34,0,3.4);
361 fh3JetTrackC4080=new TH3F("JetTrackC4080","",50,0,50,150,0.,150.,34,0,3.4);
362 if(fAngStructCloseTracks>0){
363 fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,10.);
364 fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,10.);
365 fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,10.);
366 fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,10.);
367 fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,10.);
368 fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,10.);
369 fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,10.);
370 fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,10.);
371 fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,10.);
372 fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,10.);
373 fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,10.);
374 fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,10.);
375 fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,10.);
376 fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,10.);
377 fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,10.);
378 fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
382 fh2Ntriggers=new TH2F("# of triggers","",10,0.,100.,50,0.,50.);
383 fh2JetDensity=new TH2F("Jet density vs centrality A>0.4","",100,0.,4000.,100,0.,5.);
384 fh2JetDensityA4=new TH2F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.);
385 fh2RPJets=new TH2F("RPJet","",3,0.,3.,150,0.,150.);
386 fh3spectriggeredC4080 = new TH3F("Triggered spectrumC4080","",5,0.,1.,140,-80.,200.,50,0.,50.);
387 fh3spectriggeredC20 = new TH3F("Triggered spectrumC20","",5,0.,1.,140,-80.,200.,50,0.,50.);
388 fh3spectriggeredC3060 = new TH3F("Triggered spectrumC3060","",5,0.,1.,140,-80.,200.,50,0.,50.);
392 fOutputList->Add(fHistEvtSelection);
394 fOutputList->Add(fhnDeltaR);
396 fOutputList->Add(fhnMixedEvents);
402 fOutputList->Add(fh2JetCoreMethod1C10);
403 fOutputList->Add(fh2JetCoreMethod2C10);
404 fOutputList->Add(fh2JetCoreMethod1C20);
405 fOutputList->Add(fh2JetCoreMethod2C20);
406 fOutputList->Add(fh2JetCoreMethod1C30);
407 fOutputList->Add(fh2JetCoreMethod2C30);
408 fOutputList->Add(fh2JetCoreMethod1C60);
409 fOutputList->Add(fh2JetCoreMethod2C60);}
411 fOutputList->Add(fh3JetTrackC3060);
412 fOutputList->Add(fh3JetTrackC20);
413 fOutputList->Add(fh3JetTrackC4080);
418 if(fAngStructCloseTracks>0){
419 fOutputList->Add(fh2AngStructpt1C10);
420 fOutputList->Add(fh2AngStructpt2C10);
421 fOutputList->Add(fh2AngStructpt3C10);
422 fOutputList->Add(fh2AngStructpt4C10);
423 fOutputList->Add(fh2AngStructpt1C20);
424 fOutputList->Add(fh2AngStructpt2C20);
425 fOutputList->Add(fh2AngStructpt3C20);
426 fOutputList->Add(fh2AngStructpt4C20);
427 fOutputList->Add(fh2AngStructpt1C30);
428 fOutputList->Add(fh2AngStructpt2C30);
429 fOutputList->Add(fh2AngStructpt3C30);
430 fOutputList->Add(fh2AngStructpt4C30);
431 fOutputList->Add(fh2AngStructpt1C60);
432 fOutputList->Add(fh2AngStructpt2C60);
433 fOutputList->Add(fh2AngStructpt3C60);
434 fOutputList->Add(fh2AngStructpt4C60);}
440 fOutputList->Add(fh2Ntriggers);
441 fOutputList->Add(fh2JetDensity);
442 fOutputList->Add(fh2JetDensityA4);
443 fOutputList->Add(fh2RPJets);
444 fOutputList->Add(fh3spectriggeredC4080);
445 fOutputList->Add(fh3spectriggeredC20);
446 fOutputList->Add(fh3spectriggeredC3060);
449 // =========== Switch on Sumw2 for all histos ===========
450 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
451 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
456 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
461 TH1::AddDirectory(oldStatus);
463 PostData(1, fOutputList);
466 void AliAnalysisTaskJetCore::UserExec(Option_t *)
470 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
471 AliError("Jet branch name not set.");
475 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
477 AliError("ESD not available");
478 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
480 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
482 static AliAODEvent* aod = 0;
483 // take all other information from the aod we take the tracks from
485 if(!fESD)aod = fAODIn;
490 if(fNonStdFile.Length()!=0){
491 // case that we have an AOD extension we need can fetch the jets from the extended output
492 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
493 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
495 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
502 // -- event selection --
503 fHistEvtSelection->Fill(1); // number of events before event selection
506 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
507 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
508 cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
509 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
510 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
511 fHistEvtSelection->Fill(2);
512 PostData(1, fOutputList);
518 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
519 fHistEvtSelection->Fill(3);
520 PostData(1, fOutputList);
522 AliAODVertex* primVtx = aod->GetPrimaryVertex();
525 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
526 fHistEvtSelection->Fill(3);
527 PostData(1, fOutputList);
531 Int_t nTracksPrim = primVtx->GetNContributors();
532 if ((nTracksPrim < fMinContribVtx) ||
533 (primVtx->GetZ() < fVtxZMin) ||
534 (primVtx->GetZ() > fVtxZMax) ){
535 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
536 fHistEvtSelection->Fill(3);
537 PostData(1, fOutputList);
541 // event class selection (from jet helper task)
542 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
543 if(fDebug) Printf("Event class %d", eventClass);
544 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
545 fHistEvtSelection->Fill(4);
546 PostData(1, fOutputList);
550 // centrality selection
551 AliCentrality *cent = 0x0;
552 Double_t centValue = 0.;
554 if(fESD) {cent = fESD->GetCentrality();
555 if(cent) centValue = cent->GetCentralityPercentile("V0M");}
556 else centValue=aod->GetHeader()->GetCentrality();
558 if(fDebug) printf("centrality: %f\n", centValue);
559 if (centValue < fCentMin || centValue > fCentMax){
560 fHistEvtSelection->Fill(4);
561 PostData(1, fOutputList);
566 fHistEvtSelection->Fill(0);
568 // -- end event selection --
571 AliAODJetEventBackground* externalBackground = 0;
572 if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
573 externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
574 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
576 if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
577 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
578 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
581 if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
582 externalBackground = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
583 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
589 if(externalBackground)rho = externalBackground->GetBackground(0);}
591 if(externalBackground)rho = externalBackground->GetBackground(2);}
594 TClonesArray *aodJets[2];
596 if(fAODOut&&!aodJets[0]){
597 aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
598 aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); }
599 if(fAODExtension && !aodJets[0]){
600 aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
601 aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
602 if(fAODIn&&!aodJets[0]){
603 aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
604 aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
607 //Double_t ptsub[aodJets[0]->GetEntriesFast()];
608 //Int_t inord[aodJets[0]->GetEntriesFast()];
609 //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
614 Int_t nT = GetListOfTracks(&ParticleList);
615 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
616 fListJets[iJetType]->Clear();
617 if (!aodJets[iJetType]) continue;
619 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
622 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
623 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
624 if (jet) fListJets[iJetType]->Add(jet);
626 // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
635 Double_t areasmall=0;
636 Double_t phismall=0.;
642 Int_t trigBBTrack=-1;
643 Int_t trigInTrack=-1;
644 fRPAngle = aod->GetHeader()->GetEventplane();
646 AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);
648 PostData(1, fOutputList);
650 fh2Ntriggers->Fill(centValue,partback->Pt());
651 Double_t accep=2.*TMath::Pi()*1.8;
654 for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
655 AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
656 etabig = jetbig->Eta();
657 phibig = jetbig->Phi();
658 ptbig = jetbig->Pt();
659 if(ptbig==0) continue;
660 Int_t phiBin = GetPhiBin(phibig-fRPAngle);
661 areabig = jetbig->EffectiveAreaCharged();
662 Double_t ptcorr=ptbig-rho*areabig;
663 if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
665 if(areabig>=0.2) injet=injet+1;
666 if(areabig>=0.4) injet4=injet4+1;
667 Double_t dphi=RelativePhi(partback->Phi(),phibig);
670 Double_t etadif= partback->Eta()-etabig;
671 if(TMath::Abs(etadif)<=0.5){
672 if(centValue>40. && centValue<80.) fh3JetTrackC4080->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
673 if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
674 if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}}
677 if(centValue>40. && centValue<80.) fh3JetTrackC4080->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
678 if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
679 if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}
684 if(fFlagPhiBkg!=0) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
685 if(fFlagPhiBkg==0) if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
686 if(centValue<20.) fh2RPJets->Fill(phiBin, ptcorr);
687 Double_t dismin=100.;
691 if(centValue>40. && centValue<80.) fh3spectriggeredC4080->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
692 if(centValue<20.) fh3spectriggeredC20->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
693 if(centValue>30. && centValue<60.) fh3spectriggeredC3060->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
695 if(ptcorr<=0) continue;
696 AliAODTrack* leadtrack;
699 TRefArray *genTrackList = jetbig->GetRefTracks();
700 Int_t nTracksGenJet = genTrackList->GetEntriesFast();
701 AliAODTrack* genTrack;
702 for(Int_t ir=0; ir<nTracksGenJet; ++ir){
703 genTrack = (AliAODTrack*)(genTrackList->At(ir));
704 if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
706 leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
707 if(!leadtrack) continue;
709 //store one trigger info
718 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
719 AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
720 etasmall = jetsmall->Eta();
721 phismall = jetsmall->Phi();
722 ptsmall = jetsmall->Pt();
723 areasmall = jetsmall->EffectiveAreaCharged();
724 Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
725 tmpDeltaR=TMath::Sqrt(tmpDeltaR);
726 //Fraction in the jet core
727 if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
729 if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
730 index1=j;}} //en of loop over R=0.2 jets
731 //method1:most concentric jet=core
732 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
733 if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
734 if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
735 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
736 if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
737 //method2:hardest contained jet=core
739 AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
740 if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
741 if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig);
742 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
743 if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
746 if(fDoEventMixing==0){
747 for(int it = 0;it<ParticleList.GetEntries();++it){
748 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
749 Double_t deltaR = jetbig->DeltaR(part);
750 Double_t deltaEta = etabig-part->Eta();
753 Double_t deltaPhi=phibig-part->Phi();
754 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
755 if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
756 Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,leadtrack->Pt(),partback->Pt()}; fhnDeltaR->Fill(jetEntries);}
759 //end of track loop, we only do it if EM is switched off
770 if(injet>0) fh2JetDensity->Fill(ParticleList.GetEntries(),injet/accep);
771 if(injet4>0)fh2JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep);
777 if(fDoEventMixing>0){
778 //check before if the trigger exists
779 // fTrigBuffer[i][0] = zvtx
780 // fTrigBuffer[i][1] = phi
781 // fTrigBuffer[i][2] = eta
782 // fTrigBuffer[i][3] = pt_jet
783 // fTrigBuffer[i][4] = pt_trig
784 // fTrigBuffer[i][5]= centrality
785 if(fTindex==10) fTindex=0;
786 if(fTrigBuffer[fTindex][3]>0){
787 if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
788 if (TMath::Abs(fTrigBuffer[fTindex][5]-centValue<5)){
790 for(int it = 0;it<nT;++it){
791 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
792 Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
793 Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
794 Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
795 if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
796 if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
797 Double_t triggerEntries[7] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4]};
798 fhnMixedEvents->Fill(triggerEntries);
801 if(fNevents==10) fTindex=fTindex+1;
803 if(fTindex==10&&fNevents==10) fCountAgain=0;
805 // Copy the triggers from the current event into the buffer.
806 //again, only if the trigger exists:
809 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet)); AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);
810 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
811 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
812 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
813 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
814 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
815 fTrigBuffer[fTrigBufferIndex][5] = centValue;
817 if(fTrigBufferIndex==9) {fTrigBufferIndex=0;
826 //////////////////ANGULAR STRUCTURE//////////////////////////////////////
828 //tracks up to R=0.8 distant from the jet axis
829 // if(fAngStructCloseTracks==1){
830 // TList CloseTrackList;
831 // Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
832 // Double_t difR=0.04;
833 // for(Int_t l=0;l<15;l++){
834 // Double_t rr=l*0.1+0.1;
835 // for(int it = 0;it<nn;++it){
836 // AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
837 // for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
838 // AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
839 // Double_t ptm=part1->Pt();
840 // Double_t ptn=part2->Pt();
841 // Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
842 // Rnm=TMath::Sqrt(Rnm);
843 // Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
844 // Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
845 // if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
846 // down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
847 // if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
848 // down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
849 // if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
850 // down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
851 // if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
852 // down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
855 // //only jet constituents
856 // if(fAngStructCloseTracks==2){
858 // Double_t difR=0.04;
859 // for(Int_t l=0;l<15;l++){
860 // Double_t rr=l*0.1+0.1;
863 // AliAODTrack* part1;
864 // AliAODTrack* part2;
866 // TRefArray *genTrackListb = jetbig->GetRefTracks();
867 // Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
871 // for(Int_t it=0; it<nTracksGenJetb; ++it){
872 // part1 = (AliAODTrack*)(genTrackListb->At(it));
873 // for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
874 // part2 = (AliAODTrack*)(genTrackListb->At(itu));
875 // Double_t ptm=part1->Pt();
876 // Double_t ptn=part2->Pt();
877 // Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
878 // Rnm=TMath::Sqrt(Rnm);
879 // Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
880 // Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
881 // if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
882 // down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
883 // if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
884 // down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
885 // if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
886 // down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
887 // if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
888 // down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
890 // //end loop over R=0.4 jets
891 // if(fAngStructCloseTracks>0){
892 // for(Int_t l=0;l<15;l++){
893 // Double_t rr=l*0.1+0.1;
895 // if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
896 // if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
897 // if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
898 // if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
900 // if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
901 // if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
902 // if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
903 // if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
905 // if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
906 // if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
907 // if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
908 // if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
910 // if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
911 // if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
912 // if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
913 // if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
921 PostData(1, fOutputList);
924 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
926 // Draw result to the screen
927 // Called once at the end of the query
929 if (!GetOutputData(1))
938 Int_t AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
941 AliAODEvent *aod = 0;
942 if(!fESD)aod = fAODIn;
946 for(int it = 0;it < aod->GetNumberOfTracks();++it){
947 AliAODTrack *tr = aod->GetTrack(it);
948 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
949 if(TMath::Abs(tr->Eta())>0.9)continue;
950 if(tr->Pt()<0.15)continue;
953 if(tr->Pt()>ptmax){ ptmax=tr->Pt();
963 Int_t AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
965 AliAODEvent *aod = 0;
966 if(!fESD)aod = fAODIn;
973 for(int it = 0;it < aod->GetNumberOfTracks();++it){
974 AliAODTrack *tr = aod->GetTrack(it);
975 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
976 if(TMath::Abs(tr->Eta())>0.9)continue;
977 if(tr->Pt()<0.15)continue;
979 dphi=RelativePhi(tr->Phi(),jetbig->Phi());
980 if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
981 if(tr->Pt()>ptmax){ ptmax=tr->Pt();
997 Int_t AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1000 AliAODEvent *aod = 0;
1001 if(!fESD)aod = fAODIn;
1004 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1005 AliAODTrack *tr = aod->GetTrack(it);
1006 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1007 if(TMath::Abs(tr->Eta())>0.9)continue;
1008 if(tr->Pt()<0.15)continue;
1009 Double_t disR=jetbig->DeltaR(tr);
1010 if(disR>0.8) continue;
1012 //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1031 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1034 Int_t nInputTracks = 0;
1035 AliAODEvent *aod = 0;
1036 if(!fESD)aod = fAODIn;
1038 TString jbname(fJetBranchName[1]);
1039 //needs complete event, use jets without background subtraction
1040 for(Int_t i=1; i<=3; ++i){
1041 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1043 // use only HI event
1044 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1045 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1047 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1048 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1050 Printf("Jet branch %s not found", jbname.Data());
1051 Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1055 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1056 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1058 TRefArray *trackList = jet->GetRefTracks();
1059 Int_t nTracks = trackList->GetEntriesFast();
1060 nInputTracks += nTracks;
1061 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1063 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1065 return nInputTracks;
1070 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1072 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1073 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1074 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1075 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1076 double dphi = mphi-vphi;
1077 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1078 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1079 return dphi;//dphi in [-Pi, Pi]
1082 Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1085 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1086 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1087 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1088 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1095 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1097 // generate new THnSparseF, axes are defined in GetDimParams()
1100 UInt_t tmp = entries;
1103 tmp = tmp &~ -tmp; // clear lowest bit
1106 TString hnTitle(name);
1107 const Int_t dim = count;
1114 while(c<dim && i<32){
1118 GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1119 hnTitle += Form(";%s",label.Data());
1127 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1130 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1132 // stores label and binning of axis for THnSparse
1134 const Double_t pi = TMath::Pi();
1139 label = "V0 centrality (%)";
1148 label = "corrected jet pt";
1191 label = "leading track";
1199 label = "trigger track";