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