if(fCentralityCut) {
fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
fOutputList->Add(fMultCorAfterCuts);
- fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 9, -0.5, 8.5, 100, 0, 3000);
+ fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 9, -0.5, 100.5, 101, 0, 3000);
fOutputList->Add(fMultvsCentr);
}
}
ClassImp(AliAnalysisTwoParticleResonanceFlowTask)
AliAnalysisTwoParticleResonanceFlowTask::AliAnalysisTwoParticleResonanceFlowTask() : AliAnalysisTaskSE(),
- fSpeciesA(0), fSpeciesB(0), fChargeA(0), fChargeB(0), fMassA(0), fMassB(0), fMinPtA(0), fMaxPtA(0), fMinPtB(0), fMaxPtB(0), fIsMC(0), fEventMixing(0), fPhiMinusPsiMethod(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fPermissiveMixing(0), fNPtBins(18), fNdPhiBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDp(0), fPtP(0), fPtN(0), fPtSpeciesA(0), fPtSpeciesB(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fPOICuts(0), fVertexRange(0), fPhi(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fAnalysisSummary(0)
+ fSpeciesA(0), fSpeciesB(0), fChargeA(0), fChargeB(0), fMassA(0), fMassB(0), fMinPtA(0), fMaxPtA(0), fMinPtB(0), fMaxPtB(0), fIsMC(0), fEventMixing(0), fPhiMinusPsiMethod(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fPermissiveMixing(0), fNPtBins(18), fNdPhiBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDp(0), fPtP(0), fPtN(0), fPtSpeciesA(0), fPtSpeciesB(0), fMultCorAfterCuts(0), fMultvsCentr(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fCentralityCut(0), fPOICuts(0), fVertexRange(0), fPhi(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fAnalysisSummary(0)
{
// Default constructor
for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
}
}
//_____________________________________________________________________________
-AliAnalysisTwoParticleResonanceFlowTask::AliAnalysisTwoParticleResonanceFlowTask(const char *name) : AliAnalysisTaskSE(name), fSpeciesA(0), fSpeciesB(0), fChargeA(0), fChargeB(0), fMassA(0), fMassB(0), fMinPtA(0), fMaxPtA(0), fMinPtB(0), fMaxPtB(0), fIsMC(0), fEventMixing(0), fPhiMinusPsiMethod(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0),fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fPermissiveMixing(0), fNPtBins(18), fNdPhiBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDp(0), fPtP(0), fPtN(0), fPtSpeciesA(0), fPtSpeciesB(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fPOICuts(0), fVertexRange(0), fPhi(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fAnalysisSummary(0)
+AliAnalysisTwoParticleResonanceFlowTask::AliAnalysisTwoParticleResonanceFlowTask(const char *name) : AliAnalysisTaskSE(name), fSpeciesA(0), fSpeciesB(0), fChargeA(0), fChargeB(0), fMassA(0), fMassB(0), fMinPtA(0), fMaxPtA(0), fMinPtB(0), fMaxPtB(0), fIsMC(0), fEventMixing(0), fPhiMinusPsiMethod(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0),fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fPermissiveMixing(0), fNPtBins(18), fNdPhiBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDp(0), fPtP(0), fPtN(0), fPtSpeciesA(0), fPtSpeciesB(0), fMultCorAfterCuts(0), fMultvsCentr(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fCentralityCut(0), fPOICuts(0), fVertexRange(0), fPhi(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fAnalysisSummary(0)
{
// Constructor
for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
fOutputList->Add(fDCAXYQA);
fDCAZQA = new TH1F("fDCAZQA", "fDCAZQA", 1000, -5, 5);
fOutputList->Add(fDCAZQA);
+ if(fCentralityCut) {
+ fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
+ fOutputList->Add(fMultCorAfterCuts);
+ fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 9, -0.5, 100.5, 101, 0, 3000);
+ fOutputList->Add(fMultvsCentr);
+ }
}
if(fIsMC || fQA) {
fDCAAll = new TH2F("fDCAAll", "fDCAAll", 1000, 0, 10, 1000, -5, 5);
}
Double_t centralityBins[_c];
Double_t vertexBins[_v];
- for(Int_t i(0); i < _c + 1; i++) centralityBins[i] = (Double_t)fCentralityMixingBins[i];
- for(Int_t i(0); i < _v + 1; i++) vertexBins[i] = (Double_t)fVertexMixingBins[i];
- return new AliEventPoolManager(fMixingParameters[0], fMixingParameters[1], -1+(Int_t)(sizeof(centralityBins)/sizeof(centralityBins[0])), (Double_t*)centralityBins, -1+(Int_t)(sizeof(vertexBins)/sizeof(vertexBins[0])), (Double_t*)vertexBins);
+ for(Int_t i(0); i < _c + 1; i++) centralityBins[i] = fCentralityMixingBins[i];
+ for(Int_t i(0); i < _v + 1; i++) vertexBins[i] = fVertexMixingBins[i];
+ return new AliEventPoolManager(fMixingParameters[0], fMixingParameters[1], _c, (Double_t*)centralityBins, _v, (Double_t*)vertexBins);
}
//_____________________________________________________________________________
template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::InvariantMass(const T* track1, const T* track2) const
if(fQA) fCentralityNoPass->Fill(fCentrality) ;
return kFALSE;
}
+ const Int_t nGoodTracks = event->GetNumberOfTracks();
+ if(fCentralityCut) { // cut on outliers
+ Float_t multTPC(0.); // tpc mult estimate
+ Float_t multGlob(0.); // global multiplicity
+ for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
+ AliAODTrack* trackAOD = event->GetTrack(iTracks);
+ if (!trackAOD) continue;
+ if (!(trackAOD->TestFilterBit(1))) continue;
+ if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
+ multTPC++;
+ }
+ for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
+ AliAODTrack* trackAOD = event->GetTrack(iTracks);
+ if (!trackAOD) continue;
+ if (!(trackAOD->TestFilterBit(16))) continue;
+ if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
+ Double_t b[2] = {-99., -99.};
+ Double_t bCov[3] = {-99., -99., -99.};
+ if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
+ if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
+ multGlob++;
+ } //track loop
+ // printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
+ if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))) return kFALSE;
+ if(fQA) {
+ fMultCorAfterCuts->Fill(multGlob, multTPC);
+ fMultvsCentr->Fill(fCentrality, multTPC);
+ }
+ }
fCentralityPass->Fill(fCentrality);
return kTRUE;
}
// setters
void SetPtBins(Float_t bin[19], Int_t n) { for(Int_t i = 0; i < n+1; i++) fPtBins[i] = bin[i]; fNPtBins = n; }
void SetdPhiBins(Float_t bin[19], Int_t n) { for(Int_t i = 0; i < n+1; i++) fdPhiBins[i] = bin[i]; fNdPhiBins = n;}
- void SetCentralityParameters(Double_t min, Double_t max, const char* a, const char* b) {
+ void SetCentralityParameters(Double_t min, Double_t max, const char* a, const char* b, Bool_t c) {
fCentralityMin = min;
fCentralityMax = max;
fkCentralityMethodA = a;
- fkCentralityMethodB = b; }
+ fkCentralityMethodB = b;
+ fCentralityCut = c; }
void SetPOICuts(AliFlowTrackCuts *cutsPOI) { fPOICuts = cutsPOI; }
void SetRPCuts(AliFlowTrackCuts *cutsRP) { fCutsRP = cutsRP; }
void SetPIDConfiguration(Float_t prob[7]) { for(Int_t i = 0; i < 7; i++) fPIDConfig[i] = prob[i]; }
Bool_t SetQA(Bool_t qa) {fQA = qa; return fQA;}
void SetAddTaskMacroSummary(Float_t m[12]) {for(Int_t i(0); i < 12; i++) fAddTaskMacroSummary[i] = m[i];}
void SetPOIDCAXYZ(Float_t dca[5]) { for(Int_t i = 0; i < 5; i++) fDCAConfig[i] = dca[i]; }
+ void SetMixingBins(Int_t c[20], Int_t v[20]) {for(Int_t i = 0; i < 20; i++) { fCentralityMixingBins[i] = c[i];
+ fVertexMixingBins[i] = v[i]; } }
void SetMixingParameters(Int_t p[3]) { for(Int_t i = 0; i < 3; i++) fMixingParameters[i] = p[i]; }
void SetPermissiveMixing(Bool_t p) { fPermissiveMixing = p; }
void SetupSpeciesA(Int_t species, Int_t charge, Float_t mass, Float_t minPtA, Float_t maxPtA) {fSpeciesA = species; fChargeA = charge; fMassA = mass; fMinPtA = minPtA; fMaxPtA = maxPtA;}
void SetMaxDeltaDipAngleAndPt(Float_t a, Float_t pt) { fDeltaDipAngle = a;
fDeltaDipPt = pt;
fApplyDeltaDipCut = kTRUE; };
- void SetMixingBins(Int_t c[20], Int_t v[20]) {for(Int_t i = 0; i < 20; i++) { fCentralityMixingBins[i] = c[i];
- fVertexMixingBins[i] = v[i]; } }
//getters
void GetMixingParameters(Int_t p[3]) const { for(Int_t i = 0; i < 3; i++) p[i] = fMixingParameters[i]; }
Float_t GetCenMin() const {return fCentralityMin; }
TH1F *fPtN; //! QA histogram of p_t distribution of negative particles
TH1F *fPtSpeciesA; //! QA histogram of p_t distribution of species A
TH1F *fPtSpeciesB; //! QA histogram of p_t distribution of species B
+ TH2F *fMultCorAfterCuts; //! QA profile global and tpc multiplicity after outlier cut
+ TH2F *fMultvsCentr; //! QA profile of centralty vs multiplicity
Float_t fCentralityMin; // lower bound of cenrality bin
Float_t fCentralityMax; // upper bound of centrality bin
const char *fkCentralityMethodA; // centrality determiantion (primary method)
const char *fkCentralityMethodB; // centrality determination fallback
+ Bool_t fCentralityCut; // 3 sigma cut for multiplicity outliers
AliFlowTrackCuts *fPOICuts; // cuts for particles of interest (flow package)
Float_t fVertexRange; // absolute value of maximum distance of vertex along the z-axis
TH1F *fPhi; //! QA plot of azimuthal distribution of POI daughters
}
if(task->SetVZEROSubEvents(EP3sub)) cout << " --> Setting up VZERO subevents method ... " << endl;
if(event_mixing) {
- if(debug) cout << " --> Enabeling event mixing for phi reconstruction - try at your own risk !!! <-- " << endl;
+ if(debug) cout << " --> Enabeling event mixing for reconstruction - try at your own risk !!! <-- " << endl;
// set vertex and mixing bins - arrays MUST have length 20!
Int_t c[] = {0, 2, 4, 6, 8, 10, 20, 30, 40, 50, 60, 70, 80, 90, 101, 0, 0, 0, 0, 0};
Int_t v[] = {-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0};
return 0x0;
}
else task->SetMixingBins(c, v);
- }
+ }
// set triggers
if (bCentralTrigger) task->SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral);
else task->SelectCollisionCandidates(AliVEvent::kMB);
if((deltaDip>0.005)&&(deltaDipMaxPt>0.005)) task->SetMaxDeltaDipAngleAndPt(deltaDip, deltaDipMaxPt);
else cout << " --> Disabled Delta-Dip exclusion. <-- " << endl;
task->SetCandidateEtaAndPt(POIEtaMin, POIEtaMax, 0., 15.);
- task->SetCentralityParameters(centrMin, centrMax, "TRK", "V0M");
+ task->SetCentralityParameters(centrMin, centrMax, "TRK", "V0M", kTRUE);
task->SetVertexZ(vertexZ);
if(debug) cout << " --> Set pair cuts and event cuts" << endl;
// specify the PID procedure which will be used
if((deltaDip>0.005)&&(deltaDipMaxPt>0.005)) task->SetMaxDeltaDipAngleAndPt(deltaDip, deltaDipMaxPt);
else cout << " --> Disabled Delta-Dip exclusion. <-- " << endl;
task->SetCandidateEtaAndPt(POIEtaMin, POIEtaMax, 0., 10.);
- task->SetCentralityParameters(centrMin, centrMax, "TRK", "V0M");
+ task->SetCentralityParameters(centrMin, centrMax, "TRK", "V0M", kTRUE);
task->SetVertexZ(vertexZ);
if(debug) cout << " --> Set pair cuts and event cuts" << endl;
// set the kaon cuts, and specify the PID procedure which will be used
* if the macro is called with kFALSE as argument,
* an example function of how to do analysis on the output is called
*
+ * at the end of the macro, an example function is given that can be used
+ * to connect the on the fly events to the AliAnalysistwoParticleResonanceFlowTask
+ *
* Redmer Alexander Bertens (rbertens@cern.ch)
* Carlos Eugenio Perez Lara
* Andrea Dubla
*/
class AliFlowOnTheFlyGenerator;
-
+class AliAnalysisTwoParticleResonanceFlowTask;
+class AliFlowAnalysisWithQCumulants;
// main function
//
// you can call this macro in different modes.
// doing the desired v2 analysis)
void GenerateEventsOnTheFly(Bool_t generate = kTRUE) {
if(!generate) {
+ LoadLibraries();
AnalyseEventsFromFile();
return;
}
// note that by default gamma's are 'decayed' to e+ e- pairs
Int_t decayModes[] = { TPythia6Decayer::kHadronicD }; // some decay modes
// h) number of events you want to generate.
- const int events = 1000;
+ const int events = 100;
// i) write output to file
Bool_t writeOutput = kTRUE;
// j) do qa. qa histograms will be created at the end of generation adn written to a rootfile in your working directory
// k) write analysis to an output file
TString trunkName = "OnTheFlyEvent"; // trunk of output file name
// l) specify the max number of events per file (new files will be generated automatically)
- const int maxEventsPerFile = 500; // specify the maximum number of events that is stored per file
+ const int maxEventsPerFile = 50; // specify the maximum number of events that is stored per file
// note that events are stored temporarily in RAM memory before writing to file
// so setting this number to a very large value could lead to an unresponsive system
///====== END OF GENERATOR INTERFACE===///
for(Int_t i(0); i < nDecayModes; i++) decayer->SetForceDecay(decayModes[i]);
// get an instance of the class that we'll use to generate events
AliFlowOnTheFlyEventGenerator* eventGenerator = new AliFlowOnTheFlyEventGenerator( qa, // make some QA plots
+ NULL, // introduce flow fluctuations (not implemented yet)
nMult, // set initial value for evnet multiplcity
decayer, // specify the decayer (NULL is no decay)
addMotherV2, // add v2 for mothers
//_____________________________________________________________________________
void AnalyseEventsFromFile() {
// example function: read events from file and do some analysis.
- gSystem->Load("libPWGflowBase");
Int_t nFiles(0); // file counter
// setup analysis methods
- AliFlowAnalysisWithQCumulants *qc(0x0); // for cumulants
- AliFlowAnalysisWithScalarProduct *sp(0x0); // for scalar product
- sp = new AliFlowAnalysisWithScalarProduct();
- sp->SetHarmonic(2);
- sp->SetApplyCorrectionForNUA(kFALSE);
- sp->Init();
- qc = new AliFlowAnalysisWithQCumulants();
- qc->SetHarmonic(2);
- qc->SetCalculateDiffFlow(kTRUE);
- qc->SetCalculate2DDiffFlow(kFALSE); // vs (pt,eta)
- qc->SetApplyCorrectionForNUA(kFALSE);
- qc->SetFillMultipleControlHistograms(kFALSE);
- qc->SetMultiplicityWeight("combinations"); // default (other supported options are "unit" and "multiplicity")
- qc->SetCalculateCumulantsVsM(kFALSE);
- qc->SetCalculateAllCorrelationsVsM(kFALSE); // calculate all correlations in mixed harmonics "vs M"
- qc->SetnBinsMult(10000);
- qc->SetMinMult(0);
- qc->SetMaxMult(10000);
- qc->SetBookOnlyBasicCCH(kFALSE); // book only basic common control histograms
- qc->SetCalculateDiffFlowVsEta(kTRUE); // if you set kFALSE only differential flow vs pt is calculated
- qc->SetCalculateMixedHarmonics(kFALSE); // calculate all multi-partice mixed-harmonics correlators
- qc->Init();
+ AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
+ PrepareQCumulants(qc);
while(kTRUE) {
TFile file(Form("OnTheFlyEvent_%i.root",nFiles)); // open the root file
}
// insert the poor guy into the flow anaysis method
qc->Make(event);
- sp->Make(event);
delete event;
}
nFiles++;
// prepare output for the flow package analyses
TString outputFileName = "AnalysisResults.root";
TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
- const Int_t nMethods(2);
- TString method[] = {"SP", "QC"};
+ const Int_t nMethods(1);
+ TString method[] = {"QC"};
TDirectoryFile *dirFileFinal[nMethods] = {NULL};
TString fileName[nMethods];
for(Int_t i(0); i < nMethods; i++) {
fileName[i]+="analysis";
dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
}
- sp->Finish();sp->WriteHistograms(dirFileFinal[0]);
- qc->Finish();qc->WriteHistograms(dirFileFinal[1]);
+ qc->Finish();qc->WriteHistograms(dirFileFinal[0]);
outputFile->Close();
delete outputFile;
}
+
+/* ====================================================
+ *
+ * the following two functions serve as an exmaple
+ * of how to run the phi and kstar analysis on on the
+ * fly events
+ *
+ * ====================================================
+ */
+
+//_____________________________________________________________________________
+void TwoParticleResonanceFlowOnTheFly() {
+
+ // load additional libraries
+ gSystem->Load("libGeom");
+ gSystem->Load("libVMC");
+ gSystem->Load("libXMLIO");
+ gSystem->Load("libPhysics");
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+ gSystem->Load("libPWGflowBase");
+ gSystem->Load("libPWGflowTasks");
+
+ TString dirName = "reconstruction";
+ Int_t nFiles(0);
+ // define the poi cuts, see AddTwoParticelResonanceFlowTask.C for more info
+ const Int_t mb(30); // no of massbands available for analysis
+ Float_t minMass(.99), maxMass(1.092); // upper and lower bound
+ AliFlowTrackSimpleCuts* POIfilterQC[mb]; // pointers to poi filters
+ Double_t flowBands[2][mb]; // define the invm regins
+ Double_t _inc = (maxMass-minMass)/(float)mb;
+ for (Int_t _mb = 0; _mb < mb; _mb++) { // create the poi filters
+ flowBands[0][_mb] = minMass + _mb * _inc;
+ flowBands[1][_mb] = minMass + (_mb + 1) * _inc;
+ POIfilterQC[_mb] = new AliFlowTrackSimpleCuts(Form("FilterPOIQC_MB%d", _mb));
+ POIfilterQC[_mb]->SetEtaMin(-0.8); // eta range
+ POIfilterQC[_mb]->SetEtaMax(+0.8);
+ POIfilterQC[_mb]->SetMassMin(flowBands[0][_mb]); // invm range lower bound
+ POIfilterQC[_mb]->SetMassMax(flowBands[1][_mb]); // invm rnage upper bound
+ }
+ // do the flow analysis
+ AliFlowAnalysisWithQCumulants* qc[mb];
+ for(int i(0); i < mb; i++) { // init the q-cumulant tasks, one for each invm bin
+ qc[i] = new AliFlowAnalysisWithQCumulants();
+ PrepareQCumulants(qc[i]);
+ printf(" > init qc task %i < \n", i);
+ }
+ AliAnalysisTwoParticleResonanceFlowTask* task = new AliAnalysisTwoParticleResonanceFlowTask("onthefly");
+ SetUpTask(task, minMass,maxMass);// setup the task
+ // open files
+ while(kTRUE) { // infinite loop which we will break when we've looked through all the files
+ TFile file(Form("OnTheFlyEvent_%i.root", nFiles)); // open the root file
+ if(nFiles==0&&file.IsZombie()) { // something went wrong ...
+ printf(" > cannot get input files ... exiting < \n");
+ exit(0);
+ }
+ if(file.IsZombie()) break; // break the loop on the first empty file
+ TIter iter(file.GetListOfKeys()); // get a list of keys
+ while(kTRUE) { // infinite loop over events ...
+ TKey* key = iter(); // get the next key from the file
+ if(!key) break; // ... and exit the loop if the key is empty
+ AliFlowEventSimple* event = (AliFlowEventSimple*)key->ReadObj(); // cast to a flow event
+ printf(" - read event with %i tracks from file %i \n > task ", event->NumberOfTracks(), nFiles);
+ task->DoAnalysisOnTheFly(event); // do the on the fly analysis
+ AliFlowEventSimple* flowEvent = task->GetFlowEvent(); // retrieve the flow event
+ for(int j(0); j < mb; j++) { // loop over all invm bands
+ flowEvent->TagPOI(POIfilterQC[j]); // 'offline' tagging of poi's in certain mass range
+ flowEvent->TagSubeventsInEta(-.8, 0, 0, .8); // setup subevents
+ qc[j]->Make(flowEvent); // do qc analysis
+ printf(" %i", j);
+ }
+ }
+ nFiles++;
+ }
+ // prepare output for the flow package analyses
+ TFile *outputFile = new TFile("AnalysisResults.root","RECREATE"); // common outputfile
+ TDirectoryFile* dirFileFinal[mb]; // tdirectory files
+ TString fileName[mb]; // dir names
+ for(Int_t i(0); i < mb; i++) { // loop over all bands ...
+ fileName[i]+=Form("QC_minv_%i", i);
+ dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
+ qc[i]->Finish(); // finalize the method
+ qc[i]->WriteHistograms(dirFileFinal[i]); // and write it to file
+ }
+ // write the output of the phi reconstruction to the same file
+ TDirectoryFile* dir = new TDirectoryFile(dirName.Data(), dirName.Data());
+ task->DoAnalysisOnTheFly(dir);
+ // end of analysis
+ outputFile->Close();
+ delete outputFile;
+}
+//_____________________________________________________________________________
+void SetUpTask(AliAnalysisTwoParticleResonanceFlowTask* task, Float_t minMass, Float_t maxMass)
+{
+ // some magic which is necessary to 'trick' the analysis task into thinking
+ // that we're actually doing an analysis on aod events.
+ // some of these configurations are used (e.g. setupSpecies, common constants
+ // and the binning in pt)
+ // others have no meaning (PID, DCA, RP and POI cuts)
+ // note that UserCreateOutputObjects is necessary since it initializes
+ // the output of the analysis
+ gSystem->Load("libPWGflowTasks");
+ Float_t PIDconfig[] = {0,0,0,0,0,0,0.3}; // not used
+ task->SetPIDConfiguration(PIDconfig); // not used
+ task->SetupSpeciesA(321, 1, 4.93676999e-01, 0.15, 5.); // pid and charge
+ task->SetupSpeciesB(-321, -1, 4.93676999e-01, 0.15, 5.); // pid and charge
+ task->SetRPCuts(new AliFlowTrackCuts("unused_rp_cuts")); // not used
+ task->SetPOICuts(new AliFlowTrackCuts("unused_poi_cuts")); // not used
+ Float_t POIDCA[] = {0,0,0,0,0}; // not used
+ task->SetPOIDCAXYZ(POIDCA); // not used
+ task->SetCommonConstants(30, minMass, maxMass); // necesssary
+ Float_t _pt[] = {0., 0.6, 1.2, 1.8, 2.4, 3., 4., 5., 6., 7.}; // same
+ task->SetPtBins(_pt, (Int_t)(sizeof(_pt)/sizeof(_pt[1]))-1); // same
+ task->UserCreateOutputObjects(); // necessary
+}
+//_____________________________________________________________________________
+void PrepareQCumulants(AliFlowAnalysisWithQCumulants* qc)
+{
+ // prepare the cumulants task
+ qc->SetHarmonic(2);
+ qc->SetCalculateDiffFlow(kTRUE);
+ qc->SetCalculate2DDiffFlow(kFALSE); // vs (pt,eta)
+ qc->SetApplyCorrectionForNUA(kFALSE);
+ qc->SetFillMultipleControlHistograms(kFALSE);
+ qc->SetMultiplicityWeight("combinations"); // default (other supported options are "unit" and "multiplicity")
+ qc->SetCalculateCumulantsVsM(kFALSE);
+ qc->SetCalculateAllCorrelationsVsM(kFALSE); // calculate all correlations in mixed harmonics "vs M"
+ qc->SetnBinsMult(10000);
+ qc->SetMinMult(0);
+ qc->SetMaxMult(10000);
+ qc->SetBookOnlyBasicCCH(kFALSE); // book only basic common control histograms
+ qc->SetCalculateDiffFlowVsEta(kTRUE); // if you set kFALSE only differential flow vs pt is calculated
+ qc->Init();
+}
//_____________________________________________________________________________
Bool_t LoadLibraries() {
// load libraries
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff