FLOW/Tasks/AliFlowVZEROQA.cxx
FLOW/Tasks/AliAnalysisTaskFlowEPCascade.cxx
FLOW/Tasks/AliAnalysisTaskJetFlow.cxx
+ FLOW/Tasks/AliAnalysisTaskJetFlowMC.cxx
)
string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
--- /dev/null
+// Simple class to generate toy events which can be fed into the jet finder\r
+//\r
+// Author: Redmer Alexander Bertens, Utrecht University, 2013\r
+// rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl\r
+\r
+// root includes\r
+#include "TChain.h"\r
+#include "TList.h"\r
+#include "TClonesArray.h"\r
+#include "TArrayI.h"\r
+// aliroot includes\r
+#include "AliAODEvent.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliAnalysisTaskJetFlowMC.h"\r
+#include "AliLog.h"\r
+#include "AliPicoTrack.h"\r
+\r
+class AliAnalysisTaskJetFlowMC;\r
+using namespace std;\r
+\r
+ClassImp(AliAnalysisTaskJetFlowMC)\r
+\r
+//_____________________________________________________________________________\r
+AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC() : AliAnalysisTaskSE("AliAnalysisTaskJetFlowMC"), fTracksOutName("JetFlowMC"), fTracksInName("PicoTrack"), fTracksIn(0), fTracksOut(0), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), fJetSpectrumSF(0), fNoOfSFJets(0), fHistIntV2(0), fHistIntV3(0), fFlowFluctuations(-10), fMaxNumberOfIterations(100), fPsi2(-10), fPsi3(-10), fPrecisionPhi(1e-10), fDetectorType(kVZEROC), fHistSFJetSpectrum(0), fHistSFJetEtaPhi(0) {\r
+ // default constructor for root IO\r
+ for(Int_t i(0); i < 10; i++) {\r
+ fFuncDiffV2[i] = 0x0;\r
+ fFuncDiffV3[i] = 0x0;\r
+ fHistDiffV2[i] = 0x0;\r
+ fHistDiffV3[i] = 0x0;\r
+ fHistOriginalSpectrum[i] = 0x0;\r
+ fHistOriginalEtaPhi[i] = 0x0;\r
+ fHistToySpectrum[i] = 0x0;\r
+ fHistToyEtaPhi[i] = 0x0;\r
+ fHistOriginalDeltaPhi[i] = 0x0;\r
+ fHistToyDeltaPhi[i] = 0x0;\r
+ fHistToyVn[i] = 0x0;\r
+ }\r
+}\r
+//_____________________________________________________________________________\r
+AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC(const char *name) : AliAnalysisTaskSE(name), fTracksOutName("JetFlowMC"), fTracksInName("PicoTrack"), fTracksIn(0), fTracksOut(0), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), fJetSpectrumSF(0), fNoOfSFJets(0), fHistIntV2(0), fHistIntV3(0), fFlowFluctuations(-10), fMaxNumberOfIterations(100), fPsi2(-10), fPsi3(-10), fPrecisionPhi(1e-10), fDetectorType(kVZEROC), fHistSFJetSpectrum(0), fHistSFJetEtaPhi(0) {\r
+ // constructor\r
+ DefineInput(0, TChain::Class());\r
+ DefineOutput(1, TList::Class());\r
+ for(Int_t i(0); i < 10; i++) {\r
+ fFuncDiffV2[i] = 0x0;\r
+ fFuncDiffV3[i] = 0x0;\r
+ fHistDiffV2[i] = 0x0;\r
+ fHistDiffV3[i] = 0x0;\r
+ fHistOriginalSpectrum[i] = 0x0;\r
+ fHistOriginalEtaPhi[i] = 0x0;\r
+ fHistToySpectrum[i] = 0x0;\r
+ fHistToyEtaPhi[i] = 0x0;\r
+ fHistOriginalDeltaPhi[i] = 0x0;\r
+ fHistToyDeltaPhi[i] = 0x0;\r
+ fHistToyVn[i] = 0x0;\r
+ }\r
+}\r
+//_____________________________________________________________________________\r
+AliAnalysisTaskJetFlowMC::~AliAnalysisTaskJetFlowMC()\r
+{\r
+ // desctructor, claim ownership of deleted objects by setting member pointers to NULL\r
+ if(fOutputList) {\r
+ delete fOutputList;\r
+ fOutputList= NULL;\r
+ }\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::UserCreateOutputObjects()\r
+{\r
+ // Create my user objects.\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);\r
+ fTracksOut = new TClonesArray("AliPicoTrack");\r
+ fTracksOut->SetName(fTracksOutName);\r
+ fOutputList = new TList();\r
+ fOutputList->SetOwner(kTRUE);\r
+ if(!fCentralityClasses) { // classes must be defined at this point\r
+ Int_t c[] = {0, 90};\r
+ fCentralityClasses = new TArrayI(sizeof(c)/sizeof(c[0]), c);\r
+ }\r
+ if(fHistIntV2 && !fHistIntV3) { // define function\r
+ fFuncVn = new TF1("kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());\r
+ fFuncVn->SetParameter(0, 1.); // normalization\r
+ fFuncVn->SetParameter(3, 0.2); // v2\r
+ fFuncVn->FixParameter(1, 1.); // constant\r
+ fFuncVn->FixParameter(2, 2.); // constant\r
+ } else if (!fHistIntV2 && fHistIntV3) {\r
+ fFuncVn = new TF1("kV3", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());\r
+ fFuncVn->SetParameter(0, 1.); // normalization\r
+ fFuncVn->SetParameter(3, 0.2); // v3\r
+ fFuncVn->FixParameter(1, 1.); // constant\r
+ fFuncVn->FixParameter(2, 3.); // constant\r
+ } else if (fHistIntV2 && fHistIntV3) {\r
+ fFuncVn = new TF1("kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi());\r
+ fFuncVn->SetParameter(0, 1.); // normalization\r
+ fFuncVn->SetParameter(3, 0.2); // v2\r
+ fFuncVn->FixParameter(1, 1.); // constant\r
+ fFuncVn->FixParameter(2, 2.); // constant\r
+ fFuncVn->FixParameter(5, 3.); // constant\r
+ fFuncVn->SetParameter(7, 0.2); // v3\r
+ }\r
+ // add the generator objects that have been added to the task\r
+ if(fTrackSpectrum) fOutputList->Add(fTrackSpectrum);\r
+ if(fJetSpectrumSF) fOutputList->Add(fJetSpectrumSF);\r
+ if(fHistIntV2) fOutputList->Add(fHistIntV2);\r
+ if(fHistIntV3) fOutputList->Add(fHistIntV3);\r
+ // create the QA histos\r
+ for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {\r
+ fHistOriginalSpectrum[i] = BookTH1F("fHistOriginalSpectrum", "p_{t} [GeV/c]", 100, 0, 100, i);\r
+ fHistOriginalEtaPhi[i] = BookTH2F("fHistOriginalEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);\r
+ fHistToySpectrum[i] = BookTH1F("fHistToySpectrum", "p_{t} [GeV/c]", 100, 0, 100, i);\r
+ fHistToyEtaPhi[i] = BookTH2F("fHistToyEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);\r
+ fHistOriginalDeltaPhi[i] = BookTH1F("fHistOriginalDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);\r
+ fHistToyDeltaPhi[i] = BookTH1F("fHistToyDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);\r
+ fHistToyVn[i] = BookTH2F("fHistToyVn", "p_{t} [GeV/c]", "v_{n}", 100, 0, 20, 80, 0, .8, i);\r
+ // add to outputlist\r
+ if(fFuncDiffV2[i]) fOutputList->Add(fFuncDiffV2[i]);\r
+ if(fFuncDiffV3[i]) fOutputList->Add(fFuncDiffV3[i]);\r
+ if(fHistDiffV2[i]) fOutputList->Add(fHistDiffV2[i]);\r
+ if(fHistDiffV3[i]) fOutputList->Add(fHistDiffV3[i]);\r
+ }\r
+ if(fJetSpectrumSF) {\r
+ fHistSFJetSpectrum = BookTH1F("fHistSFJetSpectrum", "p_{t} SF jets [GeV/c]", 100, 0, 200);\r
+ fHistSFJetEtaPhi = BookTH2F("fHistSFJetEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi());\r
+ }\r
+ fOutputList->Sort();\r
+ PostData(1, fOutputList);\r
+}\r
+//_____________________________________________________________________________\r
+TH1F* AliAnalysisTaskJetFlowMC::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c, Bool_t append)\r
+{\r
+ // book a TH1F and connect it to the output container\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); if(!fOutputList) return 0x0;\r
+ TString title(name);\r
+ if(c!=-1) { // format centrality dependent histograms accordingly\r
+ name = Form("%s_%i", name, c);\r
+ title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));\r
+ }\r
+ title += Form(";%s;[counts]", x);\r
+ TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);\r
+ histogram->Sumw2();\r
+ if(append) fOutputList->Add(histogram);\r
+ return histogram; \r
+}\r
+//_____________________________________________________________________________\r
+TH2F* AliAnalysisTaskJetFlowMC::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, Bool_t append)\r
+{\r
+ // book a TH2F and connect it to the output container\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); if(!fOutputList) return 0x0;\r
+ TString title(name);\r
+ if(c!=-1) { // format centrality dependent histograms accordingly\r
+ name = Form("%s_%i", name, c);\r
+ title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));\r
+ }\r
+ title += Form(";%s;%s", x, y);\r
+ TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);\r
+ histogram->Sumw2();\r
+ if(append) fOutputList->Add(histogram);\r
+ return histogram; \r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::UserExec(Option_t *) \r
+{\r
+ // user exec, called for each event.\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); if(!AliAnalysisManager::GetAnalysisManager()) return;\r
+ // retrieve tracks from input.\r
+ if (!fTracksIn) { \r
+ fTracksIn = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksInName));\r
+ if (!fTracksIn) {\r
+ AliError(Form("Could not retrieve tracks %s!", fTracksInName.Data())); \r
+ return;\r
+ }\r
+ }\r
+ // get the centrality bin for qa plots\r
+ Double_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));\r
+ fCenBin = -1;\r
+ for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {\r
+ if(cent >= fCentralityClasses->At(i) && cent <= fCentralityClasses->At(1+i)) {\r
+ fCenBin = i;\r
+ break; }\r
+ }\r
+ if(fCenBin < 0) return;\r
+ // add tracks to event if not yet there\r
+ fTracksOut->Delete();\r
+ if (!(InputEvent()->FindListObject(fTracksOutName))) InputEvent()->AddObject(fTracksOut);\r
+ fTracksOut->Clear();\r
+ // get the event plane\r
+ CalculateEventPlane();\r
+ const Int_t Ntracks = fTracksIn->GetEntriesFast();\r
+ Int_t nacc(0);\r
+ for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {\r
+ AliPicoTrack* track = static_cast<AliPicoTrack*>(fTracksIn->At(iTracks));\r
+ if(!track) continue;\r
+ Double_t phi(track->Phi()), pt((fTrackSpectrum) ? GetTrackPt() : track->Pt()), eta(track->Eta());\r
+ // fill qa histo's before applying any (possible) afterburner\r
+ FillHistogramsOriginalData(pt, eta, phi);\r
+ if(fHistDiffV2[fCenBin] || fFuncDiffV2[fCenBin]) V2AfterBurner(phi, eta, pt);\r
+ else if(fHistDiffV3[fCenBin] || fFuncDiffV3[fCenBin]) V3AfterBurner(phi, eta, pt);\r
+ else if(fHistDiffV2 || fHistDiffV3) SampleVnFromTF1(phi); \r
+ /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, track->Charge(), track->GetLabel(), 4, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal(), 1); \r
+ nacc++;\r
+ }\r
+ if(fJetSpectrumSF && fNoOfSFJets > 0) InjectSingleFragmentationJetSpectrum(nacc);\r
+ PostData(1, fOutputList);\r
+ if(fDebug > 0) PrintInfo();\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::V2AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const\r
+{\r
+ // similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations\r
+ if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); phi = PhaseShift(gRandom->Uniform(0, TMath::TwoPi()) + fPsi2);\r
+ Double_t phi0(phi), v2(GetV2(pt)), f(0.), fp(0.), phiprev(0.);\r
+ if(TMath::AreEqualAbs(v2, 0, 1e-5)) { \r
+ FillHistogramsToyData(pt, eta, phi, v2);\r
+ return;\r
+ }\r
+ // introduce flow fluctuations (gaussian)\r
+ if(fFlowFluctuations > -10) GetFlowFluctuation(v2);\r
+ for (Int_t i=0; i!=fMaxNumberOfIterations; ++i) {\r
+ phiprev=phi; //store last value for comparison\r
+ f = phi-phi0+v2*TMath::Sin(2.*(phi-fPsi2));\r
+ fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi2)); //first derivative\r
+ phi -= f/fp;\r
+ if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break; \r
+ }\r
+ FillHistogramsToyData(pt, eta, phi, v2);\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::V3AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const\r
+{\r
+ // similar to AliFlowTrackSimple::AddV3, except for the flow fluctuations\r
+ if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); phi = PhaseShift(gRandom->Uniform(0, TMath::TwoPi()) + fPsi3);\r
+ Double_t phi0(phi), v3(GetV3(pt)), f(0.), fp(0.), phiprev(0.);\r
+ if(TMath::AreEqualAbs(v3, 0, 1e-5)) {\r
+ FillHistogramsToyData(pt, eta, phi, v3);\r
+ return;\r
+ }\r
+ // introduce flow fluctuations (gaussian)\r
+ if(fFlowFluctuations > -10) GetFlowFluctuation(v3);\r
+ for (Int_t i=0; i<fMaxNumberOfIterations; i++) {\r
+ phiprev=phi; //store last value for comparison\r
+ f = phi-phi0+2./3.*v3*TMath::Sin(3.*(phi-fPsi3));\r
+ fp = 1.0+2.0*v3*TMath::Cos(3.*(phi-fPsi3)); //first derivative\r
+ phi -= f/fp;\r
+ if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;\r
+ }\r
+ FillHistogramsToyData(pt, eta, phi, v3);\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::InjectSingleFragmentationJetSpectrum(Int_t nacc)\r
+{\r
+ // inject single fragmentation jet spectrum to the tclones array, note that emcal params \r
+ // equal the barrel kinematics to pass the track and jet cuts later on\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); for(Int_t i(nacc); i < (nacc + fNoOfSFJets); i++) {\r
+ Double_t eta(gRandom->Uniform(-.5, .5)), phi(gRandom->Uniform(0, TMath::TwoPi())), pt(fJetSpectrumSF->GetRandom());\r
+ /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[i]) AliPicoTrack(pt, eta, phi, +1, 0, 0, eta, phi, pt, 0);\r
+ fHistSFJetSpectrum->Fill(pt);\r
+ fHistSFJetEtaPhi->Fill(eta, phi);\r
+ ++i;\r
+ }\r
+ nacc = 0;\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::CalculateEventPlane() {\r
+ // grab the event plane orientation from the AliVEvent header\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); Double_t a(0), b(0), e(0), f(0);\r
+ switch (fDetectorType) {\r
+ case kVZEROA : {\r
+ fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, e, f);\r
+ fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, a, b);\r
+ break;\r
+ }\r
+ case kVZEROC : {\r
+ fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, e, f);\r
+ fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, a, b);\r
+ break;\r
+ }\r
+ case kVZEROComb : {\r
+ fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, e, f) ;\r
+ fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, a, b);\r
+ break;\r
+ }\r
+ default : break;\r
+ }\r
+ // if requested, pass the event plane to the phi distribution\r
+ if(fHistIntV2 && fHistIntV3) {\r
+ fFuncVn->SetParameter(4, fPsi2); fFuncVn->SetParameter(6, fPsi3);\r
+ Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));\r
+ Double_t v3(fHistIntV2->GetBinContent(fHistIntV3->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));\r
+ if(fFlowFluctuations > -10) {\r
+ GetFlowFluctuation(v2); GetFlowFluctuation(v3);\r
+ }\r
+ fFuncVn->SetParameter(3, v2); fFuncVn->SetParameter(7, v3);\r
+ } else if (fHistIntV2) {\r
+ fFuncVn->SetParameter(4, fPsi2);\r
+ Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));\r
+ if(fFlowFluctuations > -10) GetFlowFluctuation(v2);\r
+ fFuncVn->SetParameter(3, v2);\r
+ } else if (fHistIntV3) {\r
+ fFuncVn->SetParameter(4, fPsi3);\r
+ Double_t v3(fHistIntV3->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));\r
+ if(fFlowFluctuations > -10) GetFlowFluctuation(v3);\r
+ fFuncVn->SetParameter(3, v3);\r
+ }\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::Terminate(Option_t *)\r
+{\r
+ // terminate\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::PrintInfo() const\r
+{\r
+ // print info\r
+ printf(" > No of retrieved tracks from %s \n \t %i \n", fTracksInName.Data(), fTracksIn->GetEntriesFast());\r
+ printf(" > No of created tracks in %s \n \t %i \n", fTracksOutName.Data(), fTracksOut->GetEntriesFast());\r
+\r
+}\r
+//_____________________________________________________________________________\r
--- /dev/null
+#ifndef AliAnalysisTaskJetFlowMC_H\r
+#define AliAnalysisTaskJetFlowMC_H\r
+\r
+// root includes\r
+#include "TF1.h"\r
+#include "TH1F.h"\r
+#include "TH2F.h" \r
+#include "TRandom.h"\r
+// aliroot includes\r
+#include "AliAnalysisTaskSE.h"\r
+// forward declarations\r
+class TList;\r
+class TClonesArray;\r
+class TArrayI;\r
+\r
+class AliAnalysisTaskJetFlowMC : public AliAnalysisTaskSE \r
+{\r
+ public:\r
+ // enumerators\r
+ enum detectorType {kVZEROA, kVZEROC, kVZEROComb}; // detector that was used\r
+ // constructors, destructor\r
+ AliAnalysisTaskJetFlowMC();\r
+ AliAnalysisTaskJetFlowMC(const char *name);\r
+ virtual ~AliAnalysisTaskJetFlowMC();\r
+ // mirror image of PickTrackMaker\r
+ void UserCreateOutputObjects();\r
+ 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);\r
+ 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);\r
+\r
+ void UserExec(Option_t *option);\r
+ void SetDebugMode(Bool_t d) {fDebug = d;}\r
+ void SetTracksInName(const char *name) { fTracksInName = name; }\r
+ void SetTracksOutName(const char *name) { fTracksOutName = name; }\r
+ // additional setters\r
+ void SetCentralityClasses(TArrayI* c) { fCentralityClasses = c; }\r
+ void SetReferenceDetector(detectorType type) { fDetectorType = type; }\r
+ void SetDifferentialV2(TF1* v2, Int_t c = 0) { fFuncDiffV2[c] = v2; }\r
+ void SetDifferentialV3(TF1* v3, Int_t c = 0) { fFuncDiffV3[c] = v3; }\r
+ void SetDifferentialV2(TH1* v2, Int_t c = 0) { fHistDiffV2[c] = v2; }\r
+ void SetDifferentialV3(TH1* v3, Int_t c = 0) { fHistDiffV3[c] = v3; }\r
+ void SetIntegratedV2(TH1* v2) { fHistIntV2 = v2; }\r
+ void SetIntegratedV3(TH1* v3) { fHistIntV3 = v3; }\r
+ void SetTrackSpectrum(TF1* ts) { fTrackSpectrum = ts; }\r
+ void SetSingleFragmentationJetSpectrum(TF1* js) { fJetSpectrumSF = js; }\r
+ void SetNoOfSFJets(Int_t n) { fNoOfSFJets = n; }\r
+ // additional methods\r
+ void V2AfterBurner(Double_t& phi, Double_t& eta, Double_t& pt) const;\r
+ void V3AfterBurner(Double_t& phi, Double_t& eta, Double_t& pt) const;\r
+ void InjectSingleFragmentationJetSpectrum(Int_t nacc);\r
+ void CalculateEventPlane();\r
+ // inlined helper calculations\r
+ Double_t GetTrackPt() const { return fTrackSpectrum->GetRandom();}\r
+ /* inline */ Double_t GetV2(Double_t pt) const { \r
+ return (fFuncDiffV2[fCenBin]) ? fFuncDiffV2[fCenBin]->Eval(pt) :\r
+ fHistDiffV2[fCenBin]->GetBinContent(fHistDiffV2[fCenBin]->GetXaxis()->FindBin(pt));\r
+ }\r
+ /* inline */ Double_t GetV3(Double_t pt) const { \r
+ return (fFuncDiffV3[fCenBin]) ? fFuncDiffV3[fCenBin]->Eval(pt) : \r
+ fHistDiffV3[fCenBin]->GetBinContent(fHistDiffV3[fCenBin]->GetXaxis()->FindBin(pt));\r
+ }\r
+ /* inline */ void GetFlowFluctuation(Double_t& vn) const {\r
+ vn += TMath::Sqrt(2*(vn*.25)*(vn*.25))*TMath::ErfInverse(2*(gRandom->Uniform(0, fFlowFluctuations))-1); \r
+ }\r
+ /* inline */ Double_t PhaseShift(Double_t x) const { \r
+ while (x>=TMath::TwoPi())x-=TMath::TwoPi();\r
+ while (x<0.)x+=TMath::TwoPi();\r
+ return x; }\r
+ /* inline */ Double_t PhaseShift(Double_t x, Double_t n) const {\r
+ x = PhaseShift(x);\r
+ if(TMath::Nint(n)==2) while (x>TMath::Pi()) x-=TMath::Pi();\r
+ if(TMath::Nint(n)==3) {\r
+ if(x>2.*TMath::TwoPi()/n) x = TMath::TwoPi() - x;\r
+ if(x>TMath::TwoPi()/n) x = TMath::TwoPi()-(x+TMath::TwoPi()/n);\r
+ }\r
+ return x; }\r
+ /* inline */ void SampleVnFromTF1(Double_t &phi) const {\r
+ phi = (fFuncVn) ? fFuncVn->GetRandom(0., TMath::TwoPi()) : 0; }\r
+ /* inline */ void FillHistogramsOriginalData(Double_t& pt, Double_t& eta, Double_t& phi) const {\r
+ fHistOriginalSpectrum[fCenBin]->Fill(pt); fHistOriginalEtaPhi[fCenBin]->Fill(eta, phi);\r
+ fHistOriginalDeltaPhi[fCenBin]->Fill(PhaseShift(phi-fPsi2, 2));\r
+ }\r
+ /* inline */ void FillHistogramsToyData(Double_t& pt, Double_t& eta, Double_t& phi, Double_t& vn) const {\r
+ fHistToySpectrum[fCenBin]->Fill(pt); fHistToyEtaPhi[fCenBin]->Fill(eta, phi);\r
+ fHistToyDeltaPhi[fCenBin]->Fill(PhaseShift(phi-fPsi2, 2)); fHistToyVn[fCenBin]->Fill(pt, vn);\r
+ }\r
+ void Terminate(Option_t* option);\r
+ void PrintInfo() const;\r
+ protected:\r
+ TString fTracksOutName; // name of output track array\r
+ TString fTracksInName; // name of input track array\r
+ TClonesArray *fTracksIn; //!track array in\r
+ TClonesArray *fTracksOut; //!track array out\r
+ Int_t fCenBin; //! centrality bin\r
+ TArrayI* fCentralityClasses; // centrality classes (max 10) \r
+ TF1* fFuncVn; //! vn function\r
+ TList* fOutputList; // output list\r
+ TF1* fTrackSpectrum; // track pt spectrum\r
+ TF1* fJetSpectrumSF; // single fragmentation spectrum of jets\r
+ Int_t fNoOfSFJets; // number of single fragmentation jets that will be added\r
+ TF1* fFuncDiffV2[10]; // differential v2 of charged tracks\r
+ TF1* fFuncDiffV3[10]; // differential v3 of charged tracks\r
+ TH1* fHistDiffV2[10]; // differential v2 of charged tracks\r
+ TH1* fHistDiffV3[10]; // differential v3 of charged tracks\r
+ TH1* fHistIntV2; // integrated v2 of charged tracks\r
+ TH1* fHistIntV3; // integrated v3 of charged tracks\r
+ Float_t fFlowFluctuations; // introduce gaussian flow fluctuations of this magnitude\r
+ Int_t fMaxNumberOfIterations; // max number of iterations for afterburner\r
+ Double_t fPsi2; //! 2nd order event plane orientation\r
+ Double_t fPsi3; //! 3rd order event plane orientation\r
+ Double_t fPrecisionPhi; // afterburner precision\r
+ detectorType fDetectorType; // type of detector from which the EP is taken\r
+ // output histograms\r
+ TH1F* fHistOriginalSpectrum[10]; //! original pt spectrum of accepted tracks\r
+ TH2F* fHistOriginalEtaPhi[10]; //! original eta phi spectrum of accepted tracks\r
+ TH1F* fHistToySpectrum[10]; //! spectrum of toy (generated) tracks\r
+ TH2F* fHistToyEtaPhi[10]; //! eta phi spectrum of toy (generated) tracks\r
+ TH1F* fHistOriginalDeltaPhi[10]; //! original delta phi spectrum\r
+ TH1F* fHistToyDeltaPhi[10]; //! delta phi spectrum of toy (generated) tracks\r
+ TH2F* fHistToyVn[10]; //! generated differential vn values (should equal the differential spectrum)\r
+ TH1F* fHistSFJetSpectrum; //! spectrum of generated sf jets\r
+ TH2F* fHistSFJetEtaPhi; //! eta phi of generated sf jets\r
+ private:\r
+ AliAnalysisTaskJetFlowMC(const AliAnalysisTaskJetFlowMC&); // not implemented\r
+ AliAnalysisTaskJetFlowMC &operator=(const AliAnalysisTaskJetFlowMC&); // not implemented\r
+\r
+ ClassDef(AliAnalysisTaskJetFlowMC, 1); // Task to generate toy mc PicoTracks based on real events\r
+};\r
+#endif\r
#pragma link C++ class AliFlowVZEROQA+;
#pragma link C++ class AliAnalysisTaskFlowEPCascade+;
#pragma link C++ class AliAnalysisTaskJetFlow+;
+#pragma link C++ class AliAnalysisTaskJetFlowMC+;
#endif
--- /dev/null
+///////////////////////////////////////////////////////////////////////////////\r
+// AddTaskJetFlowToyMC //\r
+// Author: Redmer A. Bertens, Utrecht University, 2013, rbertens@cern.ch //\r
+///////////////////////////////////////////////////////////////////////////////\r
+class AliAnalysisDataContainer;\r
+class AliAnalysisTaskJetFlowMC;\r
+\r
+AliAnalysisTaskJetFlowMC* AddTaskJetFlowMC(\r
+ const char *outputTracks = "JetFlowToyMC",\r
+ const char *inputTracks = "PicoTracks",\r
+ const char *name ="AliAnalysisTaskJetFlowMC"\r
+ )\r
+{ \r
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+ if (!mgr) return 0x0;\r
+ if (!mgr->GetInputEventHandler()) return 0x0;\r
+ TString fileName = AliAnalysisManager::GetCommonFileName();\r
+ fileName += ":";\r
+ fileName += name;\r
+ // create the task\r
+ AliAnalysisTaskJetFlowMC *eTask = new AliAnalysisTaskJetFlowMC(name);\r
+ eTask->SetTracksOutName(outputTracks);\r
+ eTask->SetTracksInName(inputTracks);\r
+ // connect input and output\r
+ mgr->AddTask(eTask);\r
+ mgr->ConnectInput (eTask, 0, mgr->GetCommonInputContainer());\r
+ mgr->ConnectOutput (eTask, 1, mgr->CreateContainer(Form("%s_container", fileName.Data()), TList::Class(), AliAnalysisManager::kOutputContainer, fileName.Data()));\r
+ return eTask;\r
+}\r
TDirectoryFile* dirComb(0x0);
TDirectoryFile* dirInt(0x0);
// centralities
- Double_t _c[] = {0, 10, 30, 50, 90};
+ Double_t _c[] = {0, 10, 20, 30, 40, 50, 60, 70, 90};
TArrayD* centralities = new TArrayD((int)(sizeof(_c)/sizeof(_c[0])), _c);
static const int maxCen((int)(sizeof(_c)/sizeof(_c[0])));// max number of centrality bins
// pt array for jet flow analysis
// pt array for hybrid flow analysis
Double_t ptH[] = {0., 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5};
TArrayD* _ptH = new TArrayD(sizeof(ptH)/sizeof(ptH[0]), ptH);
- Double_t jetRadius(0.3);
+ Double_t jetRadius(0.2);
AliAnalysisTaskRhoVnModulation* rho = new AliAnalysisTaskRhoVnModulation();
TH1F* r2V0A(0x0); // container for second order ep resolution
TH1F* r2V0C(0x0);
TH1F* RMSdPtkNoFit(0x0);
TH1F* RMSdPtkComb(0x0);
TH1F* RMSdPtkInt(0x0);
-
//_____________________________________________________________________________
void extractJetFlow() {
// macro to read output of v1.0 of jet flow tasks
if(lNoFit) {
w.mkdir(Form("DeltaPt_HISTO_%s", lNoFit->GetName()));
w.cd(Form("DeltaPt_HISTO_%s", lNoFit->GetName()));
- RMSdPtkNoFit = GetDeltaPtRMS(lNoFit);
+ RMSdPtkNoFit = GetDeltaPtRMS(lNoFit, TString("kNoFit"));
RMSdPtkNoFit->Write();
}
if(lComb) {
w.mkdir(Form("DeltaPt_HISTO_%s", lComb->GetName()));
w.cd(Form("DeltaPt_HISTO_%s", lComb->GetName()));
- RMSdPtkComb = GetDeltaPtRMS(lComb);
+ RMSdPtkComb = GetDeltaPtRMS(lComb, TString("kComb"));
RMSdPtkComb->Write();
}
if(lInt) {
w.mkdir(Form("DeltaPt_HISTO_%s", lInt->GetName()));
w.cd(Form("DeltaPt_HISTO_%s", lInt->GetName()));
- RMSdPtkInt = GetDeltaPtRMS(lInt);
+ RMSdPtkInt = GetDeltaPtRMS(lInt, TString("kInt"));
RMSdPtkInt->Write();
}
// Get the delta pt info from doing a iterative LHS gaus fit
if(lNoFit) {
w.mkdir(Form("DeltaPt_LHSFIT_%s", lNoFit->GetName()));
w.cd(Form("DeltaPt_LHSFIT_%s", lNoFit->GetName()));
- dPtkNoFit = GetDeltaPtSigma(lNoFit);
+ dPtkNoFit = GetDeltaPtSigma(lNoFit, TString("kNoFit"));
dPtkNoFit->Write();
- GetDeltaPtMean(lNoFit)->Write();
+ GetDeltaPtMean(lNoFit, TString("kNoFit"))->Write();
}
if(lComb) {
w.mkdir(Form("DeltaPt_LHSFIT_%s", lComb->GetName()));
w.cd(Form("DeltaPt_LHSFIT_%s", lComb->GetName()));
- dPtkComb = GetDeltaPtSigma(lComb);
+ dPtkComb = GetDeltaPtSigma(lComb, TString("kComb"));
dPtkComb->Write();
- GetDeltaPtMean(lComb)->Write();
+ GetDeltaPtMean(lComb, TString("kComb"))->Write();
}
if(lInt) {
w.mkdir(Form("DeltaPt_LHSFIT_%s", lInt->GetName()));
w.cd(Form("DeltaPt_LHSFIT_%s", lInt->GetName()));
- dPtkInt = GetDeltaPtSigma(lInt);
+ dPtkInt = GetDeltaPtSigma(lInt, TString("kInt"));
dPtkInt->Write();
- GetDeltaPtMean(lInt)->Write();
+ GetDeltaPtMean(lInt, TString("kInt"))->Write();
}
// Get the delta pt predictions
w.mkdir("DeltaPt_PREDICTION");
- GetPredictedDeltaPtSigma(lComb);
+ GetPredictedDeltaPtSigma(lComb, "");
// extract the flow
TList* listNoFit = dirNoFit->GetListOfKeys();
for(Int_t i(0); i < listNoFit->GetEntries(); i++) {
TString string = listNoFit->At(i)->GetName();
if(string.Contains("JetFlow")) {
- for(Int_t j(0); j < 5; j++) {
+ for(Int_t j(0); j < 8; j++) {
if(string.EndsWith(Form("_%i_histograms", j*10))) {
TList* op = (TList*)dirNoFit->Get(string);
GetJetTrackFlow(op, j*10-5);
}
}
}
+ // get the integrated v2 and v3 that were used for the jet
+ // background subtraction
+ if(lComb) GetIntegratedVn(lComb, TString("VZEROC"));
// get the relative improvements
GetRelativeImprovements();
GetRelativeImprovementsFromRMS();
for(Int_t i(0); i < listNoFit->GetEntries(); i++) {
TString string = listNoFit->At(i)->GetName();
if(string.Contains("HybridFlow")) {
- for(Int_t j(0); j < 5; j++) {
+ for(Int_t j(0); j < 8; j++) {
if(string.EndsWith(Form("_%i_histograms", j*10))) {
TList* op = (TList*)dirNoFit->Get(string);
GetHybridTrackFlow(op, j*10-5);
for(Int_t i(0); i < listComb->GetEntries(); i++) {
TString string = listComb->At(i)->GetName();
if(string.Contains("JetFlow")) {
- for(Int_t j(0); j < 5; j++) {
+ for(Int_t j(0); j < 8; j++) {
if(string.EndsWith(Form("_%i_histograms", j*10))) {
TList* op = (TList*)dirComb->Get(string);
GetJetTrackFlow(op, j*10-5);
for(Int_t i(0); i < listInt->GetEntries(); i++) {
TString string = listInt->At(i)->GetName();
if(string.Contains("JetFlow")) {
- for(Int_t j(0); j < 5; j++) {
+ for(Int_t j(0); j < 8; j++) {
if(string.EndsWith(Form("_%i_histograms", j*10))) {
TList* op = (TList*)dirInt->Get(string);
GetJetTrackFlow(op, j*10-5);
}
}
}
+ // save the analysis summary histogram
+ if(lNoFit) GetAnalysisSummary(lNoFit);
+ if(lComb) GetAnalysisSummary(lComb);
+ if(lInt) GetAnalysisSummary(lInt);
// lock and write the output file
w.Close();
+ // get started !
+ TBrowser* browser = new TBrowser();
+ SetStyle();
}
//_____________________________________________________________________________
void LoadLibraries() {
// note that the error propagation towards the relative
// improvement is NOT CORRECT !
if(dPtTheoryVn && dPtTheoryNoVn ) {
- TH1F* impTheory = new TH1F("relative improvement #delta p_{T} #sigma, theory ", "relative improvement #delta p_{T} #sigma, theory", centralities->GetSize()-1, centralities->GetArray());
- impTheory->GetYaxis()->SetTitle("relative improvement [ #frac{#delta p_{T} #sigma no v_{n} - #delta p_{T} #sigma v_{n}}{#delta p_{T} #sigma no v_{n}} ]");
+ TH1F* impTheory = new TH1F("theory, LHS fit ", "theory, LHS fit", centralities->GetSize()-1, centralities->GetArray());
+ impTheory->GetYaxis()->SetTitle("relative improvement");
impTheory->GetXaxis()->SetTitle("centrality percentile");
for(Int_t i (0); i < centralities->GetSize(); i++) {
Double_t a = dPtTheoryNoVn->GetBinContent(i+1);
}
}
if(dPtkNoFit && dPtkComb) {
- TH1F* impComb = new TH1F("relative improvement #delta p_{T} #sigma, kCombined ", "relative improvement #delta p_{T} #sigma, kCombined", centralities->GetSize()-1, centralities->GetArray());
- impComb->GetYaxis()->SetTitle("relative improvement [ #frac{#delta p_{T} #sigma no v_{n} - #delta p_{T} #sigma v_{n}}{#delta p_{T} #sigma no v_{n}} ]");
+ TH1F* impComb = new TH1F("measured, LHS fit", "measured, LHS fit", centralities->GetSize()-1, centralities->GetArray());
+ impComb->GetYaxis()->SetTitle("relative improvement");
impComb->GetXaxis()->SetTitle("centrality percentile");
for(Int_t i (0); i < centralities->GetSize(); i++) {
Double_t a = dPtkComb->GetBinContent(i+1);
}
if(dPtkNoFit && dPtkInt) {
TH1F* impInt = new TH1F("relative improvement #delta p_{T} #sigma, kInt ", "relative improvement #delta p_{T} #sigma, kInt", centralities->GetSize()-1, centralities->GetArray());
- impInt->GetYaxis()->SetTitle("relative improvement [ #frac{#delta p_{T} #sigma no v_{n} - #delta p_{T} #sigma v_{n}}{#delta p_{T} #sigma no v_{n}} ]");
+ impInt->GetYaxis()->SetTitle("relative improvement");
impInt->GetXaxis()->SetTitle("centrality percentile");
for(Int_t i (0); i < centralities->GetSize(); i++) {
Double_t a = dPtkInt->GetBinContent(i+1);
w.cd();
w.mkdir("Relative improvement delta pt distributions");
w.cd("Relative improvement delta pt distributions");
+ FormatMe(impTheory);
impTheory->Write();
+ FormatMe(impComb);
impComb->Write();
+ FormatMe(impInt);
impInt->Write();
}
// note that the error propagation towards the relative
// improvement is NOT CORRECT !
if(dPtTheoryVn && dPtTheoryNoVn ) {
- TH1F* impTheory = new TH1F("relative improvement #delta p_{T} #sigma, theory", "relative improvement #delta p_{T} #sigma, theory ", centralities->GetSize()-1, centralities->GetArray());
- impTheory->GetYaxis()->SetTitle("relative improvement [ #frac{#delta p_{T} #sigma no v_{n} - #delta p_{T} #sigma v_{n}}{#delta p_{T} #sigma no v_{n}} ]");
+ TH1F* impTheory = new TH1F("theory", "theory", centralities->GetSize()-1, centralities->GetArray());
+ impTheory->GetYaxis()->SetTitle("relative improvement");
impTheory->GetXaxis()->SetTitle("centrality percentile");
for(Int_t i (0); i < centralities->GetSize(); i++) {
Double_t a = RMSdPtTheoryNoVn->GetBinContent(i+1);
}
}
if(RMSdPtkNoFit && RMSdPtkComb) {
- TH1F* impComb = new TH1F("relative improvement #delta p_{T} #sigma, kCombined RMS", "relative improvement #delta p_{T} #sigma, kCombined RMS", centralities->GetSize()-1, centralities->GetArray());
- impComb->GetYaxis()->SetTitle("relative improvement [ #frac{#delta p_{T} #sigma no v_{n} - #delta p_{T} #sigma v_{n}}{#delta p_{T} #sigma no v_{n}} ]");
+ TH1F* impComb = new TH1F("measured", "measured", centralities->GetSize()-1, centralities->GetArray());
+ impComb->GetYaxis()->SetTitle("relative improvement");
impComb->GetXaxis()->SetTitle("centrality percentile");
for(Int_t i (0); i < centralities->GetSize(); i++) {
Double_t a = RMSdPtkComb->GetBinContent(i+1);
}
}
if(RMSdPtkNoFit && RMSdPtkInt) {
- TH1F* impInt = new TH1F("relative improvement #delta p_{T} #sigma, kInt RMS", "relative improvement #delta p_{T} #sigma, kInt RMS", centralities->GetSize()-1, centralities->GetArray());
- impInt->GetYaxis()->SetTitle("relative improvement [ #frac{#delta p_{T} #sigma no v_{n} - #delta p_{T} #sigma v_{n}}{#delta p_{T} #sigma no v_{n}} ]");
+ TH1F* impInt = new TH1F("measured ", "measured ", centralities->GetSize()-1, centralities->GetArray());
+ impInt->GetYaxis()->SetTitle("relative improvement");
impInt->GetXaxis()->SetTitle("centrality percentile");
for(Int_t i (0); i < centralities->GetSize(); i++) {
Double_t a = RMSdPtkInt->GetBinContent(i+1);
w.cd();
w.mkdir("Relative improvement delta pt distributions from RMS");
w.cd("Relative improvement delta pt distributions from RMS");
+ FormatMe(impTheory);
impTheory->Write();
+ FormatMe(impComb);
impComb->Write();
+ FormatMe(impInt);
impInt->Write();
}
TProfile* rc = (TProfile*)jf->FindObject("Reference cumulants");
if(qc2 && rc && v0a ) {
TH1F* result = rho->GetDifferentialQC(rc, qc2, _ptH, 2);
+ FormatMe(result);
TString t = "qc2_";
t+=jf->GetName();
result->SetNameTitle(t.Data(), t.Data());
}
if(v0a) {
TH1F* result = rho->CorrectForResolutionDiff((TH1F*)v0a, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, centralities, c, 2);
+ FormatMe(result);
+ result->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ result->GetYaxis()->SetTitle("v_{2}");
result->Write();
}
if(v0c) {
TH1F* result = rho->CorrectForResolutionDiff((TH1F*)v0c, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, centralities, c, 2);
- result->Write();
+ FormatMe(result);
+ result->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ result->GetYaxis()->SetTitle("v_{2}");
+ result->Write();
}
// attempt to get the flow from the qc analysis
TDirectoryFile* qc = (TDirectoryFile*)f.Get("QC");
TH1F* result = rho->GetDifferentialQC(rc, qc2, _ptJ, 2);
TString t = "qc2_";
t+=jf->GetName();
- result->SetNameTitle(t.Data(), t.Data());
+ result->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ result->GetYaxis()->SetTitle("v_{2}");
+ result->SetNameTitle(t.Data(), t.Data());
+ FormatMe(result);
result->Write();
}
if(v0a) {
TH1F* result = rho->CorrectForResolutionDiff((TH1F*)v0a, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, centralities, c, 2);
- result->Write();
+ result->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ result->GetYaxis()->SetTitle("v_{2}");
+ FormatMe(result);
+ result->Write();
}
if(v0c) {
TH1F* result = rho->CorrectForResolutionDiff((TH1F*)v0c, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, centralities, c, 2);
+ result->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ result->GetYaxis()->SetTitle("v_{2}");
+ FormatMe(result);
result->Write();
}
}
//_____________________________________________________________________________
-TH1F* GetDeltaPtRMS(TList* l) {
+TH1F* GetDeltaPtRMS(TList* l, TString suffix) {
// get the RMS value of delta pt
- TH1F* deltaPtRMS = new TH1F("#delta p_{T} RMS", "#delta p_{T} RMS", centralities->GetSize()-1, centralities->GetArray());
+ TH1F* deltaPtRMS = new TH1F(Form("#delta p_{T} RMS, %s", suffix.Data()), Form("#delta p_{T} RMS, %s", suffix.Data()), centralities->GetSize()-1, centralities->GetArray());
deltaPtRMS->GetXaxis()->SetTitle("centrality percentile");
deltaPtRMS->GetYaxis()->SetTitle("RMS [GeV/c]");
for(Int_t i(0); i < maxCen; i++) {
deltaPtRMS->SetBinContent(i+1, dpt->GetRMS(2));
deltaPtRMS->SetBinError(i+1, dpt->GetRMSError(2));
}
+ FormatMe(deltaPtRMS);
return deltaPtRMS;
}
//_____________________________________________________________________________
-TH1F* GetDeltaPtSigma(TList* l) {
+TH1F* GetDeltaPtSigma(TList* l, TString suffix) {
// get the sigma of the delta pt distribution from a recursive LHS gauss fit
- TH1F* deltaPtMean = new TH1F("#delta p_{T} #sigma", "#delta p_{T} #sigma", centralities->GetSize()-1, centralities->GetArray());
+ TH1F* deltaPtMean = new TH1F(Form("#delta p_{T} #sigma, %s", suffix.Data()), Form("#delta p_{T} #sigma, %s",suffix.Data()), centralities->GetSize()-1, centralities->GetArray());
deltaPtMean->GetYaxis()->SetTitle("#sigma [GeV/c]");
deltaPtMean->GetXaxis()->SetTitle("centrality percentile");
for(Int_t i(0); i < maxCen; i++) {
deltaPtMean->SetBinContent(1+i, fit->GetParameter(2));
deltaPtMean->SetBinError(1+i, fit->GetParError(2));
}
+ FormatMe(deltaPtMean);
return deltaPtMean;
}
//_____________________________________________________________________________
-TH1F* GetDeltaPtMean(TList* l) {
+TH1F* GetDeltaPtMean(TList* l, TString suffix) {
// get the mean of the delta pt distribution from a recursive LHS gauss fit
- TH1F* deltaPtMean = new TH1F("#delta p_{T} mean", "#delta p_{T} meam", centralities->GetSize()-1, centralities->GetArray());
+ TH1F* deltaPtMean = new TH1F(Form("#delta p_{T} mean %s", suffix.Data()), Form("#delta p_{T} mean %s", suffix.Data()), centralities->GetSize()-1, centralities->GetArray());
deltaPtMean->GetYaxis()->SetTitle("mean [GeV/c]");
deltaPtMean->GetXaxis()->SetTitle("centrality percentile");
for(Int_t i(0); i < maxCen; i++) {
Double_t s = temp->GetRMS();
Double_t m = temp->GetMean();
TF1* fit = new TF1(Form("mean_%s", temp->GetName()), "gaus", m-3*s, m+0.5*s);
- TH1F* qam = new TH1F(Form("QA_M_mean_%s", temp->GetName()), Form("QA_M_mean_%s", temp->GetName()), 10, 0, 10);
- TH1F* qas = new TH1F(Form("QA_S_mean_%s", temp->GetName()), Form("QA_S_mean_%s", temp->GetName()), 10, 0, 10);
+ TH1F* qam = new TH1F(Form("mu_%s", temp->GetName()), "#mu_{i} / #mu_{i-1}", 10, 0, 10);
+ TH1F* qas = new TH1F(Form("sigma_%s", temp->GetName()), "#sigma_{i} / #sigma_{i-1}", 10, 0, 10);
fit->SetParLimits(2, s/2., s*2.);
for(Int_t j(0); j < 10; j++) {
Double_t _m(m), _s(s);
- temp->Fit(fit, "QILR");
fit->SetRange(m-3*s, m+0.5*s);
+ temp->Fit(fit, "QILR");
m = fit->GetParameter(1);
s = fit->GetParameter(2);
if(!m == 0) qam->SetBinContent(j+1, _m/m);
qas->Write();
qam->Write();
}
+ FormatMe(deltaPtMean);
return deltaPtMean;
}
//_____________________________________________________________________________
-void GetPredictedDeltaPtSigma(TList* l) {
+void GetPredictedDeltaPtSigma(TList* l, TString suffix) {
// get predicted delta pt sigma
- TH1F* deltaPtSigma = new TH1F("predicted #delta p_{T} #sigma ", "predicted #delta p_{T} #sigma", centralities->GetSize()-1, centralities->GetArray());
+ TH1F* deltaPtSigma = new TH1F(Form("predicted #delta p_{T} #sigma %s", suffix.Data()), Form("predicted #delta p_{T} #sigma %s", suffix.Data()), centralities->GetSize()-1, centralities->GetArray());
deltaPtSigma->GetYaxis()->SetTitle("predicted #sigma [GeV/c]");
deltaPtSigma->GetXaxis()->SetTitle("centrality percentile");
TH1F* deltaPtSigmaNoV = new TH1F("predicted #delta p_{T} #sigma no vn", "predicted #delta p_{T} #sigma no vn", centralities->GetSize()-1, centralities->GetArray());
rho->SetOutputList((TList*)l->Clone());
// get the resolution for the desired detector
r2V0A = rho->GetResolutionFromOuptutFile(AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, 2, centralities);
+ r2V0A->SetNameTitle("VZEROA resolution for #Psi_{2}", "VZEROA resolution for #Psi_{2}");
r3V0A = rho->GetResolutionFromOuptutFile(AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, 3, centralities);
+ r3V0A->SetNameTitle("VZEROA resoltuion for #Psi_{3}", "VZEROA resolution for #Psi_{3}");
r2V0C = rho->GetResolutionFromOuptutFile(AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, 2, centralities);
+ r2V0C->SetNameTitle("VZEROC resolution for #Psi_{2}", "VZEROC resolution for #Psi_{2}");
r3V0C = rho->GetResolutionFromOuptutFile(AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, 3, centralities);
+ r3V0C->SetNameTitle("VZEROC resolution for #Psi_{3}", "VZEROC resolution for #Psi_{3}");
// grab the v2 and v3 values and do a resolution correction
TH1F* v2 = new TH1F("v2", "v2", centralities->GetSize()-1, centralities->GetArray());
TH1F* v3 = new TH1F("v3", "v3", centralities->GetSize()-1, centralities->GetArray());
v3->SetBinContent(1+i, pv3->GetBinContent(1+i));
v3->SetBinError(1+i, pv3->GetBinError(1+i));
}
- TH1F* cv2 = new TH1F("v2int","v2int",10,0,100);
- Double_t c_v2[] = {0, 0.036416,0.064765,0.084340,0.096771,0.104257,0.105902,0.104897,0.104811,0.104811,0.104811};
+ // from
+ TH1F* cv2 = new TH1F("v2 in from literaturet","v2 int from literature",10,0,100);
+ Double_t c_v2[] = {0, 0.03565825, 0.06394614, 0.08306863, 0.09470311, 0.09927855, 0.09630484, 0.08708335, 0.07051519, 0, 0};
cv2->SetContent(c_v2);
- TH1F* cv3 = new TH1F("v3int","v3int", 10, 0, 100);
- Double_t c_v3[] = {0, 0.0236149685, 0.02875255, 0.03241459, 0.03507416, 0.03730817, 0.03889757,0.04285879,0.05030896, 0, 0};
+ TH1F* cv3 = new TH1F("v3 int from literature","v3 int from literature", 10, 0, 100);
+ Double_t c_v3[] = {0, 0.02159083, 0.02642751, 0.02928424, 0.03052121, 0.03004316, 0, 0, 0, 0, 0};
cv3->SetContent(c_v3);
for(Int_t i(0); i < centralities->GetSize()-1; i++) {
TH1F* h = (TH1F*)l->FindObject(Form("fHistPicoTrackMult_%i", i));
}
w.mkdir("RhoTaskVnEstimates");
w.cd("RhoTaskVnEstimates");
+ FormatMe(r2V0A);
r2V0A->Write();
+ FormatMe(r3V0A);
r3V0A->Write();
+ FormatMe(r2V0C);
+ r2V0C->Write();
+ FormatMe(r3V0C);
+ r3V0C->Write();
+ FormatMe(cv2);
cv2->Write();
+ FormatMe(cv3);
cv3->Write();
w.cd("DeltaPt_PREDICTION");
dPtTheoryVn = deltaPtSigma;
RMSdPtTheoryVn = deltaPtSigma;
+ FormatMe(dPtTheoryVn);
dPtTheoryVn->Write();
dPtTheoryNoVn = deltaPtSigmaNoV;
RMSdPtTheoryNoVn = deltaPtSigmaNoV;
+ FormatMe(dPtTheoryNoVn);
dPtTheoryNoVn->Write();
}
//_____________________________________________________________________________
+void GetIntegratedVn(TList* l, TString det = "VZEROC") {
+ // get the v2 and v3 values that were used to estimate local energy denstity
+ w.cd("RhoTaskVnEstimates");
+ TH1F* v2 = new TH1F("v2obs", "v2obs", centralities->GetSize()-1, centralities->GetArray());
+ TH1F* v3 = new TH1F("v3obs", "v3obs", centralities->GetSize()-1, centralities->GetArray());
+ TH1F* cv2(0x0);
+ TH1F* cv3(0x0);
+ TProfile* pv2 = (TProfile*)l->FindObject("fProfV2");
+ TProfile* pv3 = (TProfile*)l->FindObject("fProfV3");
+ for(Int_t i(0); i < maxCen; i++) {
+ v2->SetBinContent(1+i, pv2->GetBinContent(1+i));
+ v2->SetBinError(1+i, pv2->GetBinError(1+i));
+ v3->SetBinContent(1+i, pv3->GetBinContent(1+i));
+ v3->SetBinError(1+i, pv3->GetBinError(1+i));
+ }
+ if(det.EqualTo("VZEROA")) {
+ cv2 = rho->CorrectForResolutionInt(v2, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, centralities, 2);
+ cv3 = rho->CorrectForResolutionInt(v3, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROA, centralities, 3);
+ } else if (det.EqualTo("VZEROC")) {
+ cv2 = rho->CorrectForResolutionInt(v2, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, centralities, 2);
+ cv3 = rho->CorrectForResolutionInt(v3, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROC, centralities, 3);
+ } else if (det.EqualTo("VZEROComb")) {
+ cv2 = rho->CorrectForResolutionInt(v2, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROComb, centralities, 2);
+ cv3 = rho->CorrectForResolutionInt(v3, AliAnalysisTaskRhoVnModulation::detectorType::kVZEROComb, centralities, 3);
+ }
+ TString nt2 = Form("v2, %s", det.Data());
+ TString nt3 = Form("v3, %s", det.Data());
+ cv2->SetNameTitle(nt2.Data(), nt2.Data());
+ cv3->SetNameTitle(nt3.Data(), nt3.Data());
+ FormatMe(cv2);
+ cv2->Write();
+ FormatMe(cv3);
+ cv3->Write();
+}
+//_____________________________________________________________________________
+void GetAnalysisSummary(TList* l) {
+ // get and format the analyis summary histogram
+ TH1F* h = (TH1F*)l->FindObject("fHistAnalysisSummary");
+ if(h) {
+ Double_t iter = h->GetBinContent(37);
+ if(iter <= 0) return; // zero events in sample ...
+ Int_t type = TMath::Nint(h->GetBinContent(34)/iter);
+ TString name = "";
+ if(type==0) name+="kNoFit";
+ if(type==1) name+="kV2";
+ if(type==2) name+="kV3";
+ if(type==3) name+="kCombined";
+ if(type==4) name+="kFourierSeries";
+ if(type==5) name+="kIntegratedFlow";
+ if(type==6) name+="kQC2";
+ if(type==7) name+="kQC4";
+ for(Int_t i(0); i < h->GetXaxis()->GetNbins(); i++) h->SetBinContent(i+1, h->GetBinContent(i+1)/iter);
+ w.mkdir(Form("Summary_%s", name.Data()));
+ w.cd(Form("Summary_%s", name.Data()));
+ h->Write();
+ }
+}
+//_____________________________________________________________________________
+void SetStyle() {
+ // set global style
+ gStyle->SetOptStat(0);
+ gPad->SetGrid(1,1);
+ gPad->SetTicks(1,1);
+ gStyle->ToggleEditor();
+}
+//_____________________________________________________________________________
+TH1* FormatMe(TObject* object, Int_t color = -1) {
+ if(color=-1) color = TMath::Nint(gRandom->Uniform(1, 9));
+ TH1* dud = (dynamic_cast<TH1*>object);
+ if (dud) return FormatHistogram(dud, color);
+}
+//_____________________________________________________________________________
+TH1F* FormatHistogram(TH1* hist, Int_t color) {
+ // return a more readable TH1F
+ hist->SetLineWidth(3);
+ hist->SetLineColor(color);
+ hist->SetMarkerStyle(20);
+ hist->SetMarkerColor(color);
+ hist->SetMarkerColor(color);
+ TString name = Form("%s R = %.2f", hist->GetTitle(), jetRadius);
+ hist->SetNameTitle(name.Data(), name.Data());
+
+}
+//_____________________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff