#include "AliCentrality.h"
#include "AliESDEvent.h"
#include "AliVTrack.h"
-#include "AliESDtrackCuts.h"
+#include "AliESDtrack.h"
#include "AliAODTrack.h"
+#include "AliAnalysisFilter.h"
ClassImp(AliForwardFlowTaskQC)
#if 0
fSumList(0), // Event sum list
fOutputList(0), // Result output list
fAOD(0), // AOD input event
- fESDTrackCuts(0), // ESD track cuts
+ fTrackCuts(0), // ESD track cuts
fMaxMoment(0), // Max flow moment
fVtx(1111), // Z vertex coordinate
fCent(-1), // Centrality
fSumList(0), // Event sum list
fOutputList(0), // Result output list
fAOD(0), // AOD input event
- fESDTrackCuts(0), // ESD track cuts
+ fTrackCuts(0), // ESD track cuts
fMaxMoment(4), // Max flow moment
fVtx(1111), // Z vertex coordinate
fCent(-1), // Centrality
fSumList(o.fSumList), // Event sum list
fOutputList(o.fOutputList), // Result output list
fAOD(o.fAOD), // AOD input event
- fESDTrackCuts(o.fESDTrackCuts), // ESD track cuts
+ fTrackCuts(o.fTrackCuts), // ESD track cuts
fMaxMoment(o.fMaxMoment), // Flow moments
fVtx(o.fVtx), // Z vertex coordinate
fCent(o.fCent), // Centrality
fSumList = o.fSumList;
fOutputList = o.fOutputList;
fAOD = o.fAOD;
- fESDTrackCuts = o.fESDTrackCuts;
+ fTrackCuts = o.fTrackCuts;
fMaxMoment = o.fMaxMoment;
fVtx = o.fVtx;
fCent = o.fCent;
//
InitVertexBins();
InitHists();
- if (fFlowFlags & kTPC) {
- fESDTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
- fESDTrackCuts->SetPtRange(0.2, 5.0);
- fESDTrackCuts->SetEtaRange(-0.8, 0.8);
- fESDTrackCuts->SetMinNClustersTPC(70);
- }
+ if ((fFlowFlags & kTracks) && !fTrackCuts) AliFatal("No track cuts set!");
PrintFlowSetup();
PostData(1, fSumList);
}
else if ((fFlowFlags & kVZERO)) {
fBinsForward.Add(new VertexBin(vL, vH, fMaxMoment, "VZERO", fFlowFlags, 0, fEtaGap));
- if ((fFlowFlags & kEtaGap) && !(fFlowFlags & kTPC))
+ if ((fFlowFlags & kEtaGap) && !(fFlowFlags & kTracks))
fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-VZERO", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
}
}
Int_t nVtxBins = fVtxAxis->GetNbins();
while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
+ i++;
// If no tracks do things normally
- if (!(fFlowFlags & kTPC) && !bin->FillHists(h, fCent, kFillBoth|flags|kReset)) return;
+ if (!(fFlowFlags & kTracks) || (flags & kMC)) {
+ if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) continue;
+ }
// if tracks things are more complicated
- else if ((fFlowFlags & kTPC)) {
- TObjArray* trList = GetTracks();
- if (!trList) return;
- Bool_t useEvent = bin->FillTracks(trList, kFillRef|kReset|flags);
- // If esd input trList is a new object owned by this task and should be cleaned up
- if (AliForwardUtil::CheckForAOD() == 2) delete trList;
- if (!useEvent) return;
- if (!bin->FillHists(h, fCent, kFillDiff|kReset|flags)) return;
+ else if ((fFlowFlags & kTracks)) {
+ if (!FillTracks(bin, kFillRef|kReset|flags)) continue;
+ if (!bin->FillHists(h, fCent, kFillDiff|kReset|flags)) continue;
}
bin->CumulantsAccumulate(fCent);
- i++;
}
return;
Int_t nVtxBins = fVtxAxis->GetNbins();
while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
- if (!(fFlowFlags & kTPC) && !bin->FillHists(href, fCent, kFillRef|flags|kReset)) return;
- else if ((fFlowFlags & kTPC)) {
- TObjArray* trList = GetTracks();
- if (!trList) return;
- Bool_t useEvent = bin->FillTracks(trList, kFillRef|kReset|flags);
- // If esd input trList is a new object owned by this task and should be cleaned up
- if (AliForwardUtil::CheckForAOD() == 2) delete trList;
- if (!useEvent) return;
+ i++;
+ if (!(fFlowFlags & kTracks) || (flags & kMC)) {
+ if(!bin->FillHists(href, fCent, kFillRef|flags|kReset)) continue;
}
- bin->FillHists(hdiff, fCent, kFillDiff|kReset);
+ else if ((fFlowFlags & kTracks)) {
+ if (!FillTracks(bin, kFillRef|kReset|flags)) continue;
+ }
+ if (!bin->FillHists(hdiff, fCent, kFillDiff|kReset|flags)) continue;
bin->CumulantsAccumulate(fCent);
- i++;
}
return;
TH2D& h = CombineHists(hcent, hfwd);
while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
- if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) return;
- bin->CumulantsAccumulate3Cor(fCent);
i++;
+ if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) continue;
+ bin->CumulantsAccumulate3Cor(fCent);
}
return;
return fHistdNdedp3Cor;
}
//_____________________________________________________________________
-TObjArray* AliForwardFlowTaskQC::GetTracks() const
+Bool_t AliForwardFlowTaskQC::FillTracks(VertexBin* bin, UShort_t mode) const
{
//
// Get TPC tracks to use for reference flow.
// Return: TObjArray with tracks
//
TObjArray* trList = 0;
- // Get input type
- UShort_t input = AliForwardUtil::CheckForAOD();
- switch (input) {
- // If AOD input, simply get the track array from the event
- case 1: trList = static_cast<TObjArray*>(fAOD->GetTracks());
- break;
- case 2: {
- // If ESD input get event, apply track cuts
- AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
- if (!esd) return 0;
- // Warning! trList is now a new array, we need to delete it after use
- // this is not a very good implementation!
- trList = fESDTrackCuts->GetAcceptedTracks(esd, kTRUE);
- break;
- }
- default: AliFatal("Neither ESD or AOD input. This should never happen");
- break;
- }
- return trList;
+ AliESDEvent* esdEv = 0;
+ if (AliForwardUtil::CheckForAOD() == 1) // AOD tracks
+ trList = static_cast<TObjArray*>(fAOD->GetTracks());
+ else
+ esdEv = dynamic_cast<AliESDEvent*>(InputEvent());
+
+ Bool_t useEvent = bin->FillTracks(trList, esdEv, fTrackCuts, mode);
+ return useEvent;
}
//_____________________________________________________________________
void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
fOutputList->SetName("Results");
fOutputList->SetOwner();
- if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTPC)) {
+ if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks)) {
TParameter<Double_t>* etaGap = new TParameter<Double_t>("EtaGap", fEtaGap);
fOutputList->Add(etaGap);
}
0.521922, 0.785915, 0.703658, 0.832479, 0.77461, 0.73129, 0.778697, 0.710265,
0.89686, 0.967688, 0.974225, 0.873445, 0.811096, 0.828493, 0.889609, 0.586056,
1.15877, 0.954656, 0.914557, 0.979028, 1.04907, 0.748518, 0.928043, 0.98175 };
+
for (Int_t i = 0; i < 64; i++) {
if (i % 8 == 0) {
ring++;
- bin = (ring < 5 ? ring+1 : 15-ring);
+ bin = (ring < 5 ? 11-ring : ring-3);
eta = fHistdNdedpV0.GetXaxis()->GetBinCenter(bin);
fHistdNdedpV0.SetBinContent(bin, 0, 1);
}
// Get bin numbers and binning defined
Int_t nHBins = GetBinNumberSin();
- Int_t nEtaBins = 48;
+ Int_t nEtaBins = 24;
if ((fFlags & k3Cor)) {
if ((fFlags & kFMD)) nEtaBins = 24;
else if ((fFlags & kVZERO)) nEtaBins = 19;
Int_t nRefBins = nEtaBins; // needs to be something as default
if ((fFlags & kStdQC)) {
- if ((fFlags & kSymEta) && !((fFlags & kTPC) && (fFlags & kSPD))) nRefBins = 1;
+ if ((fFlags & kSymEta) && !((fFlags & kTracks) && (fFlags & kSPD))) nRefBins = 1;
else nRefBins = 2;
} else if ((fFlags & kEtaGap )) {
nRefBins = 2;
if ((mode & kFillRef)) fCumuRef->Reset();
if ((mode & kFillDiff)) fCumuDiff->Reset();
}
-
// Then we loop over the input and calculate sum cos(k*n*phi)
// and fill it in the reference and differential histograms
Int_t nBadBins = 0;
TMath::Abs(eta) < fEtaGap) break;
// Backward and forward eta gap break for reference flow
if ((fFlags & kEtaGap) && (mode & kFillRef) && TMath::Abs(eta) > TMath::Abs(limit)) break;
- if ((fFlags & kStdQC) && (fFlags & kMC)) {
+ if ((fFlags & kStdQC) && (fFlags & kMC) && !(fFlags & kTracks)) {
if (!(fFlags & kSPD) && TMath::Abs(eta) < 1.75) break;
if ((fFlags & kSPD) && TMath::Abs(eta) > 2.00) break;
}
nInAvg++;
if (weight == 0) continue;
if (weight > max) max = weight;
-
// Fill into Cos() and Sin() hists
- if ((mode & kFillRef)) {
+ if ((mode & kFillRef) && !((fFlags & kTracks) && (fFlags & kMC) && TMath::Abs(eta) > 0.75)) {
fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
fdNdedpRefAcc->Fill(eta, phi, weight);
}
Double_t cosnPhi = weight*TMath::Cos(n*phi);
Double_t sinnPhi = weight*TMath::Sin(n*phi);
// fill ref
- if ((mode & kFillRef)) {
+ if ((mode & kFillRef) && !((fFlags & kTracks) && (fFlags & kMC) && TMath::Abs(eta) > 0.75)) {
fCumuRef->Fill(eta, cosBin, cosnPhi);
fCumuRef->Fill(eta, sinBin, sinnPhi);
}
return useEvent;
}
//_____________________________________________________________________
-Bool_t AliForwardFlowTaskQC::VertexBin::FillTracks(TObjArray* trList, UShort_t mode)
+Bool_t AliForwardFlowTaskQC::VertexBin::FillTracks(TObjArray* trList, AliESDEvent* esd,
+ AliAnalysisFilter* trFilter, UShort_t mode)
{
//
// Fill reference and differential eta-histograms
// mode: filling mode: kFillRef/kFillDiff/kFillBoth
//
if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
+ if (!trList && !esd) {
+ AliError("FillTracks: No AOD track list or ESD event - something might be wrong!");
+ return kFALSE;
+ }
// Fist we reset histograms
if ((mode & kReset)) {
if ((mode & kFillDiff)) fCumuDiff->Reset();
}
- UShort_t input = AliForwardUtil::CheckForAOD();
// Then we loop over the input and calculate sum cos(k*n*phi)
// and fill it in the reference and differential histograms
- Int_t nTr = trList->GetEntries();
+ Int_t nTr = 0;
+ if (trList) nTr = trList->GetEntries();
+ if (esd) nTr = esd->GetNumberOfTracks();
if (nTr == 0) return kFALSE;
AliVTrack* tr = 0;
- AliAODTrack* aodTr = 0;
// Cuts for AOD tracks (have already been applied to ESD tracks) - except dEdx
- const Double_t pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, tpcdEdx = 10.;
+// const tpcdEdx = 10;
for (Int_t i = 0; i < nTr; i++) { // track loop
- tr = (AliVTrack*)trList->At(i);
+ tr = (trList ? (AliVTrack*)trList->At(i) : (AliVTrack*)esd->GetTrack(i));
if (!tr) continue;
- if (input == 1) { // If AOD input
- // A dynamic cast would be more safe here, but this is faster...
- aodTr = (AliAODTrack*)tr;
+ if (esd) {
+ AliESDtrack* esdTr = (AliESDtrack*)tr;
+ if (!trFilter->IsSelected(esdTr)) continue;
+ }
+ else if (trList) { // If AOD input
+ Double_t pTMin = 0, pTMax = 0, etaMin = 0, etaMax = 0, minNCl = 0;
+ UInt_t bit = 0;
+ if ((fFlags & kTPC) == kTPC) pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 128;
+ if ((fFlags & kHybrid) == kHybrid) pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 272;
+
+ AliAODTrack* aodTr = (AliAODTrack*)tr;
if (aodTr->GetID() > -1) continue;
- if (!aodTr->TestFilterBit(128) || !aodTr->Pt() > pTMax || aodTr->Pt() < pTMin ||
+ if (!aodTr->TestFilterBit(bit) || aodTr->Pt() > pTMax || aodTr->Pt() < pTMin ||
aodTr->Eta() > etaMax || aodTr->Eta() < etaMin || aodTr->GetTPCNcls() < minNCl) continue;
}
- if (tr->GetTPCsignal() < tpcdEdx) continue;
+
+// if (tr->GetTPCsignal() < tpcdEdx) continue;
// Track accepted
Double_t eta = tr->Eta();
if (((fFlags & kSPD) || (fFlags & kEtaGap)) && TMath::Abs(eta) < fEtaGap) continue;
Double_t phi = tr->Phi();
+ Double_t weight = 1.;
+
if ((mode & kFillRef)) {
- fCumuRef->Fill(eta, 0.);// mult goes in underflowbin - no visual, but not needed?
- fdNdedpRefAcc->Fill(eta, phi);
+ fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
+ fdNdedpRefAcc->Fill(eta, phi, weight);
}
if ((mode & kFillDiff)) {
- fCumuDiff->Fill(eta, 0);
- fdNdedpDiffAcc->Fill(eta, phi);
+ fCumuDiff->Fill(eta, 0., weight);
+ fdNdedpDiffAcc->Fill(eta, phi, weight);
}
for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
- Double_t cosnPhi = TMath::Cos(n*phi);
- Double_t sinnPhi = TMath::Sin(n*phi);
+ Double_t cosnPhi = weight*TMath::Cos(n*phi);
+ Double_t sinnPhi = weight*TMath::Sin(n*phi);
// fill ref
if ((mode & kFillRef)) {
fCumuRef->Fill(eta, cosBin, cosnPhi);
// Fill out NUA hists
for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
- if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
- if ((fFlags & kTPC) && (fFlags && kSPD) && !(fFlags & kEtaGap)) eta = -eta;
+ if (fCumuRef->GetBinContent(etaBin, 0) <= 3) continue;
+ if ((fFlags & kTracks) && (fFlags && kSPD) && !(fFlags & kEtaGap)) eta = -eta;
for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
}
}
for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
+ Double_t refetaBin = fCumuRef->GetXaxis()->FindBin(eta);
+ if (fCumuRef->GetBinContent(refetaBin, 0) <= 3) continue;
if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
Double_t refEta = eta;
- if ((fFlags & kTPC) && (fFlags && kSPD) && !(fFlags & kEtaGap)) refEta = -eta;
+ if ((fFlags & kTracks) && (fFlags && kSPD) && !(fFlags & kEtaGap)) refEta = -eta;
Int_t refEtaBinA = fCumuRef->GetXaxis()->FindBin(refEta);
if ((fFlags & kEtaGap)) refEta = -eta;
Int_t refEtaBinB = fCumuRef->GetXaxis()->FindBin(refEta);
if (refEtaBinA != prevRefEtaBin) {
- prevRefEtaBin = refEtaBinA;
// Reference flow
multA = fCumuRef->GetBinContent(refEtaBinA, 0);
dQnReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n));
dQnReB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberCos(n));
dQnImB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberSin(n));
- if (multA <= 3 || multB <= 3) return;
+ if (multA <= 3 || multB <= 3) continue;
// The reference flow is calculated
// 2-particle
if ((fFlags & kStdQC)) {
cumuRef->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
cumuRef->Fill(eta, cent, k3pWeight, multA*(multA-1.)*(multA-2.));
} // End of QC{4}
+ prevRefEtaBin = refEtaBinA;
} // End of reference flow
// For each etaBin bin the necessary values for differential flow is calculated
Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
// Differential flow calculations for each eta bin is done:
// 2-particle differential flow
- if ((fFlags & kStdQC) && !(fFlags & kTPC)) {
+ if ((fFlags & kStdQC) && (!(fFlags & kTracks) || ((fFlags & kTracks) && (fFlags & kMC) && !(fFlags & kSPD) && TMath::Abs(eta) < 0.75))) {
mq = mp;
qnRe = pnRe;
qnIm = pnIm;
// The next line covers both cases.
qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
// Extra NUA term from 2n cosines and sines
- Double_t den = 1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB);
+ Double_t den = 1+(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB);
if (den != 0) qc2 /= den;
else qc2 = 0;
}
// Old nua
qc2Prime -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB;
// Extra NUA term from 2n cosines and sines
- qc2Prime /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
+ qc2Prime /= (1.+(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
}
if (!TMath::IsNaN(qc2Prime)) {
quality->Fill((n-2)*qualityFactor+3, Int_t(cent));
// Old nua
qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
// Extra NUA term from 2n cosines and sines
- qc2 /= (1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB));
+ qc2 /= (1+(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB));
}
if (qc2 <= 0) {
if (fDebug > 0)
qc2PrimeA -= cosP1nPsi*cosP1nPhiA + sinP1nPsi*sinP1nPhiA;
qc2PrimeB -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB; // Is this OK?
// Extra NUA term from 2n cosines and sines
- if (cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA != 1.) qc2PrimeA /= (1.-(cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA));
- if (cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB != 1.) qc2PrimeB /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
+ if (cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA != -1.) qc2PrimeA /= (1.+(cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA));
+ if (cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB != -1.) qc2PrimeB /= (1.+(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
}
if (!TMath::IsNaN(qc2PrimeA) && !TMath::IsNaN(qc2PrimeB) && qc2 != 0) {
if (qc2PrimeA*qc2PrimeB >= 0) {
TString type = "Standard QC{2} and QC{4} calculations.";
if ((fFlowFlags & kEtaGap)) type = "QC{2} with a rapidity gap.";
if ((fFlowFlags & k3Cor)) type = "QC{2} with 3 correlators.";
- if ((fFlowFlags & kTPC)) type.ReplaceAll(".", " with TPC tracks for reference.");
+ if ((fFlowFlags & kTPC) == kTPC) type.ReplaceAll(".", " with TPC tracks for reference.");
+ if ((fFlowFlags & kHybrid) == kHybrid) type.ReplaceAll(".", " with hybrid tracks for reference.");
Printf("QC calculation type :\t%s", type.Data());
Printf("Symmetrize ref. flow wrt. eta = 0 :\t%s", ((fFlowFlags & kSymEta) ? "true" : "false"));
Printf("Apply NUA correction terms :\t%s", ((fFlowFlags & kNUAcorr) ? "true" : "false"));
Printf("Satellite vertex flag :\t%s", ((fFlowFlags & kSatVtx) ? "true" : "false"));
Printf("FMD sigma cut: :\t%f", fFMDCut);
Printf("SPD sigma cut: :\t%f", fSPDCut);
- if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTPC))
+ if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks))
Printf("Eta gap: :\t%f", fEtaGap);
Printf("=======================================================");
}
else if ((flags & k3Cor)) type = "3Cor";
else type = "UNKNOWN";
if (prependUS) type.Prepend("_");
- if ((flags & kTPC)) type.Append("TPCTr");
+ if ((flags & kTPC) == kTPC) type.Append("TPCTr");
+ if ((flags & kHybrid) == kHybrid) type.Append("HybTr");
+
return type.Data();
}
//_____________________________________________________________________