]> 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)");
09c95e16 1252 fhPsiVsR = Hist3D("hPsiVsR", fPsiVsRNbin, 0., fRmax, 100, 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 {
09c95e16 1274 fhPsiVsR_MCtr = Hist3D("hPsiVsR_MCtr", fPsiVsRNbin, 0., fRmax, 100, 0., 1., fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{Jt, frac}", "P_{J} GeV/c");
1275 // fhPsiVsR_MCtr = Hist1D("hPsiVsR_MCtr", fPsiVsRNbin, 0., fRmax, "R", 1, "#Psi(R)");
89605a8b 1276 fhPsiVsRPtJ_MCtr = Hist2D("hPsiVsRPtJ_MCtr", fPsiVsRNbin, 0., fRmax, fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{tJ} GeV/c", 1, "#Psi(R)");
1277 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 1278 const char *ch[2]={"m","p"};
89605a8b 1279 for(Int_t i=0; i<3; i++)
1280 {
1281 for(Int_t j=0; j<2; j++)
1282 {
1283 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" );
1284 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" );
1285 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}" );
1286 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 1287 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]");
1288 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 1289 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");
1290 }
1291 }
1292 }
c87f1b5e 1293
1294}
1295
1296
1297
1298
1299
1300void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddToList(TList *list)
1301{
1302 list->Add(fhEvent);
1303 list->Add(fhMult);
1304 list->Add(fhPhiEtaTrack);
1305 list->Add(fhPtJ);
1306 list->Add(fhPsiVsR);
1307 list->Add(fhPsiVsRPtJ);
1308 list->Add(fhPtJvsPtCorr);
1309
89605a8b 1310
c87f1b5e 1311 for(Int_t i=0; i<3; i++)
1312 {
1313 list->Add(fhTMA_JAA[i]);
1314 list->Add(fhTMA_JAp[i]);
1315
1316 list->Add(fhTMA_B1AA[i]);
1317 list->Add(fhTMA_B1Ap[i]);
1318
1319 list->Add(fhTMA_B2AA[i]);
1320 list->Add(fhTMA_B2Ap[i]);
1321 }
1322
89605a8b 1323 if(fkMCprod)
1324 {
1325 list->Add(fhPsiVsR_MCtr);
1326 list->Add(fhPsiVsRPtJ_MCtr);
1327 list->Add(fhJetTrPtVsPartPt);
c87f1b5e 1328
89605a8b 1329 for(Int_t i=0; i<3; i++)
1330 {
1331 for(Int_t j=0; j<2; j++)
1332 {
1333 list->Add(fhPtresVsPt[i][j]);
1334 list->Add(fhPhiresVsPhi[i][j]);
1335 list->Add(fhEtaresVsEta[i][j]);
1336 list->Add(fhRresVsPt[i][j]);
1337 list->Add(fhDCAxy[i][j]);
1338 list->Add(fhDCAz[i][j]);
1339 list->Add(fhTrackPtEtaPhi[i][j]);
1340 }
1341 }
1342 }
c87f1b5e 1343}
1344
1345
1346
1347
1348Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddEvent(AliAODEvent* aodE, AliAODJet *jet, Double_t DeltaPtJ)
1349{
1350
1351
1352
1353 Int_t size = aodE->GetNumberOfTracks();
1354
1355 TClonesArray *arrP = 0;
1356
89605a8b 1357 if(fkMC || fkMCprod)
c87f1b5e 1358 {
1359 arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
1360 if(!arrP)
1361 {
1362 printf("ERROR: no Info about particles in AOD\n");
1363 return kFALSE;
1364 }
c87f1b5e 1365 }
1366
89605a8b 1367 if(fkMC)
1368 size = arrP->GetEntriesFast();
c87f1b5e 1369
1370 Int_t IndexArray[size];
89605a8b 1371 Int_t IndexArrayMC[size];
c87f1b5e 1372
1373 TClonesArray farray("TVector3", size);
1374
1375 Int_t counter = 0;
89605a8b 1376 Int_t counterMC = 0;
c87f1b5e 1377 Int_t counterAll = 0;
1378 Double_t pJ[3] = {0, 0, 0};
89605a8b 1379 Double_t pJmc[3] = {0, 0, 0};
c87f1b5e 1380
89605a8b 1381 AliAODVertex* primVtx = aodE->GetPrimaryVertex();
1382 Double_t bfield = aodE->GetMagneticField();
1383 Double_t dca[2];
1384 Double_t cov[3];
c87f1b5e 1385
1386 TVector3 vecJ(jet->Px(),jet->Py(),jet->Pz());
1387
c87f1b5e 1388
c560b734 1389
1390
1391 for(int it = 0;it < size;++it){
c87f1b5e 1392 TVector3 vec;
89605a8b 1393 Int_t ch = -999;
c560b734 1394 Int_t label = 0;
c87f1b5e 1395
89605a8b 1396
c87f1b5e 1397 if(fkMC)
1398 {
1399 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(it));
1400 if(!part)continue;
1401 if(!part->IsPhysicalPrimary())continue;
1402 if(part->Charge()==0)continue;
1403 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1404 }
1405 else
1406 {
1407 AliAODTrack *tr = aodE->GetTrack(it);
1408 if(!tr) continue;
1409 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
89605a8b 1410 label = tr->GetLabel();
1411 AliAODTrack tmp(*tr);
1412 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
c87f1b5e 1413 vec.SetXYZ(tr->Px(), tr->Py(), tr->Pz());
89605a8b 1414 ch = tr->Charge();
c87f1b5e 1415 }
1416
1417
1418 if(TMath::Abs(vec.Eta())>fEtaTrackMax)continue;
1419 if(vec.Pt()< fPtTrackMin)continue;
1420 if(vec.Pt()> fPtTrackMax) {return kFALSE;}
1421
1422
1423 new(farray[counterAll]) TVector3(vec);
1424 counterAll++;
1425
1426
1427 Double_t R = CalcR(vecJ, vec);
1428 if(R> fRmax) continue;
1429
1430 pJ[0]+=vec.Px();
1431 pJ[1]+=vec.Py();
1432 pJ[2]+=vec.Pz();
1433
1434 IndexArray[counter] = it;
1435 counter++;
89605a8b 1436
1437
1438
1439 // effics
1440 if(fkMCprod)
1441 {
1442 // printf("A1\n");
1443 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(label)));
1444 if(!part)continue;
1445
1446 Int_t type = 0;
1447 if(!part->IsPhysicalPrimary()) type=1;
1448 if(label <0) type=2;
1449
1450 if(!(ch==-1 || ch==1)) AliFatal("ch != +/- 1!!!\n\n");
1451 // printf("A4\n");
1452 if(ch==-1) ch=0;
1453
1454 Double_t phip = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(part->Phi());
1455 Double_t dphi = AliAnalysisTaskJetShapeTool::CalcdPhiSigned(part->Phi(), vec.Phi());
c560b734 1456 Double_t phiT = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
89605a8b 1457
1458 fhPtresVsPt[type][ch]->Fill(part->Pt(), 1-vec.Pt()/part->Pt());
1459 fhPhiresVsPhi[type][ch]->Fill(phip, dphi);
1460 fhEtaresVsEta[type][ch]->Fill(part->Eta(), vec.Eta() - part->Eta());
1461 fhDCAxy[type][ch]->Fill(part->Pt(), dca[0]);
1462 fhDCAz[type][ch]->Fill( part->Pt(), dca[1]);
c560b734 1463 fhTrackPtEtaPhi[type][ch]->Fill(vec.Pt(), vec.Eta(), phiT);
89605a8b 1464
1465 TVector3 vecP(part->Px(), part->Py(), part->Pz());
1466 Double_t Rgen = CalcR(vecJ, vecP);
1467 fhRresVsPt[type][ch]->Fill(part->Pt(), Rgen - R);
1468
1469 pJmc[0]+=part->Px();
1470 pJmc[1]+=part->Py();
1471 pJmc[2]+=part->Pz();
1472
1473 IndexArrayMC[counterMC] = label;
1474 counterMC++;
1475 // printf("A5\n");
1476 }
1477
c87f1b5e 1478 }
1479
1480 fhMult->Fill(counter);
1481 if(counter<1) return kFALSE;
1482
1483 fJvec.SetXYZ(pJ[0], pJ[1], pJ[2]);
1484
89605a8b 1485
c87f1b5e 1486 fPtJ = TMath::Sqrt(pJ[0]*pJ[0] + pJ[1]*pJ[1]);
1487 if((fPtJ < fPtJmin) || (fPtJ > fPtJmax)) return kFALSE;
1488 fhPtJ->Fill(fPtJ);
1489
89605a8b 1490
1491
1492
1493
1494 fhPtJvsPtCorr->Fill(fPtJ, jet->Pt() - DeltaPtJ);
c87f1b5e 1495
1496 for(Int_t it = 0; it<counter; it++)
1497 {
1498 TVector3 vec;
1499
1500 if(fkMC)
1501 {
1502 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(IndexArray[it]));
89605a8b 1503 if(!part) continue;
1504 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
c87f1b5e 1505 }
1506 else
1507 {
1508 AliAODTrack *tr = aodE->GetTrack(IndexArray[it]);
89605a8b 1509 if(!tr) continue;
1510 AliAODTrack tmp(*tr);
1511 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1512 if(tr)vec.SetXYZ(tmp.Px(), tmp.Py(), tmp.Pz());
c87f1b5e 1513 }
1514
1515 Double_t R = CalcR(fJvec, vec);
1516 // Double_t pt = (tr->Px()*pJ[0] + tr->Py()*pJ[1])/PtJ;
c560b734 1517 Double_t probL = fJvec.Dot(vec)/fJvec.Mag2();
89605a8b 1518 // fhPsiVsR->Fill(R, probL);
1519 // fhPsiVsRPtJ->Fill(R, fPtJ, probL);
c560b734 1520 fhPsiVsR->Fill(R,probL, fJvec.Mag());
89605a8b 1521 fhPsiVsRPtJ->Fill(R, fPtJ);
c87f1b5e 1522
89605a8b 1523 Double_t phi = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
c87f1b5e 1524 fhPhiEtaTrack->Fill(phi, vec.Eta());
1525
89605a8b 1526 }
1527
1528
1529 if(fkMCprod)
1530 {
1531 TVector3 fJvecMCtr(pJmc[0], pJmc[1], pJmc[2]);
c560b734 1532 Double_t fPtJMCtr= fJvecMCtr.Perp();
89605a8b 1533
c560b734 1534 fhJetTrPtVsPartPt->Fill(fPtJMCtr, 1-fPtJ/fPtJMCtr);
89605a8b 1535 for(Int_t it = 0; it<counterMC; it++)
1536 {
1537 TVector3 vec;
c87f1b5e 1538
c560b734 1539 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(IndexArrayMC[it])));
89605a8b 1540 if(!part) continue;
1541 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1542
1543 Double_t R = CalcR(fJvecMCtr, vec);
09c95e16 1544 Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
89605a8b 1545 // fhPsiVsR->Fill(R, probL);
1546 // fhPsiVsRPtJ->Fill(R, fPtJ, probL);
09c95e16 1547 //Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1548 fhPsiVsR_MCtr->Fill(R,probL, fJvecMCtr.Mag());
89605a8b 1549 fhPsiVsRPtJ_MCtr->Fill(R, fPtJMCtr);
1550 }
1551 }
c87f1b5e 1552
1553 // ang. distr.
1554 // fMyTool->Clean();
1555 AliAnalysisTaskJetShapeTool *fMyTool = new AliAnalysisTaskJetShapeTool(&farray);
1556
1557 fMyTool->SetNEntries(counterAll);
1558 fMyTool->SetVecJ(vecJ);
1559 fMyTool->SetPtMinTr(fPtTrackMin,fPtTrackMin);
1560 fMyTool->Init();
1561 Int_t scenario = 3;
1562
1563 /*
1564 to be added!!!!!!!!!
1565 fhTMA_B1AA[i]=0;
1566 fhTMA_B2AA[i]=0;
1567 */
1568
1569 for(Int_t l=0; l<3; l++)
1570 {
1571 Double_t phiA, phi1;
1572
1573 // Double_t ptRJ0= fMyTool->GetPRecInRJ().Pt();
1574 // Double_t ptRJ = fMyTool->GetPRecJ(l,0).Pt();
1575 Int_t N1 = fMyTool->GetSizeJ(l,1);
1576 Double_t dphi = -999;
1577
1578
1579 if(!fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 0) , vecJ, phiA, scenario))
1580 continue;
1581
1582 // f2JEvent->SetPhiJL(l,0, phiA);
1583
1584 if(fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 1) , vecJ, phi1, scenario))
1585 {
1586 fhTMA_JAA[l]->Fill(fMyTool->CalcdPhi0To2pi(phiA-phi1));
1587 }
1588
1589
1590 for(Int_t ip =0; ip<N1; ip++)
1591 {
1592 phi1 = fMyTool->GetLocPhiJ(l,1,ip);
1593 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1594 fhTMA_JAp[l]->Fill(dphi);
1595 }
1596
1597
1598 Int_t NB1 = fMyTool->GetSizeB1(l, 1);
1599 for(Int_t ip =0; ip<NB1; ip++)
1600 {
1601 phi1 = fMyTool->GetLocPhiB1(l,1,ip);
1602 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1603 fhTMA_B1Ap[l]->Fill(dphi);
1604 }
1605
1606 Int_t NB2 = fMyTool->GetSizeB2(l, 1);
1607 for(Int_t ip =0; ip<NB2; ip++)
1608 {
1609 phi1 = fMyTool->GetLocPhiB2(l,1,ip);
1610 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1611 fhTMA_B2Ap[l]->Fill(dphi);
1612 }
1613
1614 }
1615
1616 fhEvent->Fill(1);
1617
1618 delete fMyTool;
1619
1620 return kTRUE;
1621}
1622
1623
1624Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcR(TVector3 v1, TVector3 v2)
1625{
1626
1627 Double_t dphi = CalcdPhi(v1.Phi(), v2.Phi());
1628 Double_t deta = v1.Eta() - v2.Eta();
1629 Double_t RB = TMath::Sqrt(dphi*dphi+deta*deta);
1630 return RB;
1631}
1632
1633
1634Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcdPhi(Double_t phi1, Double_t phi2)
1635{
1636
1637 while(phi1<0) phi1+=TMath::TwoPi();
1638 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
1639
1640 while(phi2<0) phi2+=TMath::TwoPi();
1641 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
1642
1643 Double_t dphi = phi1- phi2;
1644 if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
1645 if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
1646
1647 return dphi;
1648}
1649
1650
1651
1652
89605a8b 1653
1654
1655
c87f1b5e 1656TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax, const char* xLabel, Int_t color, const char* yLabel)
1657{
1658// create a 1D histogram
1659
1660 TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xMin, xMax);
1661 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1662 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1663 res->SetLineColor(color);
1664 return res;
1665}
1666
1667
1668TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t *xArray, const char* xLabel, Int_t color, const char* yLabel)
1669{
1670// create a 1D histogram
1671
1672 TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xArray);
1673 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1674 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1675 res->SetLineColor(color);
1676 return res;
1677}
1678
1679TH2F *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)
1680{
1681// create a 2D histogram
1682
1683 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yMin, yMax);
1684 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1685 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1686 // res->SetMarkerStyle(kFullCircle);
1687 // res->SetOption("E");
1688 res->SetLineColor(color);
1689 // fOutputList->Add(res);
1690 return res;
1691}
1692
1693
1694TH2F *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)
1695{
1696// create a 2D histogram
1697
1698 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yArray);
1699 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1700 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1701 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1702 // res->SetMarkerStyle(kFullCircle);
1703 // res->SetOption("E");
1704 res->SetLineColor(color);
1705 // fOutputList->Add(res);
1706 return res;
1707}
1708
1709TH2F *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)
1710{
1711// create a 2D histogram
1712
1713 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, yArrax, nBinsy, yMin, yMax);
1714 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1715 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1716 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1717 // res->SetMarkerStyle(kFullCircle);
1718 // res->SetOption("E");
1719 res->SetLineColor(color);
1720 // fOutputList->Add(res);
1721 return res;
1722}
1723
89605a8b 1724TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1725 Int_t nBinsy, Double_t yMin, Double_t yMax,
1726 Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1727{
1728// create a 3D histogram
1729
1730 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);
1731 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1732 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1733 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1734 // res->SetMarkerStyle(kFullCircle);
c560b734 1735 // res->SetOption("E");
1736 res->SetLineColor(color);
1737 // fOutputList->Add(res);
1738 return res;
1739}
1740
1741TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1742 Int_t nBinsy, Double_t yMin, Double_t yMax,
1743 Int_t nBinsz, Double_t *zArr, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1744{
1745// create a 3D histogram
1746
1747 Double_t xArr[nBinsx+1], yArr[nBinsy+1];
1748
1749 for(Int_t i=0; i<=nBinsx; i++) xArr[i]=xMin+i*(xMax-xMin)/nBinsx;
1750 for(Int_t i=0; i<=nBinsy; i++) yArr[i]=yMin+i*(yMax-yMin)/nBinsy;
1751
1752 TH3F *res = new TH3F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xArr, nBinsy, yArr, nBinsz, zArr);
1753 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1754 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1755 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1756 // res->SetMarkerStyle(kFullCircle);
89605a8b 1757 // res->SetOption("E");
1758 res->SetLineColor(color);
1759 // fOutputList->Add(res);
1760 return res;
1761}
c87f1b5e 1762
1763
1764//////////////////////////////////////
1765
1766
1767/*
1768//________________________________________________________________________
1769void AnalysisJetMain::ConnectInputData(Option_t *)
1770{
1771 // Connect ESD
1772 // Called once
1773
1774 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
1775 // if (!fESD) {
1776 // AliError("ESD not available");
1777 // fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
1778 // fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1779
1780 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1781 // assume that the AOD is in the general output...
1782 fAODOut = AODEvent();
1783
1784
1785}
1786*/
1787
1788void AliAnalysisTaskJetShape::SetBranchNames(const TString &branch1, const TString &branch2)
1789{
1790 fJetBranchName[0] = branch1;
1791 fJetBranchName[1] = branch2;
1792}
1793
1794//________________________________________________________________________
1795void AliAnalysisTaskJetShape::UserCreateOutputObjects()
1796{
1797 //printf("Open1\n");
1798 // const char *nameF = OpenFile(1)->GetName();
1799 //printf("Open2 %s\n",nameF);
1800
1801 // fTriggerAnalysis = new AliTriggerAnalysis();
1802 // if (fkMC) fTriggerAnalysis->SetAnalyzeMC(1);
1803
1804
1805 fOutputList = new TList();
1806 fOutputList->SetOwner();
1807
1808 fhPtJL = Hist1D("hPtJL", 100, 0, 200, "Pt_{JL}"); fOutputList->Add(fhPtJL);
1809 fhAreaJL = Hist1D("hAreaJL", 100, 0., 4, "AreaJL"); fOutputList->Add(fhAreaJL);
1810
89605a8b 1811
1812
1813 printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() \n");
1814
1815
c87f1b5e 1816 fanalJetSubStrHM->SetFilterMask(fFilterMask);
89605a8b 1817 if(fkMC) fanalJetSubStrHM->MCprod();
c87f1b5e 1818 fanalJetSubStrHM->InitHistos();
1819 fanalJetSubStrHM->AddToList(fOutputList);
1820
c87f1b5e 1821
1822 if(fkMC)
1823 {
1824
1825 printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() MC\n");
1826 fanalJetSubStrMCHM->InitHistos();
1827 fanalJetSubStrMCHM->AddToList(fOutputList);
89605a8b 1828
1829 for(Int_t i=0; i<3; i++)
1830 {
1831 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]);
1832 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]);
1833 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]);
1834 fhDCAxy[i] = Hist2D(Form("hDCAxy%d",i), 100, 0, 100, 300, -3, 3, "DCAxy of prim [cm]"); fOutputList->Add(fhDCAxy[i]);
1835 fhDCAz[i] = Hist2D(Form("hDCAz%d",i) , 100, 0, 100, 300, -3, 3, "DCAz of prim [cm]") ; fOutputList->Add(fhDCAz[i]);
1836 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]);
1837 }
c87f1b5e 1838 }
1839
1840
1841
1842
89605a8b 1843 printf(" JetBranchName[0]=%s\n JetBranchName[1]=%s\n", fJetBranchName[0].Data(), fJetBranchName[1].Data());
1844
1845
1846
c87f1b5e 1847 PostData(1, fOutputList);
1848
1849 /*
1850 OpenFile(1);
1851
1852 fOutputTree = new TTree("tree","Tree with Ali2JEvent");
1853 fOutputTree->Branch("event",&f2JEvent);
1854 */
1855
1856// PostData(1, fOutputTree);
1857}
1858
1859/*
1860Bool_t AliAnalysisTaskJetSpectrum2::Notify()
1861{
1862
1863 // Fetch the aod also from the input in,
1864 // have todo it in notify
1865
1866
1867 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1868 // assume that the AOD is in the general output...
1869 fAODOut = AODEvent();
1870
1871 if(fNonStdFile.Length()!=0){
1872 // case that we have an AOD extension we need can fetch the jets from the extended output
1873 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1874 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
1875 if(!fAODExtension){
1876 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
1877 }
1878 }
1879
1880
1881
1882 return kTRUE;
1883}
1884*/
1885
1886
1887
1888
1889//________________________________________________________________________
1890void AliAnalysisTaskJetShape::UserExec(Option_t *)
1891{
1892 // return;
1893// if(f2JEvent)
1894// delete f2JEvent;
1895// f2JEvent = new Ali2JEvent();
1896// return;
1897
1898 // return;
89605a8b 1899 if(fDebug) Printf("\n\n\nEvent #%5d", (Int_t) fEntry);
c87f1b5e 1900 if(fDebug) printf("NEW EVENT 0\n");
1901
1902 if(!IsGoodEvent()) return;
1903
c87f1b5e 1904
1905 // printf("\n\n\n NEW EVENT\n");
c87f1b5e 1906
1907 AliAODEvent* aodE = 0;
1908 if(!aodE){
1909 if(!fESD)aodE = fAODIn;
1910 else aodE = fAODOut;}
1911
89605a8b 1912 /*
1913 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1914if(!aodH){
1915 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
1916 return;
1917 }
1918 */
1919 if(fDebug) {
1920 printf("NEW EVENT 1: Number of tracks = %d\n",aodE->GetNumberOfTracks());
1921 // aodE->GetList()->Print();
1922 // printf("Look for %s\n",fJetBranchName[0].Data());
1923 }
c87f1b5e 1924
1925
1926 // centrality selection
1927 AliCentrality *cent = 0x0;
1928 Double_t centrality = 0.;
1929
1930 if(fESD) {cent = fESD->GetCentrality();
1931 if(cent) centrality = cent->GetCentralityPercentile("V0M");}
1932 else centrality=aodE->GetHeader()->GetCentrality();
1933
1934
1935 if(!fkIsPbPb) {
1936 centrality=1.;
1937 }
1938
89605a8b 1939// if(fDebug) printf("NEW EVENT 2\n");
c87f1b5e 1940
1941 Bool_t fUseAOD049 = kFALSE;// kTRUE;//
1942 if(fUseAOD049&&centrality>=0){
1943 Float_t v0=0;
1944 AliAODVZERO* aodV0 = aodE->GetVZEROData();
1945 v0+=aodV0->GetMTotV0A();
1946 v0+=aodV0->GetMTotV0C();
1947 if(centrality==0 && v0 < 19500) return ;//filtering issue
1948 Float_t tkl = (Float_t)(aodE->GetTracklets()->GetNumberOfTracklets());
1949 Float_t val= 1.30552 + 0.147931 * v0;
1950 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,
1951 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,
1952 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,
1953 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,
1954 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,
1955 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,
1956 13.7544, 13.7544, 13.7544, 13.7544, 13.7544};
1957 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)centrality] )
1958 return; //outlier
1959 }
1960
1961
89605a8b 1962 if(fDebug)
1963 {
1964 printf("centrality: %f\n", centrality);
1965 aodE->Print();
1966 aodE->GetList()->Print();
1967 }
1968
1969
c87f1b5e 1970 if (centrality < fCentMin || centrality > fCentMax){
1971 // PostData(1, fOutputList);
1972 return;
1973 }
1974
1975 // fhNEvent->Fill(0);
1976
1977 // accepted events
1978 // -- end event selection --
1979
c87f1b5e 1980
1981 // get background
1982 AliAODJetEventBackground* externalBackground = 0;
1983 AliAODJetEventBackground* externalBackgroundMC = 0;
1984
1985
1986 Float_t rho = 0;
1987 Float_t rho_MC=0;
1988
1989 if(fkIsPbPb)
1990 {
1991 if(fAODOut&&!externalBackground&&fBackgroundBranch[0].Length()){
1992 externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[0].Data()));
1993 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1994 }
1995 if(fAODExtension&&!externalBackground&&fBackgroundBranch[0].Length()){
1996 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[0].Data()));
1997 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1998 }
1999 if(fAODIn&&!externalBackground&&fBackgroundBranch[0].Length()){
2000 externalBackground = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[0].Data()));
2001 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
2002 }
2003 if(fAODOut&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2004 externalBackgroundMC = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[1].Data()));
2005 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2006 }
2007 if(fAODExtension&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2008 externalBackgroundMC = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[1].Data()));
2009 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2010 }
2011 if(fAODIn&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2012 externalBackgroundMC = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[1].Data()));
2013 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2014 }
2015 if(externalBackground)rho = externalBackground->GetBackground(0);
2016 if(externalBackgroundMC)rho_MC = externalBackgroundMC->GetBackground(0);
2017 }
2018
2019
2020
89605a8b 2021 if(fkMC)
2022 {
2023 AliAODVertex* primVtx = aodE->GetPrimaryVertex();
2024 Double_t bfield = aodE->GetMagneticField();
2025
2026 TClonesArray *arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
2027 if(!arrP)
2028 {
2029 printf("ERROR: no Info about particles in AOD\n");
2030 return;
2031 }
2032
2033 for(int it = 0;it < aodE->GetNumberOfTracks(); it++)
2034 {
2035 AliAODTrack *tr = aodE->GetTrack(it);
2036 if(!tr) continue;
2037 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))) continue;
2038 if(TMath::Abs(tr->Eta())>1.) continue;
2039 Int_t label = TMath::Abs(tr->GetLabel());
2040
2041 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(label));
2042 if(!part)continue;
2043 // if(!part->IsPhysicalPrimary())continue;
2044 // if(part->Charge()==0)continue;
2045 // vec.SetXYZ(part->Px(), part->Py(), part->Pz());
2046
2047
2048
2049 Double_t dca[2];
2050 Double_t cov[3];
2051
2052 AliAODTrack tmp(*tr);
2053 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
2054
2055 Int_t type = 0;
2056 if(!part->IsPhysicalPrimary()) type=1;
2057 if(label<0) type =2;
2058
2059 fhPtresVsPt[type]->Fill(part->Pt(), 1-tr->Pt()/part->Pt());
2060 fhPhiresVsPhi[type]->Fill(part->Phi(), tr->Phi() - part->Phi());
2061 fhEtaresVsEta[type]->Fill(part->Eta(), tr->Eta() - part->Eta());
2062 fhDCAxy[type]->Fill(part->Pt(), dca[0]);
2063 fhDCAz[type]->Fill( part->Pt(), dca[1]);
2064 fhTrackPtEtaPhi[type]->Fill(tr->Pt(), tr->Eta(), tr->Phi());
2065 }
2066
2067
2068 }
2069
2070
2071
2072
c87f1b5e 2073
2074 // printf("rho = %f\n",rho);
2075
2076 // return;
2077
2078 TClonesArray *Jets[2];
2079 Jets[0]=0;
2080 Jets[1]=0;
2081 if(fAODOut&&!Jets[0]){
2082 Jets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
2083 Jets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); }
2084 if(fAODExtension && !Jets[0]){
2085 Jets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
2086 Jets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
2087 if(fAODIn&&!Jets[0]){
2088 Jets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
2089 Jets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
2090
2091
2092 if(!Jets[0]) {
2093 Printf("no friend!!!\n");
2094 return;
2095 }
2096 Int_t nJ = Jets[0]->GetEntries();
2097 Float_t ptmax = 0.;
2098 Int_t imax = -1;
2099
89605a8b 2100 if(fDebug) {
2101 printf("NEW EVENT 3: Number of Rec. jets %d\n",nJ);
2102 }
c87f1b5e 2103
2104//
2105// Find highest pT jet with pt > 20 GeV
2106//
89605a8b 2107// fPtJcorrMin=0;
c87f1b5e 2108
2109 Float_t etaJmax = 0.4;
2110
2111 for (Int_t i = 0; i < nJ; i++) {
2112 AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[0]->At(i));
2113 if (!jet) continue;
89605a8b 2114 // jet->Print("0");
c87f1b5e 2115 Float_t ptJ = jet->Pt();
2116 Float_t etaJ = TMath::Abs(jet->Eta());
2117
2118 Float_t area = jet->EffectiveAreaCharged();
2119 Float_t ptJcorr=ptJ-rho*area;
2120
2121 if ((ptJcorr > fPtJcorrMin) && (ptJcorr > ptmax) && etaJ < etaJmax) {
2122 ptmax = ptJcorr;
2123 imax = i;
2124 }
2125
2126 }
2127
89605a8b 2128
2129 if(fDebug) {
2130 printf("NEW EVENT 4:\n");
2131 }
2132
c87f1b5e 2133 TVector3 vecJ;
2134
2135 AliAODJet* jetL = 0;
2136
2137 Float_t areaJL, ptJLcorr;
2138
2139 if (imax != -1) {
2140
0bf6381d 2141 jetL = dynamic_cast<AliAODJet*> (Jets[0]->At(imax));
2142 if(jetL){
c87f1b5e 2143
0bf6381d 2144 areaJL = jetL->EffectiveAreaCharged();
2145 ptJLcorr = jetL->Pt()-rho*areaJL;
2146
2147 fhPtJL->Fill(ptJLcorr);
2148 fhAreaJL->Fill(areaJL);
2149 vecJ.SetXYZ(jetL->Px(), jetL->Py(), jetL->Pz());
2150 fanalJetSubStrHM->AddEvent(aodE, jetL, rho*areaJL);
89605a8b 2151
2152 if(fDebug) {
2153 Printf("Leading Jet");
2154 jetL->Print("0");
2155 }
2156
0bf6381d 2157 }
c87f1b5e 2158 }
2159
2160
c87f1b5e 2161 if(!fkMC)
2162 {
2163 PostData(1, fOutputList);
89605a8b 2164 if(fDebug) Printf("End of event #%5d", (Int_t) fEntry);
c87f1b5e 2165 return;
2166 }
2167
2168
89605a8b 2169
2170
c87f1b5e 2171 nJ = Jets[1]->GetEntries();
89605a8b 2172 if(fDebug) {
2173 printf("NEW EVENT 5: Number of Rec. jets %d\n",nJ);
2174 }
2175
c87f1b5e 2176 ptmax = 0;
2177 imax = -1;
2178 Float_t areaJL_MC=0;
2179
2180 for (Int_t i = 0; i < nJ; i++) {
2181 AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[1]->At(i));
2182 if (!jet) continue;
2183 Float_t ptJ1 = jet->Pt();
2184 Float_t etaJ1 = TMath::Abs(jet->Eta());
2185
2186 areaJL_MC = jet->EffectiveAreaCharged();
2187 Float_t ptJcorr=ptJ1-rho*areaJL_MC;
2188
2189 // jet->Print();
2190
2191 if ((ptJcorr > fPtJcorrMin) && (ptJcorr > ptmax) && etaJ1 < etaJmax) {
2192 ptmax = ptJcorr;
2193 imax = i;
2194 }
2195
2196 }
2197
89605a8b 2198 if(fDebug) {
2199 printf("NEW EVENT 6:\n");
2200 }
2201
c87f1b5e 2202
2203 AliAODJet* jetMCL=0;
2204
2205 if (imax != -1) {
2206 jetMCL = dynamic_cast<AliAODJet*> (Jets[1]->At(imax));
0bf6381d 2207 if(jetMCL){
2208 fanalJetSubStrMCHM->AddEvent(aodE, jetMCL, rho_MC*areaJL_MC);
2209 }
89605a8b 2210 }
c87f1b5e 2211
2212
2213
2214
2215
2216 PostData(1, fOutputList);
89605a8b 2217
2218 if(fDebug) Printf("End of event #%5d", (Int_t) fEntry);
2219
c87f1b5e 2220 return;
2221
2222
2223
2224
2225}
2226
89605a8b 2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
c87f1b5e 2239Double_t AliAnalysisTaskJetShape::CalcdPhi(Double_t phi1, Double_t phi2)
2240{
2241 while(phi1<0) phi1+=TMath::TwoPi();
2242 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
2243
2244 while(phi2<0) phi2+=TMath::TwoPi();
2245 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
2246
2247 Double_t dphi = phi1- phi2;
2248 if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
2249 if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
2250
2251 // Double_t dphi = phi1- phi2;
2252 // while(dphi<0) dphi+=TMath::TwoPi();
2253 // while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2254 return dphi;
2255}
2256
2257Double_t AliAnalysisTaskJetShape::CalcdPhi0To2pi(Double_t phi1, Double_t phi2)
2258{
2259
2260 Double_t dphi = CalcdPhi(phi1, phi2);
2261 while(dphi<0) dphi+=TMath::TwoPi();
2262 while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2263 return dphi;
2264}
2265
2266
2267Bool_t AliAnalysisTaskJetShape::IsGoodEvent()
2268{
2269
2270 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
2271 // if (!fESD) {
2272 // AliError("ESD not available");
2273 // fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
2274 // fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
2275
2276 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
2277 // assume that the AOD is in the general output...
2278 fAODOut = AODEvent();
2279
2280
2281 static AliAODEvent* aod = 0;
2282 // take all other information from the aod we take the tracks from
2283 if(!aod){
2284 if(!fESD)aod = fAODIn;
2285 else aod = fAODOut;}
2286
2287
2288
2289 if(fNonStdFile.Length()!=0){
2290 // case that we have an AOD extension we need can fetch the jets from the extended output
2291 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2292 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
2293 }
2294
2295 if(fDebug){
2296 if(fAODIn) printf("fAODIn\n");
2297 if(fAODOut) printf("fAODOut\n");
2298 if(fAODExtension) printf("fAODExtension\n");
2299 }
2300
89605a8b 2301
2302 /*
2303 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2304if(!aodH){
2305 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
2306 return kFALSE;
2307 }
2308 */
2309
c87f1b5e 2310 // -- event selection --
2311
2312 // physics selection
2313 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
2314 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
89605a8b 2315 // cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
c87f1b5e 2316 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
2317 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
2318 return kFALSE;
2319 }
2320
2321 // vertex selection
2322 // if(!aod){
2323 // if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
2324 // fhNEvent->Fill(3);
2325 // PostData(1, fOutputList);
2326 // return kFALSE;
2327 // }
2328
2329 AliAODVertex* primVtx = aod->GetPrimaryVertex();
2330
2331 if(!primVtx){
2332 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
2333 return kFALSE;
2334 }
2335
2336 Int_t nTracksPrim = primVtx->GetNContributors();
2337 if ((nTracksPrim < fVtxMinContrib) ||
2338 (primVtx->GetZ() < fVtxZMin) ||
2339 (primVtx->GetZ() > fVtxZMax) ){
2340 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
2341 return kFALSE;
2342 }
2343
2344 /*
2345 // event class selection (from jet helper task)
2346 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
2347 if(fDebug) Printf("Event class %d", eventClass);
2348 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
2349 return kFALSE;
2350 }
2351 */
2352 return kTRUE;
2353}
2354
2355
2356
2357
2358
2359
2360//________________________________________________________________________
2361void AliAnalysisTaskJetShape::Terminate(Option_t *)
2362{
2363 // Draw result to the screen
2364 // Called once at the end of the query
2365 // fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
2366
2367 printf("MyTaskTestTrigger::Terminate()\n\n\n");
2368
2369
2370 // fOutputList = dynamic_cast<TList*> (GetOutputData(1));
2371 // fOutputTree = dynamic_cast<TTree*> (GetOutputData(1));
2372
2373}
2374
2375
2376
2377/*
2378Bool_t AnalysisJetMain::GetVertex(const AliESDEvent* esd, Double_t Vxyz[3], Int_t type)
2379{
2380 // type
2381 // 0 -> SPD vtx
2382 // 1 -> ESD vtx
2383 const AliESDVertex* vtx = 0;
2384
2385 if(type==0) {
2386 vtx = esd->GetPrimaryVertexSPD();
2387 if(!vtx) return kFALSE;
2388 if(vtx->GetNContributors() < 1) return kFALSE;
2389 TString vtxTyp = vtx->GetTitle();
2390 if (vtxTyp.Contains("vertexer: Z")) {
2391 if (vtx->GetDispersion()>0.04) return kFALSE;
2392 if (vtx->GetZRes()>0.25) return kFALSE;
2393 }
2394
2395 }
2396 if(type==1) {
2397 vtx = esd->GetPrimaryVertexTracks();
2398 if(!vtx) return kFALSE;
2399 if(vtx->GetNContributors() < 1) return kFALSE;
2400 }
2401
2402
2403 Vxyz[0] = vtx->GetXv();
2404 Vxyz[1] = vtx->GetYv();
2405 Vxyz[2] = vtx->GetZv();
2406
2407 return kTRUE;
2408}
2409*/
2410
2411
2412
2413
2414
2415
2416
2417
2418
2419TH1F* AliAnalysisTaskJetShape::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax, const char* xLabel, Int_t color)
2420{
2421// create a 1D histogram
2422
2423 TH1F *h = new TH1F(name, name, nBins, xMin, xMax);
2424 if (xLabel) h->GetXaxis()->SetTitle(xLabel);
2425 // res->SetMarkerStyle(kFullCircle);
2426 // res->SetOption("E");
2427 h->SetLineColor(color);
2428 // fOutputList->Add(h);
2429 return h;
2430}
2431
2432
2433TH2F *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)
2434{
2435// create a 2D histogram
2436
2437 TH2F *res = new TH2F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax);
2438 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2439 if (xLabel) res->GetYaxis()->SetTitle(yLabel);
2440 // res->SetMarkerStyle(kFullCircle);
2441 // res->SetOption("E");
2442 res->SetLineColor(color);
2443 // fOutputList->Add(res);
2444 return res;
2445}
2446
2447
2448TH3F *AliAnalysisTaskJetShape::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
2449 Int_t nBinsy, Double_t yMin, Double_t yMax,
2450 Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
2451{
2452// create a 3D histogram
2453
2454 TH3F *res = new TH3F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax, nBinsz, zMin, zMax);
2455 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2456 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
2457 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
2458 // res->SetMarkerStyle(kFullCircle);
2459 // res->SetOption("E");
2460 res->SetLineColor(color);
2461 // fOutputList->Add(res);
2462 return res;
2463}
2464
2465
2466
2467
2468//#endif
2469
2470
2471
2472
2473
2474
2475
2476