update to AOD analysis
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskAj.cxx
CommitLineData
6134f035 1
2// ******************************************
3// This task searches for events with large dijet imbalance
4// and then looks to the jet structure of the b-t-b jets.
5// *******************************************
6
7
8/**************************************************************************
9 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10 * *
11 * Author: The ALICE Off-line Project. *
12 * Contributors are mentioned in the code where appropriate. *
13 * *
14 * Permission to use, copy, modify and distribute this software and its *
15 * documentation strictly for non-commercial purposes is hereby granted *
16 * without fee, provided that the above copyright notice appears in all *
17 * copies and that both the copyright notice and this permission notice *
18 * appear in the supporting documentation. The authors make no claims *
19 * about the suitability of this software for any purpose. It is *
20 * provided "as is" without express or implied warranty. *
21 **************************************************************************/
22
23
24#include "TChain.h"
25#include "TTree.h"
26#include "TMath.h"
27#include "TH1F.h"
28#include "TH2F.h"
29#include "TH3F.h"
30#include "THnSparse.h"
31#include "TMatrixD.h"
32#include "TMatrixDSym.h"
33#include "TMatrixDSymEigen.h"
34
35#include "TCanvas.h"
36
37#include "AliLog.h"
38
39#include "AliAnalysisTask.h"
40#include "AliAnalysisManager.h"
41
42#include "AliVEvent.h"
43#include "AliESDEvent.h"
44#include "AliESDInputHandler.h"
45#include "AliCentrality.h"
46#include "AliAnalysisHelperJetTasks.h"
47#include "AliInputEventHandler.h"
48#include "AliAODJetEventBackground.h"
49#include "AliAnalysisTaskFastEmbedding.h"
50#include "AliAODEvent.h"
51#include "AliAODHandler.h"
52#include "AliAODJet.h"
53
54#include "AliAnalysisTaskAj.h"
55
56ClassImp(AliAnalysisTaskAj)
57
58AliAnalysisTaskAj::AliAnalysisTaskAj() :
59AliAnalysisTaskSE(),
60fESD(0x0),
61fAOD(0x0),
62fAODExtension(0x0),
63fBackgroundBranch(""),
64fNonStdFile(""),
65fIsPbPb(kTRUE),
66fOfflineTrgMask(AliVEvent::kAny),
67fMinContribVtx(1),
68fVtxZMin(-10.),
69fVtxZMax(10.),
70fEvtClassMin(0),
71fEvtClassMax(4),
72fFilterMask(0),
73fRadioFrac(0.2),
74fMinDist(0.1),
75fCentMin(0.),
76fCentMax(100.),
77fNInputTracksMin(0),
78fNInputTracksMax(-1),
79fAngStructCloseTracks(0),
80fCheckMethods(0),
81fDoEventMixing(0),
82fJetEtaMin(-.5),
83fJetEtaMax(.5),
84fNevents(0x0),
85fTindex(0x0),
86fJetPtMin(20.),
87fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
88fJetPtFractionMin(0.5),
89fNMatchJets(4),
90fMatchMaxDist(0.8),
91fKeepJets(kFALSE),
92fkNbranches(2),
93fkEvtClasses(12),
94fOutputList(0x0),
95fbEvent(kTRUE),
96fHistEvtSelection(0x0),
97fh2Pt1Pt2C10(0x0),
98fh2Pt1Pt2C40(0x0),
99fh3LocalCoordinates(0x0),
100fh2Sum2pt20(0x0),
101fh2Sum4pt20(0x0),
102fh2Sum6pt20(0x0),
103fh2Sum8pt20(0x0),
104fh2Sum12pt20(0x0),
105fh2Sum2lpt20(0x0),
106fh2Sum4lpt20(0x0),
107fh2Sum6lpt20(0x0),
108fh2Sum8lpt20(0x0),
109fh2Sum12lpt20(0x0),
110fh2Sum2pt40(0x0),
111fh2Sum4pt40(0x0),
112fh2Sum6pt40(0x0),
113fh2Sum8pt40(0x0),
114fh2Sum12pt40(0x0),
115fh2Sum2lpt40(0x0),
116fh2Sum4lpt40(0x0),
117fh2Sum6lpt40(0x0),
118fh2Sum8lpt40(0x0),
119fh2Sum12lpt40(0x0),
120fhnDeltaR(0x0)
121 {
122 // default Constructor
123
124
125 fJetBranchName[0] = "";
126 fJetBranchName[1] = "";
127
128 fListJets[0] = new TList;
129 fListJets[1] = new TList;
130}
131
132AliAnalysisTaskAj::AliAnalysisTaskAj(const char *name) :
133AliAnalysisTaskSE(name),
134fESD(0x0),
135fAOD(0x0),
136fAODExtension(0x0),
137fBackgroundBranch(""),
138fNonStdFile(""),
139fIsPbPb(kTRUE),
140fOfflineTrgMask(AliVEvent::kAny),
141fMinContribVtx(1),
142fVtxZMin(-10.),
143fVtxZMax(10.),
144fEvtClassMin(0),
145fEvtClassMax(4),
146fFilterMask(0),
147fRadioFrac(0.2),
148fMinDist(0.1),
149fCentMin(0.),
150fCentMax(100.),
151fNInputTracksMin(0),
152fNInputTracksMax(-1),
153fAngStructCloseTracks(0),
154fCheckMethods(0),
155fDoEventMixing(0),
156fJetEtaMin(-.5),
157fJetEtaMax(.5),
158fNevents(0x0),
159fTindex(0x0),
160fJetPtMin(20.),
161fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
162fJetPtFractionMin(0.5),
163fNMatchJets(4),
164fMatchMaxDist(0.8),
165fKeepJets(kFALSE),
166fkNbranches(2),
167fkEvtClasses(12),
168fOutputList(0x0),
169fbEvent(kTRUE),
170fHistEvtSelection(0x0),
171fh2Pt1Pt2C10(0x0),
172fh2Pt1Pt2C40(0x0),
173fh3LocalCoordinates(0x0),
174fh2Sum2pt20(0x0),
175fh2Sum4pt20(0x0),
176fh2Sum6pt20(0x0),
177fh2Sum8pt20(0x0),
178fh2Sum12pt20(0x0),
179fh2Sum2lpt20(0x0),
180fh2Sum4lpt20(0x0),
181fh2Sum6lpt20(0x0),
182fh2Sum8lpt20(0x0),
183fh2Sum12lpt20(0x0),
184fh2Sum2pt40(0x0),
185fh2Sum4pt40(0x0),
186fh2Sum6pt40(0x0),
187fh2Sum8pt40(0x0),
188fh2Sum12pt40(0x0),
189fh2Sum2lpt40(0x0),
190fh2Sum4lpt40(0x0),
191fh2Sum6lpt40(0x0),
192fh2Sum8lpt40(0x0),
193fh2Sum12lpt40(0x0),
194fhnDeltaR(0x0)
195 {
196 // Constructor
197
198
199
200
201
202
203
204 fJetBranchName[0] = "";
205 fJetBranchName[1] = "";
206
207 fListJets[0] = new TList;
208 fListJets[1] = new TList;
209
210 DefineOutput(1, TList::Class());
211}
212
213AliAnalysisTaskAj::~AliAnalysisTaskAj()
214{
215 delete fListJets[0];
216 delete fListJets[1];
217}
218
219void AliAnalysisTaskAj::SetBranchNames(const TString &branch1, const TString &branch2)
220{
221 fJetBranchName[0] = branch1;
222 fJetBranchName[1] = branch2;
223}
224
225void AliAnalysisTaskAj::Init()
226{
227
228 // check for jet branches
229 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
230 AliError("Jet branch name not set.");
231 }
232
233}
234
235void AliAnalysisTaskAj::UserCreateOutputObjects()
236{
237 // Create histograms
238 // Called once
239 OpenFile(1);
240 if(!fOutputList) fOutputList = new TList;
241 fOutputList->SetOwner(kTRUE);
242
243 Bool_t oldStatus = TH1::AddDirectoryStatus();
244 TH1::AddDirectory(kFALSE);
245
246
247 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
248 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
249 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
250 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
251 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
252 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
253 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
254
255 UInt_t entries = 0; // bit coded, see GetDimParams() below
256 entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7|1<<8;
257 fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
258
259
260 //change binning in centrality
261 Double_t *xPt0 = new Double_t[4];
262 xPt0[0] = 0.;
263 for(int i = 1; i<=3;i++){
264 if(xPt0[i-1]<10)xPt0[i] = xPt0[i-1] + 10; // 1 - 5
265 else if(xPt0[i-1]<40)xPt0[i] = xPt0[i-1] + 30; // 5 - 12
266 else xPt0[i] = xPt0[i-1] + 60.; // 18
267 }
268 fhnDeltaR->SetBinEdges(0,xPt0);
269 delete [] xPt0;
270
271
272
273 //change binning in jTtrack
274 Double_t *xPt3 = new Double_t[3];
275 xPt3[0] = 0.;
276 for(int i = 1; i<=2;i++){
277 if(xPt3[i-1]<1)xPt3[i] = xPt3[i-1] + 1.; // 1 - 5
278 else xPt3[i] = xPt3[i-1] + 149.; // 18
279 }
280 fhnDeltaR->SetBinEdges(3,xPt3);
281 delete [] xPt3;
282
283
284
285 //change binning in pTtrack
286 Double_t *xPt4 = new Double_t[5];
287 xPt4[0] = 0.;
288 for(int i = 1; i<=4;i++){
289 if(xPt4[i-1]<0.4)xPt4[i] = xPt4[i-1] + 0.4; // 1 - 5
290 else if(xPt4[i-1]<3)xPt4[i] = xPt4[i-1] + 2.6; // 5 - 12
291 else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 7.4; // 5 - 12
292 else xPt4[i] = xPt4[i-1] + 150.; // 18
293 }
294 fhnDeltaR->SetBinEdges(4,xPt4);
295 delete [] xPt4;
296
297
298
299
300
301
302 //change binning in HTI
303 Double_t *xPt5 = new Double_t[4];
304 xPt5[0] = 0.;
305 for(int i = 1; i<=3;i++){
306 if(xPt5[i-1]<10)xPt5[i] = xPt5[i-1] + 5; // 10 - 12
307 else xPt5[i] = xPt5[i-1] + 40.; // 13
308 }
309 fhnDeltaR->SetBinEdges(8,xPt5);
310 delete [] xPt5;
311
312
313
314 fh2Pt1Pt2C10 = new TH2F("Dijet spectra central","",20,0.,200.,20,0.,200.);
315 fh2Pt1Pt2C40 = new TH2F("Dijet spectra peripheral","",20,0.,200.,20,0.,200.);
316 fh3LocalCoordinates = new TH3F("Local coordinates","",10,-2,2,10,-2,2,10,0,100);
317 fh2Sum2pt20 = new TH2F("pL R<0.2 pt20","",10,0.,1.,100,0.,100.);
318 fh2Sum4pt20 = new TH2F("pL R<0.4 pt20","",10,0.,1.,100,0.,100.);
319 fh2Sum6pt20 = new TH2F("pL R<0.6 pt20","",10,0.,1.,100,0.,100.);
320 fh2Sum8pt20 = new TH2F("pL R<0.8 pt20","",10,0.,1.,100,0.,100.);
321 fh2Sum12pt20 = new TH2F("pL R<1.2 pt20","",10,0.,1.,100,0.,100.);
322 fh2Sum2lpt20 = new TH2F("pL R<0.2 low pt pt20","",10,0.,1.,100,0,100);
323 fh2Sum4lpt20 = new TH2F("pL R<0.4 low pt pt20","",10,0.,1.,100,0,100);
324 fh2Sum6lpt20 = new TH2F("pL R<0.6 low pt pt20","",10,0.,1.,100,0,100);
325 fh2Sum8lpt20 = new TH2F("pL R<0.8 low pt pt20","",10,0.,1.,100,0,100);
326 fh2Sum12lpt20 = new TH2F("pL R<1.2 low pt pt20","",10,0.,1.,100,0,100);
327 fh2Sum2pt40 = new TH2F("pL R<0.2 pt40","",10,0.,1.,100,0.,100.);
328 fh2Sum4pt40 = new TH2F("pL R<0.4 pt40","",10,0.,1.,100,0.,100.);
329 fh2Sum6pt40 = new TH2F("pL R<0.6 pt40","",10,0.,1.,100,0.,100.);
330 fh2Sum8pt40 = new TH2F("pL R<0.8 pt40","",10,0.,1.,100,0.,100.);
331 fh2Sum12pt40 = new TH2F("pL R<1.2 pt40","",10,0.,1.,100,0.,100.);
332 fh2Sum2lpt40 = new TH2F("pL R<0.2 low pt pt40","",10,0.,1.,100,0,100);
333 fh2Sum4lpt40 = new TH2F("pL R<0.4 low pt pt40","",10,0.,1.,100,0,100);
334 fh2Sum6lpt40 = new TH2F("pL R<0.6 low pt pt40","",10,0.,1.,100,0,100);
335 fh2Sum8lpt40 = new TH2F("pL R<0.8 low pt pt40","",10,0.,1.,100,0,100);
336 fh2Sum12lpt40 = new TH2F("pL R<1.2 low pt pt40","",10,0.,1.,100,0,100);
337
338
339 fOutputList->Add(fHistEvtSelection);
340 fOutputList->Add(fh2Pt1Pt2C10);
341 fOutputList->Add(fh2Pt1Pt2C40);
342 fOutputList->Add(fh3LocalCoordinates);
343
344 fOutputList->Add(fh2Sum2pt20);
345 fOutputList->Add(fh2Sum4pt20);
346 fOutputList->Add(fh2Sum6pt20);
347 fOutputList->Add(fh2Sum8pt20);
348 fOutputList->Add(fh2Sum12pt20);
349 fOutputList->Add(fh2Sum2lpt20);
350 fOutputList->Add(fh2Sum4lpt20);
351 fOutputList->Add(fh2Sum6lpt20);
352 fOutputList->Add(fh2Sum8lpt20);
353 fOutputList->Add(fh2Sum12lpt20);
354
355 fOutputList->Add(fh2Sum2pt40);
356 fOutputList->Add(fh2Sum4pt40);
357 fOutputList->Add(fh2Sum6pt40);
358 fOutputList->Add(fh2Sum8pt40);
359 fOutputList->Add(fh2Sum12pt40);
360 fOutputList->Add(fh2Sum2lpt40);
361 fOutputList->Add(fh2Sum4lpt40);
362 fOutputList->Add(fh2Sum6lpt40);
363 fOutputList->Add(fh2Sum8lpt40);
364 fOutputList->Add(fh2Sum12lpt40);
365
366 fOutputList->Add(fhnDeltaR);
367
368
369
370 // fOutputList->Add(fh3specbiased);
371
372 // =========== Switch on Sumw2 for all histos ===========
373 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
374 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
375 if (h1){
376 h1->Sumw2();
377 continue;
378 }
379 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
380 if (hn){
381 hn->Sumw2();
382 }
383 }
384 TH1::AddDirectory(oldStatus);
385
386 PostData(1, fOutputList);
387}
388
389void AliAnalysisTaskAj::UserExec(Option_t *)
390{
391
392
393 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
394 AliError("Jet branch name not set.");
395 return;
396 }
397
398 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
399 if (!fESD) {
400 AliError("ESD not available");
401 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
402 } else {
403 fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
404 }
405
406 if(fNonStdFile.Length()!=0){
407 // case that we have an AOD extension we need can fetch the jets from the extended output
408 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
409 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
410 if(!fAODExtension){
411 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
412 }}
413
414
415
416
417
418 // -- event selection --
419 fHistEvtSelection->Fill(1); // number of events before event selection
420
421 // physics selection
422 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
423 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
424 cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
425 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
426 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
427 fHistEvtSelection->Fill(2);
428 PostData(1, fOutputList);
429 return;
430 }
431
432 // vertex selection
433 if(!fAOD){
434 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
435 fHistEvtSelection->Fill(3);
436 PostData(1, fOutputList);
437 }
438 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
439
440 if(!primVtx){
441 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
442 fHistEvtSelection->Fill(3);
443 PostData(1, fOutputList);
444 return;
445 }
446
447 Int_t nTracksPrim = primVtx->GetNContributors();
448 if ((nTracksPrim < fMinContribVtx) ||
449 (primVtx->GetZ() < fVtxZMin) ||
450 (primVtx->GetZ() > fVtxZMax) ){
451 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
452 fHistEvtSelection->Fill(3);
453 PostData(1, fOutputList);
454 return;
455 }
456
457 // event class selection (from jet helper task)
458 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
459 if(fDebug) Printf("Event class %d", eventClass);
460 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
461 fHistEvtSelection->Fill(4);
462 PostData(1, fOutputList);
463 return;
464 }
465
466 // centrality selection
467 AliCentrality *cent = 0x0;
468 Double_t centValue = 0.;
469 if(fESD) {cent = fESD->GetCentrality();
470 if(cent) centValue = cent->GetCentralityPercentile("V0M");}
471 else centValue=fAOD->GetHeader()->GetCentrality();
472
473 if(fDebug) printf("centrality: %f\n", centValue);
474 if (centValue < fCentMin || centValue > fCentMax){
475 fHistEvtSelection->Fill(4);
476 PostData(1, fOutputList);
477 return;
478 }
479
480
481 fHistEvtSelection->Fill(0);
482 // accepted events
483 // -- end event selection --
484
485 // get background
486 AliAODJetEventBackground* externalBackground = 0;
487 if(fAOD&&!externalBackground&&fBackgroundBranch.Length()){
488 externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
489 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
490 }
491 if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
492 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
493 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
494 }
495
496 Float_t rho = 0;
497 if(externalBackground)rho = externalBackground->GetBackground(0);
498
499
500 // fetch jets
501 TClonesArray *aodJets[2];
502 aodJets[0]=0;
503 if(fAOD&&!aodJets[0]){
504 aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data()));
505 aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); }
506 if(fAODExtension && !aodJets[0]){
507 aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
508 aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
509
510 //Double_t ptsub[aodJets[0]->GetEntriesFast()];
511 //Int_t inord[aodJets[0]->GetEntriesFast()];
512 //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
513 // ptsub[n]=0;
514 // inord[n]=0;}
515
516 TList ParticleList;
517 Int_t nT = GetListOfTracks(&ParticleList);
518 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
519 fListJets[iJetType]->Clear();
520 if (!aodJets[iJetType]) continue;
521
522 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
523
524
525 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
526 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
527 if (jet) fListJets[iJetType]->Add(jet);
528 // if(iJetType==0){
529 // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
530 }}
531
532
533 Double_t eta1=0;
534 Int_t selec=-1;
535 Double_t ptmax=-10;
536 Double_t areaj=0;
537 Double_t phij=0;
538 Double_t etaj=0;
539 Double_t ptj=0;
540 Double_t ptcorrj=0;
541 for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
542 AliAODJet* jetj = (AliAODJet*)(fListJets[0]->At(i));
543 etaj = jetj->Eta();
544 phij = jetj->Phi();
545 ptj = jetj->Pt();
546 if(ptj==0) continue;
547 areaj = jetj->EffectiveAreaCharged();
548 ptcorrj=ptj-rho*areaj;
549 if(ptcorrj<=0) continue;
550 if((etaj<fJetEtaMin)||(eta1>fJetEtaMax)) continue;
551 if(ptcorrj>ptmax){ptmax=ptcorrj;
552 selec=i;}}
553 ///hardest jet selected
554 if(selec<0){PostData(1, fOutputList);
555 return;}
556 AliAODJet* jet1 = (AliAODJet*)(fListJets[0]->At(selec));
557 //What is the hardest constituent track?
558 AliAODTrack* leadtrack1;
559 Int_t ippt=0;
560 Double_t ppt=-10;
561 TRefArray *genTrackList = jet1->GetRefTracks();
562 Int_t nTracksGenJet = genTrackList->GetEntriesFast();
563 AliAODTrack* genTrack;
564 for(Int_t ir=0; ir<nTracksGenJet; ++ir){
565 genTrack = (AliAODTrack*)(genTrackList->At(ir));
566 if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
567 ippt=ir;}}
568 leadtrack1=(AliAODTrack*)(genTrackList->At(ippt));
569 //If it doesn't exist or if it is greater that 100 GeV, discard.
570 if(!leadtrack1) {PostData(1, fOutputList);
571 return;}
572 if(leadtrack1->Pt()>=100.){ PostData(1, fOutputList);
573 return;}
574
575 //Look to the back-to-back jet
576 Int_t btb=-1;
577 for(Int_t j=1;j<fListJets[0]->GetEntries();j++){
578 if(j==selec) continue;
579 AliAODJet* jetb = (AliAODJet*)(fListJets[0]->At(j));
580 etaj = jetb->Eta();
581 phij = jetb->Phi();
582 ptj = jetb->Pt();
583 if(ptj<=0) continue;
584 areaj = jetb->EffectiveAreaCharged();
585 ptcorrj=ptj-rho*areaj;
586 if(ptcorrj<=0) continue;
587 if((etaj<fJetEtaMin)||(etaj>fJetEtaMax)) continue;
588 Double_t dphij=RelativePhi(jetb->Phi(),jet1->Phi());
589 if(TMath::Abs(dphij)>TMath::Pi()-0.2) { btb=j;
590 break;}}
591
592 AliAODJet* jet2 = (AliAODJet*)(fListJets[0]->At(btb));
593 //the back-to-back jet is also identified
594
595 if(btb<0){PostData(1, fOutputList);
596 return;}
597
598 Double_t ptcorr1=jet1->Pt()-rho*jet1->EffectiveAreaCharged();
599 Double_t ptcorr2=jet2->Pt()-rho*jet2->EffectiveAreaCharged();
600
601 if(centValue<10.) fh2Pt1Pt2C10->Fill(ptcorr1,ptcorr2);
602 if(centValue>40.) fh2Pt1Pt2C40->Fill(ptcorr1,ptcorr2);
603
604
605 Double_t px2=jet2->Px();
606 Double_t py2=jet2->Py();
607 Double_t pz2=jet2->Pz();
608 Double_t phi2=jet2->Phi();
609 Double_t eta2=jet2->Eta();
610 //Once we have have a dijet event,look to the structure of the back-to-back jet:
611
612 TVector3 ppJ1(px2, py2, pz2);
613 TVector3 ppJ3(- px2 * pz2, - py2 * pz2, px2 * px2 + py2 * py2);
614 ppJ3.SetMag(1.);
615 TVector3 ppJ2(-py2, px2, 0);
616 ppJ2.SetMag(1.);
617 Float_t mxx = 0.;
618 Float_t myy = 0.;
619 Float_t mxy = 0.;
620 Int_t nc = 0;
621 Float_t sump2 = 0.;
622 Float_t ptMax = 0.;
623 Float_t etaMax = 0.;
624 Float_t phiMax = 0.;
625 Int_t iMax = -1;
626 //1st loop over all tracks
627 for(int it = 0;it<nT;++it){
628 AliVParticle *track = (AliVParticle*)ParticleList.At(it);
629 TVector3 pp(track->Px(), track->Py(), track->Pz());
630 Float_t phi = track->Phi();
631 Float_t eta = track->Eta();
632 Float_t pt = track->Pt();
633 Float_t jT = pp.Perp(ppJ1);
634 Float_t deta = eta - eta2;
635 Float_t dphi = phi - phi2;
636 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
637 if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
638 Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
639 /////////To compute the TM axis, we use particles with large jT
640
641 //cout<<"Jt spectrum............"<<ptbig<<" "<<jT<<endl;
642 /////////We compute the TM with large jT particles only
643 if((jT >1.)&&(r<1.2)){
644 //longitudinal and perpendicular component of the track pT in the
645 //local frame
646 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
647 TVector3 pPerp = pp - pLong;
648 //projection onto the two perpendicular vectors defined above
649 Float_t ppjX = pPerp.Dot(ppJ2);
650 Float_t ppjY = pPerp.Dot(ppJ3);
651 Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
652 //components of the sphericity matrix
653 mxx += (ppjX * ppjX / ppjT);
654 myy += (ppjY * ppjY / ppjT);
655 mxy += (ppjX * ppjY / ppjT);
656 nc++;
657 sump2 += ppjT;
658 // max pt
659 if (pt > ptMax) {
660 ptMax = pt;
661 iMax = it;
662 etaMax = deta;
663 phiMax = dphi;
664 }
665 } // R < 0.4
666 } // 1st Track Loop
667
668 // At this point we have mxx, myy, mxy
669 if (nc == 0) return;
670// Shericity Matrix
671 const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
672 TMatrixDSym m0(2,ele);
673// Find eigenvectors
674 TMatrixDSymEigen m(m0);
675 TVectorD eval(2);
676 TMatrixD evecm = m.GetEigenVectors();
677 eval = m.GetEigenValues();
678// Largest eigenvector
679 Int_t jev = 0;
680 if (eval[0] < eval[1]) jev = 1;
681 TVectorD evec0(2);
682// Principle axis
683 evec0 = TMatrixDColumn(evecm, jev);
684 TVector2 evec(evec0[0], evec0[1]);
685// Principle axis from leading partice
686 Float_t phiM = TMath::ATan2(phiMax, etaMax);
687 TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM));
688 Float_t phistM = evecM.DeltaPhi(evec);
689 if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
690
691 //////we have now the direction
692 /////along which the sum of the projections of the particle
693 ///momentum is higher.
694 Double_t sum2lpt20=0;
695 Double_t sum4lpt20=0;
696 Double_t sum6lpt20=0;
697 Double_t sum8lpt20=0;
698 Double_t sum12lpt20=0;
699 Double_t sum2pt20=0;
700 Double_t sum4pt20=0;
701 Double_t sum6pt20=0;
702 Double_t sum8pt20=0;
703 Double_t sum12pt20=0;
704
705 Double_t sum2lpt40=0;
706 Double_t sum4lpt40=0;
707 Double_t sum6lpt40=0;
708 Double_t sum8lpt40=0;
709 Double_t sum12lpt40=0;
710 Double_t sum2pt40=0;
711 Double_t sum4pt40=0;
712 Double_t sum6pt40=0;
713 Double_t sum8pt40=0;
714 Double_t sum12pt40=0;
715
716 for (Int_t ip = 0; ip < nT; ip++) {
717 AliVParticle *track = (AliVParticle*)ParticleList.At(ip);
718 TVector3 pp(track->Px(), track->Py(), track->Pz());
719 Float_t phi = track->Phi();
720 Float_t eta = track->Eta();
721 Float_t pt = track->Pt();
722 Float_t jT = pp.Perp(ppJ1);
723 Double_t deta = eta - eta2;
724 Double_t deltaPhi = phi - phi2;
725 Double_t r = TMath::Sqrt(deltaPhi * deltaPhi + deta * deta);
726
727 if(ptcorr2>20. && ptcorr2<40.){
728 if(pt<0.4){
729 if(r<0.2) sum2lpt20=sum2lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
730 if(r<0.4) sum4lpt20=sum4lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
731 if(r<0.6) sum6lpt20=sum6lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
732 if(r<0.8) sum8lpt20=sum8lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
733 if(r<1.2) sum12lpt20=sum12lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());}
734
735 if(r<0.2) sum2pt20=sum2pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
736 if(r<0.4) sum4pt20=sum4pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
737 if(r<0.6) sum6pt20=sum6pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
738 if(r<0.8) sum8pt20=sum8pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
739 if(r<1.2) sum12pt20=sum12pt20-1.*pt*TMath::Cos(phi-jet1->Phi());}
740
741
742 if(ptcorr2>40. && ptcorr2<60.){
743 if(pt<0.4){
744 if(r<0.2) sum2lpt40=sum2lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
745 if(r<0.4) sum4lpt40=sum4lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
746 if(r<0.6) sum6lpt40=sum6lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
747 if(r<0.8) sum8lpt40=sum8lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
748 if(r<1.2) sum12lpt40=sum12lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());}
749
750 if(r<0.2) sum2pt40=sum2pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
751 if(r<0.4) sum4pt40=sum4pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
752 if(r<0.6) sum6pt40=sum6pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
753 if(r<0.8) sum8pt40=sum8pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
754 if(r<1.2) sum12pt40=sum12pt40-1.*pt*TMath::Cos(phi-jet1->Phi());}
755
756
757
758
759 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
760 TVector3 pPerp = pp - pLong;
761 Float_t ppjX = pPerp.Dot(ppJ2);
762 Float_t ppjY = pPerp.Dot(ppJ3);
763 TVector2 vr(ppjX, ppjY) ;
764 //and this is the angle between the particle and the TM axis.
765 // Float_t phistr = evec.DeltaPhi(vr);
766
767 Double_t phistr=vr.Phi()-evec.Phi();
768
769 if(centValue<10.) fh3LocalCoordinates->Fill(ppjX,ppjY,ptcorr2);
770 Double_t deltaEta = eta2-track->Eta();
771 if(phistr<-0.5*TMath::Pi()) phistr+=2.*TMath::Pi();
772 if(phistr>3./2.*TMath::Pi()) phistr-=2.*TMath::Pi();
773
774 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
775 if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
776 Double_t jetEntries[9] = {centValue,ptcorr1,ptcorr2,jT,pt,deltaEta,deltaPhi,phistr,ptMax};
777 fhnDeltaR->Fill(jetEntries);
778 }
779 if(centValue<10.){
780 if(ptcorr2>20.){
781 if(ptcorr1>80.){
782 Double_t aj=(ptcorr1-ptcorr2)/(ptcorr1+ptcorr2);
783
784 fh2Sum2pt20->Fill(aj,sum2pt20);
785 fh2Sum4pt20->Fill(aj,sum4pt20);
786 fh2Sum6pt20->Fill(aj,sum6pt20);
787 fh2Sum8pt20->Fill(aj,sum8pt20);
788 fh2Sum12pt20->Fill(aj,sum12pt20);
789 fh2Sum2lpt20->Fill(aj,sum2lpt20);
790 fh2Sum4lpt20->Fill(aj,sum4lpt20);
791 fh2Sum6lpt20->Fill(aj,sum6lpt20);
792 fh2Sum8lpt20->Fill(aj,sum8lpt20);
793 fh2Sum12lpt20->Fill(aj,sum12lpt20);
794
795 fh2Sum2pt40->Fill(aj,sum2pt40);
796 fh2Sum4pt40->Fill(aj,sum4pt40);
797 fh2Sum6pt40->Fill(aj,sum6pt40);
798 fh2Sum8pt40->Fill(aj,sum8pt40);
799 fh2Sum12pt40->Fill(aj,sum12pt40);
800 fh2Sum2lpt40->Fill(aj,sum2lpt40);
801 fh2Sum4lpt40->Fill(aj,sum4lpt40);
802 fh2Sum6lpt40->Fill(aj,sum6lpt40);
803 fh2Sum8lpt40->Fill(aj,sum8lpt40);
804 fh2Sum12lpt40->Fill(aj,sum12lpt40);
805 }}}
806
807
808
809
810
811
812
813 PostData(1, fOutputList);
814}
815
816void AliAnalysisTaskAj::Terminate(const Option_t *)
817{
818 // Draw result to the screen
819 // Called once at the end of the query
820
821 if (!GetOutputData(1))
822 return;
823}
824
825
826
827
828
829
830
831
832
833
834
835Int_t AliAnalysisTaskAj::GetListOfTracks(TList *list){
836
837 Int_t iCount = 0;
4f9d90a5 838 if(!fAOD)return iCount;
6134f035 839
840 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
841 AliAODTrack *tr = fAOD->GetTrack(it);
842 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
843 if(TMath::Abs(tr->Eta())>0.9)continue;
844 if(tr->Pt()<0.15)continue;
845 list->Add(tr);
846 //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
847 iCount++;
848 }
849
850
851 return iCount;
852
853}
854
855 Int_t AliAnalysisTaskAj::GetHardestTrackBackToJet(AliAODJet *jetbig){
856
857
858 Int_t index=-1;
859 Double_t ptmax=-10;
860 Double_t dphi=0;
861 Double_t dif=0;
862 Int_t iCount=0;
863 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
864 AliAODTrack *tr = fAOD->GetTrack(it);
865 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
866 if(TMath::Abs(tr->Eta())>0.9)continue;
867 if(tr->Pt()<0.15)continue;
868 iCount=iCount+1;
869 dphi=RelativePhi(tr->Phi(),jetbig->Phi());
870 if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
871 if(tr->Pt()>ptmax){ ptmax=tr->Pt();
872 index=iCount-1;
873 dif=dphi; }}
874
875 return index;
876
877 }
878
879
880
881
882
883
884
885
886
887 Int_t AliAnalysisTaskAj::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
888
889 Int_t iCount = 0;
890
891
892 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
893 AliAODTrack *tr = fAOD->GetTrack(it);
894 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
895 if(TMath::Abs(tr->Eta())>0.9)continue;
896 if(tr->Pt()<0.15)continue;
897 Double_t disR=jetbig->DeltaR(tr);
898 if(disR>0.8) continue;
899 list->Add(tr);
900 //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
901 iCount++;
902 }
903
904 list->Sort();
905 return iCount;
906
907}
908
909
910
911
912
913
914
915
916
917
918
919Int_t AliAnalysisTaskAj::GetNInputTracks()
920{
921
922 Int_t nInputTracks = 0;
923
924 TString jbname(fJetBranchName[1]);
925 //needs complete event, use jets without background subtraction
926 for(Int_t i=1; i<=3; ++i){
927 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
928 }
929 // use only HI event
930 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
931 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
932
933 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
934 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
935 if(!tmpAODjets){
936 Printf("Jet branch %s not found", jbname.Data());
937 Printf("AliAnalysisTaskAj::GetNInputTracks FAILED");
938 return -1;
939 }
940
941 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
942 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
943 if(!jet) continue;
944 TRefArray *trackList = jet->GetRefTracks();
945 Int_t nTracks = trackList->GetEntriesFast();
946 nInputTracks += nTracks;
947 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
948 }
949 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
950
951 return nInputTracks;
952}
953
954
955
956Double_t AliAnalysisTaskAj::RelativePhi(Double_t mphi,Double_t vphi){
957
958 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
959 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
960 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
961 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
962 double dphi = mphi-vphi;
963 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
964 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
965 return dphi;//dphi in [-Pi, Pi]
966}
967
968
969
970THnSparse* AliAnalysisTaskAj::NewTHnSparseF(const char* name, UInt_t entries)
971{
972 // generate new THnSparseF, axes are defined in GetDimParams()
973
974 Int_t count = 0;
975 UInt_t tmp = entries;
976 while(tmp!=0){
977 count++;
978 tmp = tmp &~ -tmp; // clear lowest bit
979 }
980
981 TString hnTitle(name);
982 const Int_t dim = count;
983 Int_t nbins[dim];
984 Double_t xmin[dim];
985 Double_t xmax[dim];
986
987 Int_t i=0;
988 Int_t c=0;
989 while(c<dim && i<32){
990 if(entries&(1<<i)){
991
992 TString label("");
993 GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
994 hnTitle += Form(";%s",label.Data());
995 c++;
996 }
997
998 i++;
999 }
1000 hnTitle += ";";
1001
1002 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1003}
1004
1005void AliAnalysisTaskAj::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1006{
1007 // stores label and binning of axis for THnSparse
1008
1009 const Double_t pi = TMath::Pi();
1010
1011 switch(iEntry){
1012
1013 case 0:
1014 label = "V0 centrality (%)";
1015
1016 nbins = 3;
1017 xmin = 0.;
1018 xmax = 100.;
1019 break;
1020
1021
1022 case 1:
1023 label = "corrected jet pt1";
1024 nbins = 10;
1025 xmin = 0.;
1026 xmax = 200.;
1027 break;
1028
1029 case 2:
1030 label = "corrected jet pt2";
1031 nbins = 10;
1032 xmin = 0.;
1033 xmax = 200.;
1034 break;
1035
1036
1037
1038 case 3:
1039 label = "track jT";
1040
1041 nbins = 2;
1042 xmin = 0.;
1043 xmax = 150.;
1044 break;
1045
1046 case 4:
1047 label = "track pT";
1048
1049 nbins = 4;
1050 xmin = 0.;
1051 xmax = 150.;
1052 break;
1053
1054
1055 case 5:
1056 label = "deltaEta";
1057 nbins = 8;
1058 xmin = -1.6;
1059 xmax = 1.6;
1060 break;
1061
1062
1063 case 6:
1064 label = "deltaPhi";
1065 nbins = 60;
1066 xmin = -0.5*pi;
1067 xmax = 1.5*pi;
1068 break;
1069
1070 case 7:
1071 label = "deltaPhiTM";
1072 nbins = 60;
1073 xmin = -0.5*pi;
1074 xmax = 1.5*pi;
1075 break;
1076
1077
1078
1079 case 8:
1080 label = "leading track";
1081 nbins = 3;
1082 xmin = 0;
1083 xmax = 50;
1084 break;
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095 }
1096
1097}
1098