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 = dynamic_cast<AliAODTrack*>(aodE->GetTrack(it));
1408 if(!tr) AliFatal("Not a standard AOD");
1410 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1411 label = tr->GetLabel();
1412 AliAODTrack tmp(*tr);
1413 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1414 vec.SetXYZ(tr->Px(), tr->Py(), tr->Pz());
1419 if(TMath::Abs(vec.Eta())>fEtaTrackMax)continue;
1420 if(vec.Pt()< fPtTrackMin)continue;
1421 if(vec.Pt()> fPtTrackMax) {return kFALSE;}
1424 new(farray[counterAll]) TVector3(vec);
1428 Double_t R = CalcR(vecJ, vec);
1429 if(R> fRmax) continue;
1435 IndexArray[counter] = it;
1444 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(label)));
1448 if(!part->IsPhysicalPrimary()) type=1;
1449 if(label <0) type=2;
1451 if(!(ch==-1 || ch==1)) AliFatal("ch != +/- 1!!!\n\n");
1455 Double_t phip = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(part->Phi());
1456 Double_t dphi = AliAnalysisTaskJetShapeTool::CalcdPhiSigned(part->Phi(), vec.Phi());
1457 Double_t phiT = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
1459 fhPtresVsPt[type][ch]->Fill(part->Pt(), 1-vec.Pt()/part->Pt());
1460 fhPhiresVsPhi[type][ch]->Fill(phip, dphi);
1461 fhEtaresVsEta[type][ch]->Fill(part->Eta(), vec.Eta() - part->Eta());
1462 fhDCAxy[type][ch]->Fill(part->Pt(), dca[0]);
1463 fhDCAz[type][ch]->Fill( part->Pt(), dca[1]);
1464 fhTrackPtEtaPhi[type][ch]->Fill(vec.Pt(), vec.Eta(), phiT);
1466 TVector3 vecP(part->Px(), part->Py(), part->Pz());
1467 Double_t Rgen = CalcR(vecJ, vecP);
1468 fhRresVsPt[type][ch]->Fill(part->Pt(), Rgen - R);
1470 pJmc[0]+=part->Px();
1471 pJmc[1]+=part->Py();
1472 pJmc[2]+=part->Pz();
1474 IndexArrayMC[counterMC] = label;
1481 fhMult->Fill(counter);
1482 if(counter<1) return kFALSE;
1484 fJvec.SetXYZ(pJ[0], pJ[1], pJ[2]);
1487 fPtJ = TMath::Sqrt(pJ[0]*pJ[0] + pJ[1]*pJ[1]);
1488 if((fPtJ < fPtJmin) || (fPtJ > fPtJmax)) return kFALSE;
1495 fhPtJvsPtCorr->Fill(fPtJ, jet->Pt() - DeltaPtJ);
1497 for(Int_t it = 0; it<counter; it++)
1503 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(IndexArray[it]));
1505 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1509 AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodE->GetTrack(IndexArray[it]));
1510 if(!tr) AliFatal("Not a standard AOD");
1512 AliAODTrack tmp(*tr);
1513 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1514 if(tr)vec.SetXYZ(tmp.Px(), tmp.Py(), tmp.Pz());
1517 Double_t R = CalcR(fJvec, vec);
1518 // Double_t pt = (tr->Px()*pJ[0] + tr->Py()*pJ[1])/PtJ;
1519 Double_t probL = fJvec.Dot(vec)/fJvec.Mag2();
1520 // fhPsiVsR->Fill(R, probL);
1521 // fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1522 fhPsiVsR->Fill(R,probL, fJvec.Mag());
1523 fhPsiVsRPtJ->Fill(R, fPtJ);
1525 Double_t phi = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
1526 fhPhiEtaTrack->Fill(phi, vec.Eta());
1533 TVector3 fJvecMCtr(pJmc[0], pJmc[1], pJmc[2]);
1534 Double_t fPtJMCtr= fJvecMCtr.Perp();
1536 fhJetTrPtVsPartPt->Fill(fPtJMCtr, 1-fPtJ/fPtJMCtr);
1537 for(Int_t it = 0; it<counterMC; it++)
1541 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(IndexArrayMC[it])));
1543 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1545 Double_t R = CalcR(fJvecMCtr, vec);
1546 Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1547 // fhPsiVsR->Fill(R, probL);
1548 // fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1549 //Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1550 fhPsiVsR_MCtr->Fill(R,probL, fJvecMCtr.Mag());
1551 fhPsiVsRPtJ_MCtr->Fill(R, fPtJMCtr);
1556 // fMyTool->Clean();
1557 AliAnalysisTaskJetShapeTool *fMyTool = new AliAnalysisTaskJetShapeTool(&farray);
1559 fMyTool->SetNEntries(counterAll);
1560 fMyTool->SetVecJ(vecJ);
1561 fMyTool->SetPtMinTr(fPtTrackMin,fPtTrackMin);
1566 to be added!!!!!!!!!
1571 for(Int_t l=0; l<3; l++)
1573 Double_t phiA, phi1;
1575 // Double_t ptRJ0= fMyTool->GetPRecInRJ().Pt();
1576 // Double_t ptRJ = fMyTool->GetPRecJ(l,0).Pt();
1577 Int_t N1 = fMyTool->GetSizeJ(l,1);
1578 Double_t dphi = -999;
1581 if(!fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 0) , vecJ, phiA, scenario))
1584 // f2JEvent->SetPhiJL(l,0, phiA);
1586 if(fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 1) , vecJ, phi1, scenario))
1588 fhTMA_JAA[l]->Fill(fMyTool->CalcdPhi0To2pi(phiA-phi1));
1592 for(Int_t ip =0; ip<N1; ip++)
1594 phi1 = fMyTool->GetLocPhiJ(l,1,ip);
1595 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1596 fhTMA_JAp[l]->Fill(dphi);
1600 Int_t NB1 = fMyTool->GetSizeB1(l, 1);
1601 for(Int_t ip =0; ip<NB1; ip++)
1603 phi1 = fMyTool->GetLocPhiB1(l,1,ip);
1604 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1605 fhTMA_B1Ap[l]->Fill(dphi);
1608 Int_t NB2 = fMyTool->GetSizeB2(l, 1);
1609 for(Int_t ip =0; ip<NB2; ip++)
1611 phi1 = fMyTool->GetLocPhiB2(l,1,ip);
1612 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1613 fhTMA_B2Ap[l]->Fill(dphi);
1626 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcR(TVector3 v1, TVector3 v2)
1629 Double_t dphi = CalcdPhi(v1.Phi(), v2.Phi());
1630 Double_t deta = v1.Eta() - v2.Eta();
1631 Double_t RB = TMath::Sqrt(dphi*dphi+deta*deta);
1636 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcdPhi(Double_t phi1, Double_t phi2)
1639 while(phi1<0) phi1+=TMath::TwoPi();
1640 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
1642 while(phi2<0) phi2+=TMath::TwoPi();
1643 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
1645 Double_t dphi = phi1- phi2;
1646 if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
1647 if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
1658 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)
1660 // create a 1D histogram
1662 TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xMin, xMax);
1663 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1664 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1665 res->SetLineColor(color);
1670 TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t *xArray, const char* xLabel, Int_t color, const char* yLabel)
1672 // create a 1D histogram
1674 TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xArray);
1675 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1676 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1677 res->SetLineColor(color);
1681 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)
1683 // create a 2D histogram
1685 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yMin, yMax);
1686 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1687 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1688 // res->SetMarkerStyle(kFullCircle);
1689 // res->SetOption("E");
1690 res->SetLineColor(color);
1691 // fOutputList->Add(res);
1696 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)
1698 // create a 2D histogram
1700 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yArray);
1701 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1702 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1703 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1704 // res->SetMarkerStyle(kFullCircle);
1705 // res->SetOption("E");
1706 res->SetLineColor(color);
1707 // fOutputList->Add(res);
1711 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)
1713 // create a 2D histogram
1715 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, yArrax, nBinsy, yMin, yMax);
1716 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1717 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1718 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1719 // res->SetMarkerStyle(kFullCircle);
1720 // res->SetOption("E");
1721 res->SetLineColor(color);
1722 // fOutputList->Add(res);
1726 TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1727 Int_t nBinsy, Double_t yMin, Double_t yMax,
1728 Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1730 // create a 3D histogram
1732 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);
1733 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1734 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1735 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1736 // res->SetMarkerStyle(kFullCircle);
1737 // res->SetOption("E");
1738 res->SetLineColor(color);
1739 // fOutputList->Add(res);
1743 TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1744 Int_t nBinsy, Double_t yMin, Double_t yMax,
1745 Int_t nBinsz, Double_t *zArr, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1747 // create a 3D histogram
1749 Double_t xArr[nBinsx+1], yArr[nBinsy+1];
1751 for(Int_t i=0; i<=nBinsx; i++) xArr[i]=xMin+i*(xMax-xMin)/nBinsx;
1752 for(Int_t i=0; i<=nBinsy; i++) yArr[i]=yMin+i*(yMax-yMin)/nBinsy;
1754 TH3F *res = new TH3F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xArr, nBinsy, yArr, nBinsz, zArr);
1755 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1756 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1757 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1758 // res->SetMarkerStyle(kFullCircle);
1759 // res->SetOption("E");
1760 res->SetLineColor(color);
1761 // fOutputList->Add(res);
1766 //////////////////////////////////////
1770 //________________________________________________________________________
1771 void AnalysisJetMain::ConnectInputData(Option_t *)
1776 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
1778 // AliError("ESD not available");
1779 // fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
1780 // fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1782 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1783 // assume that the AOD is in the general output...
1784 fAODOut = AODEvent();
1790 void AliAnalysisTaskJetShape::SetBranchNames(const TString &branch1, const TString &branch2)
1792 fJetBranchName[0] = branch1;
1793 fJetBranchName[1] = branch2;
1796 //________________________________________________________________________
1797 void AliAnalysisTaskJetShape::UserCreateOutputObjects()
1799 //printf("Open1\n");
1800 // const char *nameF = OpenFile(1)->GetName();
1801 //printf("Open2 %s\n",nameF);
1803 // fTriggerAnalysis = new AliTriggerAnalysis();
1804 // if (fkMC) fTriggerAnalysis->SetAnalyzeMC(1);
1807 fOutputList = new TList();
1808 fOutputList->SetOwner();
1810 fhPtJL = Hist1D("hPtJL", 100, 0, 200, "Pt_{JL}"); fOutputList->Add(fhPtJL);
1811 fhAreaJL = Hist1D("hAreaJL", 100, 0., 4, "AreaJL"); fOutputList->Add(fhAreaJL);
1815 printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() \n");
1818 fanalJetSubStrHM->SetFilterMask(fFilterMask);
1819 if(fkMC) fanalJetSubStrHM->MCprod();
1820 fanalJetSubStrHM->InitHistos();
1821 fanalJetSubStrHM->AddToList(fOutputList);
1827 printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() MC\n");
1828 fanalJetSubStrMCHM->InitHistos();
1829 fanalJetSubStrMCHM->AddToList(fOutputList);
1831 for(Int_t i=0; i<3; i++)
1833 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]);
1834 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]);
1835 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]);
1836 fhDCAxy[i] = Hist2D(Form("hDCAxy%d",i), 100, 0, 100, 300, -3, 3, "DCAxy of prim [cm]"); fOutputList->Add(fhDCAxy[i]);
1837 fhDCAz[i] = Hist2D(Form("hDCAz%d",i) , 100, 0, 100, 300, -3, 3, "DCAz of prim [cm]") ; fOutputList->Add(fhDCAz[i]);
1838 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]);
1845 printf(" JetBranchName[0]=%s\n JetBranchName[1]=%s\n", fJetBranchName[0].Data(), fJetBranchName[1].Data());
1849 PostData(1, fOutputList);
1854 fOutputTree = new TTree("tree","Tree with Ali2JEvent");
1855 fOutputTree->Branch("event",&f2JEvent);
1858 // PostData(1, fOutputTree);
1862 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
1865 // Fetch the aod also from the input in,
1866 // have todo it in notify
1869 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1870 // assume that the AOD is in the general output...
1871 fAODOut = AODEvent();
1873 if(fNonStdFile.Length()!=0){
1874 // case that we have an AOD extension we need can fetch the jets from the extended output
1875 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1876 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
1878 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
1891 //________________________________________________________________________
1892 void AliAnalysisTaskJetShape::UserExec(Option_t *)
1897 // f2JEvent = new Ali2JEvent();
1901 if(fDebug) Printf("\n\n\nEvent #%5d", (Int_t) fEntry);
1902 if(fDebug) printf("NEW EVENT 0\n");
1904 if(!IsGoodEvent()) return;
1907 // printf("\n\n\n NEW EVENT\n");
1909 AliAODEvent* aodE = 0;
1911 if(!fESD)aodE = fAODIn;
1912 else aodE = fAODOut;}
1915 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1917 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
1922 printf("NEW EVENT 1: Number of tracks = %d\n",aodE->GetNumberOfTracks());
1923 // aodE->GetList()->Print();
1924 // printf("Look for %s\n",fJetBranchName[0].Data());
1928 // centrality selection
1929 AliCentrality *cent = 0x0;
1930 Double_t centrality = 0.;
1932 if(fESD) {cent = fESD->GetCentrality();
1933 if(cent) centrality = cent->GetCentralityPercentile("V0M");}
1934 else centrality=((AliVAODHeader*)aodE->GetHeader())->GetCentrality();
1941 // if(fDebug) printf("NEW EVENT 2\n");
1943 Bool_t fUseAOD049 = kFALSE;// kTRUE;//
1944 if(fUseAOD049&¢rality>=0){
1946 AliAODVZERO* aodV0 = aodE->GetVZEROData();
1947 v0+=aodV0->GetMTotV0A();
1948 v0+=aodV0->GetMTotV0C();
1949 if(centrality==0 && v0 < 19500) return ;//filtering issue
1950 Float_t tkl = (Float_t)(aodE->GetTracklets()->GetNumberOfTracklets());
1951 Float_t val= 1.30552 + 0.147931 * v0;
1952 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,
1953 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,
1954 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,
1955 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,
1956 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,
1957 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,
1958 13.7544, 13.7544, 13.7544, 13.7544, 13.7544};
1959 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)centrality] )
1966 printf("centrality: %f\n", centrality);
1968 aodE->GetList()->Print();
1972 if (centrality < fCentMin || centrality > fCentMax){
1973 // PostData(1, fOutputList);
1977 // fhNEvent->Fill(0);
1980 // -- end event selection --
1984 AliAODJetEventBackground* externalBackground = 0;
1985 AliAODJetEventBackground* externalBackgroundMC = 0;
1993 if(fAODOut&&!externalBackground&&fBackgroundBranch[0].Length()){
1994 externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[0].Data()));
1995 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1997 if(fAODExtension&&!externalBackground&&fBackgroundBranch[0].Length()){
1998 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[0].Data()));
1999 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
2001 if(fAODIn&&!externalBackground&&fBackgroundBranch[0].Length()){
2002 externalBackground = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[0].Data()));
2003 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
2005 if(fAODOut&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2006 externalBackgroundMC = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[1].Data()));
2007 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2009 if(fAODExtension&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2010 externalBackgroundMC = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[1].Data()));
2011 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2013 if(fAODIn&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
2014 externalBackgroundMC = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[1].Data()));
2015 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
2017 if(externalBackground)rho = externalBackground->GetBackground(0);
2018 if(externalBackgroundMC)rho_MC = externalBackgroundMC->GetBackground(0);
2025 AliAODVertex* primVtx = aodE->GetPrimaryVertex();
2026 Double_t bfield = aodE->GetMagneticField();
2028 TClonesArray *arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
2031 printf("ERROR: no Info about particles in AOD\n");
2035 for(int it = 0;it < aodE->GetNumberOfTracks(); it++)
2037 AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodE->GetTrack(it));
2038 if(!tr) AliFatal("Not a standard AOD");
2040 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))) continue;
2041 if(TMath::Abs(tr->Eta())>1.) continue;
2042 Int_t label = TMath::Abs(tr->GetLabel());
2044 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(label));
2046 // if(!part->IsPhysicalPrimary())continue;
2047 // if(part->Charge()==0)continue;
2048 // vec.SetXYZ(part->Px(), part->Py(), part->Pz());
2055 AliAODTrack tmp(*tr);
2056 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
2059 if(!part->IsPhysicalPrimary()) type=1;
2060 if(label<0) type =2;
2062 fhPtresVsPt[type]->Fill(part->Pt(), 1-tr->Pt()/part->Pt());
2063 fhPhiresVsPhi[type]->Fill(part->Phi(), tr->Phi() - part->Phi());
2064 fhEtaresVsEta[type]->Fill(part->Eta(), tr->Eta() - part->Eta());
2065 fhDCAxy[type]->Fill(part->Pt(), dca[0]);
2066 fhDCAz[type]->Fill( part->Pt(), dca[1]);
2067 fhTrackPtEtaPhi[type]->Fill(tr->Pt(), tr->Eta(), tr->Phi());
2077 // printf("rho = %f\n",rho);
2081 TClonesArray *Jets[2];
2084 if(fAODOut&&!Jets[0]){
2085 Jets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
2086 Jets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); }
2087 if(fAODExtension && !Jets[0]){
2088 Jets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
2089 Jets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
2090 if(fAODIn&&!Jets[0]){
2091 Jets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
2092 Jets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
2096 Printf("no friend!!!\n");
2099 Int_t nJ = Jets[0]->GetEntries();
2104 printf("NEW EVENT 3: Number of Rec. jets %d\n",nJ);
2108 // Find highest pT jet with pt > 20 GeV
2112 Float_t etaJmax = 0.4;
2114 for (Int_t i = 0; i < nJ; i++) {
2115 AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[0]->At(i));
2118 Float_t ptJ = jet->Pt();
2119 Float_t etaJ = TMath::Abs(jet->Eta());
2121 Float_t area = jet->EffectiveAreaCharged();
2122 Float_t ptJcorr=ptJ-rho*area;
2124 if ((ptJcorr > fPtJcorrMin) && (ptJcorr > ptmax) && etaJ < etaJmax) {
2133 printf("NEW EVENT 4:\n");
2138 AliAODJet* jetL = 0;
2140 Float_t areaJL, ptJLcorr;
2144 jetL = dynamic_cast<AliAODJet*> (Jets[0]->At(imax));
2147 areaJL = jetL->EffectiveAreaCharged();
2148 ptJLcorr = jetL->Pt()-rho*areaJL;
2150 fhPtJL->Fill(ptJLcorr);
2151 fhAreaJL->Fill(areaJL);
2152 vecJ.SetXYZ(jetL->Px(), jetL->Py(), jetL->Pz());
2153 fanalJetSubStrHM->AddEvent(aodE, jetL, rho*areaJL);
2156 Printf("Leading Jet");
2166 PostData(1, fOutputList);
2167 if(fDebug) Printf("End of event #%5d", (Int_t) fEntry);
2174 nJ = Jets[1]->GetEntries();
2176 printf("NEW EVENT 5: Number of Rec. jets %d\n",nJ);
2181 Float_t areaJL_MC=0;
2183 for (Int_t i = 0; i < nJ; i++) {
2184 AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[1]->At(i));
2186 Float_t ptJ1 = jet->Pt();
2187 Float_t etaJ1 = TMath::Abs(jet->Eta());
2189 areaJL_MC = jet->EffectiveAreaCharged();
2190 Float_t ptJcorr=ptJ1-rho*areaJL_MC;
2194 if ((ptJcorr > fPtJcorrMin) && (ptJcorr > ptmax) && etaJ1 < etaJmax) {
2202 printf("NEW EVENT 6:\n");
2206 AliAODJet* jetMCL=0;
2209 jetMCL = dynamic_cast<AliAODJet*> (Jets[1]->At(imax));
2211 fanalJetSubStrMCHM->AddEvent(aodE, jetMCL, rho_MC*areaJL_MC);
2219 PostData(1, fOutputList);
2221 if(fDebug) Printf("End of event #%5d", (Int_t) fEntry);
2242 Double_t AliAnalysisTaskJetShape::CalcdPhi(Double_t phi1, Double_t phi2)
2244 while(phi1<0) phi1+=TMath::TwoPi();
2245 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
2247 while(phi2<0) phi2+=TMath::TwoPi();
2248 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
2250 Double_t dphi = phi1- phi2;
2251 if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
2252 if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
2254 // Double_t dphi = phi1- phi2;
2255 // while(dphi<0) dphi+=TMath::TwoPi();
2256 // while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2260 Double_t AliAnalysisTaskJetShape::CalcdPhi0To2pi(Double_t phi1, Double_t phi2)
2263 Double_t dphi = CalcdPhi(phi1, phi2);
2264 while(dphi<0) dphi+=TMath::TwoPi();
2265 while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2270 Bool_t AliAnalysisTaskJetShape::IsGoodEvent()
2273 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
2275 // AliError("ESD not available");
2276 // fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
2277 // fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
2279 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
2280 // assume that the AOD is in the general output...
2281 fAODOut = AODEvent();
2284 static AliAODEvent* aod = 0;
2285 // take all other information from the aod we take the tracks from
2287 if(!fESD)aod = fAODIn;
2288 else aod = fAODOut;}
2292 if(fNonStdFile.Length()!=0){
2293 // case that we have an AOD extension we need can fetch the jets from the extended output
2294 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2295 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
2299 if(fAODIn) printf("fAODIn\n");
2300 if(fAODOut) printf("fAODOut\n");
2301 if(fAODExtension) printf("fAODExtension\n");
2306 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2308 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
2313 // -- event selection --
2315 // physics selection
2316 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
2317 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2318 // cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
2319 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
2320 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
2326 // if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
2327 // fhNEvent->Fill(3);
2328 // PostData(1, fOutputList);
2332 AliAODVertex* primVtx = aod->GetPrimaryVertex();
2335 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
2339 Int_t nTracksPrim = primVtx->GetNContributors();
2340 if ((nTracksPrim < fVtxMinContrib) ||
2341 (primVtx->GetZ() < fVtxZMin) ||
2342 (primVtx->GetZ() > fVtxZMax) ){
2343 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
2348 // event class selection (from jet helper task)
2349 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
2350 if(fDebug) Printf("Event class %d", eventClass);
2351 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
2363 //________________________________________________________________________
2364 void AliAnalysisTaskJetShape::Terminate(Option_t *)
2366 // Draw result to the screen
2367 // Called once at the end of the query
2368 // fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
2370 printf("MyTaskTestTrigger::Terminate()\n\n\n");
2373 // fOutputList = dynamic_cast<TList*> (GetOutputData(1));
2374 // fOutputTree = dynamic_cast<TTree*> (GetOutputData(1));
2381 Bool_t AnalysisJetMain::GetVertex(const AliESDEvent* esd, Double_t Vxyz[3], Int_t type)
2386 const AliESDVertex* vtx = 0;
2389 vtx = esd->GetPrimaryVertexSPD();
2390 if(!vtx) return kFALSE;
2391 if(vtx->GetNContributors() < 1) return kFALSE;
2392 TString vtxTyp = vtx->GetTitle();
2393 if (vtxTyp.Contains("vertexer: Z")) {
2394 if (vtx->GetDispersion()>0.04) return kFALSE;
2395 if (vtx->GetZRes()>0.25) return kFALSE;
2400 vtx = esd->GetPrimaryVertexTracks();
2401 if(!vtx) return kFALSE;
2402 if(vtx->GetNContributors() < 1) return kFALSE;
2406 Vxyz[0] = vtx->GetX();
2407 Vxyz[1] = vtx->GetY();
2408 Vxyz[2] = vtx->GetZ();
2422 TH1F* AliAnalysisTaskJetShape::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax, const char* xLabel, Int_t color)
2424 // create a 1D histogram
2426 TH1F *h = new TH1F(name, name, nBins, xMin, xMax);
2427 if (xLabel) h->GetXaxis()->SetTitle(xLabel);
2428 // res->SetMarkerStyle(kFullCircle);
2429 // res->SetOption("E");
2430 h->SetLineColor(color);
2431 // fOutputList->Add(h);
2436 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)
2438 // create a 2D histogram
2440 TH2F *res = new TH2F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax);
2441 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2442 if (xLabel) res->GetYaxis()->SetTitle(yLabel);
2443 // res->SetMarkerStyle(kFullCircle);
2444 // res->SetOption("E");
2445 res->SetLineColor(color);
2446 // fOutputList->Add(res);
2451 TH3F *AliAnalysisTaskJetShape::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
2452 Int_t nBinsy, Double_t yMin, Double_t yMax,
2453 Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
2455 // create a 3D histogram
2457 TH3F *res = new TH3F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax, nBinsz, zMin, zMax);
2458 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2459 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
2460 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
2461 // res->SetMarkerStyle(kFullCircle);
2462 // res->SetOption("E");
2463 res->SetLineColor(color);
2464 // fOutputList->Add(res);