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