]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/UserTasks/AliAnalysisTaskJetShape.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskJetShape.cxx
CommitLineData
c87f1b5e 1// Class AliAnalysisTaskJetShape
2// Author martin.poghosyan@cern.ch
3//
4
5#include "TChain.h"
6#include "TTree.h"
7#include "TH1.h"
8#include "TH1F.h"
9#include "TH2F.h"
10#include "TH3F.h"
11#include "TCanvas.h"
12#include "TParticlePDG.h"
13#include "TParticle.h"
14#include "TRandom.h"
15#include "TDatabasePDG.h"
16#include <TPDGCode.h>
17#include "THnSparse.h"
18#include "TGraph.h"
19#include "TArrayF.h"
20#include "TArrayD.h"
21#include <TVector3.h>
22#include <TClonesArray.h>
23#include "AliLog.h"
24
25#include "AliCDBManager.h"
26#include "AliGeomManager.h"
27#include "AliCDBEntry.h"
28#include "AliCDBPath.h"
29#include "AliAlignObjParams.h"
30#include "AliAnalysisTaskSE.h"
31#include "AliAnalysisManager.h"
32#include "AliESDEvent.h"
33#include "AliESDInputHandler.h"
34#include "AliMultiplicity.h"
35#include "AliESDFMD.h"
36#include "AliESDVZERO.h"
37#include "AliESDInputHandler.h"
38#include "AliVEvent.h"
39#include "AliCentrality.h"
40#include "AliAnalysisHelperJetTasks.h"
41#include "AliInputEventHandler.h"
42#include "AliAODJetEventBackground.h"
43#include "AliAODEvent.h"
44#include "AliAODHandler.h"
45#include "AliAODJet.h"
46#include "AliAODEvent.h"
47#include "AliAODTrack.h"
48#include "AliAODMCParticle.h"
49#include "AliGenEventHeader.h"
50#include "AliHeader.h"
51#include "AliStack.h"
52#include "AliMCEvent.h"
53#include "AliMCEventHandler.h"
54#include "AliInputEventHandler.h"
55#include "AliLHCData.h"
56#include "AliGenPythiaEventHeader.h"
57#include "AliGenDPMjetEventHeader.h"
58#include "AliESDtrackCuts.h"
59#include "AliTriggerAnalysis.h"
60
61
62
63#include "AliAnalysisTaskJetShape.h"
64
65/////////////////////////////////////
66
67ClassImp(AliAnalysisTaskJetShape)
68
69
70 AliAnalysisTaskJetShape::AliAnalysisTaskJetShape(const char *name)
71: AliAnalysisTaskSE(name),
72 // fOutputTree(0),
73 fOutputList(0),
74 fESD(0x0),
75 fAODIn(0x0),
76 fAODOut(0x0),
77 fAODExtension(0x0),
78 fkMC(0),
79 fCMSE(0),
80 fRunNb(0),
81 fkIsPhysSel(0),
82 fNonStdFile(""),
83 fFilterMask(0),
84 fOfflineTrgMask(AliVEvent::kAny),
85 fCentMin(0.),
86 fCentMax(80.),
87 fEvtClassMin(0),
88 fEvtClassMax(4),
89 fPtJcorrMin(20),
90 fPtJBcorrMin(20),
91 fJpPtmin(1),
92 fJaPtmin(0.5),
93 fVtxMinContrib(1),
94 fVtxZMin(-10),
95 fVtxZMax(10),
96 fkIsPbPb(0),
97 fhPtJL(0),
98 fhAreaJL(0),
99 fanalJetSubStrHM(0),
100 fanalJetSubStrMCHM(0)
101{
102
103
104 fgkbinCent[0] = -0.099;
105 fgkbinCent[1] = 5;
106 fgkbinCent[2] = 10;
107 fgkbinCent[3] = 20;
108 fgkbinCent[4] = 40;
109 fgkbinCent[5] = 60;
110 fgkbinCent[6] = 80;
111 fgkbinCent[7] =100;
112
113 fBackgroundBranch[0] = "",
114 fBackgroundBranch[1] = "",
115 fJetBranchName[0] = "";
116 fJetBranchName[1] = "";
117
89605a8b 118 for(Int_t i=0; i<3; i++)
119 {
120 fhPtresVsPt[i]=0;
121 fhPhiresVsPhi[i]=0;
122 fhEtaresVsEta[i]=0;
123 fhDCAxy[i]=0;
124 fhDCAz[i]=0;
125 fhTrackPtEtaPhi[i]=0;
126 }
c87f1b5e 127 // fkIsPbPb = kFALSE;
128 // fDebug=0;
129
130
131
89605a8b 132 fanalJetSubStrHM = new AliAnalysisTaskJetShapeHM("rec", kFALSE);
133 fanalJetSubStrMCHM = new AliAnalysisTaskJetShapeHM("truth", kTRUE);
c87f1b5e 134
135 DefineOutput(1, TList::Class());
136
137 // DefineOutput(1, TTree::Class());
138}
139
140
141AliAnalysisTaskJetShape::~AliAnalysisTaskJetShape()
142{
143 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
144 printf("Deleteing output\n");
145
146 if(fOutputList){
147 delete fOutputList;
148 fOutputList = 0;
149 }
150
151 // if(fTriggerAnalysis)
152 // delete fTriggerAnalysis;
153
154 }
155}
156
157
158AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::AliAnalysisTaskJetShapeTool():TObject(),
159 fvecJ(0,0,0),
160 fvecB1(0,0,0),
161 fvecB2(0,0,0),
162 fRmax(0),
163 fPRecInRJ(0,0,0),
164 fList(0),
165 fEntries(0)
166{
167 fList = new TClonesArray("TVector3",10000);
168
169 TVector3 v(0,0,0);
170 SetVecJ(v);
171 fRmax = -0.5;
172 fR[0] =0.1;
173 fR[1] =0.2;
174 fR[2] =0.3;
175 fEntries=0;
176
177 for(Int_t i1=0; i1<fgkbinR; i1++)
178 {
179 for(Int_t i2=0; i2<2; i2++)
180 {
181 fListJ[i1][i2].Set(1000);
182 fListB1[i1][i2].Set(1000);
183 fListB2[i1][i2].Set(1000);
184 fListJc[i1][i2] = 0;
185 fListB1c[i1][i2] = 0;
186 fListB2c[i1][i2] = 0;
187 // fListJ[i1][i2]("TVector3",10000);
188 // fListB1[i1][i2] = TClonesArray("TVector3",10000);
189 // fListB2[i1][i2] = TClonesArray("TVector3",10000);
190
191 fListJ[i1][i2].Reset(-1);
192 fListB1[i1][i2].Reset(-1);
193 fListB2[i1][i2].Reset(-1);
194
195 fPRecJ[i1][i2].SetXYZ(0,0,0);
196
197 }
198 }
199
200 fPRecInRJ.SetXYZ(0,0,0);
201
202 fPtMinTr[0] = 2;
203 fPtMinTr[1] = 0.5;
204
205
206}
207
208AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::AliAnalysisTaskJetShapeTool(TClonesArray *list):TObject(),
209 fvecJ(0,0,0),
210 fvecB1(0,0,0),
211 fvecB2(0,0,0),
212 fRmax(0),
213 fPRecInRJ(0,0,0),
214 fList(list),
215 fEntries(0)
216{
217
218 TVector3 v(0,0,0);
219 SetVecJ(v);
220 fRmax = 0.4;
221 fR[0] =0.1;
222 fR[1] =0.2;
223 fR[2] =0.3;
224 fEntries=0;
225
226 for(Int_t i1=0; i1<fgkbinR; i1++)
227 {
228 for(Int_t i2=0; i2<2; i2++)
229 {
230 fListJ[i1][i2].Set(1000);
231 fListB1[i1][i2].Set(1000);
232 fListB2[i1][i2].Set(1000);
233 fListJc[i1][i2] = 0;
234 fListB1c[i1][i2] = 0;
235 fListB2c[i1][i2] = 0;
236 // fListJ[i1][i2]("TVector3",10000);
237 // fListB1[i1][i2] = TClonesArray("TVector3",10000);
238 // fListB2[i1][i2] = TClonesArray("TVector3",10000);
239 fListJ[i1][i2].Reset(-1);
240 fListB1[i1][i2].Reset(-1);
241 fListB2[i1][i2].Reset(-1);
242 fPRecJ[i1][i2].SetXYZ(0,0,0);
243
244 }
245 }
246
247 fPRecInRJ.SetXYZ(0,0,0);
248
249 fPtMinTr[0] = 2;
250 fPtMinTr[1] = 0.5;
251
252}
253
254
255
256
257AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::~AliAnalysisTaskJetShapeTool()
258{
259 fList->Delete();
260 fEntries=0;
261
262 // if(fList)
263 // delete fList;
264
265}
266
267
268void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::SetVecJ(TVector3 v)
269{
270 fvecJ.SetXYZ(v.X(), v.Y(), v.Z());
271 fvecB1.SetXYZ(-v.Y(), v.X(), v.Z());
272 fvecB2.SetXYZ(v.Y(),-v.X(), v.Z());
273}
274
275
276void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::Clean()
277{
278 // printf("AnalJetSubStrTool::Clean()\n");
279 /*
280 fList->Delete();
281 for(Int_t i1=0; i1<fgkbinR; i1++)
282 {
283 for(Int_t i2=0; i2<2; i2++)
284 {
285 fListJ[i1][i2]->Delete();
286 fListB1[i1][i2]->Delete();
287 fListB2[i1][i2]->Delete();
288 }
289 }
290
291 */
292
293 // fList.SetOwner();
294 // fList->Delete();
295 // fEntries=0;
296 fPRecInRJ.SetXYZ(0,0,0);
297
298 for(Int_t i1=0; i1<fgkbinR; i1++)
299 {
300 for(Int_t i2=0; i2<2; i2++)
301 {
302 fListJ[i1][i2].Reset(-1);
303 fListB1[i1][i2].Reset(-1);
304 fListB2[i1][i2].Reset(-1);
305 // fListJ[i1][i2].Set(0);
306 // fListB1[i1][i2].Set(0);
307 // fListB2[i1][i2].Set(0);
308 fListJc[i1][i2] = 0;
309 fListB1c[i1][i2] = 0;
310 fListB2c[i1][i2] = 0;
311
312 fPRecJ[i1][i2].SetXYZ(0,0,0);
313 }
314 }
315
316}
317
318
319//void AnalJetSubStrTool::Print(Option_t* /*option*/) const
320void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::PRINT() const
321{
322
323 // PRINT(fList, "all");
324
325 // fList->Print();
326 for(Int_t i1=0; i1<fgkbinR; i1++)
327 {
328
329 if(i1!=0) continue;
330 for(Int_t i2=0; i2<2; i2++)
331 {
332
333 // printf("# %d %d\n",i1, i2);
334 PRINT(fListJ[i1][i2], fListJc[i1][i2], Form("J_%d_%d",i1,i2));
335 // PRINT(fListB1[i1][i2], Form("B1_%d_%d",i1,i2));
336 // PRINT(fListB2[i1][i2], Form("B2_%d_%d",i1,i2));
337 // fListJ[i1][i2].Print("",1);
338 // fListB1[i1][i2]->Print();
339 // fListB2[i1][i2]->Print();
340 }
341 }
342
343}
344
345void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::PRINT(TArrayI a, Int_t n, const char *txt) const
346{
347 printf("%s :%d \n",txt, n);
348 Int_t count = 0;
349 for(Int_t i=0; i<n; i++){
350 printf("%4d ", a.At(i));
351
352 if(count==20)
353 {
354 printf("\n");
355 count=0;
356 }
357 else count++;
358 }
359 if(count!=20) printf("\n");
360
361}
362
363
364
365 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::Init()
366{
367 Int_t Nev = fEntries;
368
369 // PRINT();
370
371 // printf("Nev = %d\n",Nev);
372
373 for(Int_t iP=0; iP<Nev; iP++)
374 {
375 TVector3 *v = (TVector3*) fList->At(iP);
376 if(!v) {
377 printf("ERROR : entry #%d not found\n",iP);
378 continue;
379 }
380
381 // printf("#%3d ",iP); v->Print();
382
383 Double_t R = CalcR(fvecJ, *v);
384 // printf("R = %f\n",R);
385 if(R<fRmax)
386 {
387 for(Int_t iT = 0; iT < fgkbinR; iT++)
388 {
389 Int_t type = 0;
390 if(R>fR[iT]) type = 1;
391
392 if(v->Pt() < fPtMinTr[type]) continue;
393 fPRecJ[iT][type]+=*v;
394 AddToJ(iT, type, iP);
395 }
396
397 if(v->Pt() < fPtMinTr[0]) continue;
398 fPRecInRJ+=*v;
399
400 continue;
401 }
402
403 R = CalcR(fvecB1, *v);
404 if(R<fRmax)
405 {
406 for(Int_t iT = 0; iT < fgkbinR; iT++)
407 {
408 Int_t type = 0;
409 if(R>fR[iT]) type = 1;
410
411 if(v->Pt() < fPtMinTr[type]) continue;
412
413 AddToB1(iT, type, iP);
414 }
415 continue;
416 }
417
418 R = CalcR(fvecB2, *v);
419 if(R<fRmax)
420 {
421 for(Int_t iT = 0; iT < fgkbinR; iT++)
422 {
423 Int_t type = 0;
424 if(R>fR[iT]) type = 1;
425
426 if(v->Pt() < fPtMinTr[type]) continue;
427
428 AddToB2(iT, type, iP);
429 }
430 continue;
431 }
432 }
433
434
435
436 /*
437 for(Int_t i1=0; i1<fgkbinR; i1++)
438 {
439 for(Int_t i2=0; i2<2; i2++)
440 {
441 // if(fListJc[i1][i2]<2) fListJc[i1][i2]=0;
442 // if(fListB1c[i1][i2]<2) fListB1c[i1][i2]=0;
443 // if(fListB2c[i1][i2]<2) fListB2c[i1][i2]=0;
444 fListJ[i1][i2].Set(fListJc[i1][i2]);
445 fListB1[i1][i2].Set(fListB1c[i1][i2]);
446 fListB2[i1][i2].Set(fListB2c[i1][i2]);
447 }
448 }
449 */
450
451}
452
453
454
455Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::CalcR(TVector3 v1, TVector3 v2)
456{
457
458 Double_t dphi = CalcdPhi(v1.Phi(), v2.Phi());
459 // dphi*=(0.9/TMath::Pi());
460 Double_t deta = v1.Eta() - v2.Eta();
461 Double_t RB = TMath::Sqrt(dphi*dphi+deta*deta);
462 return RB;
463}
464
465
466Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::CalcdPhi(Double_t phi1, Double_t phi2)
467{
468
469 while(phi1<0) phi1+=TMath::TwoPi();
470 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
471
472 while(phi2<0) phi2+=TMath::TwoPi();
473 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
474
475 Double_t dphi = phi1- phi2;
476 if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
477 if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
478
479 return dphi;
480}
481
482
483void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::GetLocalMom(TVector3 vecJ, TVector3 *q)
484{
485 // theta and phi of a particle in loc.syst of fvecJ
486
487
488 Double_t PT = vecJ.Perp();
489 Double_t P0 = vecJ.Mag();
490
491 Double_t q0=q->X();
492 Double_t q1=q->Y();
493 Double_t q2=q->Z();
494
495 Double_t qx1 = vecJ(2)/P0/PT*(vecJ(0)*q0+vecJ(1)*q1) - PT/P0*q2;
496 Double_t qy1 = -vecJ(1)/PT*q0 + vecJ(0)/PT*q1;
497 Double_t qz1 = 1/P0*(vecJ(0)*q0+vecJ(1)*q1+vecJ(2)*q2);
498
499 // Double_t q00 = TMath::Sqrt(qx1*qx1 + qy1*qy1 + qz1*qz1);
500 // phi = TMath::ATan2(qy1, qx1);
501 // if(phi<0) phi+=fTwoPi;
502 // theta = TMath::ACos(qz1/q00);
503
504 q->SetXYZ(qx1, qy1, qz1);
505
506 return;
507
508}
509
510/*
511
512Bool_t AnalJetSubStrTool::FindCorrelationAxesAnd(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario)
513{
514 Double_t TwoPi = TMath::TwoPi();
515
516//
517// 1st track loop to determine the sphericity matrix
518//
519 // Initialize Shericity matrix
520 Float_t mxx = 0.;
521 Float_t myy = 0.;
522 Float_t mxy = 0.;
523 Int_t nc = 0;
524 Float_t sump2 = 0.;
525 Float_t ptMax = 0.;
526 Float_t etaMax = 0.;
527 Float_t phiMax = 0.;
528 Int_t iMax = -1;
529 Float_t ptsum = 0.;
530
531
532 Int_t Nev = list.GetSize();
533 if( Nev <2) return kFALSE;
534
535 Int_t indexpmax = -1;
536 Double_t pmax = -1;
537 Double_t phipmax = 0;
538
539 Int_t indexpzmax = -1;
540 Double_t pzmax = -1;
541 Double_t phipzmax = 0;
542
543
544 Int_t indexthetamin = -1;
545 Double_t thetamin = 1000;
546 Double_t phithetamin = 0;
547
548 for(Int_t ip=0; ip< Nev; ip++) {
549
550 TVector3* part = (TVector3*) GetAt(list.At(ip));
551 if(!part) continue;
552
553 TVector3 vtmp(*part);
554 GetLocalMom(vec, &vtmp);
555
556 Float_t ppjX = vtmp.X();
557 Float_t ppjY = vtmp.Y();
558 Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
559
560
561 Float_t pt = ppjT;//part->Z();
562 Float_t deta = part->Eta();
563 Float_t dphi = part->Phi();
564
565 mxx += (ppjX * ppjX / ppjT);
566 myy += (ppjY * ppjY / ppjT);
567 mxy += (ppjX * ppjY / ppjT);
568 nc++;
569 sump2 += ppjT;
570
571
572 if(vtmp.Mag()>pmax)
573 {
574 pmax=vtmp.Mag();
575 indexpmax = ip;
576 phipmax=vtmp.Phi();
577 }
578
579 if(vtmp.Z()>pzmax)
580 {
581 pzmax=vtmp.Z();
582 indexpzmax = ip;
583 phipzmax=vtmp.Phi();
584 }
585 if(vtmp.Theta()<thetamin)
586 {
587 thetamin=vtmp.Theta();
588 indexthetamin = ip;
589 phithetamin=vtmp.Phi();
590 }
591
592
593 }
594
595
596
597
598//
599// At this point we have mxx, myy, mxy
600 if (nc == 0) return kFALSE;
601// Shericity Matrix
602 const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
603 TMatrixDSym m0(2, ele);
604// Find eigenvectors
605 TMatrixDSymEigen m(m0);
606 TVectorD eval(2);
607 TMatrixD evecm = m.GetEigenVectors();
608 eval = m.GetEigenValues();
609// Largest eigenvector
610 Int_t jev = 0;
611 if (eval[0] < eval[1]) jev = 1;
612 TVectorD evec0(2);
613// Principle axis
614 evec0 = TMatrixDColumn(evecm, jev);
615 TVector2 evec(evec0[0], evec0[1]);
616// Principle axis from leading partice
617
618// Float_t phiM = TMath::ATan2(phiMax, etaMax);
619// TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM));
620// Float_t phistM = evecM.DeltaPhi(evec);
621// if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
622
623
624 Double_t phiTA = evec.Phi();
625
626
627 Double_t phiDir = 0;
628 if(scenario ==0)
629 {
630 phiDir = phipmax;
631 }
632 else if(scenario ==1)
633 {
634 phiDir = phipzmax;
635 }
636
637 else if(scenario ==2)
638 {
639 phiDir = phithetamin;
640 }
641 else
642 return kFALSE;
643
644
645 Phi = phiTA;
646
647 if( TMath::Abs(CalcdPhi(phiDir, phiTA)) <TMath::Pi()/2)
648 return kTRUE;
649 else {
650 Phi+=TMath::Pi();
651 if(Phi>TwoPi) Phi-=TwoPi;
652 return kTRUE;
653 }
654
655
656 return kTRUE;
657}
658*/
659
660Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::FindCorrelationAxes(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario)
661{
662 // scenario 0- Pmax, 1-Pzmax_Local, 2->cosMin, 3-PTmac, 4->PTmax_Local
663
664 Double_t xmean , ymean;
665
666 Double_t TwoPi = TMath::TwoPi();
667
668 xmean = 0;
669 ymean = 0;
670 Phi= -999;
671 //find axes
672
673 Double_t x=0;
674 Double_t x2=0;
675 Double_t y=0;
676 Double_t y2=0;
677 Double_t xy=0;
678
679 // Double_t phiLoc, thetaLoc;
680
681 Int_t N=0;
682
683 Int_t Nev = list.GetSize();
684 if( Nev <2) return kFALSE;
685
686
687
688
689 Int_t Index = -1;
690 Double_t val = -1;
691 if(scenario==2) val = 100;
692 Double_t phiDir = -99;
693
694
695
696 for(Int_t ip=0; ip< Nev; ip++) {
697 if(list.At(ip)<0) break;
698
699 TVector3* part = (TVector3*) GetAt(list.At(ip));
700 if(!part) continue;
701
702
703 Double_t xx, yy;
704 /*
705 TVector3 vtmp(*part);
706 GetLocalMom(vec, &vtmp);
707 Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
708 Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
709 */
710
711 xx = CalcdPhi(part->Phi(), vec.Phi());
712 yy = part->Eta() - vec.Eta();
713
714 x+=xx;
715 y+=yy;
716 x2+=(xx*xx);
717 y2+=(yy*yy);
718 xy+=(xx*yy);
719 N++;
720
721 /*
722 if(scenario ==0)
723 {
724 if(vtmp.Mag()>val)
725 {
726 val=vtmp.Mag();
727 Index = ip;
728 phiDir=vtmp.Phi();
729 }
730 }
731 else if(scenario ==1)
732 {
733 if(vtmp.Z()>val)
734 {
735 val=vtmp.Z();
736 Index = ip;
737 phiDir=vtmp.Phi();
738 }
739 }
740 else if(scenario ==2)
741 {
742 if(vtmp.Theta()<val)
743 {
744 val=vtmp.Theta();
745 Index = ip;
746 phiDir=vtmp.Phi();
747 }
748 }
749 else
750*/
751if(scenario ==3)
752 {
753 if(part->Pt()> val)
754 {
755 val =part->Pt();
756 Index = ip;
757 // phiDir=vtmp.Phi();
758 phiDir=TMath::ATan2(yy, xx); //vtmp.Phi();
759 }
760 }
761
762
763 }
764
765 if(N<2) return kFALSE;
766
767 // return kFALSE;
768
769 x/=N;
770 y/=N;
771 x2/=N;
772 y2/=N;
773 xy/=N;
774
775
776 Double_t A = xy-x*y;
777 Double_t B = x2-x*x+y*y-y2;
778 Double_t D = TMath::Sqrt(B*B+4*A*A);
779
780 // printf("N = %d A = %f\n",N, A);
781 // list.Print();
782
783 Double_t a1 = (-B+D)/2/A;
784 // Double_t a2 = (-B-D)/2/A;
785 // Double_t b1 = y-a1*x;
786 // Double_t b2 = y-a2*x;
787 // Double_t phiDir = TMath::ATan2(y, x);
788 // while(phiDir<0) phiDir+=TwoPi;
789 Double_t phi = TMath::ATan(a1);
790
791 while(phi>TwoPi) phi-=TwoPi;
792 while(phi< 0 ) phi+=TwoPi;
793
794 Phi=CalcdPhi0To2pi(phi, 0);
795
796 if(scenario!=4 && Index<0) return kFALSE;
797
798 if(Index>=0)
799 {
800 if( TMath::Abs(CalcdPhi(phiDir, phi)) <TMath::Pi()/2)
801 return kTRUE;
802 else
803 {
804 phi+=TMath::Pi();
805 Phi = CalcdPhi0To2pi(phi);
806 return kTRUE;
807 }
808 }
809
810
811 Double_t xmax = -100;
812 Double_t xmin = 100;
813
814 xmean = x;
815 ymean = y;
816
817
818 for(Int_t ip=0; ip< Nev; ip++) {
819 if(list.At(ip)<0) break;
820
821 TVector3* part = (TVector3*) GetAt(list.At(ip));
822 if(!part) continue;
823
824 TVector3 vtmp(*part);
825 GetLocalMom(vec, &vtmp);
826
827 Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
828 Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
829
830 Double_t x1 = (xx-xmean)*TMath::Cos(Phi) +(yy-ymean)*TMath::Sin(Phi);
831 if(x1 > xmax) xmax = x1;
832 if(x1 < xmin) xmin = x1;
833 // printf("c3\n");
834 }
835
836 if(TMath::Abs(xmin) > xmax) Phi+=TMath::Pi();
837
838
839
840
841 while(Phi>TwoPi) Phi-=TwoPi;
842
843 return kTRUE;
844}
845
846
847
848
849Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::FindCorrelationAxesCorr(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario, Double_t R)
850{
851 // scenario 0- Pmax, 1-Pzmax_Local, 2->cosMin, 3-PTmac, 4->PTmax_Local
852
853 Double_t xmean , ymean;
854
855 Double_t TwoPi = TMath::TwoPi();
856
857 xmean = 0;
858 ymean = 0;
859 Phi= -999;
860 //find axes
861
862 Double_t x=0;
863 Double_t x2=0;
864 Double_t y=0;
865 Double_t y2=0;
866 Double_t xy=0;
867
868 // Double_t phiLoc, thetaLoc;
869
870 Int_t N=0;
871
872 Int_t Nev = list.GetSize();
873 if( Nev <2) return kFALSE;
874
875
876
877
878 Int_t Index = -1;
879 Double_t val = -1;
880 if(scenario==2) val = 100;
881 Double_t phiDir = -99;
882
883
884
885 for(Int_t ip=0; ip< Nev; ip++) {
886 if(list.At(ip)<0) break;
887
888 TVector3* part = (TVector3*) GetAt(list.At(ip));
889 if(!part) continue;
890
891 Double_t xx, yy;
892 /*
893 TVector3 vtmp(*part);
894 GetLocalMom(vec, &vtmp);
895 Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
896 Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
897 */
898
899 xx = CalcdPhi(part->Phi(), vec.Phi());
900 yy = part->Eta() - vec.Eta();
901
902 x+=xx;
903 y+=yy;
904 x2+=(xx*xx);
905 y2+=(yy*yy);
906 xy+=(xx*yy);
907 N++;
908
909 /*
910 if(scenario ==0)
911 {
912 if(vtmp.Mag()>val)
913 {
914 val=vtmp.Mag();
915 Index = ip;
916 phiDir=vtmp.Phi();
917 }
918 }
919 else if(scenario ==1)
920 {
921 if(vtmp.Z()>val)
922 {
923 val=vtmp.Z();
924 Index = ip;
925 phiDir=vtmp.Phi();
926 }
927 }
928 else if(scenario ==2)
929 {
930 if(vtmp.Theta()<val)
931 {
932 val=vtmp.Theta();
933 Index = ip;
934 phiDir=vtmp.Phi();
935 }
936 }
937 else
938*/
939if(scenario ==3)
940 {
941 if(part->Pt()> val)
942 {
943 val =part->Pt();
944 Index = ip;
945 // phiDir=vtmp.Phi();
946 phiDir=TMath::ATan2(yy, xx); //vtmp.Phi();
947 }
948 }
949
950
951 }
952
953 if(N<2) return kFALSE;
954 if(scenario!=4 && Index<0) return kFALSE;
955
956 // return kFALSE;
957
958 x/=N;
959 y/=N;
960 x2/=N;
961 y2/=N;
962 xy/=N;
963
964
965 Double_t A = xy-x*y;
966 Double_t B = x2-x*x+y*y-y2;
967 Double_t D = TMath::Sqrt(B*B+4*A*A);
968
969 // printf("N = %d A = %f\n",N, A);
970 // list.Print();
971
972 Double_t a1 = (-B+D)/2/A;
973 // Double_t a2 = (-B-D)/2/A;
974 // Double_t b1 = y-a1*x;
975 // Double_t b2 = y-a2*x;
976 // Double_t phiDir = TMath::ATan2(y, x);
977 // while(phiDir<0) phiDir+=TwoPi;
978 Double_t phi = TMath::ATan(a1);
979
980 Phi=CalcdPhi0To2pi(phi, 0);
981 if( TMath::Abs(CalcdPhi(phiDir, phi)) > TMath::Pi()/2)
982 Phi=CalcdPhi0To2pi(phi+TMath::Pi(), 0);
983
984
985
986
987 if(Index>=0)
988 {
989 if(N==2)
990 return kTRUE;
991
992
993
994 Double_t phiMinus = 0;
995 Int_t Nminus = 0;
996 for(Int_t ip=0; ip< Nev; ip++)
997 {
998 if(list.At(ip)<0) break;
999
1000 TVector3* part = (TVector3*) GetAt(list.At(ip));
1001 if(!part) continue;
1002 Double_t xx0 = CalcdPhi(part->Phi(), vec.Phi());
1003 Double_t yy0 = part->Eta() - vec.Eta();
1004
1005 Double_t xx = (x*N - xx0)/(N-1);
1006 Double_t yy = (y*N - yy0)/(N-1);
1007 Double_t xx2 = (x2*N - xx0*xx0)/(N-1);
1008 Double_t yy2 = (y2*N - yy0*yy0)/(N-1);
1009 Double_t xxyy= (xy*N - xx0*yy0)/(N-1);
1010
1011
1012 Double_t AA = xxyy-xx*yy;
1013 Double_t BB = xx2-xx*xx+yy*yy-yy2;
1014 Double_t DD = TMath::Sqrt(BB*BB+4*AA*AA);
1015 Double_t aa1 = (-BB+DD)/2/AA;
1016 Double_t phi1 = TMath::ATan(aa1);
1017
1018 if( TMath::Abs(CalcdPhi(phiDir, phi1)) > TMath::Pi()/2)
1019 phi1+=TMath::Pi();
1020
1021 phiMinus+=CalcdPhi(Phi, phi1);
1022 Nminus++;
1023
1024 }
1025
1026
1027 Double_t phiPlus = 0;
1028 Int_t Nplus = 0;
1029 for(Int_t ip=0; ip< -Nev; ip++)
1030 {
1031 if(list.At(ip)<0) break;
1032
1033 TVector3* part = (TVector3*) GetAt(list.At(ip));
1034 if(!part) continue;
1035 Double_t rtmp = gRandom->Uniform(0, R);
1036 Double_t phitmp = gRandom->Uniform(0, TMath::TwoPi());
1037 Double_t xx0 = rtmp*TMath::Cos(phitmp);
1038 Double_t yy0 = rtmp*TMath::Sin(phitmp);
1039
1040 Double_t xx = (x*N + xx0)/(N+1);
1041 Double_t yy = (y*N + yy0)/(N+1);
1042 Double_t xx2 = (x2*N + xx0*xx0)/(N+1);
1043 Double_t yy2 = (y2*N + yy0*yy0)/(N+1);
1044 Double_t xxyy= (xy*N + xx0*yy0)/(N+1);
1045
1046
1047 Double_t AA = xxyy-xx*yy;
1048 Double_t BB = xx2-xx*xx+yy*yy-yy2;
1049 Double_t DD = TMath::Sqrt(BB*BB+4*AA*AA);
1050 Double_t aa1 = (-BB+DD)/2/AA;
1051 Double_t phi1 = TMath::ATan(aa1);
1052
1053 if( TMath::Abs(CalcdPhi(phiDir, phi1)) > TMath::Pi()/2)
1054 phi1+=TMath::Pi();
1055
1056 phiPlus+=CalcdPhi(Phi, phi1);
1057 Nplus++;
1058
1059 }
1060
1061
1062 Phi+=((phiMinus+phiPlus)/(Nminus+Nplus+1));
1063
1064 if( TMath::Abs(CalcdPhi(phiDir, Phi)) > TMath::Pi()/2)
1065 Phi+=TMath::Pi();
1066
1067 Phi=CalcdPhi0To2pi(Phi, 0);
1068
1069 return kTRUE;
1070 }
1071
1072
1073 Double_t xmax = -100;
1074 Double_t xmin = 100;
1075
1076 xmean = x;
1077 ymean = y;
1078
1079
1080 for(Int_t ip=0; ip< Nev; ip++) {
1081 if(list.At(ip)<0) break;
1082
1083 TVector3* part = (TVector3*) GetAt(list.At(ip));
1084 if(!part) continue;
1085
1086 TVector3 vtmp(*part);
1087 GetLocalMom(vec, &vtmp);
1088
1089 Double_t xx=vtmp.X(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Cos(phiLoc);
1090 Double_t yy=vtmp.Y(); //part->Mag()*TMath::Sin(thetaLoc)*TMath::Sin(phiLoc);
1091
1092 Double_t x1 = (xx-xmean)*TMath::Cos(Phi) +(yy-ymean)*TMath::Sin(Phi);
1093 if(x1 > xmax) xmax = x1;
1094 if(x1 < xmin) xmin = x1;
1095 // printf("c3\n");
1096 }
1097
1098 if(TMath::Abs(xmin) > xmax) Phi+=TMath::Pi();
1099
1100
1101
1102
1103 while(Phi>TwoPi) Phi-=TwoPi;
1104
1105 return kTRUE;
1106}
1107
89605a8b 1108/*
c87f1b5e 1109/////////////////////////////////////
1110
1111 TH2F *fhPhiEtaTrack;//!
1112
1113 TH1F *fhTMA_JAA[3];//!
1114 TH1F *fhTMA_JAp[3];//!
1115 TH1F *fhTMA_B1AA[3];//!
1116 TH1F *fhTMA_B1Ap[3];//!
1117 TH1F *fhTMA_B2AA[3];//!
1118 TH1F *fhTMA_B2Ap[3];//!
1119
1120
1121 Int_t fPtJetNbin;
1122 TArrayD fPtJetArray;
1123 Int_t fPsiVsRNbin;
1124 Double_t fRmax;
1125 Int_t fPsiVsPhiNbin;
1126
1127 UInt_t fFilterMask;
1128 Double_t fEtaTrackMax;
1129 Double_t fPtTrackMin;
1130 Double_t fPtTrackMax;
1131 Double_t fPtJmin;
1132 Double_t fPtJmax;
1133 Double_t fPtJ;
1134
1135 TVector3 fJvec;
89605a8b 1136*/
c87f1b5e 1137
89605a8b 1138AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AliAnalysisTaskJetShapeHM(TString comment, Bool_t kMC):TObject(),
c87f1b5e 1139fComment(comment),
1140fkMC(kMC),
89605a8b 1141fkMCprod(0),
c87f1b5e 1142fhEvent(0),
1143fhMult(0),
1144fhPtJ(0),
1145fhPtJvsPtCorr(0),
1146fhPsiVsR(0),
1147fhPsiVsRPtJ(0),
1148fhPhiEtaTrack(0),
89605a8b 1149fhPsiVsR_MCtr(0),
1150fhPsiVsRPtJ_MCtr(0),
1151fhJetTrPtVsPartPt(0),
c87f1b5e 1152fPtJetNbin(0),
1153fPtJetArray(0),
1154fPsiVsRNbin(0),
1155fRmax(0),
1156fPsiVsPhiNbin(0),
1157fFilterMask(0),
1158fEtaTrackMax(0),
1159fPtTrackMin(0),
1160fPtTrackMax(0),
1161fPtJmin(0),
1162fPtJmax(0),
1163fPtJ(0),
1164fJvec(0,0,0)
1165{
1166
1167 for(Int_t i=0; i<3; i++)
1168 {
1169 fhTMA_JAA[i]=0;
1170 fhTMA_JAp[i]=0;
1171 fhTMA_B1AA[i]=0;
1172 fhTMA_B1Ap[i]=0;
1173 fhTMA_B2AA[i]=0;
1174 fhTMA_B2Ap[i]=0;
1175 }
1176
89605a8b 1177 for(Int_t i=0; i<3; i++)
1178 {
1179 for(Int_t j=0; j<2; j++)
1180 {
1181 fhPtresVsPt[i][j]=0;
1182 fhPhiresVsPhi[i][j]=0;
1183 fhEtaresVsEta[i][j]=0;
1184 fhRresVsPt[i][j]=0;
1185 fhDCAxy[i][j]=0;
1186 fhDCAz[i][j]=0;
1187 fhTrackPtEtaPhi[i][j]=0;
1188 }
1189 }
1190
c87f1b5e 1191
1192 fPtJetNbin = 0;
1193 fPtJetArray = 0;
1194
1195 SetRBins();
1196 SetPhiNbins();
1197 SetFilterMask();
1198 SetEtaTrackMax();
1199 SetPtTrackRange();
1200 SetPtJetRange();
1201
1202 fJvec.SetXYZ(0,0,0);
1203 fPtJ = 0;
1204}
1205
1206
1207AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::~AliAnalysisTaskJetShapeHM()
1208{
1209 /*
1210 if(fPtJetArray)
1211 delete [] fPtJetArray;
1212 fPtJetArray=0;
1213 */
1214}
1215
1216
1217void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::SetPtJetBins(Int_t Nbin, Double_t *array)
1218{
1219
1220 /*
1221 if(fPtJetArray)
1222 delete [] fPtJetArray;
1223 fPtJetArray=0;
1224
1225 fPtJetNbin = Nbin;
1226 fPtJetArray = new Double_t[fPtJetNbin+1];
1227 for(Int_t i=0; i<= fPtJetNbin; i++) fPtJetArray[i]= array[i];
1228 */
1229
1230
1231 fPtJetNbin = Nbin;
1232 fPtJetArray.Set(fPtJetNbin+1, array);
1233}
1234
1235
1236void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::InitHistos()
1237{
1238
1239 fhEvent = Hist1D("hEvent" , 3 , -0.5, 2.5, "Event");
1240 fhMult = Hist1D("hMult" , 101, -0.5, 100.5, "N_{ch}");
1241
89605a8b 1242 fhPhiEtaTrack = Hist2D("hPhiEtaTrack", 100, 0, TMath::TwoPi(), 20, -1, 1., "#phi", "#eta");
c87f1b5e 1243
1244 if(fPtJetNbin<1) {
1245 Int_t Nbin = 13;
1246 Double_t array[]= {0., 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 130, 160, 200};
1247 SetPtJetBins(Nbin, array);
1248 }
1249
1250 fhPtJ = Hist1D("hPtJ" , fPtJetNbin, fPtJetArray.GetArray(), "Pt_{J} GeV/c");
c560b734 1251 // fhPsiVsR = Hist1D("hPsiVsR", fPsiVsRNbin, 0., fRmax, "R", 1, "#Psi(R)");
1252 fhPsiVsR = Hist3D("hPsiVsR", fPsiVsRNbin, 0., fRmax, 10, 0., 1., fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{Jt, frac}", "P_{J} GeV/c");
c87f1b5e 1253 fhPsiVsRPtJ = Hist2D("hPsiVsRPtJ", fPsiVsRNbin, 0., fRmax, fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{tJ} GeV/c", 1, "#Psi(R)");
1254
1255 fhPtJvsPtCorr = Hist2D("fhPtJvsPtCorr", fPtJetNbin, fPtJetArray.GetArray(), 100, -100, 100, "P_{tJ} GeV/c", "P_{tJ} - P_{tJ,corr} GeV/c" , 1);
1256
89605a8b 1257
1258
1259
c87f1b5e 1260 for(Int_t i=0; i<3; i++)
1261 {
1262 fhTMA_JAA[i] = Hist1D(Form("fhTMA_JAA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1263 fhTMA_JAp[i] = Hist1D(Form("fhTMA_JAp_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1264
1265 fhTMA_B1AA[i] = Hist1D(Form("fhTMA_B1AA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1266 fhTMA_B1Ap[i] = Hist1D(Form("fhTMA_B1Ap_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1267
1268 fhTMA_B2AA[i] = Hist1D(Form("fhTMA_B2AA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1269 fhTMA_B2Ap[i] = Hist1D(Form("fhTMA_B2Ap_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1270 }
1271
89605a8b 1272 if(fkMCprod)
1273 {
1274 fhPsiVsR_MCtr = Hist1D("hPsiVsR_MCtr", fPsiVsRNbin, 0., fRmax, "R", 1, "#Psi(R)");
1275 fhPsiVsRPtJ_MCtr = Hist2D("hPsiVsRPtJ_MCtr", fPsiVsRNbin, 0., fRmax, fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{tJ} GeV/c", 1, "#Psi(R)");
1276 fhJetTrPtVsPartPt = Hist2D("fhJetTrPtVsPartPt", fPtJetNbin, fPtJetArray.GetArray(), 100, -1, 1, "P_{tJ,part} GeV/c", "1-P_{tJ,tr}/P_{tJ,part} GeV/c" , 1);
c560b734 1277 const char *ch[2]={"m","p"};
89605a8b 1278 for(Int_t i=0; i<3; i++)
1279 {
1280 for(Int_t j=0; j<2; j++)
1281 {
1282 fhPtresVsPt[i][j] = Hist2D(Form("hPtresVsPt%d%s",i,ch[j]), 100, 0, 100, 100, -0.5, 0.5, "P_{t}^{gen} GeV/c", "1-P_{t}^{rec}/P_{t}^{gen} GeV/c" );
1283 fhPhiresVsPhi[i][j] = Hist2D(Form("hPhiresVsPhi%d%s",i,ch[j]), 600, 0, TMath::TwoPi(), 128, -0.2, 0.2, "#phi^{gen} rad", "#phi^{rec} - #phi^{gen} rad" );
1284 fhEtaresVsEta[i][j] = Hist2D(Form("hEtaresVsEta%d%s",i,ch[j]), 200, -1, 1, 40, -0.2, 0.2, "#eta^{gen}", "#eta^{rec} - #eta^{gen}" );
1285 fhRresVsPt[i][j] = Hist2D(Form("hRresVsPt%d%s",i,ch[j]) , 100, 0, 100, 500, -fRmax, fRmax, "P_{t, track}^{gen} GeV/c", "R^{gen}-R^{rec}" );
c560b734 1286 fhDCAxy[i][j] = Hist2D(Form("hDCAxy%d%s",i,ch[j]), 100, 0, 100, 300, -3, 3, Form("p_{t} of part. type %d%s [GeV/c]",i,ch[j]), "DCAxy [cm]");
1287 fhDCAz[i][j] = Hist2D(Form("hDCAz%d%s",i,ch[j]) , 100, 0, 100, 300, -3, 3, Form("p_{t} of part. type %d%s [GeV/c]",i,ch[j]), "DCAz [cm]") ;
89605a8b 1288 fhTrackPtEtaPhi[i][j]= Hist3D(Form("hTrackPtEtaPhi%d%s",i,ch[j]), 100, 0, 100, 100, -1, 1, 100, 0, TMath::TwoPi(),"P_{t} GeV/c ","#eta","#phi rad");
1289 }
1290 }
1291 }
c87f1b5e 1292
1293}
1294
1295
1296
1297
1298
1299void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddToList(TList *list)
1300{
1301 list->Add(fhEvent);
1302 list->Add(fhMult);
1303 list->Add(fhPhiEtaTrack);
1304 list->Add(fhPtJ);
1305 list->Add(fhPsiVsR);
1306 list->Add(fhPsiVsRPtJ);
1307 list->Add(fhPtJvsPtCorr);
1308
89605a8b 1309
c87f1b5e 1310 for(Int_t i=0; i<3; i++)
1311 {
1312 list->Add(fhTMA_JAA[i]);
1313 list->Add(fhTMA_JAp[i]);
1314
1315 list->Add(fhTMA_B1AA[i]);
1316 list->Add(fhTMA_B1Ap[i]);
1317
1318 list->Add(fhTMA_B2AA[i]);
1319 list->Add(fhTMA_B2Ap[i]);
1320 }
1321
89605a8b 1322 if(fkMCprod)
1323 {
1324 list->Add(fhPsiVsR_MCtr);
1325 list->Add(fhPsiVsRPtJ_MCtr);
1326 list->Add(fhJetTrPtVsPartPt);
c87f1b5e 1327
89605a8b 1328 for(Int_t i=0; i<3; i++)
1329 {
1330 for(Int_t j=0; j<2; j++)
1331 {
1332 list->Add(fhPtresVsPt[i][j]);
1333 list->Add(fhPhiresVsPhi[i][j]);
1334 list->Add(fhEtaresVsEta[i][j]);
1335 list->Add(fhRresVsPt[i][j]);
1336 list->Add(fhDCAxy[i][j]);
1337 list->Add(fhDCAz[i][j]);
1338 list->Add(fhTrackPtEtaPhi[i][j]);
1339 }
1340 }
1341 }
c87f1b5e 1342}
1343
1344
1345
1346
1347Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddEvent(AliAODEvent* aodE, AliAODJet *jet, Double_t DeltaPtJ)
1348{
1349
1350
1351
1352 Int_t size = aodE->GetNumberOfTracks();
1353
1354 TClonesArray *arrP = 0;
1355
89605a8b 1356 if(fkMC || fkMCprod)
c87f1b5e 1357 {
1358 arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
1359 if(!arrP)
1360 {
1361 printf("ERROR: no Info about particles in AOD\n");
1362 return kFALSE;
1363 }
c87f1b5e 1364 }
1365
89605a8b 1366 if(fkMC)
1367 size = arrP->GetEntriesFast();
c87f1b5e 1368
1369 Int_t IndexArray[size];
89605a8b 1370 Int_t IndexArrayMC[size];
c87f1b5e 1371
1372 TClonesArray farray("TVector3", size);
1373
1374 Int_t counter = 0;
89605a8b 1375 Int_t counterMC = 0;
c87f1b5e 1376 Int_t counterAll = 0;
1377 Double_t pJ[3] = {0, 0, 0};
89605a8b 1378 Double_t pJmc[3] = {0, 0, 0};
c87f1b5e 1379
89605a8b 1380 AliAODVertex* primVtx = aodE->GetPrimaryVertex();
1381 Double_t bfield = aodE->GetMagneticField();
1382 Double_t dca[2];
1383 Double_t cov[3];
c87f1b5e 1384
1385 TVector3 vecJ(jet->Px(),jet->Py(),jet->Pz());
1386
c87f1b5e 1387
c560b734 1388
1389
1390 for(int it = 0;it < size;++it){
c87f1b5e 1391 TVector3 vec;
89605a8b 1392 Int_t ch = -999;
c560b734 1393 Int_t label = 0;
c87f1b5e 1394
89605a8b 1395
c87f1b5e 1396 if(fkMC)
1397 {
1398 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(it));
1399 if(!part)continue;
1400 if(!part->IsPhysicalPrimary())continue;
1401 if(part->Charge()==0)continue;
1402 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1403 }
1404 else
1405 {
1406 AliAODTrack *tr = aodE->GetTrack(it);
1407 if(!tr) continue;
1408 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
89605a8b 1409 label = tr->GetLabel();
1410 AliAODTrack tmp(*tr);
1411 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
c87f1b5e 1412 vec.SetXYZ(tr->Px(), tr->Py(), tr->Pz());
89605a8b 1413 ch = tr->Charge();
c87f1b5e 1414 }
1415
1416
1417 if(TMath::Abs(vec.Eta())>fEtaTrackMax)continue;
1418 if(vec.Pt()< fPtTrackMin)continue;
1419 if(vec.Pt()> fPtTrackMax) {return kFALSE;}
1420
1421
1422 new(farray[counterAll]) TVector3(vec);
1423 counterAll++;
1424
1425
1426 Double_t R = CalcR(vecJ, vec);
1427 if(R> fRmax) continue;
1428
1429 pJ[0]+=vec.Px();
1430 pJ[1]+=vec.Py();
1431 pJ[2]+=vec.Pz();
1432
1433 IndexArray[counter] = it;
1434 counter++;
89605a8b 1435
1436
1437
1438 // effics
1439 if(fkMCprod)
1440 {
1441 // printf("A1\n");
1442 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(label)));
1443 if(!part)continue;
1444
1445 Int_t type = 0;
1446 if(!part->IsPhysicalPrimary()) type=1;
1447 if(label <0) type=2;
1448
1449 if(!(ch==-1 || ch==1)) AliFatal("ch != +/- 1!!!\n\n");
1450 // printf("A4\n");
1451 if(ch==-1) ch=0;
1452
1453 Double_t phip = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(part->Phi());
1454 Double_t dphi = AliAnalysisTaskJetShapeTool::CalcdPhiSigned(part->Phi(), vec.Phi());
c560b734 1455 Double_t phiT = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
89605a8b 1456
1457 fhPtresVsPt[type][ch]->Fill(part->Pt(), 1-vec.Pt()/part->Pt());
1458 fhPhiresVsPhi[type][ch]->Fill(phip, dphi);
1459 fhEtaresVsEta[type][ch]->Fill(part->Eta(), vec.Eta() - part->Eta());
1460 fhDCAxy[type][ch]->Fill(part->Pt(), dca[0]);
1461 fhDCAz[type][ch]->Fill( part->Pt(), dca[1]);
c560b734 1462 fhTrackPtEtaPhi[type][ch]->Fill(vec.Pt(), vec.Eta(), phiT);
89605a8b 1463
1464 TVector3 vecP(part->Px(), part->Py(), part->Pz());
1465 Double_t Rgen = CalcR(vecJ, vecP);
1466 fhRresVsPt[type][ch]->Fill(part->Pt(), Rgen - R);
1467
1468 pJmc[0]+=part->Px();
1469 pJmc[1]+=part->Py();
1470 pJmc[2]+=part->Pz();
1471
1472 IndexArrayMC[counterMC] = label;
1473 counterMC++;
1474 // printf("A5\n");
1475 }
1476
c87f1b5e 1477 }
1478
1479 fhMult->Fill(counter);
1480 if(counter<1) return kFALSE;
1481
1482 fJvec.SetXYZ(pJ[0], pJ[1], pJ[2]);
1483
89605a8b 1484
c87f1b5e 1485 fPtJ = TMath::Sqrt(pJ[0]*pJ[0] + pJ[1]*pJ[1]);
1486 if((fPtJ < fPtJmin) || (fPtJ > fPtJmax)) return kFALSE;
1487 fhPtJ->Fill(fPtJ);
1488
89605a8b 1489
1490
1491
1492
1493 fhPtJvsPtCorr->Fill(fPtJ, jet->Pt() - DeltaPtJ);
c87f1b5e 1494
1495 for(Int_t it = 0; it<counter; it++)
1496 {
1497 TVector3 vec;
1498
1499 if(fkMC)
1500 {
1501 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(IndexArray[it]));
89605a8b 1502 if(!part) continue;
1503 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
c87f1b5e 1504 }
1505 else
1506 {
1507 AliAODTrack *tr = aodE->GetTrack(IndexArray[it]);
89605a8b 1508 if(!tr) continue;
1509 AliAODTrack tmp(*tr);
1510 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1511 if(tr)vec.SetXYZ(tmp.Px(), tmp.Py(), tmp.Pz());
c87f1b5e 1512 }
1513
1514 Double_t R = CalcR(fJvec, vec);
1515 // Double_t pt = (tr->Px()*pJ[0] + tr->Py()*pJ[1])/PtJ;
c560b734 1516 Double_t probL = fJvec.Dot(vec)/fJvec.Mag2();
89605a8b 1517 // fhPsiVsR->Fill(R, probL);
1518 // fhPsiVsRPtJ->Fill(R, fPtJ, probL);
c560b734 1519 fhPsiVsR->Fill(R,probL, fJvec.Mag());
89605a8b 1520 fhPsiVsRPtJ->Fill(R, fPtJ);
c87f1b5e 1521
89605a8b 1522 Double_t phi = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
c87f1b5e 1523 fhPhiEtaTrack->Fill(phi, vec.Eta());
1524
89605a8b 1525 }
1526
1527
1528 if(fkMCprod)
1529 {
1530 TVector3 fJvecMCtr(pJmc[0], pJmc[1], pJmc[2]);
c560b734 1531 Double_t fPtJMCtr= fJvecMCtr.Perp();
89605a8b 1532
c560b734 1533 fhJetTrPtVsPartPt->Fill(fPtJMCtr, 1-fPtJ/fPtJMCtr);
89605a8b 1534 for(Int_t it = 0; it<counterMC; it++)
1535 {
1536 TVector3 vec;
c87f1b5e 1537
c560b734 1538 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(IndexArrayMC[it])));
89605a8b 1539 if(!part) continue;
1540 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1541
1542 Double_t R = CalcR(fJvecMCtr, vec);
1543 // Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1544 // fhPsiVsR->Fill(R, probL);
1545 // fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1546 // Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1547 fhPsiVsR_MCtr->Fill(R);
1548 fhPsiVsRPtJ_MCtr->Fill(R, fPtJMCtr);
1549 }
1550 }
c87f1b5e 1551
1552 // ang. distr.
1553 // fMyTool->Clean();
1554 AliAnalysisTaskJetShapeTool *fMyTool = new AliAnalysisTaskJetShapeTool(&farray);
1555
1556 fMyTool->SetNEntries(counterAll);
1557 fMyTool->SetVecJ(vecJ);
1558 fMyTool->SetPtMinTr(fPtTrackMin,fPtTrackMin);
1559 fMyTool->Init();
1560 Int_t scenario = 3;
1561
1562 /*
1563 to be added!!!!!!!!!
1564 fhTMA_B1AA[i]=0;
1565 fhTMA_B2AA[i]=0;
1566 */
1567
1568 for(Int_t l=0; l<3; l++)
1569 {
1570 Double_t phiA, phi1;
1571
1572 // Double_t ptRJ0= fMyTool->GetPRecInRJ().Pt();
1573 // Double_t ptRJ = fMyTool->GetPRecJ(l,0).Pt();
1574 Int_t N1 = fMyTool->GetSizeJ(l,1);
1575 Double_t dphi = -999;
1576
1577
1578 if(!fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 0) , vecJ, phiA, scenario))
1579 continue;
1580
1581 // f2JEvent->SetPhiJL(l,0, phiA);
1582
1583 if(fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 1) , vecJ, phi1, scenario))
1584 {
1585 fhTMA_JAA[l]->Fill(fMyTool->CalcdPhi0To2pi(phiA-phi1));
1586 }
1587
1588
1589 for(Int_t ip =0; ip<N1; ip++)
1590 {
1591 phi1 = fMyTool->GetLocPhiJ(l,1,ip);
1592 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1593 fhTMA_JAp[l]->Fill(dphi);
1594 }
1595
1596
1597 Int_t NB1 = fMyTool->GetSizeB1(l, 1);
1598 for(Int_t ip =0; ip<NB1; ip++)
1599 {
1600 phi1 = fMyTool->GetLocPhiB1(l,1,ip);
1601 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1602 fhTMA_B1Ap[l]->Fill(dphi);
1603 }
1604
1605 Int_t NB2 = fMyTool->GetSizeB2(l, 1);
1606 for(Int_t ip =0; ip<NB2; ip++)
1607 {
1608 phi1 = fMyTool->GetLocPhiB2(l,1,ip);
1609 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1610 fhTMA_B2Ap[l]->Fill(dphi);
1611 }
1612
1613 }
1614
1615 fhEvent->Fill(1);
1616
1617 delete fMyTool;
1618
1619 return kTRUE;
1620}
1621
1622
1623Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcR(TVector3 v1, TVector3 v2)
1624{
1625
1626 Double_t dphi = CalcdPhi(v1.Phi(), v2.Phi());
1627 Double_t deta = v1.Eta() - v2.Eta();
1628 Double_t RB = TMath::Sqrt(dphi*dphi+deta*deta);
1629 return RB;
1630}
1631
1632
1633Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcdPhi(Double_t phi1, Double_t phi2)
1634{
1635
1636 while(phi1<0) phi1+=TMath::TwoPi();
1637 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
1638
1639 while(phi2<0) phi2+=TMath::TwoPi();
1640 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
1641
1642 Double_t dphi = phi1- phi2;
1643 if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
1644 if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
1645
1646 return dphi;
1647}
1648
1649
1650
1651
89605a8b 1652
1653
1654
c87f1b5e 1655TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax, const char* xLabel, Int_t color, const char* yLabel)
1656{
1657// create a 1D histogram
1658
1659 TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xMin, xMax);
1660 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1661 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1662 res->SetLineColor(color);
1663 return res;
1664}
1665
1666
1667TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t *xArray, const char* xLabel, Int_t color, const char* yLabel)
1668{
1669// create a 1D histogram
1670
1671 TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xArray);
1672 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1673 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1674 res->SetLineColor(color);
1675 return res;
1676}
1677
1678TH2F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist2D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax, Int_t nBinsy, Double_t yMin, Double_t yMax, const char* xLabel, const char* yLabel, Int_t color)
1679{
1680// create a 2D histogram
1681
1682 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yMin, yMax);
1683 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1684 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1685 // res->SetMarkerStyle(kFullCircle);
1686 // res->SetOption("E");
1687 res->SetLineColor(color);
1688 // fOutputList->Add(res);
1689 return res;
1690}
1691
1692
1693TH2F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist2D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax, Int_t nBinsy, Double_t *yArray, const char* xLabel, const char* yLabel, Int_t color, const char* zLabel)
1694{
1695// create a 2D histogram
1696
1697 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yArray);
1698 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1699 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1700 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1701 // res->SetMarkerStyle(kFullCircle);
1702 // res->SetOption("E");
1703 res->SetLineColor(color);
1704 // fOutputList->Add(res);
1705 return res;
1706}
1707
1708TH2F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist2D(const char* name, Int_t nBinsx, Double_t *yArrax, Int_t nBinsy, Double_t yMin, Double_t yMax, const char* xLabel, const char* yLabel, Int_t color, const char* zLabel)
1709{
1710// create a 2D histogram
1711
1712 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, yArrax, nBinsy, yMin, yMax);
1713 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1714 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1715 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1716 // res->SetMarkerStyle(kFullCircle);
1717 // res->SetOption("E");
1718 res->SetLineColor(color);
1719 // fOutputList->Add(res);
1720 return res;
1721}
1722
89605a8b 1723TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1724 Int_t nBinsy, Double_t yMin, Double_t yMax,
1725 Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1726{
1727// create a 3D histogram
1728
1729 TH3F *res = new TH3F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yMin, yMax, nBinsz, zMin, zMax);
1730 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1731 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1732 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1733 // res->SetMarkerStyle(kFullCircle);
c560b734 1734 // res->SetOption("E");
1735 res->SetLineColor(color);
1736 // fOutputList->Add(res);
1737 return res;
1738}
1739
1740TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1741 Int_t nBinsy, Double_t yMin, Double_t yMax,
1742 Int_t nBinsz, Double_t *zArr, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1743{
1744// create a 3D histogram
1745
1746 Double_t xArr[nBinsx+1], yArr[nBinsy+1];
1747
1748 for(Int_t i=0; i<=nBinsx; i++) xArr[i]=xMin+i*(xMax-xMin)/nBinsx;
1749 for(Int_t i=0; i<=nBinsy; i++) yArr[i]=yMin+i*(yMax-yMin)/nBinsy;
1750
1751 TH3F *res = new TH3F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xArr, nBinsy, yArr, nBinsz, zArr);
1752 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1753 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1754 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1755 // res->SetMarkerStyle(kFullCircle);
89605a8b 1756 // res->SetOption("E");
1757 res->SetLineColor(color);
1758 // fOutputList->Add(res);
1759 return res;
1760}
c87f1b5e 1761
1762
1763//////////////////////////////////////
1764
1765
1766/*
1767//________________________________________________________________________
1768void AnalysisJetMain::ConnectInputData(Option_t *)
1769{
1770 // Connect ESD
1771 // Called once
1772
1773 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
1774 // if (!fESD) {
1775 // AliError("ESD not available");
1776 // fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
1777 // fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1778
1779 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1780 // assume that the AOD is in the general output...
1781 fAODOut = AODEvent();
1782
1783
1784}
1785*/
1786
1787void AliAnalysisTaskJetShape::SetBranchNames(const TString &branch1, const TString &branch2)
1788{
1789 fJetBranchName[0] = branch1;
1790 fJetBranchName[1] = branch2;
1791}
1792
1793//________________________________________________________________________
1794void AliAnalysisTaskJetShape::UserCreateOutputObjects()
1795{
1796 //printf("Open1\n");
1797 // const char *nameF = OpenFile(1)->GetName();
1798 //printf("Open2 %s\n",nameF);
1799
1800 // fTriggerAnalysis = new AliTriggerAnalysis();
1801 // if (fkMC) fTriggerAnalysis->SetAnalyzeMC(1);
1802
1803
1804 fOutputList = new TList();
1805 fOutputList->SetOwner();
1806
1807 fhPtJL = Hist1D("hPtJL", 100, 0, 200, "Pt_{JL}"); fOutputList->Add(fhPtJL);
1808 fhAreaJL = Hist1D("hAreaJL", 100, 0., 4, "AreaJL"); fOutputList->Add(fhAreaJL);
1809
89605a8b 1810
1811
1812 printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() \n");
1813
1814
c87f1b5e 1815 fanalJetSubStrHM->SetFilterMask(fFilterMask);
89605a8b 1816 if(fkMC) fanalJetSubStrHM->MCprod();
c87f1b5e 1817 fanalJetSubStrHM->InitHistos();
1818 fanalJetSubStrHM->AddToList(fOutputList);
1819
c87f1b5e 1820
1821 if(fkMC)
1822 {
1823
1824 printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() MC\n");
1825 fanalJetSubStrMCHM->InitHistos();
1826 fanalJetSubStrMCHM->AddToList(fOutputList);
89605a8b 1827
1828 for(Int_t i=0; i<3; i++)
1829 {
1830 fhPtresVsPt[i] = Hist2D(Form("hPtresVsPt%d",i), 100, 0, 100, 100, -0.5, 0.5, "P_{t}^{gen} GeV/c", "1-P_{t}^{rec}/P_{t}^{gen} GeV/c" ); fOutputList->Add(fhPtresVsPt[i]);
1831 fhPhiresVsPhi[i] = Hist2D(Form("hPhiresVsPhi%d",i), 600, 0, TMath::TwoPi(), 128, -0.2, 0.2, "#phi^{gen} rad", "#phi^{rec} - #phi^{gen} rad" ); fOutputList->Add(fhPhiresVsPhi[i]);
1832 fhEtaresVsEta[i] = Hist2D(Form("hEtaresVsEta%d",i), 200, -1, 1, 40, -0.2, 0.2, "#eta^{gen}", "#eta^{rec} - #eta^{gen}" ); fOutputList->Add(fhEtaresVsEta[i]);
1833 fhDCAxy[i] = Hist2D(Form("hDCAxy%d",i), 100, 0, 100, 300, -3, 3, "DCAxy of prim [cm]"); fOutputList->Add(fhDCAxy[i]);
1834 fhDCAz[i] = Hist2D(Form("hDCAz%d",i) , 100, 0, 100, 300, -3, 3, "DCAz of prim [cm]") ; fOutputList->Add(fhDCAz[i]);
1835 fhTrackPtEtaPhi[i]= Hist3D(Form("hTrackPtEtaPhi%d",i), 100, 0, 100, 100, -1, 1, 100, 0, TMath::TwoPi(),"P_{t} GeV/c ","#eta","#phi rad"); fOutputList->Add(fhTrackPtEtaPhi[i]);
1836 }
c87f1b5e 1837 }
1838
1839
1840
1841
89605a8b 1842 printf(" JetBranchName[0]=%s\n JetBranchName[1]=%s\n", fJetBranchName[0].Data(), fJetBranchName[1].Data());
1843
1844
1845
c87f1b5e 1846 PostData(1, fOutputList);
1847
1848 /*
1849 OpenFile(1);
1850
1851 fOutputTree = new TTree("tree","Tree with Ali2JEvent");
1852 fOutputTree->Branch("event",&f2JEvent);
1853 */
1854
1855// PostData(1, fOutputTree);
1856}
1857
1858/*
1859Bool_t AliAnalysisTaskJetSpectrum2::Notify()
1860{
1861
1862 // Fetch the aod also from the input in,
1863 // have todo it in notify
1864
1865
1866 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1867 // assume that the AOD is in the general output...
1868 fAODOut = AODEvent();
1869
1870 if(fNonStdFile.Length()!=0){
1871 // case that we have an AOD extension we need can fetch the jets from the extended output
1872 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1873 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
1874 if(!fAODExtension){
1875 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
1876 }
1877 }
1878
1879
1880
1881 return kTRUE;
1882}
1883*/
1884
1885
1886
1887
1888//________________________________________________________________________
1889void AliAnalysisTaskJetShape::UserExec(Option_t *)
1890{
1891 // return;
1892// if(f2JEvent)
1893// delete f2JEvent;
1894// f2JEvent = new Ali2JEvent();
1895// return;
1896
1897 // return;
89605a8b 1898 if(fDebug) Printf("\n\n\nEvent #%5d", (Int_t) fEntry);
c87f1b5e 1899 if(fDebug) printf("NEW EVENT 0\n");
1900
1901 if(!IsGoodEvent()) return;
1902
c87f1b5e 1903
1904 // printf("\n\n\n NEW EVENT\n");
c87f1b5e 1905
1906 AliAODEvent* aodE = 0;
1907 if(!aodE){
1908 if(!fESD)aodE = fAODIn;
1909 else aodE = fAODOut;}
1910
89605a8b 1911 /*
1912 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1913if(!aodH){
1914 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
1915 return;
1916 }
1917 */
1918 if(fDebug) {
1919 printf("NEW EVENT 1: Number of tracks = %d\n",aodE->GetNumberOfTracks());
1920 // aodE->GetList()->Print();
1921 // printf("Look for %s\n",fJetBranchName[0].Data());
1922 }
c87f1b5e 1923
1924
1925 // centrality selection
1926 AliCentrality *cent = 0x0;
1927 Double_t centrality = 0.;
1928
1929 if(fESD) {cent = fESD->GetCentrality();
1930 if(cent) centrality = cent->GetCentralityPercentile("V0M");}
1931 else centrality=aodE->GetHeader()->GetCentrality();
1932
1933
1934 if(!fkIsPbPb) {
1935 centrality=1.;
1936 }
1937
89605a8b 1938// if(fDebug) printf("NEW EVENT 2\n");
c87f1b5e 1939
1940 Bool_t fUseAOD049 = kFALSE;// kTRUE;//
1941 if(fUseAOD049&&centrality>=0){
1942 Float_t v0=0;
1943 AliAODVZERO* aodV0 = aodE->GetVZEROData();
1944 v0+=aodV0->GetMTotV0A();
1945 v0+=aodV0->GetMTotV0C();
1946 if(centrality==0 && v0 < 19500) return ;//filtering issue
1947 Float_t tkl = (Float_t)(aodE->GetTracklets()->GetNumberOfTracklets());
1948 Float_t val= 1.30552 + 0.147931 * v0;
1949 Float_t tklSigma[101]={176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86, 120.788, 115.611, 113.172, 110.496,
1950 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654, 92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
1951 68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224, 51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114,
1952 43.2083, 41.3065, 40.1863, 38.5255, 37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398, 26.6488, 25.0183,
1953 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235, 19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173,
1954 12.9504, 12.9504, 12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544,
1955 13.7544, 13.7544, 13.7544, 13.7544, 13.7544};
1956 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)centrality] )
1957 return; //outlier
1958 }
1959
1960
89605a8b 1961 if(fDebug)
1962 {
1963 printf("centrality: %f\n", centrality);
1964 aodE->Print();
1965 aodE->GetList()->Print();
1966 }
1967
1968
c87f1b5e 1969 if (centrality < fCentMin || centrality > fCentMax){
1970 // PostData(1, fOutputList);
1971 return;
1972 }
1973
1974 // fhNEvent->Fill(0);
1975
1976 // accepted events
1977 // -- end event selection --
1978
c87f1b5e 1979
1980 // get background
1981 AliAODJetEventBackground* externalBackground = 0;
1982 AliAODJetEventBackground* externalBackgroundMC = 0;
1983
1984
1985 Float_t rho = 0;
1986 Float_t rho_MC=0;
1987
1988 if(fkIsPbPb)
1989 {
1990 if(fAODOut&&!externalBackground&&fBackgroundBranch[0].Length()){
1991 externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[0].Data()));
1992 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1993 }
1994 if(fAODExtension&&!externalBackground&&fBackgroundBranch[0].Length()){
1995 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[0].Data()));
1996 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1997 }
1998 if(fAODIn&&!externalBackground&&fBackgroundBranch[0].Length()){
1999 externalBackground = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[0].Data()));
2000 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
2001 }
2002 if(fAODOut&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2003 externalBackgroundMC = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[1].Data()));
2004 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2005 }
2006 if(fAODExtension&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2007 externalBackgroundMC = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[1].Data()));
2008 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2009 }
2010 if(fAODIn&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2011 externalBackgroundMC = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[1].Data()));
2012 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2013 }
2014 if(externalBackground)rho = externalBackground->GetBackground(0);
2015 if(externalBackgroundMC)rho_MC = externalBackgroundMC->GetBackground(0);
2016 }
2017
2018
2019
89605a8b 2020 if(fkMC)
2021 {
2022 AliAODVertex* primVtx = aodE->GetPrimaryVertex();
2023 Double_t bfield = aodE->GetMagneticField();
2024
2025 TClonesArray *arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
2026 if(!arrP)
2027 {
2028 printf("ERROR: no Info about particles in AOD\n");
2029 return;
2030 }
2031
2032 for(int it = 0;it < aodE->GetNumberOfTracks(); it++)
2033 {
2034 AliAODTrack *tr = aodE->GetTrack(it);
2035 if(!tr) continue;
2036 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))) continue;
2037 if(TMath::Abs(tr->Eta())>1.) continue;
2038 Int_t label = TMath::Abs(tr->GetLabel());
2039
2040 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(label));
2041 if(!part)continue;
2042 // if(!part->IsPhysicalPrimary())continue;
2043 // if(part->Charge()==0)continue;
2044 // vec.SetXYZ(part->Px(), part->Py(), part->Pz());
2045
2046
2047
2048 Double_t dca[2];
2049 Double_t cov[3];
2050
2051 AliAODTrack tmp(*tr);
2052 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
2053
2054 Int_t type = 0;
2055 if(!part->IsPhysicalPrimary()) type=1;
2056 if(label<0) type =2;
2057
2058 fhPtresVsPt[type]->Fill(part->Pt(), 1-tr->Pt()/part->Pt());
2059 fhPhiresVsPhi[type]->Fill(part->Phi(), tr->Phi() - part->Phi());
2060 fhEtaresVsEta[type]->Fill(part->Eta(), tr->Eta() - part->Eta());
2061 fhDCAxy[type]->Fill(part->Pt(), dca[0]);
2062 fhDCAz[type]->Fill( part->Pt(), dca[1]);
2063 fhTrackPtEtaPhi[type]->Fill(tr->Pt(), tr->Eta(), tr->Phi());
2064 }
2065
2066
2067 }
2068
2069
2070
2071
c87f1b5e 2072
2073 // printf("rho = %f\n",rho);
2074
2075 // return;
2076
2077 TClonesArray *Jets[2];
2078 Jets[0]=0;
2079 Jets[1]=0;
2080 if(fAODOut&&!Jets[0]){
2081 Jets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
2082 Jets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); }
2083 if(fAODExtension && !Jets[0]){
2084 Jets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
2085 Jets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
2086 if(fAODIn&&!Jets[0]){
2087 Jets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
2088 Jets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
2089
2090
2091 if(!Jets[0]) {
2092 Printf("no friend!!!\n");
2093 return;
2094 }
2095 Int_t nJ = Jets[0]->GetEntries();
2096 Float_t ptmax = 0.;
2097 Int_t imax = -1;
2098
89605a8b 2099 if(fDebug) {
2100 printf("NEW EVENT 3: Number of Rec. jets %d\n",nJ);
2101 }
c87f1b5e 2102
2103//
2104// Find highest pT jet with pt > 20 GeV
2105//
89605a8b 2106// fPtJcorrMin=0;
c87f1b5e 2107
2108 Float_t etaJmax = 0.4;
2109
2110 for (Int_t i = 0; i < nJ; i++) {
2111 AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[0]->At(i));
2112 if (!jet) continue;
89605a8b 2113 // jet->Print("0");
c87f1b5e 2114 Float_t ptJ = jet->Pt();
2115 Float_t etaJ = TMath::Abs(jet->Eta());
2116
2117 Float_t area = jet->EffectiveAreaCharged();
2118 Float_t ptJcorr=ptJ-rho*area;
2119
2120 if ((ptJcorr > fPtJcorrMin) && (ptJcorr > ptmax) && etaJ < etaJmax) {
2121 ptmax = ptJcorr;
2122 imax = i;
2123 }
2124
2125 }
2126
89605a8b 2127
2128 if(fDebug) {
2129 printf("NEW EVENT 4:\n");
2130 }
2131
c87f1b5e 2132 TVector3 vecJ;
2133
2134 AliAODJet* jetL = 0;
2135
2136 Float_t areaJL, ptJLcorr;
2137
2138 if (imax != -1) {
2139
0bf6381d 2140 jetL = dynamic_cast<AliAODJet*> (Jets[0]->At(imax));
2141 if(jetL){
c87f1b5e 2142
0bf6381d 2143 areaJL = jetL->EffectiveAreaCharged();
2144 ptJLcorr = jetL->Pt()-rho*areaJL;
2145
2146 fhPtJL->Fill(ptJLcorr);
2147 fhAreaJL->Fill(areaJL);
2148 vecJ.SetXYZ(jetL->Px(), jetL->Py(), jetL->Pz());
2149 fanalJetSubStrHM->AddEvent(aodE, jetL, rho*areaJL);
89605a8b 2150
2151 if(fDebug) {
2152 Printf("Leading Jet");
2153 jetL->Print("0");
2154 }
2155
0bf6381d 2156 }
c87f1b5e 2157 }
2158
2159
c87f1b5e 2160 if(!fkMC)
2161 {
2162 PostData(1, fOutputList);
89605a8b 2163 if(fDebug) Printf("End of event #%5d", (Int_t) fEntry);
c87f1b5e 2164 return;
2165 }
2166
2167
89605a8b 2168
2169
c87f1b5e 2170 nJ = Jets[1]->GetEntries();
89605a8b 2171 if(fDebug) {
2172 printf("NEW EVENT 5: Number of Rec. jets %d\n",nJ);
2173 }
2174
c87f1b5e 2175 ptmax = 0;
2176 imax = -1;
2177 Float_t areaJL_MC=0;
2178
2179 for (Int_t i = 0; i < nJ; i++) {
2180 AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[1]->At(i));
2181 if (!jet) continue;
2182 Float_t ptJ1 = jet->Pt();
2183 Float_t etaJ1 = TMath::Abs(jet->Eta());
2184
2185 areaJL_MC = jet->EffectiveAreaCharged();
2186 Float_t ptJcorr=ptJ1-rho*areaJL_MC;
2187
2188 // jet->Print();
2189
2190 if ((ptJcorr > fPtJcorrMin) && (ptJcorr > ptmax) && etaJ1 < etaJmax) {
2191 ptmax = ptJcorr;
2192 imax = i;
2193 }
2194
2195 }
2196
89605a8b 2197 if(fDebug) {
2198 printf("NEW EVENT 6:\n");
2199 }
2200
c87f1b5e 2201
2202 AliAODJet* jetMCL=0;
2203
2204 if (imax != -1) {
2205 jetMCL = dynamic_cast<AliAODJet*> (Jets[1]->At(imax));
0bf6381d 2206 if(jetMCL){
2207 fanalJetSubStrMCHM->AddEvent(aodE, jetMCL, rho_MC*areaJL_MC);
2208 }
89605a8b 2209 }
c87f1b5e 2210
2211
2212
2213
2214
2215 PostData(1, fOutputList);
89605a8b 2216
2217 if(fDebug) Printf("End of event #%5d", (Int_t) fEntry);
2218
c87f1b5e 2219 return;
2220
2221
2222
2223
2224}
2225
89605a8b 2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
c87f1b5e 2238Double_t AliAnalysisTaskJetShape::CalcdPhi(Double_t phi1, Double_t phi2)
2239{
2240 while(phi1<0) phi1+=TMath::TwoPi();
2241 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
2242
2243 while(phi2<0) phi2+=TMath::TwoPi();
2244 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
2245
2246 Double_t dphi = phi1- phi2;
2247 if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
2248 if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
2249
2250 // Double_t dphi = phi1- phi2;
2251 // while(dphi<0) dphi+=TMath::TwoPi();
2252 // while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2253 return dphi;
2254}
2255
2256Double_t AliAnalysisTaskJetShape::CalcdPhi0To2pi(Double_t phi1, Double_t phi2)
2257{
2258
2259 Double_t dphi = CalcdPhi(phi1, phi2);
2260 while(dphi<0) dphi+=TMath::TwoPi();
2261 while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2262 return dphi;
2263}
2264
2265
2266Bool_t AliAnalysisTaskJetShape::IsGoodEvent()
2267{
2268
2269 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
2270 // if (!fESD) {
2271 // AliError("ESD not available");
2272 // fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
2273 // fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
2274
2275 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
2276 // assume that the AOD is in the general output...
2277 fAODOut = AODEvent();
2278
2279
2280 static AliAODEvent* aod = 0;
2281 // take all other information from the aod we take the tracks from
2282 if(!aod){
2283 if(!fESD)aod = fAODIn;
2284 else aod = fAODOut;}
2285
2286
2287
2288 if(fNonStdFile.Length()!=0){
2289 // case that we have an AOD extension we need can fetch the jets from the extended output
2290 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2291 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
2292 }
2293
2294 if(fDebug){
2295 if(fAODIn) printf("fAODIn\n");
2296 if(fAODOut) printf("fAODOut\n");
2297 if(fAODExtension) printf("fAODExtension\n");
2298 }
2299
89605a8b 2300
2301 /*
2302 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2303if(!aodH){
2304 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
2305 return kFALSE;
2306 }
2307 */
2308
c87f1b5e 2309 // -- event selection --
2310
2311 // physics selection
2312 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
2313 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
89605a8b 2314 // cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
c87f1b5e 2315 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
2316 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
2317 return kFALSE;
2318 }
2319
2320 // vertex selection
2321 // if(!aod){
2322 // if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
2323 // fhNEvent->Fill(3);
2324 // PostData(1, fOutputList);
2325 // return kFALSE;
2326 // }
2327
2328 AliAODVertex* primVtx = aod->GetPrimaryVertex();
2329
2330 if(!primVtx){
2331 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
2332 return kFALSE;
2333 }
2334
2335 Int_t nTracksPrim = primVtx->GetNContributors();
2336 if ((nTracksPrim < fVtxMinContrib) ||
2337 (primVtx->GetZ() < fVtxZMin) ||
2338 (primVtx->GetZ() > fVtxZMax) ){
2339 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
2340 return kFALSE;
2341 }
2342
2343 /*
2344 // event class selection (from jet helper task)
2345 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
2346 if(fDebug) Printf("Event class %d", eventClass);
2347 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
2348 return kFALSE;
2349 }
2350 */
2351 return kTRUE;
2352}
2353
2354
2355
2356
2357
2358
2359//________________________________________________________________________
2360void AliAnalysisTaskJetShape::Terminate(Option_t *)
2361{
2362 // Draw result to the screen
2363 // Called once at the end of the query
2364 // fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
2365
2366 printf("MyTaskTestTrigger::Terminate()\n\n\n");
2367
2368
2369 // fOutputList = dynamic_cast<TList*> (GetOutputData(1));
2370 // fOutputTree = dynamic_cast<TTree*> (GetOutputData(1));
2371
2372}
2373
2374
2375
2376/*
2377Bool_t AnalysisJetMain::GetVertex(const AliESDEvent* esd, Double_t Vxyz[3], Int_t type)
2378{
2379 // type
2380 // 0 -> SPD vtx
2381 // 1 -> ESD vtx
2382 const AliESDVertex* vtx = 0;
2383
2384 if(type==0) {
2385 vtx = esd->GetPrimaryVertexSPD();
2386 if(!vtx) return kFALSE;
2387 if(vtx->GetNContributors() < 1) return kFALSE;
2388 TString vtxTyp = vtx->GetTitle();
2389 if (vtxTyp.Contains("vertexer: Z")) {
2390 if (vtx->GetDispersion()>0.04) return kFALSE;
2391 if (vtx->GetZRes()>0.25) return kFALSE;
2392 }
2393
2394 }
2395 if(type==1) {
2396 vtx = esd->GetPrimaryVertexTracks();
2397 if(!vtx) return kFALSE;
2398 if(vtx->GetNContributors() < 1) return kFALSE;
2399 }
2400
2401
2402 Vxyz[0] = vtx->GetXv();
2403 Vxyz[1] = vtx->GetYv();
2404 Vxyz[2] = vtx->GetZv();
2405
2406 return kTRUE;
2407}
2408*/
2409
2410
2411
2412
2413
2414
2415
2416
2417
2418TH1F* AliAnalysisTaskJetShape::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax, const char* xLabel, Int_t color)
2419{
2420// create a 1D histogram
2421
2422 TH1F *h = new TH1F(name, name, nBins, xMin, xMax);
2423 if (xLabel) h->GetXaxis()->SetTitle(xLabel);
2424 // res->SetMarkerStyle(kFullCircle);
2425 // res->SetOption("E");
2426 h->SetLineColor(color);
2427 // fOutputList->Add(h);
2428 return h;
2429}
2430
2431
2432TH2F *AliAnalysisTaskJetShape::Hist2D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax, Int_t nBinsy, Double_t yMin, Double_t yMax, const char* xLabel, const char* yLabel, Int_t color)
2433{
2434// create a 2D histogram
2435
2436 TH2F *res = new TH2F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax);
2437 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2438 if (xLabel) res->GetYaxis()->SetTitle(yLabel);
2439 // res->SetMarkerStyle(kFullCircle);
2440 // res->SetOption("E");
2441 res->SetLineColor(color);
2442 // fOutputList->Add(res);
2443 return res;
2444}
2445
2446
2447TH3F *AliAnalysisTaskJetShape::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
2448 Int_t nBinsy, Double_t yMin, Double_t yMax,
2449 Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
2450{
2451// create a 3D histogram
2452
2453 TH3F *res = new TH3F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax, nBinsz, zMin, zMax);
2454 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2455 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
2456 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
2457 // res->SetMarkerStyle(kFullCircle);
2458 // res->SetOption("E");
2459 res->SetLineColor(color);
2460 // fOutputList->Add(res);
2461 return res;
2462}
2463
2464
2465
2466
2467//#endif
2468
2469
2470
2471
2472
2473
2474
2475