28 #include <AliAnalysisTask.h>
29 #include <AliAnalysisManager.h>
30 #include <AliCentrality.h>
31 #include <AliVVertex.h>
32 #include <AliESDEvent.h>
33 #include <AliAODEvent.h>
34 #include <AliAODTrack.h>
63 for(Int_t i(0); i < 10; i++) {
84 for(Int_t i(0); i < 10; i++) {
89 DefineInput(0, TChain::Class());
90 DefineOutput(1, TList::Class());
94 DefineOutput(2, TList::Class());
95 DefineOutput(3, TList::Class());
115 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
116 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
123 if(!(InputEvent()->FindListObject(
fLocalRho->GetName()))) {
126 AliFatal(Form(
"%s: Container with same name %s already present. Aborting", GetName(),
fLocalRho->GetName()));
134 AliFatal(Form(
"%s: Couldn't find container for scaled rho. Aborting !", GetName()));
137 if(!
GetJetContainer()) AliFatal(Form(
"%s: Couldn't get jet container. Aborting !", GetName()));
144 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
145 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
152 SetModulationFit(
new TF1(
"fit_kV2",
"[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
159 SetModulationFit(
new TF1(
"fit_kV3",
"[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
166 SetModulationFit(
new TF1(
"fit_kCombined",
"[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi()));
187 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
188 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
190 fHistSwap =
new TH1F(
"fHistSwap",
"fHistSwap", 20, 0, TMath::TwoPi());
192 Int_t c[] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100};
204 fHistPsi2[i] =
BookTH1F(
"fHistPsi2",
"#Psi_{2}", 100, -.5*TMath::Pi(), .5*TMath::Pi(), i);
205 fHistPsi3[i] =
BookTH1F(
"fHistPsi3",
"#Psi_{3}", 100, -1.*TMath::Pi()/3., TMath::Pi()/3., i);
257 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
258 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
263 name = Form(
"%s_%i", name, c);
266 title += Form(
";%s;[counts]", x);
267 TH1F* histogram =
new TH1F(name, title.Data(), bins, min, max);
275 Int_t binsy, Double_t miny, Double_t maxy, Int_t c, Bool_t append)
278 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
279 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
284 name = Form(
"%s_%i", name, c);
287 title += Form(
";%s;%s", x, y);
288 TH2F* histogram =
new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
298 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
299 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
305 Double_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile(
"V0M"));
315 Double_t psi2(-1), psi3(-1);
321 psi2 = tpc[0]; psi3 = tpc[1];
326 Double_t vzero[2][2];
328 psi2 = vzero[0][0]; psi3 = vzero[0][1];
333 Double_t vzero[2][2];
335 psi2 = vzero[1][0]; psi3 = vzero[1][1];
341 Double_t vzeroComb[2];
343 psi2 = vzeroComb[0]; psi3 = vzeroComb[1];
438 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
439 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
448 Double_t a(0), b(0), c(0), d(0), e(0), f(0), g(0), h(0);
449 vzero[0][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, a, b);
450 vzero[1][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, c, d);
451 vzero[0][1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, e, f);
452 vzero[1][1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, g, h);
456 if(fDebug > 0) printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
457 Double_t qxa2(0), qya2(0), qxc2(0), qyc2(0);
458 Double_t qxa3(0), qya3(0), qxc3(0), qyc3(0);
459 for(Int_t iVZERO(0); iVZERO < 64; iVZERO++) {
460 Double_t phi(TMath::PiOver4()*(.5+iVZERO%8)), weight(InputEvent()->GetVZEROEqMultiplicity(iVZERO));
463 qxa2 += weight*TMath::Cos(2.*phi);
464 qya2 += weight*TMath::Sin(2.*phi);
465 qxa3 += weight*TMath::Cos(3.*phi);
466 qya3 += weight*TMath::Sin(3.*phi);
469 qxc2 += weight*TMath::Cos(2.*phi);
470 qyc2 += weight*TMath::Sin(2.*phi);
471 qxc3 += weight*TMath::Cos(3.*phi);
472 qyc3 += weight*TMath::Sin(3.*phi);
475 vzero[0][0] = .5*TMath::ATan2(qya2, qxa2);
476 vzero[1][0] = .5*TMath::ATan2(qyc2, qxc2);
477 vzero[0][1] = (1./3.)*TMath::ATan2(qya3, qxa3);
478 vzero[1][1] = (1./3.)*TMath::ATan2(qyc3, qxc3);
487 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
488 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
491 Double_t qx2(0), qy2(0);
492 Double_t qx3(0), qy3(0);
494 Float_t excludeInEta = -999;
497 if(leadingJet) leadingJet->
Eta();
499 Int_t iTracks(
fTracks->GetEntriesFast());
500 for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
501 AliVTrack* track =
static_cast<AliVTrack*
>(
fTracks->At(iTPC));
505 qx2+= TMath::Cos(2.*track->Phi());
506 qy2+= TMath::Sin(2.*track->Phi());
507 qx3+= TMath::Cos(3.*track->Phi());
508 qy3+= TMath::Sin(3.*track->Phi());
511 tpc[0] = .5*TMath::ATan2(qy2, qx2);
512 tpc[1] = (1./3.)*TMath::ATan2(qy3, qx3);
519 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
520 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
522 Double_t a(0), b(0), c(0), d(0);
523 comb[0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, a, b);
524 comb[1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, c, d);
531 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
532 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
534 Double_t reQ(0), imQ(0), modQ(0), M11(0), M(0);
536 QCnQnk(harm, 1, reQ, imQ);
537 modQ = reQ*reQ+imQ*imQ;
539 return (M11 > 0) ? ((modQ -
QCnS(1,2))/M11) : -999;
541 QCnQnk(harm, 0, reQ, imQ);
542 modQ = reQ*reQ+imQ*imQ;
544 return (M > 1) ? (modQ - M)/(M*(M-1)) : -999;
551 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
552 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
554 Double_t reQn1(0), imQn1(0), reQ2n2(0), imQ2n2(0), reQn3(0), imQn3(0), M1111(0), M(0);
555 Double_t a(0), b(0), c(0), d(0), e(0), f(0), g(0);
557 QCnQnk(harm, 1, reQn1, imQn1);
558 QCnQnk(harm*2, 2, reQ2n2, imQ2n2);
559 QCnQnk(harm, 3, reQn3, imQn3);
561 a = (reQn1*reQn1+imQn1*imQn1)*(reQn1*reQn1+imQn1*imQn1);
562 b = reQ2n2*reQ2n2 + imQ2n2*imQ2n2;
563 c = -2.*(reQ2n2*reQn1*reQn1-reQ2n2*imQn1*imQn1+2.*imQ2n2*reQn1*imQn1);
564 d = 8.*(reQn3*reQn1+imQn3*imQn1);
565 e = -4.*
QCnS(1,2)*(reQn1*reQn1+imQn1*imQn1);
569 return (M1111 > 0) ? (a+b+c+d+e+f+g)/M1111 : -999;
571 Double_t reQn(0), imQn(0), reQ2n(0), imQ2n(0);
572 QCnQnk(harm, 0, reQn, imQn);
573 QCnQnk(harm*2, 0, reQ2n, imQ2n);
576 if(M < 4)
return -999;
577 a = (reQn*reQn+imQn*imQn)*(reQn*reQn+imQn*imQn);
578 b = reQ2n*reQ2n + imQ2n*imQ2n;
579 c = -2.*(reQ2n*reQn*reQn-reQ2n*imQn*imQn+2.*imQ2n*reQn*imQn);
580 e = -4.*(M-2)*(reQn*reQn+imQn*imQn);
582 return (a+b+c+e+f)/(M*(M-1)*(M-2)*(M-3));
589 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
590 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
594 Int_t iTracks(
fTracks->GetEntriesFast());
595 for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
596 AliVTrack* track =
static_cast<AliVTrack*
>(
fTracks->At(iTPC));
600 reQ += TMath::Power(track->Pt(), k) * TMath::Cos(((
double)n)*track->Phi());
601 imQ += TMath::Power(track->Pt(), k) * TMath::Sin(((
double)n)*track->Phi());
609 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
610 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
612 if(!
fTracks || i <= 0 || j <= 0)
return -999;
613 Int_t iTracks(
fTracks->GetEntriesFast());
615 for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
616 AliVTrack* track =
static_cast<AliVTrack*
>(
fTracks->At(iTPC));
618 Sij+=TMath::Power(track->Pt(), j);
620 return TMath::Power(Sij, i);
627 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
628 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
637 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
638 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
647 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
648 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
658 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
659 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
677 if( c2 < 0 ) c2 = -1.*TMath::Sqrt(-1.*c2);
678 if( c3 < 0 ) c3 = -1.*TMath::Sqrt(-1.*c3);
692 default :
return kFALSE;
715 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
716 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
742 }
else if (!
QCnRecovery(psi2, psi3))
return kFALSE;
772 }
else if (!
QCnRecovery(psi2, psi3))
return kFALSE;
796 TString detector(
"");
798 case kTPC : detector+=
"TPC";
800 case kVZEROA : detector+=
"VZEROA";
802 case kVZEROC : detector+=
"VZEROC";
808 Int_t iTracks(
fTracks->GetEntriesFast());
809 Double_t excludeInEta = -999;
810 Double_t excludeInPhi = -999;
811 Double_t excludeInPt = -999;
812 if(iTracks <= 0 || fLocalRho->GetVal() <= 0 )
return kFALSE;
816 excludeInEta = leadingJet->
Eta();
817 excludeInPhi = leadingJet->
Phi();
818 excludeInPt = leadingJet->
Pt();
826 _tempSwap = TH1F(
"_tempSwap",
"_tempSwap", TMath::CeilNint(TMath::Sqrt(
fNAcceptedTracks)), 0, TMath::TwoPi());
832 Double_t totalpts(0.), totalptsquares(0.), totalns(0.);
833 for(Int_t i(0); i < iTracks; i++) {
834 AliVTrack* track =
static_cast<AliVTrack*
>(
fTracks->At(i));
838 _tempSwap.Fill(track->Phi(), track->Pt());
840 totalpts += track->Pt();
841 totalptsquares += track->Pt()*track->Pt();
843 _tempSwapN.Fill(track->Phi());
846 else _tempSwap.Fill(track->Phi());
853 if(totalns < 1)
return kFALSE;
854 for(Int_t l = 0; l < _tempSwap.GetNbinsX(); l++) {
855 if(_tempSwapN.GetBinContent(l+1) == 0) {
856 _tempSwap.SetBinContent(l+1,0);
857 _tempSwap.SetBinError(l+1,0);
860 Double_t vartimesnsq = totalptsquares*totalns - totalpts*totalpts;
861 Double_t variance = vartimesnsq/(totalns*(totalns-1.));
862 Double_t SDOMSq = variance / _tempSwapN.GetBinContent(l+1);
863 Double_t SDOMSqOverMeanSq = SDOMSq * _tempSwapN.GetBinContent(l+1) * _tempSwapN.GetBinContent(l+1) / (_tempSwapN.GetBinContent(l+1) * _tempSwapN.GetBinContent(l+1));
864 Double_t poissonfrac = 1./_tempSwapN.GetBinContent(l+1);
865 Double_t vartotalfrac = SDOMSqOverMeanSq + poissonfrac;
866 Double_t vartotal = vartotalfrac * _tempSwap.GetBinContent(l+1) * _tempSwap.GetBinContent(l+1);
867 if(vartotal > 0.0001) _tempSwap.SetBinError(l+1,TMath::Sqrt(vartotal));
869 _tempSwap.SetBinContent(l+1,0);
870 _tempSwap.SetBinError(l+1,0);
893 Double_t cos2(0), sin2(0), cos3(0), sin3(0), sumPt(0);
894 for(Int_t i(0); i < iTracks; i++) {
895 AliVTrack* track =
static_cast<AliVTrack*
>(
fTracks->At(i));
897 sumPt += track->Pt();
898 cos2 += track->Pt()*TMath::Cos(2*
PhaseShift(track->Phi()-psi2));
899 sin2 += track->Pt()*TMath::Sin(2*
PhaseShift(track->Phi()-psi2));
900 cos3 += track->Pt()*TMath::Cos(3*
PhaseShift(track->Phi()-psi3));
901 sin3 += track->Pt()*TMath::Sin(3*
PhaseShift(track->Phi()-psi3));
922 static Int_t didacticCounterBest(0);
923 TProfile* didacticProfile = (TProfile*)_tempSwap.Clone(Form(
"Fit_%i_1-CDF_%.3f_cen_%i_%s", didacticCounterBest, CDF,
fInCentralitySelection, detector.Data()));
925 didacticProfile->GetListOfFunctions()->Add(didactifFit);
927 didacticCounterBest++;
928 TH2F* didacticSurface =
BookTH2F(Form(
"surface_%s", didacticProfile->GetName()),
"#phi",
"#eta", 50, 0, TMath::TwoPi(), 50, -1, 1, -1, kFALSE);
929 for(Int_t i(0); i < iTracks; i++) {
930 AliVTrack* track =
static_cast<AliVTrack*
>(
fTracks->At(i));
932 if(
fUsePtWeight) didacticSurface->Fill(track->Phi(), track->Eta(), track->Pt());
933 else didacticSurface->Fill(track->Phi(), track->Eta());
937 TF2 *f2 =
new TF2(Form(
"%s_LJ", didacticSurface->GetName()),
"[0]*TMath::Gaus(x,[1],[2])*TMath::Gaus(y,[3],[4])", 0, TMath::TwoPi(), -1, 1);
938 f2->SetParameters(excludeInPt/3.,excludeInPhi,.1,excludeInEta,.1);
939 didacticSurface->GetListOfFunctions()->Add(f2);
949 static Int_t didacticCounterWorst(0);
951 TProfile* didacticProfile = (TProfile*)_tempSwap.Clone(Form(
"Fit_%i_1-CDF_%.3f_cen_%i_%s", didacticCounterWorst, CDF,
fInCentralitySelection, detector.Data() ));
953 didacticProfile->GetListOfFunctions()->Add(didactifFit);
955 didacticCounterWorst++;
977 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
978 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1026 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
1027 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1036 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
1037 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1046 #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
1047 printf(
"__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
Double_t QCnS(Int_t i, Int_t j)
void CalculateEventPlaneCombinedVZERO(Double_t *comb) const
detectorType fDetectorType
Float_t fPercentageOfFits
void CalculateEventPlaneTPC(Double_t *tpc)
virtual void Terminate(Option_t *option)
virtual void UserCreateOutputObjects()
Bool_t PassesCuts(AliVTrack *track) const
AliJetContainer * GetJetContainer(Int_t i=0) const
Bool_t QCnRecovery(Double_t psi2, Double_t psi3)
Bool_t fAbsVnHarmonics
status of rho vs centrality
Double_t ChiSquareCDF(Int_t ndf, Double_t x) const
Double_t GetJetEtaMax() const
Bool_t fRebinSwapHistoOnTheFly
TH1F * fHistAnalysisSummary
swap histogram
Float_t fExcludeLeadingJetsFromFit
virtual ~AliAnalysisTaskLocalRho()
TH1F * fHistPsi2[10]
v3 cumulant
ClassImp(AliAnalysisTaskLocalRho) AliAnalysisTaskLocalRho
Double_t GetJetRadius(Int_t i=0) const
void CalculateEventPlaneVZERO(Double_t vzero[2][2]) const
TString fLocalRhoName
name for local rho
Bool_t CorrectRho(Double_t psi2, Double_t psi3)
TArrayI * fCentralityClasses
Int_t fNAcceptedTracksQCn
number of accepted tracks
AliEmcalJet * GetLeadingJet(const char *opt="")
void SetJetEtaLimits(Float_t min, Float_t max, Int_t c=0)
AliRhoParameter * fRho
! event rho
TList * fOutputListBad
output list for local analysis
Double_t CalculateQC4(Int_t harm)
void QCnQnk(Int_t n, Int_t k, Double_t &reQ, Double_t &imQ)
Double_t fCent
!event centrality
AliLocalRhoParameter * fLocalRho
! local event rho
Bool_t fUseV0EventPlaneFromHeader
fitModulationType fFitModulationType
centrality bin, only for QA plots
TClonesArray * fJets
! jets
void SetJetPhiLimits(Float_t min, Float_t max, Int_t c=0)
TList * fOutputListGood
output list
TH2F * fHistRhoStatusCent
cdf value of chisquare p
TProfile * fProfV3
v2 cumulant
TString fFitModulationOptions
Bool_t InitializeAnalysis()
Float_t GetJetRadius() const
void FillEventPlaneHistograms(Double_t psi2, Double_t psi3) const
TH2F * BookTH2F(const char *name, const char *x, const char *y, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t c=-1, Bool_t append=kTRUE)
TClonesArray * fTracks
!tracks
TProfile * fProfV3Cumulant
extracted v3
void FillAnalysisSummaryHistogram() const
void SetModulationFit(TF1 *fit)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Bool_t fUsePtWeightErrorPropagation
TH1F * fHistSwap
output list for local analysis
Int_t fInCentralitySelection
accepted tracks for QCn
TProfile * fProfV2Cumulant
extracted v2
Double_t PhaseShift(Double_t x) const
TH1F * BookTH1F(const char *name, const char *x, Int_t bins, Double_t min, Double_t max, Int_t c=-1, Bool_t append=kTRUE)
Bool_t fAttachToEvent
is the analysis initialized?
AliAnalysisTaskLocalRho()
Double_t CalculateQC2(Int_t harm)
Bool_t fNoEventWeightsForQC
TH1F * fHistPsi3[10]
psi 2