]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/UserTasks/AliAnalysisTaskJetAntenna.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add extra binning for embedding index
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskJetAntenna.cxx
CommitLineData
e59c88bd 1// ******************************************
2// This task computes several jet observables like
3// the fraction of energy in inner and outer coronnas,
4// jet-track correlations,triggered jet shapes and
5// correlation strength distribution of particles inside jets.
6// Author: lcunquei@cern.ch
7// *******************************************
8
9
10/**************************************************************************
11 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
12 * *
13 * Author: The ALICE Off-line Project. *
14 * Contributors are mentioned in the code where appropriate. *
15 * *
16 * Permission to use, copy, modify and distribute this software and its *
17 * documentation strictly for non-commercial purposes is hereby granted *
18 * without fee, provided that the above copyright notice appears in all *
19 * copies and that both the copyright notice and this permission notice *
20 * appear in the supporting documentation. The authors make no claims *
21 * about the suitability of this software for any purpose. It is *
22 * provided "as is" without express or implied warranty. *
23 **************************************************************************/
24
25
26#include "TChain.h"
27#include "TTree.h"
28#include "TMath.h"
29#include "TH1F.h"
30#include "TH2F.h"
31#include "TH3F.h"
32#include "THnSparse.h"
33#include "TMatrixD.h"
34#include "TMatrixDSym.h"
35#include "TMatrixDSymEigen.h"
36#include "TCanvas.h"
37#include "TRandom3.h"
38#include "AliLog.h"
39
40#include "AliAnalysisTask.h"
41#include "AliAnalysisManager.h"
42
43#include "AliVEvent.h"
44#include "AliESDEvent.h"
45#include "AliESDInputHandler.h"
46#include "AliCentrality.h"
47#include "AliAnalysisHelperJetTasks.h"
48#include "AliInputEventHandler.h"
49#include "AliAODJetEventBackground.h"
50#include "AliAODMCParticle.h"
51#include "AliAnalysisTaskFastEmbedding.h"
52#include "AliAODEvent.h"
53#include "AliAODHandler.h"
54#include "AliAODJet.h"
55
56#include "AliAnalysisTaskJetAntenna.h"
57
58using std::cout;
59using std::endl;
60
61ClassImp(AliAnalysisTaskJetAntenna)
62
63AliAnalysisTaskJetAntenna::AliAnalysisTaskJetAntenna() :
64AliAnalysisTaskSE(),
65fESD(0x0),
66fAODIn(0x0),
67fAODOut(0x0),
68fAODExtension(0x0),
69fBackgroundBranch(""),
70fNonStdFile(""),
71fIsPbPb(kTRUE),
72fOfflineTrgMask(AliVEvent::kAny),
73fMinContribVtx(1),
74fVtxZMin(-10.),
75fVtxZMax(10.),
76fEvtClassMin(0),
77fEvtClassMax(4),
78fFilterMask(0),
79fFilterMaskBestPt(0),
80fFilterType(0),
81fCentMin(0.),
82fCentMax(100.),
83fNInputTracksMin(0),
84fNInputTracksMax(-1),
85fRequireITSRefit(0),
86fApplySharedClusterCut(0),
87fTrackTypeRec(kTrackUndef),
88fRPAngle(0),
89fNRPBins(50),
90fSemigoodCorrect(0),
54aaacb8 91fDoMatching(kFALSE),
e59c88bd 92fHolePos(4.71),
93fHoleWidth(0.2),
94fCutTM(0.15),
95fJetEtaMin(-.5),
96fJetEtaMax(.5),
97fNevents(0),
98fTindex(0),
99fJetPtMin(20.),
100fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
101fJetPtFractionMin(0.5),
102fNMatchJets(4),
103fMatchMaxDist(0.8),
104fKeepJets(kFALSE),
105fkNbranches(2),
106fkEvtClasses(12),
107fOutputList(0x0),
108fHistEvtSelection(0x0),
03c25974 109fh2JetEntries(0x0),
e59c88bd 110fh2Circularity(0x0),
9a0db301 111fh2JetAxisPhi(0x0),
e59c88bd 112fhnJetTM(0x0)
113{
114 // default Constructor
115
116 fJetBranchName[0] = "";
117 fJetBranchName[1] = "";
118
119 fListJets[0] = new TList;
120 fListJets[1] = new TList;
121}
122
123AliAnalysisTaskJetAntenna::AliAnalysisTaskJetAntenna(const char *name) :
124AliAnalysisTaskSE(name),
125fESD(0x0),
126fAODIn(0x0),
127fAODOut(0x0),
128fAODExtension(0x0),
129fBackgroundBranch(""),
130fNonStdFile(""),
131fIsPbPb(kTRUE),
132fOfflineTrgMask(AliVEvent::kAny),
133fMinContribVtx(1),
134fVtxZMin(-10.),
135fVtxZMax(10.),
136fEvtClassMin(0),
137fEvtClassMax(4),
138fFilterMask(0),
139fFilterMaskBestPt(0),
140fFilterType(0),
141fCentMin(0.),
142fCentMax(100.),
143fNInputTracksMin(0),
144fNInputTracksMax(-1),
145fRequireITSRefit(0),
146fApplySharedClusterCut(0),
147fTrackTypeRec(kTrackUndef),
148fRPAngle(0),
149fNRPBins(50),
150fSemigoodCorrect(0),
54aaacb8 151fDoMatching(kFALSE),
e59c88bd 152fHolePos(4.71),
153fHoleWidth(0.2),
154fCutTM(0.15),
155fJetEtaMin(-.5),
156fJetEtaMax(.5),
157fNevents(0),
158fTindex(0),
159fJetPtMin(20.),
160fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
161fJetPtFractionMin(0.5),
162fNMatchJets(4),
163fMatchMaxDist(0.8),
164fKeepJets(kFALSE),
165fkNbranches(2),
166fkEvtClasses(12),
167fOutputList(0x0),
168fHistEvtSelection(0x0),
03c25974 169fh2JetEntries(0x0),
e59c88bd 170fh2Circularity(0x0),
9a0db301 171fh2JetAxisPhi(0x0),
e59c88bd 172fhnJetTM(0x0)
173 {
174 // Constructor
175
176
177 fJetBranchName[0] = "";
178 fJetBranchName[1] = "";
179
180 fListJets[0] = new TList;
181 fListJets[1] = new TList;
182
183 DefineOutput(1, TList::Class());
184}
185
186AliAnalysisTaskJetAntenna::~AliAnalysisTaskJetAntenna()
187{
188 delete fListJets[0];
189 delete fListJets[1];
190}
191
192void AliAnalysisTaskJetAntenna::SetBranchNames(const TString &branch1, const TString &branch2)
193{
194 fJetBranchName[0] = branch1;
195 fJetBranchName[1] = branch2;
196}
197
198void AliAnalysisTaskJetAntenna::Init()
199{
200
201 // check for jet branches
202 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
203 AliError("Jet branch name not set.");
204 }
205
206}
207
208void AliAnalysisTaskJetAntenna::UserCreateOutputObjects()
209{
210 // Create histograms
211 // Called once
212 OpenFile(1);
213 if(!fOutputList) fOutputList = new TList;
214 fOutputList->SetOwner(kTRUE);
215
216 Bool_t oldStatus = TH1::AddDirectoryStatus();
217 TH1::AddDirectory(kFALSE);
218
219 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
220 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
221 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
222 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
223 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
224 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
225 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
226 fOutputList->Add(fHistEvtSelection);
03c25974 227 fh2JetEntries=new TH2F("JetEntries","",150,0,150,10,-0.5,9.5);
228 fOutputList->Add(fh2JetEntries);
e59c88bd 229 fh2Circularity=new TH2F("Circcularity","",10,0,1,150,0,150);
230 fOutputList->Add(fh2Circularity);
551abe2f 231 fh2JetAxisPhi=new TH2F("JetAxisSmearPhi","",9,0,TMath::Pi(),10,-0.5,9.5);
232 fOutputList->Add(fh2JetAxisPhi);
233
234
9a0db301 235 Int_t nbinsJet[10]={3,9,7,9,36,10,2,10,20,2};
236 Double_t binlowJet[10]= {0,0, 0, 0,-0.5*TMath::Pi(),0,0,-0.5,0,0};
237 Double_t binupJet[10]= {100,0.9, 150,150,1.5*TMath::Pi(),1,200,9.5,200,2};
238 fhnJetTM = new THnSparseF("fhnJetTM", "fhnJetTM; cent;dr;pt_jet;pt_track;phi;circ;nc;pthard",10,nbinsJet,binlowJet,binupJet);
551abe2f 239
240 Double_t *xPt3=new Double_t[10];
e59c88bd 241 xPt3[0] = 0.;
242 for(Int_t i = 1;i<=9;i++){
243 if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.4; // 1 - 5
244 else if(xPt3[i-1]<11)xPt3[i] = xPt3[i-1] + 3; // 5 - 12
245 else xPt3[i] = xPt3[i-1] + 150.; // 18
246 }
86f50c7c 247 fhnJetTM->SetBinEdges(3,xPt3);
248
249 Double_t *xPt2=new Double_t[10];
250 xPt2[0] = 0.;
251 xPt2[1]=20;
252 xPt2[2]=40;
253 xPt2[3]=60;
254 xPt2[4]=80;
255 xPt2[5]=100;
256 xPt2[6]=120;
257 xPt2[7]=150;
258
259 fhnJetTM->SetBinEdges(2,xPt2);
551abe2f 260
261
262 Double_t *xPt4=new Double_t[4];
263 xPt4[0] = 0.;
264 xPt4[1]=10;
265 xPt4[2]=30;
266 xPt4[3]=50;
267
268
269 fhnJetTM->SetBinEdges(0,xPt4);
270
271
272 fOutputList->Add(fhnJetTM);
86f50c7c 273 delete [] xPt3;
274 delete [] xPt2;
551abe2f 275 delete [] xPt4;
e59c88bd 276 // =========== Switch on Sumw2 for all histos ===========
277 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
278 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
279 if (h1){
280 h1->Sumw2();
281 continue;
282 }
283 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
284 if (hn){
285 hn->Sumw2();
286 }
287 }
288 TH1::AddDirectory(oldStatus);
289
290 PostData(1, fOutputList);
291}
292
293void AliAnalysisTaskJetAntenna::UserExec(Option_t *)
294{
295
296
297 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
298 AliError("Jet branch name not set.");
299 return;
300 }
301
302 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
303 if (!fESD) {
304 AliError("ESD not available");
305 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
306 }
307 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
308
309 static AliAODEvent* aod = 0;
310 // take all other information from the aod we take the tracks from
311 if(!aod){
312 if(!fESD)aod = fAODIn;
313 else aod = fAODOut;}
314
315
316
317 if(fNonStdFile.Length()!=0){
318 // case that we have an AOD extension we need can fetch the jets from the extended output
319 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
320 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
321 if(!fAODExtension){
322 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
323 }
324 }
325
326
327 // -- event selection --
328 fHistEvtSelection->Fill(1); // number of events before event selection
329
330
331 Bool_t selected=kTRUE;
332 selected = AliAnalysisHelperJetTasks::Selected();
333 if(!selected){
334 // no selection by the service task, we continue
335 PostData(1,fOutputList);
336 return;}
337
338
339
340 // physics selection: this is now redundant, all should appear as accepted after service task selection
341 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
342 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
343 std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
344 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
345 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
346 fHistEvtSelection->Fill(2);
347 PostData(1, fOutputList);
348 return;
349 }
350
351
352
353 // vertex selection
354 if(!aod){
355 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
356 fHistEvtSelection->Fill(3);
357 PostData(1, fOutputList);
358 }
359 AliAODVertex* primVtx = aod->GetPrimaryVertex();
360
361 if(!primVtx){
362 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
363 fHistEvtSelection->Fill(3);
364 PostData(1, fOutputList);
365 return;
366 }
367
368 Int_t nTracksPrim = primVtx->GetNContributors();
369 if ((nTracksPrim < fMinContribVtx) ||
370 (primVtx->GetZ() < fVtxZMin) ||
371 (primVtx->GetZ() > fVtxZMax) ){
372 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
373 fHistEvtSelection->Fill(3);
374 PostData(1, fOutputList);
375 return;
376 }
377
378
379
380 // centrality selection
381 AliCentrality *cent = 0x0;
382 Double_t centValue = 0.;
383 if(fIsPbPb){
384 if(fESD) {cent = fESD->GetCentrality();
385 if(cent) centValue = cent->GetCentralityPercentile("V0M");}
386 else centValue=aod->GetHeader()->GetCentrality();
387
388 if(fDebug) printf("centrality: %f\n", centValue);
389 if (centValue < fCentMin || centValue > fCentMax){
390 fHistEvtSelection->Fill(4);
391 PostData(1, fOutputList);
392 return;
393 }}
394
395
396 fHistEvtSelection->Fill(0);
397 // accepted events
398 // -- end event selection --
86f50c7c 399 // pt hard
400 Double_t pthard=0;
401 Double_t pthardbin=0;
402 if(fDoMatching){
403 pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
404 pthardbin = GetPtHardBin(pthard);}
e59c88bd 405 // get background
406 AliAODJetEventBackground* externalBackground = 0;
407 if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
408 externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
409 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
410 }
411 if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
412 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
413 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
414 }
415
416 if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
417 externalBackground = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
418 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
419 }
420
421 Float_t rho = 0;
422
423
424 if(externalBackground)rho = externalBackground->GetBackground(0);
425
426 // fetch jets
427 TClonesArray *aodJets[2];
428 aodJets[0]=0;
429 if(fAODOut&&!aodJets[0]){
430 aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
431 aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));
432 }
433 if(fAODExtension && !aodJets[0]){
434 aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
435 aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
436 }
437 if(fAODIn&&!aodJets[0]){
438 aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
439 aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data()));
440 }
441
442
443
444 Int_t nT=0;
445 TList ParticleList;
54aaacb8 446 if(!fDoMatching) nT = GetListOfTracks(&ParticleList);
447 if(fDoMatching)nT=GetListOfTracksExtra(&ParticleList);
e59c88bd 448 if(nT<0){
449 PostData(1, fOutputList);
450 return;
451 }
452
453 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
454 fListJets[iJetType]->Clear();
455 if (!aodJets[iJetType]) continue;
456 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
457 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
458 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
459 if (jet) fListJets[iJetType]->Add(jet);
460 }
461 }
462
54aaacb8 463// jet matching
464 static TArrayI aMatchIndex(fListJets[0]->GetEntries());
465 static TArrayF aPtFraction(fListJets[0]->GetEntries());
466 if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
467 if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
e59c88bd 468
54aaacb8 469 if(fDoMatching){
470
e59c88bd 471
54aaacb8 472 // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
473 AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
474 fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
475 aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);}
476
477
478
479
480
481
482 // loop over matched jets
483 Int_t ir = -1; // index of matched reconstruced jet
484 Float_t fraction = -1.;
485 AliAODJet *jetmatched =0x0;
486
487
488
489
490
491 for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
e59c88bd 492
493 Double_t etabig=0;
494 Double_t ptbig=0;
495 Double_t areabig=0;
496 Double_t phibig=0.;
497 Double_t pxbig,pybig,pzbig;
551abe2f 498 Double_t phitrue=0;
499 Double_t etatrue=0;
500 Double_t smearphi=0;
e59c88bd 501
502 AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
503 etabig = jetbig->Eta();
504 phibig = jetbig->Phi();
551abe2f 505 phitrue=phibig;
506 etatrue=etabig;
e59c88bd 507 ptbig = jetbig->Pt();
508 if(ptbig==0) continue;
509 areabig = jetbig->EffectiveAreaCharged();
54aaacb8 510 if(!fDoMatching) ptbig=ptbig-rho*areabig;
e59c88bd 511 if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
512
513 if(fSemigoodCorrect){
514 if((phibig>fHolePos-fHoleWidth) && (phibig<fHolePos+fHoleWidth)) continue;
515 }
e59c88bd 516 pxbig=jetbig->Px();
517 pybig=jetbig->Py();
518 pzbig=jetbig->Pz();
54aaacb8 519
520 if(fDoMatching){
521 Bool_t jetAccepted=kTRUE;
522 ir = aMatchIndex[i];
523 if(ir>=0) jetmatched = (AliAODJet*)(fListJets[1]->At(ir));
524 else continue;
525 fraction = aPtFraction[i];
526 // minimum fraction required
527
528 if(fraction<fJetPtFractionMin) jetAccepted = kFALSE;
529 if(jetAccepted){
530 // jet acceptance + minimum pT check
531 if(jetmatched->Eta()>fJetEtaMax || jetmatched->Eta()<fJetEtaMin) jetAccepted = kFALSE;
532
533 }
534 if(!jetAccepted) continue;
551abe2f 535 etabig=jetmatched->Eta();
536 phibig=jetmatched->Phi();
54aaacb8 537 pxbig=jetmatched->Px();
538 pybig=jetmatched->Py();
539 pzbig=jetmatched->Pz();
551abe2f 540 ptbig=jetmatched->Pt()-rho*jetmatched->EffectiveAreaCharged();
541
542
543
544
545
546 smearphi=RelativePhi(phitrue,jetmatched->Phi());
547 smearphi=TMath::Abs(smearphi);
548
549
550}
54aaacb8 551
552
553
554
e59c88bd 555 TVector3 ppJ1(pxbig, pybig, pzbig);
556 TVector3 ppJ3(- pxbig * pzbig, - pybig * pzbig, pxbig * pxbig + pybig * pybig);
557 ppJ3.SetMag(1.);
558 TVector3 ppJ2(-pybig, pxbig, 0);
559 ppJ2.SetMag(1.);
560
54aaacb8 561
562
563
e59c88bd 564 Float_t mxx = 0.;
565 Float_t myy = 0.;
566 Float_t mxy = 0.;
567 Int_t nc = 0;
568 Float_t sump2 = 0.;
03c25974 569 Float_t ptmax=-10;
54aaacb8 570 for(int it = 0;it<ParticleList.GetEntries();++it){
e59c88bd 571 AliVParticle *track = (AliVParticle*)ParticleList.At(it);
572 TVector3 pp(track->Px(), track->Py(), track->Pz());
573 Float_t phi = track->Phi();
574 Float_t eta = track->Eta();
575 Float_t pt = track->Pt();
576 Float_t deta = eta - etabig;
577 Float_t dphi = RelativePhi(phi,phibig);
578 if(TMath::Abs(dphi)>=0.5*TMath::Pi()) continue;
579 Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
03c25974 580
581 if (r < 0.4 && pt>fCutTM) {
e59c88bd 582 //longitudinal and perpendicular component of the track pT in the
583 //local frame
03c25974 584 if(pt>ptmax) ptmax=pt;
e59c88bd 585 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
586 TVector3 pPerp = pp - pLong;
587 //projection onto the two perpendicular vectors defined above
588 Float_t ppjX = pPerp.Dot(ppJ2);
589 Float_t ppjY = pPerp.Dot(ppJ3);
590 Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
591 //components of the 2D symmetrical sphericity matrix
592 mxx += (ppjX * ppjX / ppjT);
593 myy += (ppjY * ppjY / ppjT);
594 mxy += (ppjX * ppjY / ppjT);
595 nc++;
596 sump2 += ppjT;}
03c25974 597
e59c88bd 598 if(nc<2) continue;
599
600 } // 1st Track Loop
601
602
603 // Sphericity Matrix
604 const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
605 TMatrixDSym m0(2,ele);
606 // Find eigenvectors
607 TMatrixDSymEigen m(m0);
608 TVectorD eval(2);
609 TMatrixD evecm = m.GetEigenVectors();
610 eval = m.GetEigenValues();
611 // Largest eigenvector
612 Int_t jev = 0;
613 if (eval[0] < eval[1]) jev = 1;
614 TVectorD evec0(2);
615 // Principle axis
616 evec0 = TMatrixDColumn(evecm, jev);
617 TVector2 evec(evec0[0], evec0[1]);
618 Float_t circ=0;
619 if(jev==1) circ=2*eval[0];
620 if(jev==0) circ=2*eval[1];
621 fh2Circularity->Fill(circ,ptbig);
03c25974 622 fh2JetEntries->Fill(ptbig,pthardbin);
623
551abe2f 624
625 if(fDoMatching) fh2JetAxisPhi->Fill(smearphi,pthardbin);
626
627
03c25974 628
54aaacb8 629 for (Int_t ip = 0; ip < ParticleList.GetEntries(); ip++) {
e59c88bd 630 AliVParticle *track = (AliVParticle*)ParticleList.At(ip);
9a0db301 631 Float_t isembed=0.5;
632 if(fDoMatching) if(ip<nT) isembed=1.5;
e59c88bd 633 TVector3 pp(track->Px(), track->Py(), track->Pz());
634 Float_t phi = track->Phi();
635 Float_t eta = track->Eta();
636 Float_t pt = track->Pt();
03c25974 637
e59c88bd 638 Float_t deta = eta - etabig;
639 Float_t dphi = RelativePhi(phi,phibig);
640 if(TMath::Abs(dphi)>=0.5*TMath::Pi()) continue;
641
642 Float_t dRR = TMath::Sqrt(dphi * dphi + deta * deta);
643 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
644 TVector3 pPerp = pp - pLong;
645 Float_t ppjX = pPerp.Dot(ppJ2);
646 Float_t ppjY = pPerp.Dot(ppJ3);
647 TVector2 vr(ppjX, ppjY) ;
648 //and this is the angle between the particle and the TM axis.
649 float phistr=evec.Phi()-vr.Phi();
650 if(phistr>2*TMath::Pi()) phistr -= 2*TMath::Pi();
651 if(phistr<-2*TMath::Pi()) phistr += 2*TMath::Pi();
652 if(phistr<-0.5*TMath::Pi()) phistr += 2*TMath::Pi();
653 if(phistr>1.5*TMath::Pi()) phistr -= 2*TMath::Pi();
654
9a0db301 655 double jetEntries[10] = {centValue,dRR,ptbig,pt,phistr,circ,static_cast<double>(nc),pthardbin,ptmax,isembed};
e59c88bd 656 fhnJetTM->Fill(jetEntries);
657
658 } // 2nd Track loop
659 }//jet loop
660
661 PostData(1, fOutputList);
662}
663
664void AliAnalysisTaskJetAntenna::Terminate(const Option_t *)
665{
666 // Draw result to the screen
667 // Called once at the end of the query
668
669 if (!GetOutputData(1))
670 return;
671}
672
673
674Int_t AliAnalysisTaskJetAntenna::GetListOfTracks(TList *list){
675
676 Int_t iCount = 0;
677 AliAODEvent *aod = 0;
678
679 if(!fESD)aod = fAODIn;
680 else aod = fAODOut;
e59c88bd 681 if(!aod)return 0;
682
e59c88bd 683
54aaacb8 684 for(int it = 0;it < aod->GetNumberOfTracks();++it){
e59c88bd 685 AliAODTrack *tr = aod->GetTrack(it);
686 Bool_t bGood = false;
687 if(fFilterType == 0)bGood = true;
688 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
689 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
690 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
691 if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
692 if(bGood==false) continue;
693 if (fApplySharedClusterCut) {
694 Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
695 if (frac > 0.4) continue;
696 }
697 if(TMath::Abs(tr->Eta())>0.9)continue;
698 if(tr->Pt()<0.15)continue;
699 list->Add(tr);
700 iCount++;
e59c88bd 701 }
54aaacb8 702 return iCount;
e59c88bd 703}
704
705
54aaacb8 706
707
708Int_t AliAnalysisTaskJetAntenna::GetListOfTracksExtra(TList *list){
709
710 Int_t iCount = 0;
9a0db301 711 Int_t nEmbed=0;
54aaacb8 712 AliAODEvent *aod = 0;
713
714 if(!fESD)aod = fAODIn;
715 else aod = fAODOut;
716 if(!aod)return 0;
717
c7fdb504 718
9a0db301 719 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
54aaacb8 720 if(!aodExtraTracks)return iCount;
721 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
722 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
723 if (!track) continue;
724 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
725 if(!trackAOD)continue;
726 Bool_t bGood = false;
727 if(fFilterType == 0)bGood = true;
728 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
729 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
730 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
731 if(fRequireITSRefit){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
732 if (fApplySharedClusterCut) {
733 Double_t frac = Double_t(trackAOD->GetTPCnclsS()) /Double_t(trackAOD->GetTPCncls());
734 if (frac > 0.4) continue;
735 }
736
54aaacb8 737 if(TMath::Abs(trackAOD->Eta())>0.9) continue;
738 if(trackAOD->Pt()<0.15) continue;
739 if(fDebug) printf("pt extra track %.2f \n", trackAOD->Pt());
740 list->Add(trackAOD);
741 iCount++;
9a0db301 742 }
743
744 nEmbed=iCount-1;
745
746
747
748
749 for(int it = 0;it < aod->GetNumberOfTracks();++it){
750 AliAODTrack *tr = aod->GetTrack(it);
751 Bool_t bGood = false;
752 if(fFilterType == 0)bGood = true;
753 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
754 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
755 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
756 if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
757 if(bGood==false) continue;
758 if (fApplySharedClusterCut) {
759 Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
760 if (frac > 0.4) continue;
761 }
762 if(TMath::Abs(tr->Eta())>0.9)continue;
763 if(tr->Pt()<0.15)continue;
764 list->Add(tr);
765 iCount++;
766 }
767
768 return nEmbed;
54aaacb8 769}
770
e59c88bd 771Double_t AliAnalysisTaskJetAntenna::RelativePhi(Double_t mphi,Double_t vphi){
772
773 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
774 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
775 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
776 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
777 double dphi = mphi-vphi;
778 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
779 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
780 return dphi;//dphi in [-Pi, Pi]
781}
782
783Int_t AliAnalysisTaskJetAntenna::GetPhiBin(Double_t phi)
784{
785 Int_t phibin=-1;
786 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
787 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
788 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
789 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
790 return phibin;
791}
86f50c7c 792
793Int_t AliAnalysisTaskJetAntenna::GetPtHardBin(Double_t ptHard){
794
795 const Int_t nBins = 10;
796 Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
797
798 Int_t bin = -1;
799 while(bin<nBins-1 && binLimits[bin+1]<ptHard){
800 bin++;
801 }
802
803 return bin;
804}