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 fhPsiVsRPtJ = Hist2D("hPsiVsRPtJ", fPsiVsRNbin, 0., fRmax, fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{tJ} GeV/c", 1, "#Psi(R)");
1254 fhPtJvsPtCorr = Hist2D("fhPtJvsPtCorr", fPtJetNbin, fPtJetArray.GetArray(), 100, -100, 100, "P_{tJ} GeV/c", "P_{tJ} - P_{tJ,corr} GeV/c" , 1);
1259 for(Int_t i=0; i<3; i++)
1261 fhTMA_JAA[i] = Hist1D(Form("fhTMA_JAA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1262 fhTMA_JAp[i] = Hist1D(Form("fhTMA_JAp_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1264 fhTMA_B1AA[i] = Hist1D(Form("fhTMA_B1AA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1265 fhTMA_B1Ap[i] = Hist1D(Form("fhTMA_B1Ap_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1267 fhTMA_B2AA[i] = Hist1D(Form("fhTMA_B2AA_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1268 fhTMA_B2Ap[i] = Hist1D(Form("fhTMA_B2Ap_%d",i), fPsiVsPhiNbin, 0, TMath::TwoPi(), "#phi_{R}", 1, "#Psi(#phi_{R})");
1273 fhPsiVsR_MCtr = Hist1D("hPsiVsR_MCtr", fPsiVsRNbin, 0., fRmax, "R", 1, "#Psi(R)");
1274 fhPsiVsRPtJ_MCtr = Hist2D("hPsiVsRPtJ_MCtr", fPsiVsRNbin, 0., fRmax, fPtJetNbin, fPtJetArray.GetArray(), "R", "P_{tJ} GeV/c", 1, "#Psi(R)");
1275 fhJetTrPtVsPartPt = Hist2D("fhJetTrPtVsPartPt", fPtJetNbin, fPtJetArray.GetArray(), 100, -1, 1, "P_{tJ,part} GeV/c", "1-P_{tJ,tr}/P_{tJ,part} GeV/c" , 1);
1276 const char *ch[2]={"p","m"};
1277 for(Int_t i=0; i<3; i++)
1279 for(Int_t j=0; j<2; j++)
1281 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" );
1282 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" );
1283 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}" );
1284 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}" );
1285 fhDCAxy[i][j] = Hist2D(Form("hDCAxy%d%s",i,ch[j]), 100, 0, 100, 300, -3, 3, "DCAxy of prim [cm]");
1286 fhDCAz[i][j] = Hist2D(Form("hDCAz%d%s",i,ch[j]) , 100, 0, 100, 300, -3, 3, "DCAz of prim [cm]") ;
1287 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");
1298 void AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddToList(TList *list)
1302 list->Add(fhPhiEtaTrack);
1304 list->Add(fhPsiVsR);
1305 list->Add(fhPsiVsRPtJ);
1306 list->Add(fhPtJvsPtCorr);
1309 for(Int_t i=0; i<3; i++)
1311 list->Add(fhTMA_JAA[i]);
1312 list->Add(fhTMA_JAp[i]);
1314 list->Add(fhTMA_B1AA[i]);
1315 list->Add(fhTMA_B1Ap[i]);
1317 list->Add(fhTMA_B2AA[i]);
1318 list->Add(fhTMA_B2Ap[i]);
1323 list->Add(fhPsiVsR_MCtr);
1324 list->Add(fhPsiVsRPtJ_MCtr);
1325 list->Add(fhJetTrPtVsPartPt);
1327 for(Int_t i=0; i<3; i++)
1329 for(Int_t j=0; j<2; j++)
1331 list->Add(fhPtresVsPt[i][j]);
1332 list->Add(fhPhiresVsPhi[i][j]);
1333 list->Add(fhEtaresVsEta[i][j]);
1334 list->Add(fhRresVsPt[i][j]);
1335 list->Add(fhDCAxy[i][j]);
1336 list->Add(fhDCAz[i][j]);
1337 list->Add(fhTrackPtEtaPhi[i][j]);
1346 Bool_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::AddEvent(AliAODEvent* aodE, AliAODJet *jet, Double_t DeltaPtJ)
1351 Int_t size = aodE->GetNumberOfTracks();
1353 TClonesArray *arrP = 0;
1355 if(fkMC || fkMCprod)
1357 arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
1360 printf("ERROR: no Info about particles in AOD\n");
1366 size = arrP->GetEntriesFast();
1368 Int_t IndexArray[size];
1369 Int_t IndexArrayMC[size];
1371 TClonesArray farray("TVector3", size);
1374 Int_t counterMC = 0;
1375 Int_t counterAll = 0;
1376 Double_t pJ[3] = {0, 0, 0};
1377 Double_t pJmc[3] = {0, 0, 0};
1379 AliAODVertex* primVtx = aodE->GetPrimaryVertex();
1380 Double_t bfield = aodE->GetMagneticField();
1384 TVector3 vecJ(jet->Px(),jet->Py(),jet->Pz());
1386 for(int it = 0;it < size;++it){
1395 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(it));
1397 if(!part->IsPhysicalPrimary())continue;
1398 if(part->Charge()==0)continue;
1399 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1403 AliAODTrack *tr = aodE->GetTrack(it);
1405 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1406 label = tr->GetLabel();
1407 AliAODTrack tmp(*tr);
1408 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1409 vec.SetXYZ(tr->Px(), tr->Py(), tr->Pz());
1414 if(TMath::Abs(vec.Eta())>fEtaTrackMax)continue;
1415 if(vec.Pt()< fPtTrackMin)continue;
1416 if(vec.Pt()> fPtTrackMax) {return kFALSE;}
1419 new(farray[counterAll]) TVector3(vec);
1423 Double_t R = CalcR(vecJ, vec);
1424 if(R> fRmax) continue;
1430 IndexArray[counter] = it;
1439 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(TMath::Abs(label)));
1443 if(!part->IsPhysicalPrimary()) type=1;
1444 if(label <0) type=2;
1446 if(!(ch==-1 || ch==1)) AliFatal("ch != +/- 1!!!\n\n");
1450 Double_t phip = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(part->Phi());
1451 Double_t dphi = AliAnalysisTaskJetShapeTool::CalcdPhiSigned(part->Phi(), vec.Phi());
1453 fhPtresVsPt[type][ch]->Fill(part->Pt(), 1-vec.Pt()/part->Pt());
1454 fhPhiresVsPhi[type][ch]->Fill(phip, dphi);
1455 fhEtaresVsEta[type][ch]->Fill(part->Eta(), vec.Eta() - part->Eta());
1456 fhDCAxy[type][ch]->Fill(part->Pt(), dca[0]);
1457 fhDCAz[type][ch]->Fill( part->Pt(), dca[1]);
1458 fhTrackPtEtaPhi[type][ch]->Fill(vec.Pt(), vec.Eta(), vec.Phi());
1460 TVector3 vecP(part->Px(), part->Py(), part->Pz());
1461 Double_t Rgen = CalcR(vecJ, vecP);
1462 fhRresVsPt[type][ch]->Fill(part->Pt(), Rgen - R);
1464 pJmc[0]+=part->Px();
1465 pJmc[1]+=part->Py();
1466 pJmc[2]+=part->Pz();
1468 IndexArrayMC[counterMC] = label;
1475 fhMult->Fill(counter);
1476 if(counter<1) return kFALSE;
1478 fJvec.SetXYZ(pJ[0], pJ[1], pJ[2]);
1481 fPtJ = TMath::Sqrt(pJ[0]*pJ[0] + pJ[1]*pJ[1]);
1482 if((fPtJ < fPtJmin) || (fPtJ > fPtJmax)) return kFALSE;
1489 fhPtJvsPtCorr->Fill(fPtJ, jet->Pt() - DeltaPtJ);
1491 for(Int_t it = 0; it<counter; it++)
1497 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(IndexArray[it]));
1499 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1503 AliAODTrack *tr = aodE->GetTrack(IndexArray[it]);
1505 AliAODTrack tmp(*tr);
1506 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
1507 if(tr)vec.SetXYZ(tmp.Px(), tmp.Py(), tmp.Pz());
1510 Double_t R = CalcR(fJvec, vec);
1511 // Double_t pt = (tr->Px()*pJ[0] + tr->Py()*pJ[1])/PtJ;
1512 // Double_t probL = fJvec.Dot(vec)/fJvec.Mag2();
1513 // fhPsiVsR->Fill(R, probL);
1514 // fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1516 fhPsiVsRPtJ->Fill(R, fPtJ);
1518 Double_t phi = AliAnalysisTaskJetShapeTool::CalcdPhi0To2pi(vec.Phi());
1519 fhPhiEtaTrack->Fill(phi, vec.Eta());
1526 TVector3 fJvecMCtr(pJmc[0], pJmc[1], pJmc[2]);
1527 Double_t fPtJMCtr= fJvecMCtr.Perp();
1529 fhJetTrPtVsPartPt->Fill(fPtJMCtr, 1-fPtJ/fPtJMCtr);
1530 for(Int_t it = 0; it<counterMC; it++)
1534 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(IndexArrayMC[it]));
1536 vec.SetXYZ(part->Px(), part->Py(), part->Pz());
1538 Double_t R = CalcR(fJvecMCtr, vec);
1539 // Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1540 // fhPsiVsR->Fill(R, probL);
1541 // fhPsiVsRPtJ->Fill(R, fPtJ, probL);
1542 // Double_t probL = fJvecMCtr.Dot(vec)/fJvecMCtr.Mag2();
1543 fhPsiVsR_MCtr->Fill(R);
1544 fhPsiVsRPtJ_MCtr->Fill(R, fPtJMCtr);
1549 // fMyTool->Clean();
1550 AliAnalysisTaskJetShapeTool *fMyTool = new AliAnalysisTaskJetShapeTool(&farray);
1552 fMyTool->SetNEntries(counterAll);
1553 fMyTool->SetVecJ(vecJ);
1554 fMyTool->SetPtMinTr(fPtTrackMin,fPtTrackMin);
1559 to be added!!!!!!!!!
1564 for(Int_t l=0; l<3; l++)
1566 Double_t phiA, phi1;
1568 // Double_t ptRJ0= fMyTool->GetPRecInRJ().Pt();
1569 // Double_t ptRJ = fMyTool->GetPRecJ(l,0).Pt();
1570 Int_t N1 = fMyTool->GetSizeJ(l,1);
1571 Double_t dphi = -999;
1574 if(!fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 0) , vecJ, phiA, scenario))
1577 // f2JEvent->SetPhiJL(l,0, phiA);
1579 if(fMyTool->FindCorrelationAxes(fMyTool->GetListJ(l, 1) , vecJ, phi1, scenario))
1581 fhTMA_JAA[l]->Fill(fMyTool->CalcdPhi0To2pi(phiA-phi1));
1585 for(Int_t ip =0; ip<N1; ip++)
1587 phi1 = fMyTool->GetLocPhiJ(l,1,ip);
1588 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1589 fhTMA_JAp[l]->Fill(dphi);
1593 Int_t NB1 = fMyTool->GetSizeB1(l, 1);
1594 for(Int_t ip =0; ip<NB1; ip++)
1596 phi1 = fMyTool->GetLocPhiB1(l,1,ip);
1597 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1598 fhTMA_B1Ap[l]->Fill(dphi);
1601 Int_t NB2 = fMyTool->GetSizeB2(l, 1);
1602 for(Int_t ip =0; ip<NB2; ip++)
1604 phi1 = fMyTool->GetLocPhiB2(l,1,ip);
1605 dphi = fMyTool->CalcdPhi0To2pi(phiA, phi1);
1606 fhTMA_B2Ap[l]->Fill(dphi);
1619 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcR(TVector3 v1, TVector3 v2)
1622 Double_t dphi = CalcdPhi(v1.Phi(), v2.Phi());
1623 Double_t deta = v1.Eta() - v2.Eta();
1624 Double_t RB = TMath::Sqrt(dphi*dphi+deta*deta);
1629 Double_t AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::CalcdPhi(Double_t phi1, Double_t phi2)
1632 while(phi1<0) phi1+=TMath::TwoPi();
1633 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
1635 while(phi2<0) phi2+=TMath::TwoPi();
1636 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
1638 Double_t dphi = phi1- phi2;
1639 if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
1640 if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
1651 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)
1653 // create a 1D histogram
1655 TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xMin, xMax);
1656 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1657 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1658 res->SetLineColor(color);
1663 TH1F* AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist1D(const char* name, Int_t nBins, Double_t *xArray, const char* xLabel, Int_t color, const char* yLabel)
1665 // create a 1D histogram
1667 TH1F* res = new TH1F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBins, xArray);
1668 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1669 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1670 res->SetLineColor(color);
1674 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)
1676 // create a 2D histogram
1678 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yMin, yMax);
1679 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1680 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1681 // res->SetMarkerStyle(kFullCircle);
1682 // res->SetOption("E");
1683 res->SetLineColor(color);
1684 // fOutputList->Add(res);
1689 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)
1691 // create a 2D histogram
1693 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, xMin, xMax, nBinsy, yArray);
1694 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1695 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1696 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1697 // res->SetMarkerStyle(kFullCircle);
1698 // res->SetOption("E");
1699 res->SetLineColor(color);
1700 // fOutputList->Add(res);
1704 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)
1706 // create a 2D histogram
1708 TH2F *res = new TH2F(Form("%s_%s",fComment.Data(), name), Form("%s_%s",fComment.Data(), name), nBinsx, yArrax, nBinsy, yMin, yMax);
1709 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1710 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1711 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1712 // res->SetMarkerStyle(kFullCircle);
1713 // res->SetOption("E");
1714 res->SetLineColor(color);
1715 // fOutputList->Add(res);
1719 TH3F *AliAnalysisTaskJetShape::AliAnalysisTaskJetShapeHM::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
1720 Int_t nBinsy, Double_t yMin, Double_t yMax,
1721 Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
1723 // create a 3D histogram
1725 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);
1726 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
1727 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
1728 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
1729 // res->SetMarkerStyle(kFullCircle);
1730 // res->SetOption("E");
1731 res->SetLineColor(color);
1732 // fOutputList->Add(res);
1737 //////////////////////////////////////
1741 //________________________________________________________________________
1742 void AnalysisJetMain::ConnectInputData(Option_t *)
1747 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
1749 // AliError("ESD not available");
1750 // fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
1751 // fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1753 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1754 // assume that the AOD is in the general output...
1755 fAODOut = AODEvent();
1761 void AliAnalysisTaskJetShape::SetBranchNames(const TString &branch1, const TString &branch2)
1763 fJetBranchName[0] = branch1;
1764 fJetBranchName[1] = branch2;
1767 //________________________________________________________________________
1768 void AliAnalysisTaskJetShape::UserCreateOutputObjects()
1770 //printf("Open1\n");
1771 // const char *nameF = OpenFile(1)->GetName();
1772 //printf("Open2 %s\n",nameF);
1774 // fTriggerAnalysis = new AliTriggerAnalysis();
1775 // if (fkMC) fTriggerAnalysis->SetAnalyzeMC(1);
1778 fOutputList = new TList();
1779 fOutputList->SetOwner();
1781 fhPtJL = Hist1D("hPtJL", 100, 0, 200, "Pt_{JL}"); fOutputList->Add(fhPtJL);
1782 fhAreaJL = Hist1D("hAreaJL", 100, 0., 4, "AreaJL"); fOutputList->Add(fhAreaJL);
1786 printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() \n");
1789 fanalJetSubStrHM->SetFilterMask(fFilterMask);
1790 if(fkMC) fanalJetSubStrHM->MCprod();
1791 fanalJetSubStrHM->InitHistos();
1792 fanalJetSubStrHM->AddToList(fOutputList);
1798 printf("AliAnalysisTaskJetShape::UserCreateOutputObjects() MC\n");
1799 fanalJetSubStrMCHM->InitHistos();
1800 fanalJetSubStrMCHM->AddToList(fOutputList);
1802 for(Int_t i=0; i<3; i++)
1804 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]);
1805 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]);
1806 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]);
1807 fhDCAxy[i] = Hist2D(Form("hDCAxy%d",i), 100, 0, 100, 300, -3, 3, "DCAxy of prim [cm]"); fOutputList->Add(fhDCAxy[i]);
1808 fhDCAz[i] = Hist2D(Form("hDCAz%d",i) , 100, 0, 100, 300, -3, 3, "DCAz of prim [cm]") ; fOutputList->Add(fhDCAz[i]);
1809 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]);
1816 printf(" JetBranchName[0]=%s\n JetBranchName[1]=%s\n", fJetBranchName[0].Data(), fJetBranchName[1].Data());
1820 PostData(1, fOutputList);
1825 fOutputTree = new TTree("tree","Tree with Ali2JEvent");
1826 fOutputTree->Branch("event",&f2JEvent);
1829 // PostData(1, fOutputTree);
1833 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
1836 // Fetch the aod also from the input in,
1837 // have todo it in notify
1840 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1841 // assume that the AOD is in the general output...
1842 fAODOut = AODEvent();
1844 if(fNonStdFile.Length()!=0){
1845 // case that we have an AOD extension we need can fetch the jets from the extended output
1846 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1847 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
1849 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
1862 //________________________________________________________________________
1863 void AliAnalysisTaskJetShape::UserExec(Option_t *)
1868 // f2JEvent = new Ali2JEvent();
1872 if(fDebug) Printf("\n\n\nEvent #%5d", (Int_t) fEntry);
1873 if(fDebug) printf("NEW EVENT 0\n");
1875 if(!IsGoodEvent()) return;
1878 // printf("\n\n\n NEW EVENT\n");
1880 AliAODEvent* aodE = 0;
1882 if(!fESD)aodE = fAODIn;
1883 else aodE = fAODOut;}
1886 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1888 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
1893 printf("NEW EVENT 1: Number of tracks = %d\n",aodE->GetNumberOfTracks());
1894 // aodE->GetList()->Print();
1895 // printf("Look for %s\n",fJetBranchName[0].Data());
1899 // centrality selection
1900 AliCentrality *cent = 0x0;
1901 Double_t centrality = 0.;
1903 if(fESD) {cent = fESD->GetCentrality();
1904 if(cent) centrality = cent->GetCentralityPercentile("V0M");}
1905 else centrality=aodE->GetHeader()->GetCentrality();
1912 // if(fDebug) printf("NEW EVENT 2\n");
1914 Bool_t fUseAOD049 = kFALSE;// kTRUE;//
1915 if(fUseAOD049&¢rality>=0){
1917 AliAODVZERO* aodV0 = aodE->GetVZEROData();
1918 v0+=aodV0->GetMTotV0A();
1919 v0+=aodV0->GetMTotV0C();
1920 if(centrality==0 && v0 < 19500) return ;//filtering issue
1921 Float_t tkl = (Float_t)(aodE->GetTracklets()->GetNumberOfTracklets());
1922 Float_t val= 1.30552 + 0.147931 * v0;
1923 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,
1924 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,
1925 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,
1926 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,
1927 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,
1928 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,
1929 13.7544, 13.7544, 13.7544, 13.7544, 13.7544};
1930 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)centrality] )
1937 printf("centrality: %f\n", centrality);
1939 aodE->GetList()->Print();
1943 if (centrality < fCentMin || centrality > fCentMax){
1944 // PostData(1, fOutputList);
1948 // fhNEvent->Fill(0);
1951 // -- end event selection --
1955 AliAODJetEventBackground* externalBackground = 0;
1956 AliAODJetEventBackground* externalBackgroundMC = 0;
1964 if(fAODOut&&!externalBackground&&fBackgroundBranch[0].Length()){
1965 externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[0].Data()));
1966 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1968 if(fAODExtension&&!externalBackground&&fBackgroundBranch[0].Length()){
1969 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[0].Data()));
1970 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1972 if(fAODIn&&!externalBackground&&fBackgroundBranch[0].Length()){
1973 externalBackground = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[0].Data()));
1974 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[0].Data());
1976 if(fAODOut&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
1977 externalBackgroundMC = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch[1].Data()));
1978 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
1980 if(fAODExtension&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
1981 externalBackgroundMC = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch[1].Data()));
1982 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
1984 if(fAODIn&&!externalBackgroundMC&&fBackgroundBranch[1].Length()){
1985 externalBackgroundMC = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch[1].Data()));
1986 if(!externalBackgroundMC)Printf("%s:%d MC Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch[1].Data());
1988 if(externalBackground)rho = externalBackground->GetBackground(0);
1989 if(externalBackgroundMC)rho_MC = externalBackgroundMC->GetBackground(0);
1996 AliAODVertex* primVtx = aodE->GetPrimaryVertex();
1997 Double_t bfield = aodE->GetMagneticField();
1999 TClonesArray *arrP = dynamic_cast<TClonesArray*>(aodE->FindListObject(AliAODMCParticle::StdBranchName()));
2002 printf("ERROR: no Info about particles in AOD\n");
2006 for(int it = 0;it < aodE->GetNumberOfTracks(); it++)
2008 AliAODTrack *tr = aodE->GetTrack(it);
2010 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))) continue;
2011 if(TMath::Abs(tr->Eta())>1.) continue;
2012 Int_t label = TMath::Abs(tr->GetLabel());
2014 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(arrP->At(label));
2016 // if(!part->IsPhysicalPrimary())continue;
2017 // if(part->Charge()==0)continue;
2018 // vec.SetXYZ(part->Px(), part->Py(), part->Pz());
2025 AliAODTrack tmp(*tr);
2026 tmp.PropagateToDCA(primVtx, bfield, 5., dca, cov);
2029 if(!part->IsPhysicalPrimary()) type=1;
2030 if(label<0) type =2;
2032 fhPtresVsPt[type]->Fill(part->Pt(), 1-tr->Pt()/part->Pt());
2033 fhPhiresVsPhi[type]->Fill(part->Phi(), tr->Phi() - part->Phi());
2034 fhEtaresVsEta[type]->Fill(part->Eta(), tr->Eta() - part->Eta());
2035 fhDCAxy[type]->Fill(part->Pt(), dca[0]);
2036 fhDCAz[type]->Fill( part->Pt(), dca[1]);
2037 fhTrackPtEtaPhi[type]->Fill(tr->Pt(), tr->Eta(), tr->Phi());
2047 // printf("rho = %f\n",rho);
2051 TClonesArray *Jets[2];
2054 if(fAODOut&&!Jets[0]){
2055 Jets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
2056 Jets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); }
2057 if(fAODExtension && !Jets[0]){
2058 Jets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
2059 Jets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
2060 if(fAODIn&&!Jets[0]){
2061 Jets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
2062 Jets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
2066 Printf("no friend!!!\n");
2069 Int_t nJ = Jets[0]->GetEntries();
2074 printf("NEW EVENT 3: Number of Rec. jets %d\n",nJ);
2078 // Find highest pT jet with pt > 20 GeV
2082 Float_t etaJmax = 0.4;
2084 for (Int_t i = 0; i < nJ; i++) {
2085 AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[0]->At(i));
2088 Float_t ptJ = jet->Pt();
2089 Float_t etaJ = TMath::Abs(jet->Eta());
2091 Float_t area = jet->EffectiveAreaCharged();
2092 Float_t ptJcorr=ptJ-rho*area;
2094 if ((ptJcorr > fPtJcorrMin) && (ptJcorr > ptmax) && etaJ < etaJmax) {
2103 printf("NEW EVENT 4:\n");
2108 AliAODJet* jetL = 0;
2110 Float_t areaJL, ptJLcorr;
2114 jetL = dynamic_cast<AliAODJet*> (Jets[0]->At(imax));
2117 areaJL = jetL->EffectiveAreaCharged();
2118 ptJLcorr = jetL->Pt()-rho*areaJL;
2120 fhPtJL->Fill(ptJLcorr);
2121 fhAreaJL->Fill(areaJL);
2122 vecJ.SetXYZ(jetL->Px(), jetL->Py(), jetL->Pz());
2123 fanalJetSubStrHM->AddEvent(aodE, jetL, rho*areaJL);
2126 Printf("Leading Jet");
2136 PostData(1, fOutputList);
2137 if(fDebug) Printf("End of event #%5d", (Int_t) fEntry);
2144 nJ = Jets[1]->GetEntries();
2146 printf("NEW EVENT 5: Number of Rec. jets %d\n",nJ);
2151 Float_t areaJL_MC=0;
2153 for (Int_t i = 0; i < nJ; i++) {
2154 AliAODJet* jet = dynamic_cast<AliAODJet*> (Jets[1]->At(i));
2156 Float_t ptJ1 = jet->Pt();
2157 Float_t etaJ1 = TMath::Abs(jet->Eta());
2159 areaJL_MC = jet->EffectiveAreaCharged();
2160 Float_t ptJcorr=ptJ1-rho*areaJL_MC;
2164 if ((ptJcorr > fPtJcorrMin) && (ptJcorr > ptmax) && etaJ1 < etaJmax) {
2172 printf("NEW EVENT 6:\n");
2176 AliAODJet* jetMCL=0;
2179 jetMCL = dynamic_cast<AliAODJet*> (Jets[1]->At(imax));
2181 fanalJetSubStrMCHM->AddEvent(aodE, jetMCL, rho_MC*areaJL_MC);
2189 PostData(1, fOutputList);
2191 if(fDebug) Printf("End of event #%5d", (Int_t) fEntry);
2212 Double_t AliAnalysisTaskJetShape::CalcdPhi(Double_t phi1, Double_t phi2)
2214 while(phi1<0) phi1+=TMath::TwoPi();
2215 while(phi1>TMath::TwoPi()) phi1-=TMath::TwoPi();
2217 while(phi2<0) phi2+=TMath::TwoPi();
2218 while(phi2>TMath::TwoPi()) phi2-=TMath::TwoPi();
2220 Double_t dphi = phi1- phi2;
2221 if(dphi>TMath::Pi())dphi = dphi - 2.*TMath::Pi();
2222 if(dphi<(-1.*TMath::Pi()))dphi = dphi + 2.*TMath::Pi();
2224 // Double_t dphi = phi1- phi2;
2225 // while(dphi<0) dphi+=TMath::TwoPi();
2226 // while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2230 Double_t AliAnalysisTaskJetShape::CalcdPhi0To2pi(Double_t phi1, Double_t phi2)
2233 Double_t dphi = CalcdPhi(phi1, phi2);
2234 while(dphi<0) dphi+=TMath::TwoPi();
2235 while(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
2240 Bool_t AliAnalysisTaskJetShape::IsGoodEvent()
2243 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
2245 // AliError("ESD not available");
2246 // fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
2247 // fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
2249 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
2250 // assume that the AOD is in the general output...
2251 fAODOut = AODEvent();
2254 static AliAODEvent* aod = 0;
2255 // take all other information from the aod we take the tracks from
2257 if(!fESD)aod = fAODIn;
2258 else aod = fAODOut;}
2262 if(fNonStdFile.Length()!=0){
2263 // case that we have an AOD extension we need can fetch the jets from the extended output
2264 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2265 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
2269 if(fAODIn) printf("fAODIn\n");
2270 if(fAODOut) printf("fAODOut\n");
2271 if(fAODExtension) printf("fAODExtension\n");
2276 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2278 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
2283 // -- event selection --
2285 // physics selection
2286 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
2287 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2288 // cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
2289 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
2290 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
2296 // if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
2297 // fhNEvent->Fill(3);
2298 // PostData(1, fOutputList);
2302 AliAODVertex* primVtx = aod->GetPrimaryVertex();
2305 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
2309 Int_t nTracksPrim = primVtx->GetNContributors();
2310 if ((nTracksPrim < fVtxMinContrib) ||
2311 (primVtx->GetZ() < fVtxZMin) ||
2312 (primVtx->GetZ() > fVtxZMax) ){
2313 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
2318 // event class selection (from jet helper task)
2319 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
2320 if(fDebug) Printf("Event class %d", eventClass);
2321 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
2333 //________________________________________________________________________
2334 void AliAnalysisTaskJetShape::Terminate(Option_t *)
2336 // Draw result to the screen
2337 // Called once at the end of the query
2338 // fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
2340 printf("MyTaskTestTrigger::Terminate()\n\n\n");
2343 // fOutputList = dynamic_cast<TList*> (GetOutputData(1));
2344 // fOutputTree = dynamic_cast<TTree*> (GetOutputData(1));
2351 Bool_t AnalysisJetMain::GetVertex(const AliESDEvent* esd, Double_t Vxyz[3], Int_t type)
2356 const AliESDVertex* vtx = 0;
2359 vtx = esd->GetPrimaryVertexSPD();
2360 if(!vtx) return kFALSE;
2361 if(vtx->GetNContributors() < 1) return kFALSE;
2362 TString vtxTyp = vtx->GetTitle();
2363 if (vtxTyp.Contains("vertexer: Z")) {
2364 if (vtx->GetDispersion()>0.04) return kFALSE;
2365 if (vtx->GetZRes()>0.25) return kFALSE;
2370 vtx = esd->GetPrimaryVertexTracks();
2371 if(!vtx) return kFALSE;
2372 if(vtx->GetNContributors() < 1) return kFALSE;
2376 Vxyz[0] = vtx->GetXv();
2377 Vxyz[1] = vtx->GetYv();
2378 Vxyz[2] = vtx->GetZv();
2392 TH1F* AliAnalysisTaskJetShape::Hist1D(const char* name, Int_t nBins, Double_t xMin, Double_t xMax, const char* xLabel, Int_t color)
2394 // create a 1D histogram
2396 TH1F *h = new TH1F(name, name, nBins, xMin, xMax);
2397 if (xLabel) h->GetXaxis()->SetTitle(xLabel);
2398 // res->SetMarkerStyle(kFullCircle);
2399 // res->SetOption("E");
2400 h->SetLineColor(color);
2401 // fOutputList->Add(h);
2406 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)
2408 // create a 2D histogram
2410 TH2F *res = new TH2F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax);
2411 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2412 if (xLabel) res->GetYaxis()->SetTitle(yLabel);
2413 // res->SetMarkerStyle(kFullCircle);
2414 // res->SetOption("E");
2415 res->SetLineColor(color);
2416 // fOutputList->Add(res);
2421 TH3F *AliAnalysisTaskJetShape::Hist3D(const char* name, Int_t nBinsx, Double_t xMin, Double_t xMax,
2422 Int_t nBinsy, Double_t yMin, Double_t yMax,
2423 Int_t nBinsz, Double_t zMin, Double_t zMax, const char* xLabel, const char* yLabel, const char* zLabel, Int_t color)
2425 // create a 3D histogram
2427 TH3F *res = new TH3F(name, name, nBinsx, xMin, xMax, nBinsy, yMin, yMax, nBinsz, zMin, zMax);
2428 if (xLabel) res->GetXaxis()->SetTitle(xLabel);
2429 if (yLabel) res->GetYaxis()->SetTitle(yLabel);
2430 if (zLabel) res->GetZaxis()->SetTitle(zLabel);
2431 // res->SetMarkerStyle(kFullCircle);
2432 // res->SetOption("E");
2433 res->SetLineColor(color);
2434 // fOutputList->Add(res);