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 "AliAODMCParticle.h"
49 #include "AliAnalysisTaskFastEmbedding.h"
50 #include "AliAODEvent.h"
51 #include "AliAODHandler.h"
52 #include "AliAODJet.h"
54 #include "AliAnalysisTaskJetCore.h"
59 ClassImp(AliAnalysisTaskJetCore)
61 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
67 fBackgroundBranch(""),
70 fOfflineTrgMask(AliVEvent::kAny),
85 fAngStructCloseTracks(0),
94 fTrackTypeRec(kTrackUndef),
104 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
105 fJetPtFractionMin(0.5),
113 fHistEvtSelection(0x0),
116 fh2JetCoreMethod1C10(0x0),
117 fh2JetCoreMethod2C10(0x0),
118 fh2JetCoreMethod1C20(0x0),
119 fh2JetCoreMethod2C20(0x0),
120 fh2JetCoreMethod1C30(0x0),
121 fh2JetCoreMethod2C30(0x0),
122 fh2JetCoreMethod1C60(0x0),
123 fh2JetCoreMethod2C60(0x0),
124 fh3JetTrackC3060(0x0),
126 fh2AngStructpt1C10(0x0),
127 fh2AngStructpt2C10(0x0),
128 fh2AngStructpt3C10(0x0),
129 fh2AngStructpt4C10(0x0),
130 fh2AngStructpt1C20(0x0),
131 fh2AngStructpt2C20(0x0),
132 fh2AngStructpt3C20(0x0),
133 fh2AngStructpt4C20(0x0),
134 fh2AngStructpt1C30(0x0),
135 fh2AngStructpt2C30(0x0),
136 fh2AngStructpt3C30(0x0),
137 fh2AngStructpt4C30(0x0),
138 fh2AngStructpt1C60(0x0),
139 fh2AngStructpt2C60(0x0),
140 fh2AngStructpt3C60(0x0),
141 fh2AngStructpt4C60(0x0),
143 fh2Ntriggers2C10(0x0),
144 fh2Ntriggers2C20(0x0),
146 fh3JetDensityA4(0x0),
151 fh3spectriggeredC10(0x0),
152 fh3spectriggeredC20(0x0),
153 fh3spectriggeredC3060(0x0)
157 // default Constructor
161 for(Int_t i=0; i<10; i++) {
162 for(Int_t j=0; j<6; j++) {
171 fJetBranchName[0] = "";
172 fJetBranchName[1] = "";
174 fListJets[0] = new TList;
175 fListJets[1] = new TList;
178 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
179 AliAnalysisTaskSE(name),
184 fBackgroundBranch(""),
187 fOfflineTrgMask(AliVEvent::kAny),
194 fFilterMaskBestPt(0),
201 fNInputTracksMax(-1),
202 fAngStructCloseTracks(0),
211 fTrackTypeRec(kTrackUndef),
221 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
222 fJetPtFractionMin(0.5),
230 fHistEvtSelection(0x0),
233 fh2JetCoreMethod1C10(0x0),
234 fh2JetCoreMethod2C10(0x0),
235 fh2JetCoreMethod1C20(0x0),
236 fh2JetCoreMethod2C20(0x0),
237 fh2JetCoreMethod1C30(0x0),
238 fh2JetCoreMethod2C30(0x0),
239 fh2JetCoreMethod1C60(0x0),
240 fh2JetCoreMethod2C60(0x0),
241 fh3JetTrackC3060(0x0),
243 fh2AngStructpt1C10(0x0),
244 fh2AngStructpt2C10(0x0),
245 fh2AngStructpt3C10(0x0),
246 fh2AngStructpt4C10(0x0),
247 fh2AngStructpt1C20(0x0),
248 fh2AngStructpt2C20(0x0),
249 fh2AngStructpt3C20(0x0),
250 fh2AngStructpt4C20(0x0),
251 fh2AngStructpt1C30(0x0),
252 fh2AngStructpt2C30(0x0),
253 fh2AngStructpt3C30(0x0),
254 fh2AngStructpt4C30(0x0),
255 fh2AngStructpt1C60(0x0),
256 fh2AngStructpt2C60(0x0),
257 fh2AngStructpt3C60(0x0),
258 fh2AngStructpt4C60(0x0),
260 fh2Ntriggers2C10(0x0),
261 fh2Ntriggers2C20(0x0),
263 fh3JetDensityA4(0x0),
268 fh3spectriggeredC10(0x0),
269 fh3spectriggeredC20(0x0),
270 fh3spectriggeredC3060(0x0)
276 for(Int_t i=0; i<10; i++) {
277 for(Int_t j=0; j<6; j++) {
284 fJetBranchName[0] = "";
285 fJetBranchName[1] = "";
287 fListJets[0] = new TList;
288 fListJets[1] = new TList;
290 DefineOutput(1, TList::Class());
293 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
299 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
301 fJetBranchName[0] = branch1;
302 fJetBranchName[1] = branch2;
305 void AliAnalysisTaskJetCore::Init()
308 // check for jet branches
309 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
310 AliError("Jet branch name not set.");
315 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
320 if(!fOutputList) fOutputList = new TList;
321 fOutputList->SetOwner(kTRUE);
323 Bool_t oldStatus = TH1::AddDirectoryStatus();
324 TH1::AddDirectory(kFALSE);
327 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
328 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
329 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
330 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
331 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
332 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
333 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
335 UInt_t entries = 0; // bit coded, see GetDimParams() below
336 entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7;
337 fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
339 //change binning in pTtrack
340 Double_t *xPt3 = new Double_t[10];
342 for(int i = 1; i<=9;i++){
343 if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.4; // 1 - 5
344 else if(xPt3[i-1]<11)xPt3[i] = xPt3[i-1] + 3; // 5 - 12
345 else xPt3[i] = xPt3[i-1] + 150.; // 18
347 fhnDeltaR->SetBinEdges(2,xPt3);
350 //change binning in HTI
351 Double_t *xPt4 = new Double_t[14];
353 for(int i = 1; i<=13;i++){
354 if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 1; // 1 - 10
355 else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 5; // 10 - 12
356 else xPt4[i] = xPt4[i-1] + 30.; // 13
358 fhnDeltaR->SetBinEdges(6,xPt4);
368 UInt_t cifras = 0; // bit coded, see GetDimParams() below
369 cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7;
370 fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);}
373 fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
374 fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
375 fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
376 fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
377 fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
378 fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
379 fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
380 fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
382 fh3JetTrackC3060=new TH3F("JetTrackC3060","",50,0,50,150,0.,150.,35,0.,3.5);
383 fh3JetTrackC20=new TH3F("JetTrackC20","",50,0,50,150,0.,150.,35,0.,3.5);
384 if(fAngStructCloseTracks>0){
385 fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,10.);
386 fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,10.);
387 fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,10.);
388 fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,10.);
389 fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,10.);
390 fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,10.);
391 fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,10.);
392 fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,10.);
393 fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,10.);
394 fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,10.);
395 fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,10.);
396 fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,10.);
397 fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,10.);
398 fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,10.);
399 fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,10.);
400 fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
404 fh2Ntriggers=new TH2F("# of triggers","",10,0.,100.,50,0.,50.);
405 fh2Ntriggers2C10=new TH2F("# of triggers2C10","",50,0.,50.,50,0.,50.);
406 fh2Ntriggers2C20=new TH2F("# of triggers2C20","",50,0.,50.,50,0.,50.);
407 fh3JetDensity=new TH3F("Jet density vs mutliplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
408 fh3JetDensityA4=new TH3F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
409 fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,35,100,0.,100.);
410 fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,35,100,0.,100.);
411 fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,35,50,0.,50.);
412 fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,35,50,0.,50.);
413 fh3spectriggeredC10 = new TH3F("Triggered spectrumC10","",100,0.,1.,140,-80.,200.,50,0.,50.);
414 fh3spectriggeredC20 = new TH3F("Triggered spectrumC20","",100,0.,1.,140,-80.,200.,50,0.,50.);
415 fh3spectriggeredC3060 = new TH3F("Triggered spectrumC3060","",100,0.,1.,140,-80.,200.,10,0.,50.);
419 fOutputList->Add(fHistEvtSelection);
421 fOutputList->Add(fhnDeltaR);
423 fOutputList->Add(fhnMixedEvents);
429 fOutputList->Add(fh2JetCoreMethod1C10);
430 fOutputList->Add(fh2JetCoreMethod2C10);
431 fOutputList->Add(fh2JetCoreMethod1C20);
432 fOutputList->Add(fh2JetCoreMethod2C20);
433 fOutputList->Add(fh2JetCoreMethod1C30);
434 fOutputList->Add(fh2JetCoreMethod2C30);
435 fOutputList->Add(fh2JetCoreMethod1C60);
436 fOutputList->Add(fh2JetCoreMethod2C60);}
438 fOutputList->Add(fh3JetTrackC3060);
439 fOutputList->Add(fh3JetTrackC20);
445 if(fAngStructCloseTracks>0){
446 fOutputList->Add(fh2AngStructpt1C10);
447 fOutputList->Add(fh2AngStructpt2C10);
448 fOutputList->Add(fh2AngStructpt3C10);
449 fOutputList->Add(fh2AngStructpt4C10);
450 fOutputList->Add(fh2AngStructpt1C20);
451 fOutputList->Add(fh2AngStructpt2C20);
452 fOutputList->Add(fh2AngStructpt3C20);
453 fOutputList->Add(fh2AngStructpt4C20);
454 fOutputList->Add(fh2AngStructpt1C30);
455 fOutputList->Add(fh2AngStructpt2C30);
456 fOutputList->Add(fh2AngStructpt3C30);
457 fOutputList->Add(fh2AngStructpt4C30);
458 fOutputList->Add(fh2AngStructpt1C60);
459 fOutputList->Add(fh2AngStructpt2C60);
460 fOutputList->Add(fh2AngStructpt3C60);
461 fOutputList->Add(fh2AngStructpt4C60);}
467 fOutputList->Add(fh2Ntriggers);
468 fOutputList->Add(fh2Ntriggers2C10);
469 fOutputList->Add(fh2Ntriggers2C20);
470 fOutputList->Add(fh3JetDensity);
471 fOutputList->Add(fh3JetDensityA4);
472 fOutputList->Add(fh2RPJetsC10);
473 fOutputList->Add(fh2RPJetsC20);
474 fOutputList->Add(fh2RPTC10);
475 fOutputList->Add(fh2RPTC20);
476 fOutputList->Add(fh3spectriggeredC10);
477 fOutputList->Add(fh3spectriggeredC20);
478 fOutputList->Add(fh3spectriggeredC3060);
481 // =========== Switch on Sumw2 for all histos ===========
482 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
483 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
488 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
493 TH1::AddDirectory(oldStatus);
495 PostData(1, fOutputList);
498 void AliAnalysisTaskJetCore::UserExec(Option_t *)
502 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
503 AliError("Jet branch name not set.");
507 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
509 AliError("ESD not available");
510 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
512 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
514 static AliAODEvent* aod = 0;
515 // take all other information from the aod we take the tracks from
517 if(!fESD)aod = fAODIn;
522 if(fNonStdFile.Length()!=0){
523 // case that we have an AOD extension we need can fetch the jets from the extended output
524 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
525 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
527 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
531 // -- event selection --
532 fHistEvtSelection->Fill(1); // number of events before event selection
535 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
536 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
537 std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
538 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
539 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
540 fHistEvtSelection->Fill(2);
541 PostData(1, fOutputList);
547 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
548 fHistEvtSelection->Fill(3);
549 PostData(1, fOutputList);
551 AliAODVertex* primVtx = aod->GetPrimaryVertex();
554 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
555 fHistEvtSelection->Fill(3);
556 PostData(1, fOutputList);
560 Int_t nTracksPrim = primVtx->GetNContributors();
561 if ((nTracksPrim < fMinContribVtx) ||
562 (primVtx->GetZ() < fVtxZMin) ||
563 (primVtx->GetZ() > fVtxZMax) ){
564 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
565 fHistEvtSelection->Fill(3);
566 PostData(1, fOutputList);
570 // event class selection (from jet helper task)
571 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
572 if(fDebug) Printf("Event class %d", eventClass);
573 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
574 fHistEvtSelection->Fill(4);
575 PostData(1, fOutputList);
579 // centrality selection
580 AliCentrality *cent = 0x0;
581 Double_t centValue = 0.;
583 if(fESD) {cent = fESD->GetCentrality();
584 if(cent) centValue = cent->GetCentralityPercentile("V0M");}
585 else centValue=aod->GetHeader()->GetCentrality();
587 if(fDebug) printf("centrality: %f\n", centValue);
588 if (centValue < fCentMin || centValue > fCentMax){
589 fHistEvtSelection->Fill(4);
590 PostData(1, fOutputList);
595 fHistEvtSelection->Fill(0);
597 // -- end event selection --
600 AliAODJetEventBackground* externalBackground = 0;
601 if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
602 externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
603 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
605 if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
606 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
607 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
610 if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
611 externalBackground = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
612 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
618 if(externalBackground)rho = externalBackground->GetBackground(0);}
620 if(externalBackground)rho = externalBackground->GetBackground(2);}
623 TClonesArray *aodJets[2];
625 if(fAODOut&&!aodJets[0]){
626 aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
627 aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); }
628 if(fAODExtension && !aodJets[0]){
629 aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
630 aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
631 if(fAODIn&&!aodJets[0]){
632 aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
633 aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
636 //Double_t ptsub[aodJets[0]->GetEntriesFast()];
637 //Int_t inord[aodJets[0]->GetEntriesFast()];
638 //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
643 Int_t nT = GetListOfTracks(&ParticleList);
644 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
645 fListJets[iJetType]->Clear();
646 if (!aodJets[iJetType]) continue;
647 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
648 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
649 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
650 if (jet) fListJets[iJetType]->Add(jet);
652 // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
661 Double_t areasmall=0;
662 Double_t phismall=0.;
668 Int_t trigBBTrack=-1;
669 Int_t trigInTrack=-1;
670 fRPAngle = aod->GetHeader()->GetEventplane();
672 AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);
674 PostData(1, fOutputList);
678 //for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
679 //if(fFlagOnlyHardest!=0){if(tt!=nT) continue;}
680 //AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);
682 fh2Ntriggers->Fill(centValue,partback->Pt());
683 Double_t phiBinT = RelativePhi(partback->Phi(),fRPAngle);
684 if(centValue<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
685 if(centValue<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
686 Double_t accep=2.*TMath::Pi()*1.8;
690 for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
691 AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
692 etabig = jetbig->Eta();
693 phibig = jetbig->Phi();
694 ptbig = jetbig->Pt();
695 if(ptbig==0) continue;
696 Double_t phiBin = RelativePhi(phibig,fRPAngle);
697 areabig = jetbig->EffectiveAreaCharged();
698 Double_t ptcorr=ptbig-rho*areabig;
699 if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
700 if(areabig>=0.07) injet=injet+1;
701 if(areabig>=0.4) injet4=injet4+1;
702 Double_t dphi=RelativePhi(partback->Phi(),phibig);
705 Double_t etadif= partback->Eta()-etabig;
706 if(TMath::Abs(etadif)<=0.5){
708 if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
709 if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}}
711 if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
712 if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}
715 if(fFlagJetHadron==0){
716 if(fFlagPhiBkg!=0) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
717 if(fFlagPhiBkg==0) if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;}
719 if(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
722 if(centValue<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
723 if(centValue<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
724 Double_t dismin=100.;
729 if(centValue<10.) fh3spectriggeredC10->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
730 if(centValue<20.) fh3spectriggeredC20->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
731 if(centValue>30. && centValue<60.) fh3spectriggeredC3060->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
733 if(ptcorr<=0) continue;
735 AliAODTrack* leadtrack=0;
738 if(fFlagJetHadron==0){
739 TRefArray *genTrackList = jetbig->GetRefTracks();
740 Int_t nTracksGenJet = genTrackList->GetEntriesFast();
741 AliAODTrack* genTrack;
742 for(Int_t ir=0; ir<nTracksGenJet; ++ir){
743 genTrack = (AliAODTrack*)(genTrackList->At(ir));
744 if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
746 leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
747 if(!leadtrack) continue;}
749 AliVParticle* leadtrackb=0;
750 if(fFlagJetHadron!=0){
751 Int_t nTb = GetHardestTrackBackToJet(jetbig);
752 leadtrackb = (AliVParticle*)ParticleList.At(nTb);
753 if(!leadtrackb) continue;
760 //store one trigger info
769 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
770 AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
771 etasmall = jetsmall->Eta();
772 phismall = jetsmall->Phi();
773 ptsmall = jetsmall->Pt();
774 areasmall = jetsmall->EffectiveAreaCharged();
775 Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
776 tmpDeltaR=TMath::Sqrt(tmpDeltaR);
777 //Fraction in the jet core
778 if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
780 if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
781 index1=j;}} //en of loop over R=0.2 jets
782 //method1:most concentric jet=core
783 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
784 if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
785 if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
786 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
787 if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
788 //method2:hardest contained jet=core
790 AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
791 if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
792 if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig);
793 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
794 if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
795 if(centValue<10) fh2Ntriggers2C10->Fill(leadtrack->Pt(),partback->Pt());
796 if(centValue<20) fh2Ntriggers2C20->Fill(leadtrack->Pt(),partback->Pt());
797 if(fDoEventMixing==0 && fFlagOnlyRecoil==0){
798 for(int it = 0;it<ParticleList.GetEntries();++it){
799 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
800 Double_t deltaR = jetbig->DeltaR(part);
801 Double_t deltaEta = etabig-part->Eta();
803 Double_t deltaPhi=phibig-part->Phi();
804 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
805 if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
807 if(fFlagJetHadron==0) pTcont=leadtrack->Pt();
808 if(fFlagJetHadron!=0) pTcont=leadtrackb->Pt();
809 Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,pTcont,partback->Pt()};
810 fhnDeltaR->Fill(jetEntries);}
816 //end of track loop, we only do it if EM is switched off
827 if(injet>0) fh3JetDensity->Fill(ParticleList.GetEntries(),injet/accep,partback->Pt());
828 if(injet4>0)fh3JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep,partback->Pt());
834 if(fDoEventMixing>0){
835 //check before if the trigger exists
836 // fTrigBuffer[i][0] = zvtx
837 // fTrigBuffer[i][1] = phi
838 // fTrigBuffer[i][2] = eta
839 // fTrigBuffer[i][3] = pt_jet
840 // fTrigBuffer[i][4] = pt_trig
841 // fTrigBuffer[i][5]= centrality
842 if(fTindex==10) fTindex=0;
843 if(fTrigBuffer[fTindex][3]>0){
844 if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
845 if (TMath::Abs(fTrigBuffer[fTindex][5]-centValue<5)){
847 for(int it = 0;it<nT;++it){
848 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
849 Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
850 Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
851 Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
852 if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
853 if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
854 Double_t triggerEntries[7] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4]};
855 fhnMixedEvents->Fill(triggerEntries);
858 if(fNevents==10) fTindex=fTindex+1;
861 if(fTindex==10&&fNevents==10) fCountAgain=0;
863 // Copy the triggers from the current event into the buffer.
864 //again, only if the trigger exists:
867 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet)); AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);
868 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
869 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
870 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
871 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
872 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
873 fTrigBuffer[fTrigBufferIndex][5] = centValue;
875 if(fTrigBufferIndex==9) {fTrigBufferIndex=0;
884 //////////////////ANGULAR STRUCTURE//////////////////////////////////////
886 //tracks up to R=0.8 distant from the jet axis
887 // if(fAngStructCloseTracks==1){
888 // TList CloseTrackList;
889 // Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
890 // Double_t difR=0.04;
891 // for(Int_t l=0;l<15;l++){
892 // Double_t rr=l*0.1+0.1;
893 // for(int it = 0;it<nn;++it){
894 // AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
895 // for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
896 // AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
897 // Double_t ptm=part1->Pt();
898 // Double_t ptn=part2->Pt();
899 // Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
900 // Rnm=TMath::Sqrt(Rnm);
901 // Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
902 // Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
903 // if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
904 // down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
905 // if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
906 // down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
907 // if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
908 // down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
909 // if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
910 // down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
913 // //only jet constituents
914 // if(fAngStructCloseTracks==2){
916 // Double_t difR=0.04;
917 // for(Int_t l=0;l<15;l++){
918 // Double_t rr=l*0.1+0.1;
921 // AliAODTrack* part1;
922 // AliAODTrack* part2;
924 // TRefArray *genTrackListb = jetbig->GetRefTracks();
925 // Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
929 // for(Int_t it=0; it<nTracksGenJetb; ++it){
930 // part1 = (AliAODTrack*)(genTrackListb->At(it));
931 // for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
932 // part2 = (AliAODTrack*)(genTrackListb->At(itu));
933 // Double_t ptm=part1->Pt();
934 // Double_t ptn=part2->Pt();
935 // Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
936 // Rnm=TMath::Sqrt(Rnm);
937 // Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
938 // Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
939 // if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
940 // down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
941 // if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
942 // down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
943 // if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
944 // down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
945 // if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
946 // down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
948 // //end loop over R=0.4 jets
949 // if(fAngStructCloseTracks>0){
950 // for(Int_t l=0;l<15;l++){
951 // Double_t rr=l*0.1+0.1;
953 // if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
954 // if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
955 // if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
956 // if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
958 // if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
959 // if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
960 // if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
961 // if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
963 // if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
964 // if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
965 // if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
966 // if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
968 // if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
969 // if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
970 // if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
971 // if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
979 PostData(1, fOutputList);
982 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
984 // Draw result to the screen
985 // Called once at the end of the query
987 if (!GetOutputData(1))
996 Int_t AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
999 AliAODEvent *aod = 0;
1001 if(!fESD)aod = fAODIn;
1008 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1009 AliAODTrack *tr = aod->GetTrack(it);
1010 Bool_t bGood = false;
1011 if(fFilterType == 0)bGood = true;
1012 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1013 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1014 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1015 if(bGood==false) continue;
1016 if(TMath::Abs(tr->Eta())>0.9)continue;
1017 if(tr->Pt()<0.15)continue;
1020 if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
1021 if(tr->TestFilterBit(fFilterMaskBestPt)){
1037 // else if (type == kTrackAODMCCharged) {
1038 // TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1039 // if(!tca)return iCount;
1040 // for(int it = 0;it < tca->GetEntriesFast();++it){
1041 // AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1042 // if(!part)continue;
1043 // if(part->Pt()<0.15)continue;
1044 // if(!part->IsPhysicalPrimary())continue;
1045 // if(part->Charge()==0)continue;
1046 // if(TMath::Abs(part->Eta())>0.9)continue;
1049 // if(part->Pt()>ptmax){ ptmax=part->Pt();
1050 // index=iCount-1;}}}
1055 Int_t AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
1057 AliAODEvent *aod = 0;
1058 if(!fESD)aod = fAODIn;
1065 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1066 AliAODTrack *tr = aod->GetTrack(it);
1067 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1068 if(TMath::Abs(tr->Eta())>0.9)continue;
1069 if(tr->Pt()<0.15)continue;
1071 dphi=RelativePhi(tr->Phi(),jetbig->Phi());
1072 if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
1073 if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1089 Int_t AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1092 AliAODEvent *aod = 0;
1093 if(!fESD)aod = fAODIn;
1096 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1097 AliAODTrack *tr = aod->GetTrack(it);
1098 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1099 if(TMath::Abs(tr->Eta())>0.9)continue;
1100 if(tr->Pt()<0.15)continue;
1101 Double_t disR=jetbig->DeltaR(tr);
1102 if(disR>0.8) continue;
1104 //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1123 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1126 Int_t nInputTracks = 0;
1127 AliAODEvent *aod = 0;
1128 if(!fESD)aod = fAODIn;
1130 TString jbname(fJetBranchName[1]);
1131 //needs complete event, use jets without background subtraction
1132 for(Int_t i=1; i<=3; ++i){
1133 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1135 // use only HI event
1136 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1137 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1139 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1140 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1142 Printf("Jet branch %s not found", jbname.Data());
1143 Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1147 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1148 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1150 TRefArray *trackList = jet->GetRefTracks();
1151 Int_t nTracks = trackList->GetEntriesFast();
1152 nInputTracks += nTracks;
1153 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1155 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1157 return nInputTracks;
1162 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1164 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1165 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1166 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1167 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1168 double dphi = mphi-vphi;
1169 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1170 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1171 return dphi;//dphi in [-Pi, Pi]
1174 Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1177 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1178 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1179 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1180 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1187 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1189 // generate new THnSparseF, axes are defined in GetDimParams()
1192 UInt_t tmp = entries;
1195 tmp = tmp &~ -tmp; // clear lowest bit
1198 TString hnTitle(name);
1199 const Int_t dim = count;
1206 while(c<dim && i<32){
1210 GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1211 hnTitle += Form(";%s",label.Data());
1219 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1222 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1224 // stores label and binning of axis for THnSparse
1226 const Double_t pi = TMath::Pi();
1231 label = "V0 centrality (%)";
1240 label = "corrected jet pt";
1283 label = "leading track";
1291 label = "trigger track";