1 // Class AliAnalysisTaskJetShape
2 // Author martin.poghosyan@cern.ch
12 #include "TParticlePDG.h"
13 #include "TParticle.h"
15 #include "TDatabasePDG.h"
17 #include "THnSparse.h"
22 #include <TClonesArray.h>
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"
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"
63 #include "AliAnalysisTaskJetShape.h"
65 /////////////////////////////////////
67 ClassImp(AliAnalysisTaskJetShape)
70 AliAnalysisTaskJetShape::AliAnalysisTaskJetShape(const char *name)
71 : AliAnalysisTaskSE(name),
84 fOfflineTrgMask(AliVEvent::kAny),
100 fanalJetSubStrMCHM(0)
104 fgkbinCent[0] = -0.099;
113 fBackgroundBranch[0] = "",
114 fBackgroundBranch[1] = "",
115 fJetBranchName[0] = "";
116 fJetBranchName[1] = "";
118 for(Int_t i=0; i<3; i++)
125 fhTrackPtEtaPhi[i]=0;
127 // fkIsPbPb = kFALSE;
132 fanalJetSubStrHM = new AliAnalysisTaskJetShapeHM("rec", kFALSE);
133 fanalJetSubStrMCHM = new AliAnalysisTaskJetShapeHM("truth", kTRUE);
135 DefineOutput(1, TList::Class());
137 // DefineOutput(1, TTree::Class());
141 AliAnalysisTaskJetShape::~AliAnalysisTaskJetShape()
143 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
144 printf("Deleteing output\n");
151 // if(fTriggerAnalysis)
152 // delete fTriggerAnalysis;
158 AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::AliAnalysisTaskJetShapeTool():TObject(),
167 fList = new TClonesArray("TVector3",10000);
177 for(Int_t i1=0; i1<fgkbinR; i1++)
179 for(Int_t i2=0; i2<2; i2++)
181 fListJ[i1][i2].Set(1000);
182 fListB1[i1][i2].Set(1000);
183 fListB2[i1][i2].Set(1000);
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);
191 fListJ[i1][i2].Reset(-1);
192 fListB1[i1][i2].Reset(-1);
193 fListB2[i1][i2].Reset(-1);
195 fPRecJ[i1][i2].SetXYZ(0,0,0);
200 fPRecInRJ.SetXYZ(0,0,0);
208 AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::AliAnalysisTaskJetShapeTool(TClonesArray *list):TObject(),
226 for(Int_t i1=0; i1<fgkbinR; i1++)
228 for(Int_t i2=0; i2<2; i2++)
230 fListJ[i1][i2].Set(1000);
231 fListB1[i1][i2].Set(1000);
232 fListB2[i1][i2].Set(1000);
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);
247 fPRecInRJ.SetXYZ(0,0,0);
257 AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::~AliAnalysisTaskJetShapeTool()
268 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::SetVecJ(TVector3 v)
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());
276 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::Clean()
278 // printf("AnalJetSubStrTool::Clean()\n");
281 for(Int_t i1=0; i1<fgkbinR; i1++)
283 for(Int_t i2=0; i2<2; i2++)
285 fListJ[i1][i2]->Delete();
286 fListB1[i1][i2]->Delete();
287 fListB2[i1][i2]->Delete();
296 fPRecInRJ.SetXYZ(0,0,0);
298 for(Int_t i1=0; i1<fgkbinR; i1++)
300 for(Int_t i2=0; i2<2; i2++)
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);
309 fListB1c[i1][i2] = 0;
310 fListB2c[i1][i2] = 0;
312 fPRecJ[i1][i2].SetXYZ(0,0,0);
319 //void AnalJetSubStrTool::Print(Option_t* /*option*/) const
320 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::PRINT() const
323 // PRINT(fList, "all");
326 for(Int_t i1=0; i1<fgkbinR; i1++)
330 for(Int_t i2=0; i2<2; i2++)
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();
345 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::PRINT(TArrayI a, Int_t n, const char *txt) const
347 printf("%s :%d \n",txt, n);
349 for(Int_t i=0; i<n; i++){
350 printf("%4d ", a.At(i));
359 if(count!=20) printf("\n");
365 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::Init()
367 Int_t Nev = fEntries;
371 // printf("Nev = %d\n",Nev);
373 for(Int_t iP=0; iP<Nev; iP++)
375 TVector3 *v = (TVector3*) fList->At(iP);
377 printf("ERROR : entry #%d not found\n",iP);
381 // printf("#%3d ",iP); v->Print();
383 Double_t R = CalcR(fvecJ, *v);
384 // printf("R = %f\n",R);
387 for(Int_t iT = 0; iT < fgkbinR; iT++)
390 if(R>fR[iT]) type = 1;
392 if(v->Pt() < fPtMinTr[type]) continue;
393 fPRecJ[iT][type]+=*v;
394 AddToJ(iT, type, iP);
397 if(v->Pt() < fPtMinTr[0]) continue;
403 R = CalcR(fvecB1, *v);
406 for(Int_t iT = 0; iT < fgkbinR; iT++)
409 if(R>fR[iT]) type = 1;
411 if(v->Pt() < fPtMinTr[type]) continue;
413 AddToB1(iT, type, iP);
418 R = CalcR(fvecB2, *v);
421 for(Int_t iT = 0; iT < fgkbinR; iT++)
424 if(R>fR[iT]) type = 1;
426 if(v->Pt() < fPtMinTr[type]) continue;
428 AddToB2(iT, type, iP);
437 for(Int_t i1=0; i1<fgkbinR; i1++)
439 for(Int_t i2=0; i2<2; i2++)
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]);
455 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::CalcR(TVector3 v1, TVector3 v2)
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);
466 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::CalcdPhi(Double_t phi1, Double_t phi2)
469 while(phi1<0) phi1+=TMath::TwoPi();
470 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
472 while(phi2<0) phi2+=TMath::TwoPi();
473 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
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();
483 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::GetLocalMom(TVector3 vecJ, TVector3 *q)
485 // theta and phi of a particle in loc.syst of fvecJ
488 Double_t PT = vecJ.Perp();
489 Double_t P0 = vecJ.Mag();
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);
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);
504 q->SetXYZ(qx1, qy1, qz1);
512 Bool_t AnalJetSubStrTool::FindCorrelationAxesAnd(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario)
514 Double_t TwoPi = TMath::TwoPi();
517 // 1st track loop to determine the sphericity matrix
519 // Initialize Shericity matrix
532 Int_t Nev = list.GetSize();
533 if( Nev <2) return kFALSE;
535 Int_t indexpmax = -1;
537 Double_t phipmax = 0;
539 Int_t indexpzmax = -1;
541 Double_t phipzmax = 0;
544 Int_t indexthetamin = -1;
545 Double_t thetamin = 1000;
546 Double_t phithetamin = 0;
548 for(Int_t ip=0; ip< Nev; ip++) {
550 TVector3* part = (TVector3*) GetAt(list.At(ip));
553 TVector3 vtmp(*part);
554 GetLocalMom(vec, &vtmp);
556 Float_t ppjX = vtmp.X();
557 Float_t ppjY = vtmp.Y();
558 Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
561 Float_t pt = ppjT;//part->Z();
562 Float_t deta = part->Eta();
563 Float_t dphi = part->Phi();
565 mxx += (ppjX * ppjX / ppjT);
566 myy += (ppjY * ppjY / ppjT);
567 mxy += (ppjX * ppjY / ppjT);
585 if(vtmp.Theta()<thetamin)
587 thetamin=vtmp.Theta();
589 phithetamin=vtmp.Phi();
599 // At this point we have mxx, myy, mxy
600 if (nc == 0) return kFALSE;
602 const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
603 TMatrixDSym m0(2, ele);
605 TMatrixDSymEigen m(m0);
607 TMatrixD evecm = m.GetEigenVectors();
608 eval = m.GetEigenValues();
609 // Largest eigenvector
611 if (eval[0] < eval[1]) jev = 1;
614 evec0 = TMatrixDColumn(evecm, jev);
615 TVector2 evec(evec0[0], evec0[1]);
616 // Principle axis from leading partice
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.);
624 Double_t phiTA = evec.Phi();
632 else if(scenario ==1)
637 else if(scenario ==2)
639 phiDir = phithetamin;
647 if( TMath::Abs(CalcdPhi(phiDir, phiTA)) <TMath::Pi()/2)
651 if(Phi>TwoPi) Phi-=TwoPi;
660 Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::FindCorrelationAxes(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario)
662 // scenario 0- Pmax, 1-Pzmax_Local, 2->cosMin, 3-PTmac, 4->PTmax_Local
664 Double_t xmean , ymean;
666 Double_t TwoPi = TMath::TwoPi();
679 // Double_t phiLoc, thetaLoc;
683 Int_t Nev = list.GetSize();
684 if( Nev <2) return kFALSE;
691 if(scenario==2) val = 100;
692 Double_t phiDir = -99;
696 for(Int_t ip=0; ip< Nev; ip++) {
697 if(list.At(ip)<0) break;
699 TVector3* part = (TVector3*) GetAt(list.At(ip));
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);
711 xx = CalcdPhi(part->Phi(), vec.Phi());
712 yy = part->Eta() - vec.Eta();
731 else if(scenario ==1)
740 else if(scenario ==2)
757 // phiDir=vtmp.Phi();
758 phiDir=TMath::ATan2(yy, xx); //vtmp.Phi();
765 if(N<2) return kFALSE;
777 Double_t B = x2-x*x+y*y-y2;
778 Double_t D = TMath::Sqrt(B*B+4*A*A);
780 // printf("N = %d A = %f\n",N, A);
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);
791 while(phi>TwoPi) phi-=TwoPi;
792 while(phi< 0 ) phi+=TwoPi;
794 Phi=CalcdPhi0To2pi(phi, 0);
796 if(scenario!=4 && Index<0) return kFALSE;
800 if( TMath::Abs(CalcdPhi(phiDir, phi)) <TMath::Pi()/2)
805 Phi = CalcdPhi0To2pi(phi);
811 Double_t xmax = -100;
818 for(Int_t ip=0; ip< Nev; ip++) {
819 if(list.At(ip)<0) break;
821 TVector3* part = (TVector3*) GetAt(list.At(ip));
824 TVector3 vtmp(*part);
825 GetLocalMom(vec, &vtmp);
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);
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;
836 if(TMath::Abs(xmin) > xmax) Phi+=TMath::Pi();
841 while(Phi>TwoPi) Phi-=TwoPi;
849 Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeTool::FindCorrelationAxesCorr(TArrayI list, TVector3 vec, Double_t &Phi, Int_t scenario, Double_t R)
851 // scenario 0- Pmax, 1-Pzmax_Local, 2->cosMin, 3-PTmac, 4->PTmax_Local
853 Double_t xmean , ymean;
855 Double_t TwoPi = TMath::TwoPi();
868 // Double_t phiLoc, thetaLoc;
872 Int_t Nev = list.GetSize();
873 if( Nev <2) return kFALSE;
880 if(scenario==2) val = 100;
881 Double_t phiDir = -99;
885 for(Int_t ip=0; ip< Nev; ip++) {
886 if(list.At(ip)<0) break;
888 TVector3* part = (TVector3*) GetAt(list.At(ip));
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);
899 xx = CalcdPhi(part->Phi(), vec.Phi());
900 yy = part->Eta() - vec.Eta();
919 else if(scenario ==1)
928 else if(scenario ==2)
945 // phiDir=vtmp.Phi();
946 phiDir=TMath::ATan2(yy, xx); //vtmp.Phi();
953 if(N<2) return kFALSE;
954 if(scenario!=4 && Index<0) return kFALSE;
966 Double_t B = x2-x*x+y*y-y2;
967 Double_t D = TMath::Sqrt(B*B+4*A*A);
969 // printf("N = %d A = %f\n",N, A);
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);
980 Phi=CalcdPhi0To2pi(phi, 0);
981 if( TMath::Abs(CalcdPhi(phiDir, phi)) > TMath::Pi()/2)
982 Phi=CalcdPhi0To2pi(phi+TMath::Pi(), 0);
994 Double_t phiMinus = 0;
996 for(Int_t ip=0; ip< Nev; ip++)
998 if(list.At(ip)<0) break;
1000 TVector3* part = (TVector3*) GetAt(list.At(ip));
1002 Double_t xx0 = CalcdPhi(part->Phi(), vec.Phi());
1003 Double_t yy0 = part->Eta() - vec.Eta();
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);
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);
1018 if( TMath::Abs(CalcdPhi(phiDir, phi1)) > TMath::Pi()/2)
1021 phiMinus+=CalcdPhi(Phi, phi1);
1027 Double_t phiPlus = 0;
1029 for(Int_t ip=0; ip< -Nev; ip++)
1031 if(list.At(ip)<0) break;
1033 TVector3* part = (TVector3*) GetAt(list.At(ip));
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);
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);
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);
1053 if( TMath::Abs(CalcdPhi(phiDir, phi1)) > TMath::Pi()/2)
1056 phiPlus+=CalcdPhi(Phi, phi1);
1062 Phi+=((phiMinus+phiPlus)/(Nminus+Nplus+1));
1064 if( TMath::Abs(CalcdPhi(phiDir, Phi)) > TMath::Pi()/2)
1067 Phi=CalcdPhi0To2pi(Phi, 0);
1073 Double_t xmax = -100;
1074 Double_t xmin = 100;
1080 for(Int_t ip=0; ip< Nev; ip++) {
1081 if(list.At(ip)<0) break;
1083 TVector3* part = (TVector3*) GetAt(list.At(ip));
1086 TVector3 vtmp(*part);
1087 GetLocalMom(vec, &vtmp);
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);
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;
1098 if(TMath::Abs(xmin) > xmax) Phi+=TMath::Pi();
1103 while(Phi>TwoPi) Phi-=TwoPi;
1109 /////////////////////////////////////
1111 TH2F *fhPhiEtaTrack;//!
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];//!
1122 TArrayD fPtJetArray;
1125 Int_t fPsiVsPhiNbin;
1128 Double_t fEtaTrackMax;
1129 Double_t fPtTrackMin;
1130 Double_t fPtTrackMax;
1138 AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AliAnalysisTaskJetShapeHM(TString comment, Bool_t kMC):TObject(),
1150 fhPsiVsRPtJ_MCtr(0),
1151 fhJetTrPtVsPartPt(0),
1167 for(Int_t i=0; i<3; i++)
1177 for(Int_t i=0; i<3; i++)
1179 for(Int_t j=0; j<2; j++)
1181 fhPtresVsPt[i][j]=0;
1182 fhPhiresVsPhi[i][j]=0;
1183 fhEtaresVsEta[i][j]=0;
1187 fhTrackPtEtaPhi[i][j]=0;
1202 fJvec.SetXYZ(0,0,0);
1207 AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::~AliAnalysisTaskJetShapeHM()
1211 delete [] fPtJetArray;
1217 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::SetPtJetBins(Int_t Nbin, Double_t *array)
1222 delete [] fPtJetArray;
1226 fPtJetArray = new Double_t[fPtJetNbin+1];
1227 for(Int_t i=0; i<= fPtJetNbin; i++) fPtJetArray[i]= array[i];
1232 fPtJetArray.Set(fPtJetNbin+1, array);
1236 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::InitHistos()
1239 fhEvent = Hist1D("hEvent" , 3 , -0.5, 2.5, "Event");
1240 fhMult = Hist1D("hMult" , 101, -0.5, 100.5, "N_{ch}");
1242 fhPhiEtaTrack = Hist2D("hPhiEtaTrack", 100, 0, TMath::TwoPi(), 20, -1, 1., "#phi", "#eta");
1246 Double_t array[]= {0., 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 130, 160, 200};
1247 SetPtJetBins(Nbin, array);
1250 fhPtJ = Hist1D("hPtJ" , fPtJetNbin, fPtJetArray.GetArray(), "Pt_{J} GeV/c");
1251 // fhPsiVsR = Hist1D("hPsiVsR", fPsiVsRNbin, 0., fRmax, "R", 1, "#Psi(R)");
1252 fhPsiVsR = Hist3D("hPsiVsR", fPsiVsRNbin, 0., fRmax, 100, 0., 1., fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{Jt, frac}", "P_{J} GeV/c");
1253 fhPsiVsRPtJ = Hist2D("hPsiVsRPtJ", fPsiVsRNbin, 0., fRmax, fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{tJ} GeV/c", 1, "#Psi(R)");
1255 fhPtJvsPtCorr = Hist2D("fhPtJvsPtCorr", fPtJetNbin, fPtJetArray.GetArray(), 100, -100, 100, "P_{tJ} GeV/c", "P_{tJ} - P_{tJ,corr} GeV/c" , 1);
1260 for(Int_t i=0; i<3; i++)
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})");
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})");
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})");
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)");
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);
1278 const char *ch[2]={"m","p"};
1279 for(Int_t i=0; i<3; i++)
1281 for(Int_t j=0; j<2; j++)
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}" );
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]") ;
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");
1300 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddToList(TList *list)
1304 list->Add(fhPhiEtaTrack);
1306 list->Add(fhPsiVsR);
1307 list->Add(fhPsiVsRPtJ);
1308 list->Add(fhPtJvsPtCorr);
1311 for(Int_t i=0; i<3; i++)
1313 list->Add(fhTMA_JAA[i]);
1314 list->Add(fhTMA_JAp[i]);
1316 list->Add(fhTMA_B1AA[i]);
1317 list->Add(fhTMA_B1Ap[i]);
1319 list->Add(fhTMA_B2AA[i]);
1320 list->Add(fhTMA_B2Ap[i]);
1325 list->Add(fhPsiVsR_MCtr);
1326 list->Add(fhPsiVsRPtJ_MCtr);
1327 list->Add(fhJetTrPtVsPartPt);
1329 for(Int_t i=0; i<3; i++)
1331 for(Int_t j=0; j<2; j++)
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]);
1348 Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddEvent(AliAODEvent* aodE, AliAODJet *jet, Double_t DeltaPtJ)
1353 Int_t size = aodE->GetNumberOfTracks();
1355 TClonesArray *arrP = 0;
1357 if(fkMC || fkMCprod)
1359 arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
1362 printf("ERROR: no Info about particles in AOD\n");
1368 size = arrP->GetEntriesFast();
1370 Int_t IndexArray[size];
1371 Int_t IndexArrayMC[size];
1373 TClonesArray farray("TVector3", size);
1376 Int_t counterMC = 0;
1377 Int_t counterAll = 0;
1378 Double_t pJ[3] = {0, 0, 0};
1379 Double_t pJmc[3] = {0, 0, 0};
1381 AliAODVertex* primVtx = aodE->GetPrimaryVertex();
1382 Double_t bfield = aodE->GetMagneticField();
1383 Double_t dca[2] = {0., 0.};
1384 Double_t cov[3] = {0., 0., 0.};
1386 TVector3 vecJ(jet->Px(),jet->Py(),jet->Pz());
1391 for(int it = 0;it < size;++it){
1399 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(it));
1401 if(!part->IsPhysicalPrimary())continue;
1402 if(part->Charge()==0)continue;
1403 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1407 AliAODTrack *tr = aodE->GetTrack(it);
1409 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1410 label = tr->GetLabel();
1411 AliAODTrack tmp(*tr);
1412 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1413 vec.SetXYZ(tr->Px(), tr->Py(), tr->Pz());
1418 if(TMath::Abs(vec.Eta())>fEtaTrackMax)continue;
1419 if(vec.Pt()< fPtTrackMin)continue;
1420 if(vec.Pt()> fPtTrackMax) {return kFALSE;}
1423 new(farray[counterAll]) TVector3(vec);
1427 Double_t R = CalcR(vecJ, vec);
1428 if(R> fRmax) continue;
1434 IndexArray[counter] = it;
1443 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(label)));
1447 if(!part->IsPhysicalPrimary()) type=1;
1448 if(label <0) type=2;
1450 if(!(ch==-1 || ch==1)) AliFatal("ch != +/- 1!!!\n\n");
1454 Double_t phip = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(part->Phi());
1455 Double_t dphi = AliAnalysisTaskJetShapeTool::CalcdPhiSigned(part->Phi(), vec.Phi());
1456 Double_t phiT = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
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]);
1463 fhTrackPtEtaPhi[type][ch]->Fill(vec.Pt(), vec.Eta(), phiT);
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);
1469 pJmc[0]+=part->Px();
1470 pJmc[1]+=part->Py();
1471 pJmc[2]+=part->Pz();
1473 IndexArrayMC[counterMC] = label;
1480 fhMult->Fill(counter);
1481 if(counter<1) return kFALSE;
1483 fJvec.SetXYZ(pJ[0], pJ[1], pJ[2]);
1486 fPtJ = TMath::Sqrt(pJ[0]*pJ[0] + pJ[1]*pJ[1]);
1487 if((fPtJ < fPtJmin) || (fPtJ > fPtJmax)) return kFALSE;
1494 fhPtJvsPtCorr->Fill(fPtJ, jet->Pt() - DeltaPtJ);
1496 for(Int_t it = 0; it<counter; it++)
1502 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(IndexArray[it]));
1504 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1508 AliAODTrack *tr = aodE->GetTrack(IndexArray[it]);
1510 AliAODTrack tmp(*tr);
1511 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1512 if(tr)vec.SetXYZ(tmp.Px(), tmp.Py(), tmp.Pz());
1515 Double_t R = CalcR(fJvec, vec);
1516 // Double_t pt = (tr->Px()*pJ[0] + tr->Py()*pJ[1])/PtJ;
1517 Double_t probL = fJvec.Dot(vec)/fJvec.Mag2();
1518 // fhPsiVsR->Fill(R, probL);
1519 // fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1520 fhPsiVsR->Fill(R,probL, fJvec.Mag());
1521 fhPsiVsRPtJ->Fill(R, fPtJ);
1523 Double_t phi = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
1524 fhPhiEtaTrack->Fill(phi, vec.Eta());
1531 TVector3 fJvecMCtr(pJmc[0], pJmc[1], pJmc[2]);
1532 Double_t fPtJMCtr= fJvecMCtr.Perp();
1534 fhJetTrPtVsPartPt->Fill(fPtJMCtr, 1-fPtJ/fPtJMCtr);
1535 for(Int_t it = 0; it<counterMC; it++)
1539 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(IndexArrayMC[it])));
1541 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1543 Double_t R = CalcR(fJvecMCtr, vec);
1544 Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1545 // fhPsiVsR->Fill(R, probL);
1546 // fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1547 //Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1548 fhPsiVsR_MCtr->Fill(R,probL, fJvecMCtr.Mag());
1549 fhPsiVsRPtJ_MCtr->Fill(R, fPtJMCtr);
1554 // fMyTool->Clean();
1555 AliAnalysisTaskJetShapeTool *fMyTool = new AliAnalysisTaskJetShapeTool(&farray);
1557 fMyTool->SetNEntries(counterAll);
1558 fMyTool->SetVecJ(vecJ);
1559 fMyTool->SetPtMinTr(fPtTrackMin,fPtTrackMin);
1564 to be added!!!!!!!!!
1569 for(Int_t l=0; l<3; l++)
1571 Double_t phiA, phi1;
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;
1579 if(!fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 0) , vecJ, phiA, scenario))
1582 // f2JEvent->SetPhiJL(l,0, phiA);
1584 if(fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 1) , vecJ, phi1, scenario))
1586 fhTMA_JAA[l]->Fill(fMyTool->CalcdPhi0To2pi(phiA-phi1));
1590 for(Int_t ip =0; ip<N1; ip++)
1592 phi1 = fMyTool->GetLocPhiJ(l,1,ip);
1593 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1594 fhTMA_JAp[l]->Fill(dphi);
1598 Int_t NB1 = fMyTool->GetSizeB1(l, 1);
1599 for(Int_t ip =0; ip<NB1; ip++)
1601 phi1 = fMyTool->GetLocPhiB1(l,1,ip);
1602 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1603 fhTMA_B1Ap[l]->Fill(dphi);
1606 Int_t NB2 = fMyTool->GetSizeB2(l, 1);
1607 for(Int_t ip =0; ip<NB2; ip++)
1609 phi1 = fMyTool->GetLocPhiB2(l,1,ip);
1610 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1611 fhTMA_B2Ap[l]->Fill(dphi);
1624 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcR(TVector3 v1, TVector3 v2)
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);
1634 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcdPhi(Double_t phi1, Double_t phi2)
1637 while(phi1<0) phi1+=TMath::TwoPi();
1638 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
1640 while(phi2<0) phi2+=TMath::TwoPi();
1641 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
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();
1656 TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax, const char* xLabel, Int_t color, const char* yLabel)
1658 // create a 1D histogram
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);
1668 TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t *xArray, const char* xLabel, Int_t color, const char* yLabel)
1670 // create a 1D histogram
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);
1679 TH2F *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)
1681 // create a 2D histogram
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);
1694 TH2F *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)
1696 // create a 2D histogram
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);
1709 TH2F *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)
1711 // create a 2D histogram
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);
1724 TH3F *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)
1728 // create a 3D histogram
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);
1735 // res->SetOption("E");
1736 res->SetLineColor(color);
1737 // fOutputList->Add(res);
1741 TH3F *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)
1745 // create a 3D histogram
1747 Double_t xArr[nBinsx+1], yArr[nBinsy+1];
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;
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);
1757 // res->SetOption("E");
1758 res->SetLineColor(color);
1759 // fOutputList->Add(res);
1764 //////////////////////////////////////
1768 //________________________________________________________________________
1769 void AnalysisJetMain::ConnectInputData(Option_t *)
1774 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
1776 // AliError("ESD not available");
1777 // fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
1778 // fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1780 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1781 // assume that the AOD is in the general output...
1782 fAODOut = AODEvent();
1788 void AliAnalysisTaskJetShape::SetBranchNames(const TString &branch1, const TString &branch2)
1790 fJetBranchName[0] = branch1;
1791 fJetBranchName[1] = branch2;
1794 //________________________________________________________________________
1795 void AliAnalysisTaskJetShape::UserCreateOutputObjects()
1797 //printf("Open1\n");
1798 // const char *nameF = OpenFile(1)->GetName();
1799 //printf("Open2 %s\n",nameF);
1801 // fTriggerAnalysis = new AliTriggerAnalysis();
1802 // if (fkMC) fTriggerAnalysis->SetAnalyzeMC(1);
1805 fOutputList = new TList();
1806 fOutputList->SetOwner();
1808 fhPtJL = Hist1D("hPtJL", 100, 0, 200, "Pt_{JL}"); fOutputList->Add(fhPtJL);
1809 fhAreaJL = Hist1D("hAreaJL", 100, 0., 4, "AreaJL"); fOutputList->Add(fhAreaJL);
1813 printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() \n");
1816 fanalJetSubStrHM->SetFilterMask(fFilterMask);
1817 if(fkMC) fanalJetSubStrHM->MCprod();
1818 fanalJetSubStrHM->InitHistos();
1819 fanalJetSubStrHM->AddToList(fOutputList);
1825 printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() MC\n");
1826 fanalJetSubStrMCHM->InitHistos();
1827 fanalJetSubStrMCHM->AddToList(fOutputList);
1829 for(Int_t i=0; i<3; i++)
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]);
1843 printf(" JetBranchName[0]=%s\n JetBranchName[1]=%s\n", fJetBranchName[0].Data(), fJetBranchName[1].Data());
1847 PostData(1, fOutputList);
1852 fOutputTree = new TTree("tree","Tree with Ali2JEvent");
1853 fOutputTree->Branch("event",&f2JEvent);
1856 // PostData(1, fOutputTree);
1860 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
1863 // Fetch the aod also from the input in,
1864 // have todo it in notify
1867 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1868 // assume that the AOD is in the general output...
1869 fAODOut = AODEvent();
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);
1876 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
1889 //________________________________________________________________________
1890 void AliAnalysisTaskJetShape::UserExec(Option_t *)
1895 // f2JEvent = new Ali2JEvent();
1899 if(fDebug) Printf("\n\n\nEvent #%5d", (Int_t) fEntry);
1900 if(fDebug) printf("NEW EVENT 0\n");
1902 if(!IsGoodEvent()) return;
1905 // printf("\n\n\n NEW EVENT\n");
1907 AliAODEvent* aodE = 0;
1909 if(!fESD)aodE = fAODIn;
1910 else aodE = fAODOut;}
1913 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1915 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
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());
1926 // centrality selection
1927 AliCentrality *cent = 0x0;
1928 Double_t centrality = 0.;
1930 if(fESD) {cent = fESD->GetCentrality();
1931 if(cent) centrality = cent->GetCentralityPercentile("V0M");}
1932 else centrality=aodE->GetHeader()->GetCentrality();
1939 // if(fDebug) printf("NEW EVENT 2\n");
1941 Bool_t fUseAOD049 = kFALSE;// kTRUE;//
1942 if(fUseAOD049&¢rality>=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] )
1964 printf("centrality: %f\n", centrality);
1966 aodE->GetList()->Print();
1970 if (centrality < fCentMin || centrality > fCentMax){
1971 // PostData(1, fOutputList);
1975 // fhNEvent->Fill(0);
1978 // -- end event selection --
1982 AliAODJetEventBackground* externalBackground = 0;
1983 AliAODJetEventBackground* externalBackgroundMC = 0;
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());
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());
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());
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());
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());
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());
2015 if(externalBackground)rho = externalBackground->GetBackground(0);
2016 if(externalBackgroundMC)rho_MC = externalBackgroundMC->GetBackground(0);
2023 AliAODVertex* primVtx = aodE->GetPrimaryVertex();
2024 Double_t bfield = aodE->GetMagneticField();
2026 TClonesArray *arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
2029 printf("ERROR: no Info about particles in AOD\n");
2033 for(int it = 0;it < aodE->GetNumberOfTracks(); it++)
2035 AliAODTrack *tr = aodE->GetTrack(it);
2037 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))) continue;
2038 if(TMath::Abs(tr->Eta())>1.) continue;
2039 Int_t label = TMath::Abs(tr->GetLabel());
2041 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(label));
2043 // if(!part->IsPhysicalPrimary())continue;
2044 // if(part->Charge()==0)continue;
2045 // vec.SetXYZ(part->Px(), part->Py(), part->Pz());
2052 AliAODTrack tmp(*tr);
2053 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
2056 if(!part->IsPhysicalPrimary()) type=1;
2057 if(label<0) type =2;
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());
2074 // printf("rho = %f\n",rho);
2078 TClonesArray *Jets[2];
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())); }
2093 Printf("no friend!!!\n");
2096 Int_t nJ = Jets[0]->GetEntries();
2101 printf("NEW EVENT 3: Number of Rec. jets %d\n",nJ);
2105 // Find highest pT jet with pt > 20 GeV
2109 Float_t etaJmax = 0.4;
2111 for (Int_t i = 0; i < nJ; i++) {
2112 AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[0]->At(i));
2115 Float_t ptJ = jet->Pt();
2116 Float_t etaJ = TMath::Abs(jet->Eta());
2118 Float_t area = jet->EffectiveAreaCharged();
2119 Float_t ptJcorr=ptJ-rho*area;
2121 if ((ptJcorr > fPtJcorrMin) && (ptJcorr > ptmax) && etaJ < etaJmax) {
2130 printf("NEW EVENT 4:\n");
2135 AliAODJet* jetL = 0;
2137 Float_t areaJL, ptJLcorr;
2141 jetL = dynamic_cast<AliAODJet*> (Jets[0]->At(imax));
2144 areaJL = jetL->EffectiveAreaCharged();
2145 ptJLcorr = jetL->Pt()-rho*areaJL;
2147 fhPtJL->Fill(ptJLcorr);
2148 fhAreaJL->Fill(areaJL);
2149 vecJ.SetXYZ(jetL->Px(), jetL->Py(), jetL->Pz());
2150 fanalJetSubStrHM->AddEvent(aodE, jetL, rho*areaJL);
2153 Printf("Leading Jet");
2163 PostData(1, fOutputList);
2164 if(fDebug) Printf("End of event #%5d", (Int_t) fEntry);
2171 nJ = Jets[1]->GetEntries();
2173 printf("NEW EVENT 5: Number of Rec. jets %d\n",nJ);
2178 Float_t areaJL_MC=0;
2180 for (Int_t i = 0; i < nJ; i++) {
2181 AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[1]->At(i));
2183 Float_t ptJ1 = jet->Pt();
2184 Float_t etaJ1 = TMath::Abs(jet->Eta());
2186 areaJL_MC = jet->EffectiveAreaCharged();
2187 Float_t ptJcorr=ptJ1-rho*areaJL_MC;
2191 if ((ptJcorr > fPtJcorrMin) && (ptJcorr > ptmax) && etaJ1 < etaJmax) {
2199 printf("NEW EVENT 6:\n");
2203 AliAODJet* jetMCL=0;
2206 jetMCL = dynamic_cast<AliAODJet*> (Jets[1]->At(imax));
2208 fanalJetSubStrMCHM->AddEvent(aodE, jetMCL, rho_MC*areaJL_MC);
2216 PostData(1, fOutputList);
2218 if(fDebug) Printf("End of event #%5d", (Int_t) fEntry);
2239 Double_t AliAnalysisTaskJetShape::CalcdPhi(Double_t phi1, Double_t phi2)
2241 while(phi1<0) phi1+=TMath::TwoPi();
2242 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
2244 while(phi2<0) phi2+=TMath::TwoPi();
2245 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
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();
2251 // Double_t dphi = phi1- phi2;
2252 // while(dphi<0) dphi+=TMath::TwoPi();
2253 // while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2257 Double_t AliAnalysisTaskJetShape::CalcdPhi0To2pi(Double_t phi1, Double_t phi2)
2260 Double_t dphi = CalcdPhi(phi1, phi2);
2261 while(dphi<0) dphi+=TMath::TwoPi();
2262 while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2267 Bool_t AliAnalysisTaskJetShape::IsGoodEvent()
2270 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
2272 // AliError("ESD not available");
2273 // fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
2274 // fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
2276 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
2277 // assume that the AOD is in the general output...
2278 fAODOut = AODEvent();
2281 static AliAODEvent* aod = 0;
2282 // take all other information from the aod we take the tracks from
2284 if(!fESD)aod = fAODIn;
2285 else aod = fAODOut;}
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);
2296 if(fAODIn) printf("fAODIn\n");
2297 if(fAODOut) printf("fAODOut\n");
2298 if(fAODExtension) printf("fAODExtension\n");
2303 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2305 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
2310 // -- event selection --
2312 // physics selection
2313 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
2314 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2315 // cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
2316 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
2317 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
2323 // if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
2324 // fhNEvent->Fill(3);
2325 // PostData(1, fOutputList);
2329 AliAODVertex* primVtx = aod->GetPrimaryVertex();
2332 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
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());
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){
2360 //________________________________________________________________________
2361 void AliAnalysisTaskJetShape::Terminate(Option_t *)
2363 // Draw result to the screen
2364 // Called once at the end of the query
2365 // fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
2367 printf("MyTaskTestTrigger::Terminate()\n\n\n");
2370 // fOutputList = dynamic_cast<TList*> (GetOutputData(1));
2371 // fOutputTree = dynamic_cast<TTree*> (GetOutputData(1));
2378 Bool_t AnalysisJetMain::GetVertex(const AliESDEvent* esd, Double_t Vxyz[3], Int_t type)
2383 const AliESDVertex* vtx = 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;
2397 vtx = esd->GetPrimaryVertexTracks();
2398 if(!vtx) return kFALSE;
2399 if(vtx->GetNContributors() < 1) return kFALSE;
2403 Vxyz[0] = vtx->GetXv();
2404 Vxyz[1] = vtx->GetYv();
2405 Vxyz[2] = vtx->GetZv();
2419 TH1F* AliAnalysisTaskJetShape::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax, const char* xLabel, Int_t color)
2421 // create a 1D histogram
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);
2433 TH2F *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)
2435 // create a 2D histogram
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);
2448 TH3F *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)
2452 // create a 3D histogram
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);