// --- Make the task and add it to the manager ---------------------
AliForwardQATask* task = new AliForwardQATask("forwardQA");
- // task->SetBLow(blow);
- // task->SetBLow(bhigh);
mgr->AddTask(task);
+
+ // --- Cuts --------------------------------------------------------
+ // Old style cuts
Double_t nXi = mc ? 2 : 2; //HHD test
Bool_t includeSigma = false; //true;
-
+ // High cuts for sharing filter
+ AliFMDMultCuts cHigh;
+ cHigh.SetMPVFraction(0.7);
+ cHigh.SetMultCuts(-1);
+ // Low cuts for sharing and density calculator
+ AliFMDMultCuts cLow;
+ cLow.SetMultCuts(0.1, 0.1, 0.12, 0.1, 0.12);
+ // Density cuts
AliFMDMultCuts cDensity;
- //c2.SetNXi(mc ? 1 : 1);
- // c2.SetIncludeSigma(false);
- //c2.SetMPVFraction(0.5);
- //Double_t factor = 1.2;
- cDensity.SetMultCuts(0.3, 0.3, 0.3, 0.3, 0.3);
+ cDensity.SetMPVFraction(0.7);
+ cDensity.SetMultCuts(-1);
+
// --- Set parameters on the algorithms ----------------------------
// Set the number of SPD tracklets for which we consider the event a
task->GetEventInspector().SetLowFluxCut(1000);
// Set the maximum error on v_z [cm]
task->GetEventInspector().SetMaxVzErr(0.2);
+ // Disable use of code from 1st Physics
+ task->GetEventInspector().SetUseFirstPhysicsVtx(kFALSE);
// --- Set parameters on energy loss fitter ------------------------
// Set the eta axis to use - note, this overrides whatever is used
// Disable use of angle corrected signals in the algorithm
task->GetSharingFilter().SetZeroSharedHitsBelowThreshold(false);
// Enable use of angle corrected signals in the algorithm
- task->GetSharingFilter().GetHCuts().SetMultCuts(-1);
+ task->GetSharingFilter().SetHCuts(cHigh);
+ // Multiplicity cut
+ // task->GetSharingFilter().GetHCuts().SetMultCuts(-1);
// Set the number of xi's (width of landau peak) to stop at
- task->GetSharingFilter().GetHCuts().SetNXi(nXi);
+ // task->GetSharingFilter().GetHCuts().SetNXi(nXi);
// Set whether or not to include sigma in cut
- task->GetSharingFilter().GetHCuts().SetIncludeSigma(includeSigma);
- // Enable use of angle corrected signals in the algorithm
- task->GetSharingFilter().SetLCuts(cDensity);
+ // task->GetSharingFilter().GetHCuts().SetIncludeSigma(includeSigma);
+ // Lower cuts from object
+ task->GetSharingFilter().SetLCuts(cLow);
+ // Use simplified sharing algorithm
+ task->GetSharingFilter().SetUseSimpleSharing(kTRUE);
// --- Density calculator ------------------------------------------
// Set the maximum number of particle to try to reconstruct
- task->GetDensityCalculator().SetMaxParticles(10);
+ task->GetDensityCalculator().SetMaxParticles(3);
// Wet whether to use poisson statistics to estimate N_ch
task->GetDensityCalculator().SetUsePoisson(true);
+ // How to lump for Poisson calculator - 64 strips, 4 regions
+ task->GetDensityCalculator().SetLumping(64,4);
// Set whether or not to include sigma in cut
task->GetDensityCalculator().SetCuts(cDensity);
// Set whether or not to use the phi acceptance
// AliFMDDensityCalculator::kPhiNoCorrect
// AliFMDDensityCalculator::kPhiCorrectNch
// AliFMDDensityCalculator::kPhiCorrectELoss
- task->GetDensityCalculator()
- .SetUsePhiAcceptance(AliFMDDensityCalculator::kPhiCorrectNch);
+ task->GetDensityCalculator().
+ SetUsePhiAcceptance(AliFMDDensityCalculator::kPhiNoCorrect);
+
// --- Set limits on fits the energy -------------------------------
// Maximum relative error on parameters
AliFMDCorrELossFit::ELossFit::fgMaxRelError = .12;
AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
// Maximum value of reduced chi^2
AliFMDCorrELossFit::ELossFit::fgMaxChi2nu = 20;
+
+ // --- Debug -------------------------------------------------------
+ // Set the overall debug level (1: some output, 3: a lot of output)
+ task->SetDebug(3);
// --- Make the output container and connect it --------------------
// TString outputfile = ;
AliAnalysisDataContainer *output =
mgr->CreateContainer("ForwardResults", TList::Class(),
AliAnalysisManager::kParamContainer,
- AliAnalysisManager::GetCommonFileName());
+ "trending.root");
+ // AliAnalysisManager::GetCommonFileName());
mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
mgr->ConnectOutput(task, 1, histOut);
mgr->ConnectOutput(task, 2, output);
-
return task;
}
+//
+// EOF
+//
Int_t nDists = dists->GetEntries();
Int_t low = nDists;
Int_t high = 0;
+ Int_t nEmpty = 0;
+ Int_t nLow = 0;
+ Int_t nFitted= 0;
for (Int_t i = 0; i < nDists; i++) {
TH1D* dist = static_cast<TH1D*>(dists->At(i));
// Ignore empty histograms altoghether
- if (!dist || dist->GetEntries() <= 0) continue;
+ if (!dist || dist->GetEntries() <= 0) {
+ nEmpty++;
+ continue;
+ }
// Scale to the bin-width
dist->Scale(1., "width");
dist->Scale(1/max);
// Check that we have enough entries
- if (dist->GetEntries() <= minEntries) {
+ Int_t nEntries = Int_t(dist->GetEntries());
+ if (nEntries <= minEntries) {
AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)",
- i, dist->GetName(), Int_t(dist->GetEntries()),
- minEntries));
+ i, dist->GetName(), nEntries, minEntries));
+ nLow++;
continue;
}
TF1* res = FitHist(dist,lowCut,nParticles,minusBins,
relErrorCut,chi2nuCut);
if (!res) continue;
+ nFitted++;
// dist->GetListOfFunctions()->Add(res);
// Store eta limits
hA[j]->SetBinError(i+1, res->GetParError(kA+j));
}
}
+ printf("%s: Out of %d histograms, %d where empty, %d had too little data,"
+ "leaving %d to be fitted, of which %d succeeded\n",
+ GetName(), nDists, nEmpty, nLow, nDists-nEmpty-nLow, nFitted);
+
+ TH1* status = new TH1I("status", "Status of Fits", 5, 0, 5);
+ status->GetXaxis()->SetBinLabel(1, "Total");
+ status->GetXaxis()->SetBinLabel(2, "Empty");
+ status->GetXaxis()->SetBinLabel(3, Form("<%d", minEntries));
+ status->GetXaxis()->SetBinLabel(4, "Candidates");
+ status->GetXaxis()->SetBinLabel(5, "Fitted");
+ status->SetXTitle("Status");
+ status->SetYTitle("# of #Delta distributions");
+ status->SetBinContent(1, nDists);
+ status->SetBinContent(2, nEmpty);
+ status->SetBinContent(3, nLow);
+ status->SetBinContent(4, nDists-nLow-nEmpty);
+ status->SetBinContent(5, nFitted);
+ status->SetFillColor(Color());
+ status->SetFillStyle(3001);
+ status->SetLineColor(Color());
+ status->SetDirectory(0);
+ status->SetStats(0);
+ pars->Add(status);
// Fit the full-ring histogram
TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
//____________________________________________________________________
void
-AliFMDEventInspector::StoreInformation()
+AliFMDEventInspector::StoreInformation(Int_t runNo)
{
// Write TNamed objects to output list containing information about
// the running conditions
TNamed* sys = new TNamed("sys", "");
TNamed* sNN = new TNamed("sNN", "");
TNamed* fld = new TNamed("field", "");
+ TNamed* run = new TNamed("runNo", Form("%d", runNo));
sys->SetTitle(AliForwardUtil::CollisionSystemString(fCollisionSystem));
sNN->SetTitle(AliForwardUtil::CenterOfMassEnergyString(fEnergy));
fld->SetTitle(AliForwardUtil::MagneticFieldString(fField));
sys->SetUniqueID(fCollisionSystem);
sNN->SetUniqueID(fEnergy);
fld->SetUniqueID(fField);
+ run->SetUniqueID(runNo);
fList->Add(sys);
fList->Add(sNN);
fList->Add(fld);
+ fList->Add(run);
}
//____________________________________________________________________
fField =
AliForwardUtil::ParseMagneticField(esd->GetMagneticField());
- StoreInformation();
+ StoreInformation(esd->GetRunNumber());
if (fCollisionSystem == AliForwardUtil::kUnknown ||
fEnergy <= 0 ||
TMath::Abs(fField) > 10)
* - sys Contains the collision system string and identifier.
* - sNN Contains the center-of-mass energy per nucleon (GeV)
* - field Contains the L3 magnetic field (kG)
+ * - run Contains the run number
*
+ * @param runNo Run number - read off from ESD event
*/
- virtual void StoreInformation();
+ virtual void StoreInformation(Int_t runNo);
protected:
/**
* Read the trigger information from the ESD event
//____________________________________________________________________
void
-AliFMDMCEventInspector::StoreInformation()
+AliFMDMCEventInspector::StoreInformation(Int_t runNo)
{
// Store information about running conditions in the output list
if (!fList) return;
- AliFMDEventInspector::StoreInformation();
+ AliFMDEventInspector::StoreInformation(runNo);
TNamed* mc = new TNamed("mc", fProduction.Data());
mc->SetUniqueID(1);
fList->Add(mc);
* The presence of this indicate MC data.
*
* - mc Nothing special, and unique id is 1
+ *
+ * @param runNo Run number
*/
- virtual void StoreInformation();
+ virtual void StoreInformation(Int_t runNo);
/**
* Read the production details
*
// d detector
// r ring
//
- fBefore = new TH1D("esdEloss", "Energy loss (reconstruction)",
- 600, 0, 15);
+ fBefore = new TH1D("esdEloss", Form("Energy loss in %s (reconstruction)",
+ GetName()), 600, 0, 15);
fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
fBefore->SetFillColor(Color());
fBefore->SetDirectory(0);
fAfter = static_cast<TH1D*>(fBefore->Clone("anaEloss"));
- fAfter->SetTitle("Energy loss in %s (sharing corrected)");
+ fAfter->SetTitle(Form("Energy loss in %s (sharing corrected)", GetName()));
fAfter->SetFillColor(Color()+2);
fAfter->SetLineStyle(1);
fAfter->SetDirectory(0);
fSingle = new TH1D("singleEloss", "Energy loss (single strips)",
600, 0, 15);
- fSingle->SetXTitle("#Delta E/#Delta E_{mip}");
- fSingle->SetYTitle("P(#Delta E/#Delta E_{mip})");
- fSingle->SetFillColor(kMagenta);
+ fSingle->SetXTitle("#Delta/#Delta_{mip}");
+ fSingle->SetYTitle("P(#Delta/#Delta_{mip})");
+ fSingle->SetFillColor(Color());
fSingle->SetFillStyle(3001);
fSingle->SetLineColor(kBlack);
fSingle->SetLineStyle(2);
fSingle->SetDirectory(0);
- fDouble = new TH1D("doubleEloss", "Energy loss (two strips)",
- 600, 0, 15);
- fDouble->SetXTitle("#Delta E/#Delta E_{mip}");
- fDouble->SetYTitle("P(#Delta E/#Delta E_{mip})");
- fDouble->SetFillColor(kMagenta+1);
- fDouble->SetFillStyle(3001);
- fDouble->SetLineColor(kBlack);
- fDouble->SetLineStyle(2);
+ fDouble = static_cast<TH1D*>(fSingle->Clone("doubleEloss"));
+ fDouble->SetTitle("Energy loss (two strips)");
+ fDouble->SetFillColor(Color()+1);
fDouble->SetDirectory(0);
-
- fTriple = new TH1D("tripleEloss", "Energy loss (three strips)",
- 600, 0, 15);
- fTriple->SetXTitle("#Delta E/#Delta E_{mip}");
- fTriple->SetYTitle("P(#Delta E/#Delta E_{mip})");
- fTriple->SetFillColor(kMagenta+2);
- fTriple->SetFillStyle(3001);
- fTriple->SetLineColor(kBlack);
- fTriple->SetLineStyle(2);
+
+ fTriple = static_cast<TH1D*>(fSingle->Clone("tripleEloss"));
+ fTriple->SetTitle("Energy loss (three strips)");
+ fTriple->SetFillColor(Color()+2);
fTriple->SetDirectory(0);
//Int_t nBinsForInner = (r == 'I' ? 32 : 16);
fSinglePerStrip = new TH2D("singlePerStrip", "SinglePerStrip",
600,0,15, nBinsForInner,0,nStrips);
- fSinglePerStrip->SetXTitle("Eloss");
- fSinglePerStrip->SetYTitle("Strip");
+ fSinglePerStrip->SetXTitle("#Delta/#Delta_{mip}");
+ fSinglePerStrip->SetYTitle("Strip #");
fSinglePerStrip->SetZTitle("Counts");
fSinglePerStrip->SetDirectory(0);
fDistanceBefore = new TH1D("distanceBefore", "Distance before sharing",
- nStrips , 0,nStrips );
+ nStrips , 0,nStrips );
fDistanceBefore->SetXTitle("Distance");
fDistanceBefore->SetYTitle("Counts");
fDistanceBefore->SetFillColor(kGreen+2);
fDistanceBefore->SetLineStyle(2);
fDistanceBefore->SetDirectory(0);
- fDistanceAfter = new TH1D("distanceAfter", "Distance after sharing",
- nStrips , 0,nStrips );
- fDistanceAfter->SetXTitle("Distance");
- fDistanceAfter->SetYTitle("Counts");
+ fDistanceAfter = static_cast<TH1D*>(fDistanceBefore->Clone("distanceAfter"));
+ fDistanceAfter->SetTitle("Distance after sharing");
fDistanceAfter->SetFillColor(kGreen+1);
- fDistanceAfter->SetFillStyle(3001);
- fDistanceAfter->SetLineColor(kBlack);
- fDistanceAfter->SetLineStyle(2);
fDistanceAfter->SetDirectory(0);
fEnergyFitter(),
fSharingFilter(),
fDensityCalculator(),
- fList(0)
+ fList(0),
+ fDebug(0)
{
//
// Constructor
fEnergyFitter("energy"),
fSharingFilter("sharing"),
fDensityCalculator("density"),
- fList(0)
+ fList(0),
+ fDebug(0)
{
//
// Constructor
fEnergyFitter(o.fEnergyFitter),
fSharingFilter(o.fSharingFilter),
fDensityCalculator(o.fDensityCalculator),
- fList(o.fList)
+ fList(o.fList),
+ fDebug(o.fDebug)
{
//
// Copy constructor
fDensityCalculator = o.fDensityCalculator;
fHistos = o.fHistos;
fList = o.fList;
+ fDebug = o.fDebug;
return *this;
}
// Parameters:
// dbg Debug level
//
+ fDebug = dbg;
fEventInspector.SetDebug(dbg);
fEnergyFitter.SetDebug(dbg);
fSharingFilter.SetDebug(dbg);
// Parameters:
// option Not used
//
+ if (fDebug) AliInfo("In Forwards terminate");
TStopwatch swt;
swt.Start();
// Make a deep copy and post that as output 2
TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults",
list->GetName())));
+ if (fDebug) AliInfoF("Posting post processing results to %s",
+ list2->GetName());
list2->SetOwner();
PostData(2, list2);
AliFMDDensityCalculator fDensityCalculator; // Algorithm
TList* fList; // Output list
+ Int_t fDebug; // Debug flag
ClassDef(AliForwardQATask,1) // Forward QA class
};
//c.SetNXi(mc ? 1 : 1);
//c.SetIncludeSigma(true);
//c.SetMPVFraction(0.5);
-
Double_t factor = 1.;
//if(mc) factor = 1.15;
-
cSharing.SetMultCuts(0.3, 0.3, 0.3, 0.3, 0.3);
AliFMDMultCuts cDensity;
// Disable use of angle corrected signals in the algorithm
task->GetSharingFilter().SetZeroSharedHitsBelowThreshold(false);
// Enable use of angle corrected signals in the algorithm
-
task->GetSharingFilter().GetHCuts().SetMultCuts(-1);
// Set the number of xi's (width of landau peak) to stop at
task->GetSharingFilter().GetHCuts().SetNXi(nXi);
task->GetSharingFilter().SetLCuts(cSharing);
- // --- Density calculator
+ // --- Density calculator ------------------------------------------
// Set the maximum number of particle to try to reconstruct
task->GetDensityCalculator().SetMaxParticles(10);
// Wet whether to use poisson statistics to estimate N_ch
task->GetDensityCalculator().SetUsePoisson(true);
-
// Set whether or not to include sigma in cut
task->GetDensityCalculator().SetCuts(cDensity);
--- /dev/null
+#!/bin/bash
+
+# --------------------------------------------------------------------
+uid=`id -u`
+genv_file=/tmp/gclient_env_${uid}
+
+if test ! -f ${genv_file} ; then
+ echo "No such file: ${genv_file}, please do alien-token-init" >/dev/stderr
+ exit 1
+fi
+alien-token-info | grep -q "Token is still valid"
+if test $? -ne 0 ; then
+ echo "Token not valid, please re-new" > /dev/stderr
+ exit 1
+fi
+
+# --------------------------------------------------------------------
+verb=0
+
+mess()
+{
+ if test $verb -lt 1 ; then return ; fi
+ echo $@
+}
+
+# --------------------------------------------------------------------
+usage()
+{
+ cat <<EOF
+Usage: $0 -p PRODUCTION [OPTIONS]
+
+Options:
+ -h,--help This help
+ -v,--verbose Be verbose
+ -p,--production IDENTIFIER Production identifier [$prod]
+ -y,--year YEAR Year of production [$year]
+ -P,--pass NUMBER Reconstruction pass number [$pass]
+ -d,--destination DIRECTORY Directory to store result in [$dest]
+ -r,--run NUMBER Run number [$runn]
+ -q,--qa NUMBER QA number
+ -a,--archives Get ZIP archives
+ -n,--no-action Run dry (do not copy files)
+
+If run number is set, get the parts of next-to-last merge of that run only.
+
+EOF
+}
+
+# --------------------------------------------------------------------
+check_file()
+{
+ inp=$1
+ scr=`mktemp -u Check_XXXXXXXX`
+ cat <<EOF > ${scr}.C
+void ${scr}()
+{
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+ gSystem->Load("libTENDER");
+ // gSystem->Load("libTENDERSupplies");
+ gSystem->Load("libPWG1");
+ gSystem->Load("libPWG3base");
+ TFile* file = TFile::Open("${inp}", "READ");
+ if (!file) {
+ Error("${scr}", "No such file ${inp}");
+ exit(1);
+ return false;
+ }
+ TObject* forward1 = file->Get("Forward");
+ if (!forward1) {
+ Error("${scr}", "No Forward object found in ${inp}");
+ exit(2);
+ ret = false;
+ }
+ TObject* forward2 = file->Get("ForwardResults");
+ if (!forward2) {
+ Error("${scr}", "No ForwardResults object found in ${inp}");
+ exit(4);
+ }
+ exit(0);
+}
+EOF
+ # cat ${scr}.C
+ mess -n "aliroot -l -b -q ${scr}.C "
+ aliroot -l -b -q ${scr}.C > /dev/null 2>&1
+ ret=$?
+ mess "-> $ret (0: good, 1: no file, 2: no Forward, 4: no ForwardResults"
+ rm -f ${scr}.C
+ return $ret
+}
+
+# --------------------------------------------------------------------
+year=
+prod=
+pass=1
+qual=
+dest=.
+runn=0
+qano=0
+noac=0
+arch=0
+
+while test $# -gt 0 ; do
+ case $1 in
+ -h|--help) usage ; exit 0 ;;
+ -v|--verbose) let verb=$verb+1 ;;
+ -p|--production) prod=$2 ; shift ;;
+ -P|--pass) pass=$2 ; shift ;;
+ -Q|--prepass) qual=$2 ; shift ;;
+ -y|--year) year=$2 ; shift ;;
+ -d|--destination) dest=$2 ; shift ;;
+ -r|--run) runn=$2 ; shift ;;
+ -q|--qa) qano=$2 ; shift ;;
+ -n|--no-action) noac=1 ;;
+ -a|--archives) arch=1 ;;
+ *) echo "$0: Unknown option $1" > /dev/stderr ; exit 2 ;;
+ esac
+ shift
+done
+
+# --------------------------------------------------------------------
+if test "x$prod" = "x" ; then
+ echo "No production identifier given" > /dev/stderr
+ exit 2
+fi
+
+if test "x$year" = "x" ; then
+ year=`echo $prod | sed -e 's/LHC\(..\).*/\1/'`
+ if test "x$year" = "x" ; then
+ echo "Couldn't get year from production identifier $prod" > /dev/stderr
+ exit 2
+ fi
+fi
+
+redir="/dev/null"
+if test $verb -gt 1 ; then redir1=/dev/stderr ; fi
+
+# --------------------------------------------------------------------
+path=/alice/data/20${year}/${prod}/
+store=${dest}/${prod}/${qual}pass${pass}
+search="ESDs/${qual}pass${pass}/"
+if test $runn -gt 0 ; then
+ path=`printf "${path}%09d/ESDs/${qual}pass${pass}/" $runn`
+ store=`printf "${store}/%09d" $runn`
+ search=
+fi
+if test $qano -gt 0 ; then
+ path=`printf "%sQA%02d/" $path $qano`
+fi
+if test $arch -gt 0 ; then
+ search="${search}QA_archive.zip"
+else
+ search="${search}QAresults.root"
+fi
+
+cat <<EOF
+Settings:
+
+ Production: $prod
+ Year: $year
+ Pass: $pass
+ Pass qualifier: $qual
+ Path: $path
+ Destination: $dest
+ Store: $store
+ Run number: $runn
+ Search string: $search
+EOF
+
+# --------------------------------------------------------------------
+mkdir -p ${store}
+mess "alien_find ${path} ${search}"
+files=`alien_find ${path} ${search} | grep -v "files found" 2> ${redir}`
+j=0
+runs=
+for i in $files ; do
+ b=`echo $i | sed -e "s,${path},,"`
+ d=
+ if test $runn -gt 0 ; then
+ r=`echo $b | sed -e "s,[0-9]*${runn}\([0-9.]*\)/.*,\1," | tr '.' '_'`
+ if test $arch -lt 1 ; then
+ o=${store}/QAresults_${r}.zip
+ else
+ d=${store}/${r}
+ mkdir -p $d
+ o=${d}/QAarchive.zip
+ fi
+ else
+ r=`echo $b | sed -e "s,/.*,,"`
+ o=${store}/QAresults_${r}.root
+ fi
+ runs="$runs $r"
+ mess "$i -> $o"
+ if test $noac -lt 1 && test ! -f $o ; then
+ mess "alien_cp alien:${i} file:${o}"
+ alien_cp alien:${i} file:${o} > ${redir} 2>&1
+ fi
+ if test $noac -lt 1 && test $arch -lt 1 ; then check_file $o ; fi
+ if test $noac -lt 1 && test $arch -gt 0 ; then
+ (cd ${d} && unzip QAarchive.zip) > $redir 2>&1
+ fi
+ let j=$j+1
+done
+mess "Got a total of $j files for runs $runs"
+
+# --------------------------------------------------------------------
+#
+# EOF
+#
+
+
\ No newline at end of file
--- /dev/null
+// --- Get a single histogram ----------------------------------------
+TH1*
+GetOne(const TList* list, const char* which, bool mirror)
+{
+ if (!list) {
+ Error("GetOne", "No list passed");
+ return 0;
+ }
+ TString n(Form("dndeta%s_rebin05", which));
+ if (mirror) n.Append("_mirror");
+
+ TObject* o = list->FindObject(n);
+ if (!o) {
+ Error("GetOne", "Object %s not found in %s", n.Data(), list->GetName());
+ return 0;
+ }
+ TH1* ret = static_cast<TH1*>(o);
+ ret->SetLineColor(ret->GetMarkerColor());
+ ret->SetDirectory(0);
+ return ret;
+}
+
+// --- Make point-to-point systematic errors -------------------------
+TH1*
+MakeSysError(TH1* h, Double_t sysErr)
+{
+ TString n(h->GetName());
+ n.Append("_syserr");
+ TH1* ret = static_cast<TH1*>(h->Clone(n));
+ ret->SetMarkerStyle(1);
+ ret->SetMarkerSize(0);
+ ret->SetMarkerColor(kBlue-10);
+ ret->SetFillColor(kBlue-10);
+ ret->SetFillStyle(1001);
+ ret->SetLineColor(kBlue-10);
+ ret->SetLineWidth(0);
+ ret->SetDirectory(0);
+
+ for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
+ Double_t c = ret->GetBinContent(i);
+ if (c < 0.001) {
+ ret->SetBinContent(i,0);
+ ret->SetBinError(i,0);
+ }
+ Double_t e = c * sysErr/100;
+ ret->SetBinError(i,e);
+ }
+
+ return ret;
+}
+
+// --- Turn a TGraphAsymmErrors into a histogram ---------------------
+TH1* Graph2Hist(const TGraphAsymmErrors* g)
+{
+ Int_t nBins = g->GetN();
+ TArrayF bins(nBins+1);
+ TArrayF y(nBins);
+ TArrayF ey(nBins);
+ Double_t dx = 0;
+ Double_t xmin = 10000;
+ Double_t xmax = -10000;
+ for (Int_t i = 0; i < nBins; i++) {
+ Double_t x = g->GetX()[i];
+ Double_t exl = g->GetEXlow()[i];
+ Double_t exh = g->GetEXhigh()[i];
+ xmin = TMath::Min(x-exl, xmin);
+ xmax = TMath::Max(x+exh, xmax);
+ bins.fArray[i] = x-exl;
+ bins.fArray[i+1] = x+exh;
+ Double_t dxi = exh+exl;
+ if (dxi == 0 && i != 0) dxi = bins.fArray[i]-bins.fArray[i-1];
+ if (dx == 0) dx = dxi;
+ else if (dxi != dx) dx = 0;
+
+ y.fArray[i] = g->GetY()[i];
+ ey.fArray[i] = TMath::Max(g->GetEYlow()[i],g->GetEYhigh()[i]);
+
+ }
+ TString name(g->GetName());
+ TString title(g->GetTitle());
+ TH1D* h = 0;
+ if (dx != 0) {
+ h = new TH1D(name.Data(), title.Data(), nBins,
+ bins[0]-dx/2, bins[nBins]+dx/2);
+ }
+ else {
+ h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
+ }
+ for (Int_t i = 1; i <= nBins; i++) {
+ h->SetBinContent(i, y.fArray[i-1]);
+ h->SetBinError(i, ey.fArray[i-1]);
+ }
+ h->SetMarkerStyle(g->GetMarkerStyle());
+ h->SetMarkerColor(g->GetMarkerColor());
+ h->SetMarkerSize(g->GetMarkerSize());
+ h->SetDirectory(0);
+
+ return h;
+}
+
+// --- Set Histogram properties --------------------------------------
+void SetAttributes(TH1* h, UShort_t sNN, Int_t which, Bool_t mirror)
+{
+ Int_t marker = (sNN == 900 ? 20 : 21) + (mirror ? 4 : 0);
+ Int_t color = (which == 0 ? kRed+2 :
+ which == 1 ? kMagenta+2 :
+ kBlue + 2);
+ h->SetMarkerColor(color);
+ h->SetLineColor(color);
+ h->SetMarkerStyle(marker);
+}
+
+
+// --- Get publlished data via script --------------------------------
+TH1* GetPublished(UShort_t sNN, Bool_t isNSD)
+{
+ Long_t ptr = 0;
+ Int_t typ = (isNSD ? 4 : 1);
+
+ ptr = gROOT->ProcessLine(Form("GetSingle(2,1,%d,%d)", sNN, typ));
+ if (!ptr) return 0;
+
+ TGraphAsymmErrors* g = reinterpret_cast<TGraphAsymmErrors*>(ptr);
+ TH1* h = Graph2Hist(g);
+ h->SetLineColor(h->GetMarkerColor());
+
+ return h;
+}
+
+// --- Add points to stack -------------------------------------------
+Double_t
+AddToStack(THStack* stack, UShort_t sNN, Bool_t isNSD,
+ TList* forward, TList* central,
+ Double_t fwdSysErr, Double_t cenSysErr,
+ Double_t strangeCorr)
+{
+ TH1* fd1 = GetOne(forward, "Forward", false);
+ TH1* fd2 = GetOne(forward, "Forward", true);
+ fd1->Scale(1-strangeCorr/100);
+ fd2->Scale(1-strangeCorr/100);
+ TH1* fs1 = MakeSysError(fd1, fwdSysErr);
+ TH1* fs2 = MakeSysError(fd2, fwdSysErr);
+
+ TH1* cd = 0;
+ TH1* cs = 0;
+ if (central) {
+ cd = GetOne(central, "Central", false);
+ cs = MakeSysError(cd, cenSysErr);
+ }
+ else {
+ cd = GetPublished(sNN, isNSD);
+ }
+ SetAttributes(fd1, sNN, 0, false);
+ SetAttributes(fd2, sNN, 0, true);
+ if (cd) SetAttributes(cd, sNN, central ? 1 : 2, false);
+
+ if (cs) stack->Add(cs, "e2");
+ stack->Add(fs1, "e2");
+ stack->Add(fs2, "e2");
+ if (cd) stack->Add(cd, "ep");
+ stack->Add(fd1, "ep");
+ stack->Add(fd2, "ep");
+
+ Double_t mcd = (cd ? cd->GetMaximum() : 0);
+ Double_t mfs = fs1->GetMaximum();
+
+ return TMath::Max(mcd, mfs);
+}
+
+// --- Find a list in directory --------------------------------------
+TList*
+GetList(const TDirectory* d, const char* what)
+{
+ if (!d) {
+ Error("GetList", "No diretory passed");
+ return 0;
+ }
+ TList* p = static_cast<TList*>(d->Get(Form("%sResults", what)));
+ if (!p) {
+ Error("GetList", "%sResults not found in %s", what, d->GetName());
+ return 0;
+ }
+ TList* r = static_cast<TList*>(p->FindObject("all"));
+ if (!r) {
+ Error("GetList", "all not found in %s", p->GetName());
+ return 0;
+ }
+ return r;
+}
+
+void
+AdjustAxis(TAxis* a)
+{
+ a->SetTitleFont(132);
+ a->SetLabelFont(132);
+ a->SetTitleSize(0.08);
+ a->SetLabelSize(0.08);
+ a->SetTitleOffset(0.5);
+ a->SetNdivisions(10);
+}
+
+// --- Make a stack for a single plot --------------------------------
+THStack*
+MakeStack(Bool_t isNSD, Bool_t showClusters,
+ Double_t strangeCorr, Double_t maxFactor=1.1)
+{
+ const char* trg = isNSD ? "nsd" : "inel";
+ TFile* f0900 =TFile::Open(Form("forward_dndeta_%s%04d.root",trg,900, "READ"));
+ TFile* f7000 =TFile::Open(Form("forward_dndeta_%s%04d.root",trg,7000,"READ"));
+ if (!f0900 || !f7000) {
+ Error("MakeStack", "Failed to open one or more files (%p, %p)",
+ f0900, f7000);
+ return 0;
+ }
+ TList* forward0900 = GetList(f0900, "Forward");
+ TList* forward7000 = GetList(f7000, "Forward");
+ TList* central0900 = 0;
+ TList* central7000 = showClusters ? GetList(f7000, "Central") : 0;
+
+ THStack* stack = new THStack("stack", "Stack");
+ Double_t sysDen = 5;
+ Double_t sysPt = 2;
+ Double_t sysMix = 2;
+ Double_t sysNch = 5;
+ Double_t sysNor = 2;
+ Double_t fwdSys = TMath::Sqrt(sysDen * sysDen +
+ sysPt * sysPt +
+ sysMix * sysMix +
+ sysNch * sysNch +
+ sysNor * sysNor);
+ Info("", "Forward systematic error: %4.1f", fwdSys);
+ Double_t m0900 =
+ AddToStack(stack, 900, isNSD, forward0900, central0900,
+ fwdSys, 5, strangeCorr);
+ Double_t m7000 =
+ AddToStack(stack, 7000, isNSD, forward7000, central7000,
+ fwdSys, 5, strangeCorr);
+
+ stack->SetMaximum(maxFactor * TMath::Max(m0900, m7000));
+
+ f0900->Close();
+ f7000->Close();
+
+ return stack;
+}
+
+// --- Draw one stack ------------------------------------------------
+Double_t
+DrawStack(Bool_t isNSD, Bool_t showClusters,
+ Double_t strangeCorr, Double_t maxFactor=1.1)
+{
+ THStack* stack = MakeStack(isNSD, showClusters, strangeCorr, maxFactor);
+ stack->Draw("nostack ep");
+ stack->GetHistogram()->SetXTitle("#eta");
+ if (!isNSD)
+ stack->GetHistogram()->SetYTitle("#frac{1}{N} #frac{dN_{ch}}{d#eta}");
+ if (!isNSD)
+ stack->SetMinimum(.0001);
+ AdjustAxis(stack->GetHistogram()->GetXaxis());
+ AdjustAxis(stack->GetHistogram()->GetYaxis());
+ stack->GetHistogram()->GetXaxis()->SetTitleOffset(1);
+ gPad->Clear();
+ stack->Draw("nostack ep");
+
+ TLatex* title = new TLatex(1-gPad->GetRightMargin()-.01,
+ 1-gPad->GetTopMargin()-.01,
+ (isNSD ? "NSD" : "INEL"));
+ title->SetNDC();
+ title->SetTextSize(0.1);
+ title->SetTextFont(132);
+ title->SetTextAlign(33);
+ title->SetTextColor(kBlue+2);
+ title->Draw();
+
+ return stack->GetMaximum();
+}
+
+// --- Draw one stack ------------------------------------------------
+void
+ExtractData(Bool_t showClusters = true)
+{
+ gROOT->LoadMacro("OtherData.C");
+ gStyle->SetOptTitle(0);
+ gStyle->SetGridColor(kGray);
+ TCanvas* c = new TCanvas("c", "C", 1200, 850);
+ c->SetRightMargin(0.01);
+ c->SetLeftMargin(0.1);
+ c->SetTopMargin(0.02);
+ c->SetBottomMargin(0.15);
+ c->SetFillColor(0);
+ c->SetBorderSize(0);
+ c->SetBorderMode(0);
+ c->cd();
+
+ c->Divide(1, 2, 0, 0);
+
+ TVirtualPad* p = c->cd(1);
+ p->SetFillColor(0);
+ p->SetRightMargin(0.02);
+ p->SetGridx();
+ p->SetGridy();
+ DrawStack(false, showClusters, 1.5);
+
+ Double_t ms = 1.2;
+
+ TLegend* l = new TLegend(.33, .01, .53, .35);
+ l->SetBorderSize(0);
+ l->SetFillColor(0);
+ l->SetTextFont(132);
+ l->SetFillStyle(0);
+ TLegendEntry* e = l->AddEntry("", "900GeV", "p");
+ e->SetMarkerStyle(20);
+ e->SetMarkerSize(ms+.1);
+ e = l->AddEntry("", " 7TeV", "p");
+ e->SetMarkerStyle(21);
+ e->SetMarkerSize(ms);
+ e = l->AddEntry("", "Mirrored data", "p");
+ e->SetMarkerStyle(24);
+ e->SetMarkerSize(ms);
+ l->Draw();
+
+ TLegend* l2 = new TLegend(.5, .01, .8, .35);
+ l2->SetBorderSize(0);
+ l2->SetFillColor(0);
+ l2->SetTextFont(132);
+ l2->SetFillStyle(0);
+ l2->SetColumnSeperation(-.05);
+ e = l2->AddEntry("", "Forward", "p");
+ e->SetMarkerStyle(20);
+ e->SetMarkerColor(kRed+2);
+ e->SetLineColor(kRed+2);
+ e->SetMarkerSize(ms);
+ if (showClusters) {
+ e = l2->AddEntry("", "Central", "p");
+ e->SetMarkerStyle(20);
+ e->SetMarkerColor(kMagenta+2);
+ e->SetLineColor(kMagenta+2);
+ e->SetMarkerSize(ms);
+ }
+ e = l2->AddEntry("", "#splitline{Eur.Phys.J.#font[22]{C68}:89-108}"
+ "{Eur.Phys.J.#font[22]{C68}:345--354}", "p");
+ e->SetMarkerStyle(20);
+ e->SetMarkerColor(kBlue+2);
+ e->SetLineColor(kBlue+2);
+ e->SetMarkerSize(ms);
+ l2->Draw();
+
+
+ p = c->cd(2);
+ p->SetGridx();
+ p->SetGridy();
+ p->SetFillColor(0);
+ p->SetRightMargin(0.02);
+ DrawStack(true, showClusters, 1.5);
+
+ // TPad, x1, y1, x2, y2, ts, t1, t2, prel
+ gROOT->LoadMacro("AddLogo.C");
+ AddLogo((TPad*)p, .4, p->GetBottomMargin()+.05,
+ .5, .44, 0.08, "", "pp data", true);
+
+ c->cd();
+ TString out("dndeta_pp_forward");
+ if (!showClusters) out.Append("_noclusters");
+ c->SaveAs(Form("%s.png", out.Data()));
+ c->SaveAs(Form("%s.eps", out.Data()));
+}
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff