-
+/**
+ * @file DrawELossPoisson.C
+ * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
+ * @date Thu Jul 7 10:24:58 2011
+ *
+ * @brief A script to draw the Poisson vs Energy Loss correlation
+ *
+ *
+ */
+
+/**
+ * Draw the poisson @f$N_{ch}@f$ estimate against the @f$\Delta@f$
+ * @f$N_{ch}@f$ estimate and do a regression line between the two
+ *
+ * @param p List
+ * @param d Detector
+ * @param r Ring
+ * @param xmin Minimum
+ * @param xmax Maximum
+ *
+ * @return The regression coefficient
+ */
Double_t
DrawRingELossPoisson(TList* p, UShort_t d, Char_t r,
Double_t xmin, Double_t xmax)
return corr->GetCorrelationFactor();
}
-
+/**
+ * Draw the correlation between the Poisson @f$N_{ch}@f$ estimate
+ * and the @f$\Delta@f$ @f$N_{ch}@f$ estimate and do a regression
+ * line between the two for each ring
+ *
+ * @param filename File to read
+ * @param xmax Minimum X
+ * @param xmin Maximum X
+ */
void
DrawELossPoisson(const char* filename="forward.root",
Double_t xmax=-1,
c->cd();
c->SaveAs("elossVsPoisson.png");
}
-
-
-
-
+//
+// EOF
+//
void
LoadLibs()
{
+ gROOT->LoadClass("TVirtualMC", "libVMC");
+ gROOT->LoadClass("AliVEvent", "libSTEERBase");
+ gROOT->LoadClass("AliESDEvent", "libESD");
+ gROOT->LoadClass("AliAnalysisManager", "libANALYSIS");
+ gROOT->LoadClass("AliAnalysisTaskSE", "libANALYSISalice");
+ gROOT->LoadClass("AliAODForwardMult", "libPWG2forward2");
+
+#if 0
const char* test = gSystem->GetLibraries("PWG2forward2","D",false);
if (test && test[0] != '\0') {
// TInterpreter* inter = gROOT->GetInterpreter();
gSystem->Load("libANALYSISalice");
gSystem->Load("libPWG0base");
gSystem->Load("libPWG2forward2");
+#endif
}
//
// EOF
+++ /dev/null
-#include <TMath.h>
-#include <TGraph.h>
-#include <TMultiGraph.h>
-#include <TMarker.h>
-#include <TLine.h>
-#include <TArc.h>
-#include <TCanvas.h>
-#include <TH2F.h>
-#include <THStack.h>
-
-//_____________________________________________________________________
-void AcceptanceCorrection(Char_t r, UShort_t t, Float_t& oldm, Float_t& newm)
-{
- // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
-
- //AliFMDRing fmdring(ring);
- //fmdring.Init();
- // Float_t rad = pars->GetMaxR(r)-pars->GetMinR(r);
- // Float_t slen = pars->GetStripLength(r,t);
- // Float_t sblen = pars->GetBaseStripLength(r,t);
- const Double_t ic1[] = { 4.9895, 15.3560 };
- const Double_t ic2[] = { 1.8007, 17.2000 };
- const Double_t oc1[] = { 4.2231, 26.6638 };
- const Double_t oc2[] = { 1.8357, 27.9500 };
-
- Double_t minR = (r == 'I' ? 4.5213 : 15.4);
- Double_t maxR = (r == 'I' ? 17.2 : 28.0);
- Double_t rad = maxR - minR;
-
- Int_t nStrips = (r == 'I' ? 512 : 256);
- Int_t nSec = (r == 'I' ? 20 : 40);
- Float_t segment = rad / nStrips;
- Float_t basearc = 2 * TMath::Pi() / (0.5 * nSec);
- Float_t radius = minR + t * segment;
- Float_t baselen = 0.5 * basearc * radius;
-
- const Double_t* c1 = (r == 'I' ? ic1 : oc1);
- const Double_t* c2 = (r == 'I' ? ic2 : oc2);
-
- // If the radius of the strip is smaller than the radius corresponding
- // to the first corner we have a full strip length
- Float_t cr = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
- if (radius <= cr) newm = 1;
- else {
- // Next, we should find the end-point of the strip - that is,
- // the coordinates where the line from c1 to c2 intersects a circle
- // with radius given by the strip.
- // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
- Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
- Float_t dx = c2[0] - c1[0];
- Float_t dy = c2[1] - c1[1];
- Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
- Float_t det = radius * radius * dr * dr - D*D;
-
- if (det < 0) newm = 1; // No intersection
- else if (det == 0) newm = 1; // Exactly tangent
- else {
- Float_t x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
- Float_t y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
- Float_t th = TMath::ATan2(x, y);
-
- newm = th / (basearc/2);
- }
- }
-
- Float_t slope = (c1[1] - c2[1]) / (c1[0] - c2[0]);
- Float_t constant = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
-
- Float_t d = (TMath::Power(TMath::Abs(radius*slope),2) +
- TMath::Power(radius,2) - TMath::Power(constant,2));
-
- // If below corners return 1
- Float_t arclen = baselen;
- if (d > 0) {
- Float_t x = ((-TMath::Sqrt(d) - slope * constant) /
- (1+TMath::Power(slope, 2)));
- Float_t y = slope*x + constant;
-
- // If x is larger than corner x or y less than corner y, we have full
- // length strip
- if(x < c1[0] && y> c1[1]) {
- //One sector since theta is by definition half-hybrid
- Float_t theta = TMath::ATan2(x,y);
- arclen = radius * theta;
- }
- }
- // Calculate the area of a strip with no cut
- Float_t basearea1 = 0.5 * baselen * TMath::Power(radius,2);
- Float_t basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
- Float_t basearea = basearea1 - basearea2;
-
- // Calculate the area of a strip with cut
- Float_t area1 = 0.5 * arclen * TMath::Power(radius,2);
- Float_t area2 = 0.5 * arclen * TMath::Power((radius-segment),2);
- Float_t area = area1 - area2;
-
- // Acceptance is ratio
- oldm = area/basearea;
-}
-
-void DrawSolution(Char_t r, UShort_t dt=16, UShort_t offT=128)
-{
- TCanvas* c = new TCanvas(Form("c%c", r), r == 'I' ?
- "Inner" : "Outer", 800, 500);
-
- const Double_t ic1[] = { 4.9895, 15.3560 };
- const Double_t ic2[] = { 1.8007, 17.2000 };
- const Double_t oc1[] = { 4.2231, 26.6638 };
- const Double_t oc2[] = { 1.8357, 27.9500 };
-
- const Double_t* c1 = (r == 'I' ? ic1 : oc1);
- const Double_t* c2 = (r == 'I' ? ic2 : oc2);
- Double_t minR = (r == 'I' ? 4.5213 : 15.4);
- Double_t maxR = (r == 'I' ? 17.2 : 28.0);
- Double_t rad = maxR - minR;
- Float_t cr = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
- Double_t theta = (r == 'I' ? 18 : 9);
- TH2* frame = new TH2F("frame", "Frame", 100, -10, 10, 100, 0, 30);
- frame->Draw();
-
- TLine* line = new TLine(c1[0], c1[1], c2[0], c2[1]);
- line->SetLineColor(kBlue+1);
- line->Draw();
-
- UShort_t nT = (r == 'I' ? 512 : 256);
- for (Int_t t = offT; t < nT; t += dt) {
- Double_t radius = minR + t * rad / nT;
-
- TArc* circle = new TArc(0, 0, radius, 90-theta, 90+theta);
- circle->SetFillColor(0);
- circle->SetFillStyle(0);
- circle->Draw();
-
- // Now find the intersection
- if (radius <= cr) continue;
-
- // Next, we should find the end-point of the strip - that is,
- // the coordinates where the line from c1 to c2 intersects a circle
- // with radius given by the strip.
- // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
- Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
- Float_t dx = c2[0] - c1[0];
- Float_t dy = c2[1] - c1[1];
- Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
- Float_t det = radius * radius * dr * dr - D*D;
-
- if (det < 0) continue;
- if (det == 0) continue;
-
-
- Float_t x = (+D * dy - dx * TMath::Sqrt(det)) / dr / dr;
- Float_t y = (-D * dx - dy * TMath::Sqrt(det)) / dr / dr;
-
- TMarker* m = new TMarker(x, y, 20);
- m->SetMarkerColor(kRed+1);
- m->Draw();
-
- x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
- y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
-
- // Float_t th = TMath::ATan2(x,y);
- // Printf("theta of strip %3d: %f degrees", 180. / TMath::Pi() * th);
-
- m = new TMarker(x, y, 20);
- m->SetMarkerColor(kGreen+1);
- m->Draw();
- }
- c->cd();
-}
-
-
-void TestAcc()
-{
- TCanvas* c = new TCanvas("c", "C");
- TGraph* innerOld = new TGraph(512);
- TGraph* outerOld = new TGraph(256);
- TGraph* innerNew = new TGraph(512);
- TGraph* outerNew = new TGraph(256);
- innerOld->SetName("innerOld");
- innerOld->SetName("Inner type, old");
- innerOld->SetMarkerStyle(21);
- innerOld->SetMarkerColor(kRed+1);
- outerOld->SetName("outerOld");
- outerOld->SetName("Outer type, old");
- outerOld->SetMarkerStyle(21);
- outerOld->SetMarkerColor(kBlue+1);
- innerNew->SetName("innerNew");
- innerNew->SetName("Inner type, new");
- innerNew->SetMarkerStyle(20);
- innerNew->SetMarkerColor(kGreen+1);
- outerNew->SetName("outerNew");
- outerNew->SetName("Outer type, new");
- outerNew->SetMarkerStyle(20);
- outerNew->SetMarkerColor(kMagenta+1);
-
- TMultiGraph* all = new TMultiGraph("all", "Acceptances");
- all->Add(innerOld);
- all->Add(outerOld);
- all->Add(innerNew);
- all->Add(outerNew);
-
- TH2* innerCor = new TH2F("innerCor", "Inner correlation", 100,0,1,100,0,1);
- TH2* outerCor = new TH2F("outerCor", "Outer correlation", 100,0,1,100,0,1);
- innerCor->SetMarkerStyle(20);
- outerCor->SetMarkerStyle(20);
- innerCor->SetMarkerColor(kRed+1);
- outerCor->SetMarkerColor(kBlue+1);
- // innerCor->SetMarkerSize(1.2);
- // outerCor->SetMarkerSize(1.2);
- THStack* stack = new THStack("corr", "Correlations");
- stack->Add(innerCor);
- stack->Add(outerCor);
-
- for (Int_t i = 0; i < 512; i++) {
- Float_t oldAcc, newAcc;
-
- AcceptanceCorrection('I', i, oldAcc, newAcc);
- innerOld->SetPoint(i, i, oldAcc);
- innerNew->SetPoint(i, i, newAcc);
- // Printf("Inner strip # %3d: old=%8.6f, new=%8.6f", i, oldAcc, newAcc);
-
- innerCor->Fill(oldAcc, newAcc);
- if (i >= 256) continue;
-
- AcceptanceCorrection('O', i, oldAcc, newAcc);
- outerOld->SetPoint(i, i, oldAcc);
- outerNew->SetPoint(i, i, newAcc);
- // Printf("Outer strip # %3d: old=%8.6f, new=%8.6f", i, oldAcc, newAcc);
- outerCor->Fill(oldAcc, newAcc);
- }
-
- all->Draw("ap");
-
- DrawSolution('I',4,256);
- DrawSolution('O',4,128);
-
- TCanvas* c2 = new TCanvas("cc", "cc");
- stack->Draw("nostack p");
- TLine* l = new TLine(0,0,1,1);
- l->SetLineStyle(2);
- l->Draw();
-
- c2->cd();
-}
+++ /dev/null
-#include <iomanip>
-#ifndef __CINT__
-#include <iostream>
-#include <TH1.h>
-#include <TH2.h>
-#include <TRandom.h>
-#include <TCanvas.h>
-#include <TArrayD.h>
-#include <TStyle.h>
-// #include <TFitResult.h>
-#include <TGraphErrors.h>
-#include <TF1.h>
-#include <TMath.h>
-#include <THStack.h>
-#include <TList.h>
-#include <TLatex.h>
-#include <TProfile.h>
-#include <TLegend.h>
-#include <TLegendEntry.h>
-#else
-class TH1;
-class TH2;
-class TArrayD;
-class TRandom;
-class TCanvas;
-class TF1;
-#endif
-
-//___________________________________________________________________
-static Double_t landauGaus1(Double_t* xp, Double_t* pp);
-static Double_t landauGausN(Double_t* xp, Double_t* pp);
-static Double_t landauGausI(Double_t* xp, Double_t* pp);
-
-//====================================================================
-struct Function
-{
- /**
- * Enumeration of parameters
- *
- */
- enum {
- /// Constant
- kC,
- /// MPV of landau
- kDelta,
- /// Width of landau
- kXi,
- /// Sigma of gaussian
- kSigma,
- /// Number of particles
- kN,
- /// Weights
- kA
- };
- /// MP shift
- static const Double_t fgkMPShift;
- /// 1 / sqrt(2 * pi)
- static const Double_t fgkInvRoot2Pi;
- /// Number of sigmas to integrate over
- static const Double_t fgkConvNSigma;
- /// Number of steps in integration
- static const Double_t fgkConvNSteps;
- /**
- * Evaluate shifted landay
- * @f[
- * f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
- * @f]
- *
- * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
- * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
- * @f$\Delta=0,\xi=1@f$.
- *
- * @param x Where to evaluate @f$ f'_{L}@f$
- * @param delta Most probable value
- * @param xi The 'width' of the distribution
- *
- * @return @f$ f'_{L}(x;\Delta,\xi) @f$
- */
- static Double_t Landau(Double_t x, Double_t delta, Double_t xi)
- {
- return TMath::Landau(x, delta - xi * fgkMPShift, xi);
- }
- /**
- * Calculate the value of a Landau convolved with a Gaussian
- *
- * @f[
- * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
- * \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
- * \exp{-\frac{(\Delta-\Delta')^2}{2\sigma^2}}
- * @f]
- *
- * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
- * energy loss, @f$ \xi@f$ the width of the Landau, and @f$\sigma@f$
- * is the variance of the Gaussian.
- *
- * Note that this function uses the constants fgkConvNSteps and
- * fgkConvNSigma
- *
- * References:
- * - <a href="dx.doi.org/10.1016/0168-583X(84)90472-5">
- * Nucl.Instrum.Meth.B1:16</a>
- * - <a href="dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
- * - <a href="root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">
- * ROOT implementation</a>
- *
- * @param x where to evaluate @f$ f@f$
- * @param delta @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
- * @param xi @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
- * @param sigma @f$ \sigma@f$ of Gaussian
- *
- * @return @f$ f@f$ evaluated at @f$ x@f$.
- */
- static Double_t LandauGaus(Double_t x, Double_t delta,
- Double_t xi, Double_t sigma)
- {
- Double_t deltap = delta - xi * fgkMPShift;
- Double_t xlow = x - fgkConvNSigma * sigma;
- Double_t xhigh = x + fgkConvNSigma * sigma;
- Double_t step = (xhigh - xlow) / fgkConvNSteps;
- Double_t sum = 0;
-
- for (Int_t i = 0; i < fgkConvNSteps / 2; i++) {
- Double_t x1 = xlow + (i - .5) * step;
- Double_t x2 = xhigh - (i - .5) * step;
-
- Double_t s1 =
- TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma);
- Double_t s2 =
- TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma);
- sum += s1 + s2;
- }
- Double_t ret = step * sum * fgkInvRoot2Pi / sigma;
- return ret;
- }
- /**
- * Evaluate
- * @f[
- * f_i(x;\Delta,\xi,\sigma) = f(x;\Delta_i,\xi_i,\sigma_i)
- * @f]
- * corresponding to @f$ i@f$ particles i.e., with the substitutions
- * @f{eqnarray*}{
- * \Delta \rightarrow \Delta_i &=& i(\Delta + \xi\log(i))\\
- * \xi \rightarrow \xi_i &=& i \xi\\
- * \sigma \rightarrow \sigma_i &=& \sqrt{i}\sigma\\
- * @f}
- *
- * @param x Where to evaluate
- * @param delta @f$ \Delta@f$
- * @param xi @f$ \xi@f$
- * @param sigma @f$ \sigma@f$
- * @param i @f$ i @f$
- *
- * @return @f$ f_i @f$ evaluated *
- */
- static Double_t ILandauGaus(Int_t i, Double_t x, Double_t delta,
- Double_t xi, Double_t sigma)
- {
- if (i <= 1) return LandauGaus(x, delta, xi, sigma);
- Double_t di = i;
- Double_t deltai = i * (delta + xi * TMath::Log(di));
- Double_t xii = i * xi;
- Double_t sigmai = TMath::Sqrt(di) * sigma;
-
- if (sigmai < 1e-10) return Landau(x, deltai, xii);
- Double_t ret = LandauGaus(x, deltai, xii, sigmai);;
- // Info("ILandauGaus", "Fi(%f;%d,%f,%f,%f)->%f",
- // x, i, deltai, xii, sigmai, ret);
- return ret;
- }
- /**
- * Numerically evaluate
- * @f[
- * \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
- * @f]
- * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter. The mapping
- * of the parameters is given by
- *
- * - 0: @f$\Delta@f$
- * - 1: @f$\xi@f$
- * - 2: @f$\sigma@f$
- *
- * This is the partial derivative with respect to the parameter of
- * the response function corresponding to @f$ i@f$ particles i.e.,
- * with the substitutions
- * @f[
- * \Delta \rightarrow \Delta_i = i(\Delta + \xi\log(i))\\
- * \xi \rightarrow \xi_i = i \xi\\
- * \sigma \rightarrow \sigma_i = \sqrt{i}\sigma\\
- * @f]
- *
- * @param x Where to evaluate
- * @param ipar Parameter number
- * @param dp @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
- * @param delta @f$ \Delta@f$
- * @param xi @f$ \xi@f$
- * @param sigma @f$ \sigma@f$
- * @param i @f$ i@f$
- *
- * @return @f$ f_i@f$ evaluated
- */
- static Double_t IdLandauGausdPar(Int_t i, Double_t x,
- UShort_t ipar, Double_t dPar,
- Double_t delta, Double_t xi,
- Double_t sigma)
- {
- if (dPar == 0) return 0;
- Double_t dp = dPar;
- Double_t d2 = dPar / 2;
- Double_t deltaI = i * (delta + xi * TMath::Log(i));
- Double_t xiI = i * xi;
- Double_t si = TMath::Sqrt(Double_t(i));
- Double_t sigmaI = si*sigma;
- Double_t y1 = 0;
- Double_t y2 = 0;
- Double_t y3 = 0;
- Double_t y4 = 0;
- switch (ipar) {
- case 0:
- y1 = ILandauGaus(i, x, deltaI+i*dp, xiI, sigmaI);
- y2 = ILandauGaus(i, x, deltaI+i*d2, xiI, sigmaI);
- y3 = ILandauGaus(i, x, deltaI-i*d2, xiI, sigmaI);
- y4 = ILandauGaus(i, x, deltaI-i*dp, xiI, sigmaI);
- break;
- case 1:
- y1 = ILandauGaus(i, x, deltaI, xiI+i*dp, sigmaI);
- y2 = ILandauGaus(i, x, deltaI, xiI+i*d2, sigmaI);
- y3 = ILandauGaus(i, x, deltaI, xiI-i*d2, sigmaI);
- y4 = ILandauGaus(i, x, deltaI, xiI-i*dp, sigmaI);
- break;
- case 2:
- y1 = ILandauGaus(i, x, deltaI, xiI, sigmaI+si*dp);
- y2 = ILandauGaus(i, x, deltaI, xiI, sigmaI+si*d2);
- y3 = ILandauGaus(i, x, deltaI, xiI, sigmaI-si*d2);
- y4 = ILandauGaus(i, x, deltaI, xiI, sigmaI-si*dp);
- break;
- default:
- return 0;
- }
-
- Double_t d0 = y1 - y4;
- Double_t d1 = 2 * (y2 - y3);
-
- Double_t g = 1/(2*dp) * (4*d1 - d0) / 3;
-
- return g;
- }
- /**
- * Evaluate
- * @f[
- * f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
- * @f]
- *
- * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
- * Landau with a Gaussian (see LandauGaus). Note that
- * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,
- * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.
- *
- * References:
- * - <a href="dx.doi.org/10.1016/0168-583X(84)90472-5">
- * Nucl.Instrum.Meth.B1:16</a>
- * - <a href="dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
- * @param x Where to evaluate @f$ f_N@f$
- * @param delta @f$ \Delta_1@f$
- * @param xi @f$ \xi_1@f$
- * @param sigma @f$ \sigma_1@f$
- * @param n @f$ N@f$ in the sum above.
- * @param a Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
- * @f$ i > 1@f$
- *
- * @return @f$ f_N(x;\Delta,\xi,\sigma)@f$
- */
- static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi,
- Double_t sigma, Int_t n, Double_t* a)
- {
- Double_t res = LandauGaus(x, delta, xi, sigma);
- for (Int_t i = 2; i <= n; i++)
- res += a[i-2] * ILandauGaus(i, x, delta, xi, sigma);
- return res;
- }
- /**
- * Calculate the the estimate of number of particle contribtutions
- * at energy loss @a x
- * @f[
- * E(\Delta;\Delta_p,\xi,\sigma,\mathbf{a}) =
- * \frac{\sum_{i=1}^N i a_i f_i(x;\Delta_p,\xi,\sigma)}{
- * \sum_{i=1}^N a_i f_i(x;\Delta_p,\xi,\sigma)}
- * @f]
- * where @f$a_1=1@f$
- *
- * @param x Energy loss @f$\Delta@f$
- * @param delta @f$\Delta_{p}@f$
- * @param xi @f$\xi@f$
- * @param sigma @f$\sigma@f$
- * @param n Maximum number of particles @f$N@f$
- * @param a Weights @f$a_i@f$ for @f$i=2,...,N@f$
- *
- * @return @f$E(\Delta;\Delta_p,\xi,\sigma,\mathbf{a})@f$
- */
- Double_t NEstimate(Double_t x, Double_t delta, Double_t xi,
- Double_t sigma, Int_t n, Double_t* a)
- {
- Double_t num = LandauGaus(x, delta, xi, sigma);
- Double_t den = num;
- for (Int_t i = 2; i <= n; i++) {
- Double_t f = ILandauGaus(i, x, delta, xi, sigma);
- num += a[i-2] * f * i;
- den += a[i-2] * f;
- }
- if (den < 1e-4) return 0;
- return num / den;
- }
- /**
- * Calculate the partial derivative of the function
- * @f$E(\Delta;\Delta_p,\xi,\sigma,\mathbf{a})@f$ with respect to
- * one of the parameters @f$\Delta_{p}@f$, @f$\xi@f$, @f$\sigma@f$,
- * or @f$a_i@f$. Note for that
- * @f{eqnarray*}
- * \frac{\partial E}{\partial C} & \equiv & 0\\
- * \frac{\partial E}{\partial a_1} & \equiv & 0\\
- * \frac{\partial E}{\partial N} & \equiv & 0\\
- * @f{eqnarray*}
- *
- * @param x Where to evaluate the derivative
- * @param ipar Parameter to differentiate relative to
- * @param dPar Variation in parameter to use in calculation
- * @param delta @f$\Delta_{p}@f$
- * @param xi @f$\xi@f$
- * @param sigma @f$\sigma@f$
- * @param n @f$N@f$
- * @param a @f$\mathbf{a} = (a_2,...,a_N)@f$
- *
- * @return @f$\partial E/\partial p|_{x}@f$
- */
- Double_t DEstimatePar(Double_t x, Int_t ipar, Double_t delta, Double_t xi,
- Double_t sigma, UShort_t n, Double_t* a,
- Double_t dPar=0.001)
- {
- Double_t dp = dPar;
- Double_t d2 = dPar / 2;
- Double_t y1 = 0;
- Double_t y2 = 0;
- Double_t y3 = 0;
- Double_t y4 = 0;
- switch (ipar) {
- case Function::kC:
- case Function::kN:
- return 0;
- case Function::kDelta:
- y1 = NEstimate(x, delta+dp, xi, sigma, n, a);
- y2 = NEstimate(x, delta+d2, xi, sigma, n, a);
- y3 = NEstimate(x, delta-d2, xi, sigma, n, a);
- y4 = NEstimate(x, delta-dp, xi, sigma, n, a);
- break;
- case Function::kXi:
- y1 = NEstimate(x, delta, xi+dp, sigma, n, a);
- y2 = NEstimate(x, delta, xi+d2, sigma, n, a);
- y3 = NEstimate(x, delta, xi-d2, sigma, n, a);
- y4 = NEstimate(x, delta, xi-dp, sigma, n, a);
- break;
- case Function::kSigma:
- y1 = NEstimate(x, delta, xi, sigma+dp, n, a);
- y2 = NEstimate(x, delta, xi, sigma+d2, n, a);
- y3 = NEstimate(x, delta, xi, sigma-d2, n, a);
- y4 = NEstimate(x, delta, xi, sigma-dp, n, a);
- break;
- default: {
- Int_t j = ipar-kA;
- if (j+1 > n) return 0;
- Double_t aa = a[j];
- a[j] = aa+dp; y1 = NEstimate(x, delta, xi, sigma, n, a);
- a[j] = aa+d2; y2 = NEstimate(x, delta, xi, sigma, n, a);
- a[j] = aa-d2; y3 = NEstimate(x, delta, xi, sigma, n, a);
- a[j] = aa-dp; y4 = NEstimate(x, delta, xi, sigma, n, a);
- a[j] = aa;
- }
- break;
- }
- Double_t d0 = y1 - y4;
- Double_t d1 = 2 * (y2 - y3);
-
- Double_t g = 1/(2*dp) * (4*d1 - d0) / 3;
- return g;
- }
- /**
- * Generate a TF1 object of @f$ f_I@f$
- *
- * @param n @f$ n@f$ - the number of particles
- * @param c Constant
- * @param delta @f$ \Delta@f$
- * @param xi @f$ \xi_1@f$
- * @param sigma @f$ \sigma_1@f$
- * @param a Array of @f$a_i@f$ for @f$i > 1@f$
- * @param xmin Least value of range
- * @param xmax Largest value of range
- *
- * @return Newly allocated TF1 object
- */
- static TF1* MakeFunc(Int_t n, Double_t c, Double_t delta,
- Double_t xi, Double_t sigma, Double_t* a,
- Double_t xmin, Double_t xmax)
- {
- Int_t nPar = kN + n;
- TF1* f = new TF1(Form("landGaus%d",n), &landauGausN, xmin, xmax, nPar);
- f->SetNpx(500);
- f->SetParName(kC, "C");
- f->SetParName(kDelta, "#Delta_{p}");
- f->SetParName(kXi, "#xi");
- f->SetParName(kSigma, "#sigma");
- f->SetParName(kN, "N");
- f->SetParameter(kC, c);
- f->SetParameter(kDelta, delta);
- f->SetParameter(kXi, xi);
- f->SetParameter(kSigma, sigma);
- f->FixParameter(kN, n);
- f->SetParLimits(kDelta, xmin, 4);
- f->SetParLimits(kXi, 0.00, 4);
- f->SetParLimits(kSigma, 0.01, 0.11);
-
- for (Int_t i = 2; i <= n; i++) {
- Int_t j = i-2;
- f->SetParName(kA+j, Form("a_{%d}", i));
- f->SetParameter(kA+j, a[j]);
- f->SetParLimits(kA+j, 0, 1);
- }
- return f;
- }
- /**
- * Make a ROOT TF1 function object corresponding to a single
- * component for @f$ i@f$ particles
- *
- * @param i Number of particles
- * @param c @f$C@f$
- * @param delta @f$ \Delta@f$
- * @param xi @f$ \xi_1@f$
- * @param sigma @f$ \sigma_1@f$
- * @param xmin Minimum of range of @f$\Delta@f$
- * @param xmax Maximum of range of @f$\Delta@f$
- *
- * @return Pointer to newly allocated ROOT TF1 object
- */
- static TF1* MakeIFunc(Int_t i, Double_t c, Double_t delta,
- Double_t xi, Double_t sigma,
- Double_t xmin, Double_t xmax)
- {
- Int_t nPar = 5;
- TF1* f = new TF1(Form("landGausI%d",i), &landauGausI, xmin, xmax, nPar);
- f->SetNpx(100);
- f->SetParName(kC, "C");
- f->SetParName(kDelta, "#Delta_{p}");
- f->SetParName(kXi, "#xi");
- f->SetParName(kSigma, "#sigma");
- f->SetParName(kN, "N");
- f->SetParameter(kC, c);
- f->SetParameter(kDelta, delta);
- f->SetParameter(kXi, xi);
- f->SetParameter(kSigma, sigma);
- f->FixParameter(kN, i);
- return f;
- }
-
- // --- Object code -------------------------------------------------
- /**
- * Constructor
- *
- * @param n
- * @param c
- * @param delta
- * @param xi
- * @param sigma
- * @param a
- * @param xmin
- * @param xmax
- *
- * @return
- */
- Function(Int_t n, Double_t c, Double_t delta, Double_t xi, Double_t sigma,
- Double_t* a, Double_t xmin, Double_t xmax)
- : fF(MakeFunc(n, c, delta, xi, sigma, a, xmin, xmax))
- {}
- Function(TF1* f) : fF(f) {}
- /**
- * Get pointer to ROOT TF1 object
- *
- * @return Pointer to TF1 object
- */
- TF1* GetF1() { return fF; }
- /**
- * Evaluate the function at @a x
- * @f[
- * f_N(x;\Delta,\xi,\sigma) =
- * \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i)
- * @f]
- *
- * @param x Where to evaluate
- *
- * @return
- */
- Double_t Evaluate(Double_t x) const
- {
- return fF->Eval(x);
- }
- /**
- * Evaluate the function at @a x
- * @f[
- * f_N(x;\Delta,\xi,\sigma) =
- * \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i)
- * @f]
- *
- * This also calculates the error on the point like
- * @f[
- * \delta^2 f_N =
- * \left(\frac{\partial f_N}{\partial\Delta_p}\right)^2\delta^2\Delta_p+
- * \left(\frac{\partial f_N}{\partial\xi }\right)^2\delta^2\xi+
- * \left(\frac{\partial f_N}{\partial\sigma }\right)^2\delta^2\sigma+
- * \sum_{i=2}^N\left(\frac{\partial f_N}{\partial a_i}\right)^2\delta^2a_i
- * @f]
- *
- * @param x Where to evaluate
- * @param e On return, contains the error @f$\delta f_N@f$ on the
- * function value
- *
- * @return @f$ f_N(x;\Delta,\xi,\sigma)@f$
- */
- Double_t Evaluate(Double_t x, Double_t& e) const
- {
- Double_t delta = GetDelta();
- Double_t xi = GetXi();
- Double_t sigma = GetSigma();
- Double_t dFdDelta2 = 0;
- Double_t dFdXi2 = 0;
- Double_t dFdSigma2 = 0;
- e = 0;
- for (Int_t i = 1; i <= GetN(); i++) {
- Double_t a = GetA(i);
- Double_t dDelta = a*IdLandauGausdPar(i, x, 1, 0.001, delta, xi, sigma);
- Double_t dXi = a*IdLandauGausdPar(i, x, 2, 0.001, delta, xi, sigma);
- Double_t dSigma = a*IdLandauGausdPar(i, x, 3, 0.001, delta, xi, sigma);
- Double_t dFda = a*ILandauGaus(i, x, delta, xi, sigma);
- dFdDelta2 += dDelta * dDelta;
- dFdXi2 += dXi * dXi;
- dFdSigma2 += dSigma * dSigma;
- e += TMath::Power(dFda * GetEA(i),2);
- }
- Double_t edelta = GetEDelta();
- Double_t exi = GetEXi();
- Double_t esigma = GetESigma();
- e += (dFdDelta2 * edelta * edelta +
- dFdXi2 * exi * exi +
- dFdSigma2 * esigma * esigma);
- e = TMath::Sqrt(e);
- return fF->Eval(x);
- }
-
- /**
- * Evaluate
- * @f[
- * E(x;\Delta,\xi,\sigma,\mathbf{a}) =
- * \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma)}{
- * f_N(x;\Delta,\xi,\sigma)} =
- * \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i)}{
- * \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i)}
- * @f]
- *
- * If the denominator is less than @f$10^{-4}@f$, then 0 is returned.
- *
- * @param x Where to evaluate @f$ E@f$
- * @param maxN Maximum number of weights
- *
- * @return @f$ E(x;\Delta,\xi,\sigma,\mathbf{a})@f$
- */
- Double_t EstimateNParticles(Double_t x, Short_t maxN=-1)
- {
- UShort_t n = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
- Double_t delta = GetDelta();
- Double_t xi = GetXi();
- Double_t sigma = GetSigma();
- return NEstimate(x, delta, xi, sigma, n, GetAs());
- }
- /**
- * Evaluate
- * @f[
- * E(x;\Delta,\xi,\sigma,\mathbf{a}) =
- * \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma)}{
- * f_N(x;\Delta,\xi,\sigma,\mathbf{a})} =
- * \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i)}{
- * \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i)}
- * @f]
- *
- * If @f$f_N(x;\Delta,\xi,\sigma,\mathbf{a})<10^{-4}@f$ then 0 is
- * returned.
- *
- * This also calculatues the error @f$\delta E@f$ of the
- * function value at @f$x@f$:
- *
- * @f[
- * \delta^2 E =
- * \left(\frac{\partial E}{\partial\Delta_p}\right)^2\delta^2\Delta_p+
- * \left(\frac{\partial E}{\partial\xi }\right)^2\delta^2\xi+
- * \left(\frac{\partial E}{\partial\sigma }\right)^2\delta^2\sigma+
- * \sum_{i=2}^N\left(\frac{\partial E}{\partial a_i}\right)^2\delta^2a_i
- * @f]
- * The partial derivatives are evaluated numerically.
- *
- * @param x Where to evaluate @f$ E@f$
- * @param e On return, @f$\delta E|_{x}@f$
- * @param maxN Maximum number of weights
- *
- * @return @f$ E(x;\Delta,\xi,\sigma,\mathbf{a})@f$
- */
- Double_t EstimateNParticles(Double_t x, Double_t& e, Short_t maxN=-1)
- {
-
- UShort_t n = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
- Double_t delta = GetDelta();
- Double_t xi = GetXi();
- Double_t sigma = GetSigma();
- Double_t* a = GetAs();
- Double_t dDelta = (DEstimatePar(x,Function::kDelta,delta,xi,sigma,n,a,\
- 0.01 * GetEDelta()) * GetEDelta());
- Double_t dXi = (DEstimatePar(x,Function::kXi, delta,xi,sigma,n,a,
- 0.01 * GetEXi()) * GetEXi());
- Double_t dSigma = (DEstimatePar(x,Function::kSigma,delta,xi,sigma,n,a,
- 0.01 * GetESigma()) * GetESigma());
- e = dDelta * dDelta + dXi * dXi + dSigma * dSigma;
- for (Int_t i = 2; i <= n; i++) {
- Double_t dAi = (DEstimatePar(x, Function::kA+i-2,delta,xi,sigma,n,a,
- 0.01 * GetEA(i)) * GetEA(i));
- e += dAi * dAi;
- }
- e = TMath::Sqrt(e);
- return NEstimate(x, GetDelta(), GetXi(), GetSigma(), n, GetAs());
- }
- /**
- * Estimate the number of particles by calculating the weights
- * @f[
- * w_i = a_i f_i(x;\Delta_p,\xi,\sigma)
- * @f]
- *
- * and then draw a random number @f$r@f$ in the range @f$[0,\sum_i^N
- * w_i]@f$ and return the @f$i@f$ for which @f$w_{i-1} < @r \leq
- * w_{i}@f$ (notem @f$w_{-1}=0@f$).
- *
- * @param x Where to evaluate the component functions
- * @param maxN Maximum weight to use
- *
- * @return Estimate of the number of particles using the procedure
- * outlined above.
- */
- Double_t RandomEstimateNParticles(Double_t x, Short_t maxN=-1)
- {
- UShort_t n = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
- TArrayD p(n);
- Double_t delta = GetDelta();
- Double_t xi = GetXi();
- Double_t sigma = GetSigma();
-
- if (Evaluate(x) < 1e-4) return 0;
-
- // std::cout << "RandomEstimateNParticles(" << x << "):";
- for (Int_t i = 1; i <= n; i++) {
- Double_t a = GetA(i);
- Double_t f = ILandauGaus(i, x, delta, xi, sigma);
- p[i-1] = a * f;
- // std::cout << p[i-1] << ",";
- }
- Double_t r = gRandom->Uniform(p.GetSum());
- Double_t l = 0;
- // std::cout << "s=" << p.GetSum() << ",r=" << r << ",";
- for (Int_t i = 1; i <= n; i++) {
- if (r > l && r <= l + p[i-1]) {
- // std::cout << "l=" << l << ",l+p=" << l+p[i-1] << "-> " << i
- // << std::endl;
- return i;
- }
- l += p[i-1];
- }
- return 0;
- }
- /**
- * Draw the function
- *
- * @param option Drawing option
- */
- void Draw(Option_t* option="") { fF->Draw(option); }
- /** @return @f$C@f$ */
- Double_t GetC() const { return fF->GetParameter(kC); }
- /** @return @f$\Delta_{mp}@f$ */
- Double_t GetDelta() const { return fF->GetParameter(kDelta); }
- /** @return @f$\xi@f$ */
- Double_t GetXi() const { return fF->GetParameter(kXi); }
- /** @return @f$\sigma@f$ */
- Double_t GetSigma() const { return fF->GetParameter(kSigma); }
- /** @return @f$N@f$ */
- UShort_t GetN() const { return fF->GetParameter(kN); }
- /** @return @f$a_{i}@f$ */
- Double_t GetA(Int_t i) const { return i == 1 ? 1 : fF->GetParameter(kA+i-2);}
- /** @return @f$a_2,...,a_N@f$ */
- Double_t* GetAs() const { return &(fF->GetParameters()[kA]);}
- /** @return @f$\Delta C@f$ */
- Double_t GetEC() const { return fF->GetParError(kC); }
- /** @return @f$\Delta \Delta_{mp}@f$ */
- Double_t GetEDelta() const { return fF->GetParError(kDelta); }
- /** @return @f$\Delta \xi@f$ */
- Double_t GetEXi() const { return fF->GetParError(kXi); }
- /** @return @f$\Delta \sigma@f$ */
- Double_t GetESigma() const { return fF->GetParError(kSigma); }
- /** @return @f$\Delta N@f$ */
- UShort_t GetEN() const { return 0; }
- /** @return @f$\Delta a_{i}@f$ */
- Double_t GetEA(Int_t i) const { return i == 1 ? 0 : fF->GetParError(kA+i-2);}
- /** @return @f$\Delta a_2,...,\Delta a_N@f$ */
- Double_t* GetEAs() const { return &(fF->GetParErrors()[kA]);}
- /** The ROOT object */
- TF1* fF;
-};
-#ifndef __CINT__
-const Double_t Function::fgkMPShift = -0.22278298;
-const Double_t Function::fgkInvRoot2Pi = 1. / TMath::Sqrt(2*TMath::Pi());
-const Double_t Function::fgkConvNSigma = 5;
-const Double_t Function::fgkConvNSteps = 100;
-#endif
-
-//====================================================================
-struct Fitter
-{
- // --- Object code ------------------------------------------------
- Fitter(Int_t nMinus)
- : fNMinus(nMinus)
- {
- }
- /**
- *
- *
- * @param dist
- *
- * @return
- */
- TF1* FitFirst(TH1* dist)
- {
- fFunctions.Clear();
-
- // Get upper edge
- Int_t midBin = dist->GetMaximumBin();
- Double_t midE = dist->GetBinLowEdge(midBin);
-
- // Get low edge
- Int_t minBin = midBin - fNMinus;
- Double_t minE = dist->GetBinCenter(minBin);
- Double_t maxE = dist->GetBinCenter(midBin+2*fNMinus);
- Double_t a[] = { 0 };
- TF1* f = Function::MakeFunc(1, 0.5, midE, 0.07, 0.1, a, minE, maxE);
-
- // Do the fit
- dist->Fit(f, "RNQS", "", minE, maxE);
- fFunctions.AddAtAndExpand(f, 0);
- return f;
- }
- /**
- *
- *
- * @param dist
- * @param n
- *
- * @return
- */
- Function* FitN(TH1* dist, Int_t n)
- {
- if (n == 1) return new Function(FitFirst(dist));
- TF1* f1 = static_cast<TF1*>(fFunctions.At(0));
- if (!f1) {
- f1 = FitFirst(dist);
- if (!f1) {
- ::Warning("FitN", "First fit missing");
- return 0;
- }
- }
-
- Double_t delta1 = f1->GetParameter(Function::kDelta);
- Double_t xi1 = f1->GetParameter(Function::kXi);
- Double_t maxEi = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
- Double_t minE = f1->GetXmin();
-
- TArrayD a(n-1);
- for (UShort_t i = 2; i <= n; i++) a.fArray[i-2] = (n==2? 0.05 : 0.000001);
-
- TF1* fn = Function::MakeFunc(n, f1->GetParameter(Function::kC),
- delta1, xi1,
- f1->GetParameter(Function::kSigma),
- a.fArray, minE, maxEi);
-
- dist->Fit(fn, "RSNQ", "", minE, maxEi);
- dist->Fit(fn, "RSNQM", "", minE, maxEi);
- fFunctions.AddAtAndExpand(fn, n-1);
-
- return new Function(fn);
- }
- Int_t fNMinus;
- TObjArray fFunctions;
-};
-
-//====================================================================
-/**
- * Utility function
- *
- * @param xp Pointer to independent variables
- * @param pp Pointer to parameters
- *
- * @return Function evaluated at xp[0]
- */
-static Double_t landauGaus1(Double_t* xp, Double_t* pp)
-{
- Double_t x = xp[0];
- Double_t c = pp[Function::kC];
- Double_t delta = pp[Function::kDelta];
- Double_t xi = pp[Function::kXi];
- Double_t sigma = pp[Function::kSigma];
-
- return c * Function::LandauGaus(x, delta, xi, sigma);
-}
-/**
- * Utility function
- *
- * @param xp Pointer to independent variables
- * @param pp Pointer to parameters
- *
- * @return Function evaluated at xp[0]
- */
-static Double_t landauGausN(Double_t* xp, Double_t* pp)
-{
- Double_t x = xp[0];
- Double_t c = pp[Function::kC];
- Double_t delta = pp[Function::kDelta];
- Double_t xi = pp[Function::kXi];
- Double_t sigma = pp[Function::kSigma];
- Int_t n = Int_t(pp[Function::kN]);
- Double_t* a = &pp[Function::kA];
-
- return c * Function::NLandauGaus(x, delta, xi, sigma, n, a);
-}
-/**
- * Utility function
- *
- * @param xp Pointer to independent variables
- * @param pp Pointer to parameters
- *
- * @return Function evaluated at xp[0]
- */
-static Double_t landauGausI(Double_t* xp, Double_t* pp)
-{
- Double_t x = xp[0];
- Double_t c = pp[Function::kC];
- Double_t delta = pp[Function::kDelta];
- Double_t xi = pp[Function::kXi];
- Double_t sigma = pp[Function::kSigma];
- Int_t i = Int_t(pp[Function::kN]);
-
- return c * Function::ILandauGaus(i, x, delta, xi, sigma);
-}
-
-//====================================================================
-/**
- *
- *
- */
-struct Generator
-{
- Double_t fDelta;
- Double_t fXi;
- Double_t fSigma;
- Int_t fMaxN;
- TH1* fSingle;
- TH1* fSummed;
- TH1* fSimple;
- TH2* fNSum;
- TH1* fN;
- TArrayD fRP; // Raw probabilities
- TArrayD fP; // Normalized probabilities
- TF1* fF;
- Int_t fNObs;
-
- /**
- *
- *
- *
- * @return
- */
- Generator(Double_t delta=.55, Double_t xi=0.06, Double_t sigma=0.02,
- Int_t maxN=5, const Double_t* a=0)
- : fDelta(delta), fXi(xi), fSigma(sigma), fMaxN(maxN),
- fSingle(0), fSummed(0), fSimple(0), fNSum(0), fRP(), fP(fMaxN), fF(0),
- fNObs(0)
- {
- fMaxN = TMath::Min(fMaxN, 5);
- const Double_t ps[] = { 1, .05, .02, .005, .002 };
- fRP.Set(5, ps);
- if (a) {
- for (Int_t i = 0; i < maxN; i++) fRP[i] = a[i];
- }
-
- Double_t sum = 0;
- for (Int_t i = 0; i < fMaxN; i++) sum += fRP[i];
- for (Int_t i = 0; i < fP.fN; i++) {
- fP[i] = fRP[i] / sum;
- if (i != 0) fP[i] += fP[i-1];
- }
-
- fSingle = new TH1D("single", "Single signals", 500, 0, 10);
- fSingle->SetXTitle("Signal");
- fSingle->SetYTitle("Events");
- fSingle->SetFillColor(kRed+1);
- fSingle->SetMarkerColor(kRed+1);
- fSingle->SetLineColor(kRed+1);
- fSingle->SetFillStyle(3001);
- fSingle->SetMarkerStyle(21);
- fSingle->SetDirectory(0);
- fSingle->Sumw2();
-
- fSummed = static_cast<TH1D*>(fSingle->Clone("summed"));
- fSummed->SetTitle("Summed signals");
- fSummed->SetFillColor(kBlue+1);
- fSummed->SetMarkerColor(kBlue+1);
- fSummed->SetLineColor(kBlue+1);
- fSummed->SetDirectory(0);
-
- fSimple = static_cast<TH1D*>(fSingle->Clone("summed"));
- fSimple->SetTitle("Summed signals");
- fSimple->SetFillColor(kGreen+1);
- fSimple->SetMarkerColor(kGreen+1);
- fSimple->SetLineColor(kGreen+1);
- fSimple->SetDirectory(0);
-
- fNSum = new TH2D("single", "Observation signals",
- fMaxN, .5, fMaxN+.5, 100, 0, 10);
- fNSum->SetMarkerColor(kMagenta+1);
- fNSum->SetLineColor(kMagenta+1);
- fNSum->SetFillColor(kMagenta+1);
- fNSum->SetYTitle("Signal");
- fNSum->SetXTitle("N");
- fNSum->SetZTitle("Events");
- fNSum->SetDirectory(0);
- fNSum->Sumw2();
-
- fN = new TH1D("n", "Number of particles", fMaxN, .5, fMaxN+.5);
- fN->SetXTitle("N");
- fN->SetYTitle("Events");
- fN->SetMarkerColor(kMagenta+1);
- fN->SetLineColor(kMagenta+1);
- fN->SetFillColor(kMagenta+1);
- fN->SetMarkerStyle(22);
- fN->SetFillStyle(3001);
- fN->SetDirectory(0);
-
- std::cout << "Probabilities:" << std::setprecision(4)
- << "\n ";
- for (Int_t i = 0; i < fRP.fN; i++)
- std::cout << std::setw(7) << fRP[i] << " ";
- std::cout << " = " << fRP.GetSum() << "\n ";
- for (Int_t i = 0; i < fP.fN; i++)
- std::cout << std::setw(7) << fP[i] << " ";
- std::cout << std::endl;
-
- Double_t aa[] = { 0 };
- fF = Function::MakeFunc(1, 1, fDelta, fXi, fSigma, aa, 0, 10);
- }
- /**
- *
- *
- */
- void Clear()
- {
- fSingle->Reset();
- fSummed->Reset();
- fNSum->Reset();
- fN->Reset();
- }
- /**
- *
- *
- *
- * @return
- */
- Int_t GenerateN()
- {
- Double_t rndm = gRandom->Uniform();
- Int_t ret = 0;
- if ( rndm < fP[0]) ret = 1;
- else if (rndm >= fP[0] && rndm < fP[1]) ret = 2;
- else if (rndm >= fP[1] && rndm < fP[2]) ret = 3;
- else if (rndm >= fP[2] && rndm < fP[3]) ret = 4;
- else if (rndm >= fP[3] && rndm < fP[4]) ret = 5;
- // Printf("%f -> %d", rndm, ret);
- fN->Fill(ret);
- return ret;
- }
- /**
- *
- *
- * @param mpv
- * @param xi
- * @param sigma
- *
- * @return
- */
- Double_t Generate1Signal()
- {
- Double_t ret = 0;
- if (fF) ret = fF->GetRandom();
- else {
- Double_t rmpv = gRandom->Gaus(fDelta, fSigma);
- ret = gRandom->Landau(rmpv, fXi);
- }
- fSingle->Fill(ret);
- fSimple->Fill(gRandom->Landau(fDelta - fXi * Function::fgkMPShift, fXi));
- return ret;
- }
- /**
- *
- *
- * @param n
- * @param mpv
- * @param xi
- * @param sigma
- *
- * @return
- */
- Double_t GenerateNSignal(Int_t n)
- {
- Double_t ret = 0;
- for (Int_t i = 0; i < n; i++) ret += Generate1Signal();
- fSummed->Fill(ret);
- fNSum->Fill(n, ret);
- return ret;
- }
- /**
- *
- *
- * @param nObs
- * @param reset
- */
- void Generate(Int_t nObs=1000, Bool_t reset=true)
- {
- if (reset) Clear();
- fNObs = nObs;
-
- for (Int_t i = 0; i < nObs; i++) {
- // if (((i+1) % (nObs/10)) == 0)
- // std::cout << "Event " << std::setw(6) << i << std::endl;
- Int_t n = GenerateN();
- GenerateNSignal(n);
- }
- Double_t m = fSingle->GetMaximum();
- fSingle->Scale(1. / m);
- m = fSummed->GetMaximum();
- fSummed->Scale(1. / m);
- m = fSimple->GetMaximum();
- fSimple->Scale(1. / m);
-
- std::cout << "Resulting probabilities:\n ";
- for (Int_t i = 1; i <= fN->GetNbinsX(); i++)
- std::cout << std::setprecision(4) << std::setw(7)
- << fN->GetBinContent(i) << " ";
- std::cout << std::endl;
- }
- /**
- * Pretty print of parameters.
- *
- * @param input
- * @param f
- * @param i
- */
- void PrintParameter(Double_t input, TF1* f, Int_t i)
- {
- Double_t val = 0, err = 0;
- if (i < f->GetNpar()) {
- val = f->GetParameter(i);
- err = f->GetParError(i);
- }
- std::cout << " " << std::setw(16) << f->GetParName(i) << ": "
- << std::fixed
- << std::setprecision(4) << std::setw(6) << input << " -> "
- << std::setprecision(4) << std::setw(6) << val << " +/- "
- << std::setprecision(5) << std::setw(7) << err << " [delta:"
- << std::setw(3) << int(100 * (input-val) / input) << "%,err:"
- << std::setw(3) << int(100 * err / val) << "%]"
- << std::endl;
- }
- /**
- * Print the result
- *
- * @param f
- */
- void PrintResult(TF1* f)
- {
- Double_t chi2 = f->GetChisquare();
- Int_t ndf = f->GetNDF();
- std::cout << "Chi^2/NDF = " << chi2 << "/" << ndf << "="
- << chi2/ndf << std::endl;
- PrintParameter(1, f, Function::kC);
- PrintParameter(fDelta, f, Function::kDelta);
- PrintParameter(fXi, f, Function::kXi);
- PrintParameter(fSigma, f, Function::kSigma);
- PrintParameter(fMaxN, f, Function::kN);
- for (Int_t i = 0; i < fMaxN-1; i++)
- PrintParameter(fRP[i+1], f, Function::kA+i);
- }
- /**
- *
- *
- * @param x
- * @param y
- * @param intput
- * @param f
- * @param i
- */
- void DrawParameter(Double_t x, Double_t y, Double_t input, TF1* f, Int_t i)
- {
- Double_t val = 0, err = 0;
- if (i < f->GetNpar()) {
- val = f->GetParameter(i);
- err = f->GetParError(i);
- }
-
- TLatex* ltx = new TLatex(x, y, f->GetParName(i));
- ltx->SetTextSize(.04);
- ltx->SetTextFont(132);
- ltx->SetTextAlign(11);
- ltx->SetNDC();
- ltx->Draw();
-
- ltx->DrawLatex(x+.05, y, Form("%5.3f", input));
- ltx->DrawLatex(x+.14, y, "#rightarrow");
- ltx->SetTextAlign(21);
- ltx->DrawLatex(x+.3, y, Form("%5.3f #pm %6.4f", val, err));
- ltx->SetTextAlign(11);
- ltx->DrawLatex(x+.48, y, Form("[#Delta=%3d%%, #delta=%3d%%]",
- int(100 * (input-val)/input),
- int(100*err/val)));
- }
- /**
- *
- *
- * @param x
- * @param y
- * @param f
- */
- void DrawResult(Double_t x, Double_t y, TF1* f)
- {
- Double_t chi2 = f->GetChisquare();
- Int_t ndf = f->GetNDF();
- TLatex* ltx = new TLatex(x, y, Form("#chi^{2}/NDF=%5.2f/%3d=%6.3f",
- chi2, ndf, chi2/ndf));
- ltx->SetTextSize(.04);
- ltx->SetNDC();
- ltx->SetTextFont(132);
- ltx->Draw();
-
- x += .05;
- y -= .035; DrawParameter(x, y, 1, f, Function::kC);
- y -= .035; DrawParameter(x, y, fDelta, f, Function::kDelta);
- y -= .035; DrawParameter(x, y, fXi, f, Function::kXi);
- y -= .035; DrawParameter(x, y, fSigma, f, Function::kSigma);
- y -= .035; DrawParameter(x, y, fMaxN, f, Function::kN);
- for (Int_t i = 0; i < fMaxN-1; i++) {
- y -= .035; DrawParameter(x, y, fRP[i+1], f, Function::kA+i);
- }
- }
- /**
- *
- *
- */
- void Draw(const char* type="png")
- {
- gStyle->SetOptTitle(0);
- gStyle->SetOptStat(0);
- gStyle->SetPalette(1);
-
- TCanvas* c = new TCanvas("c", "Generator distributions", 1200, 1000);
- c->SetFillColor(0);
- c->SetBorderSize(0);
- c->SetBorderMode(0);
- c->SetTopMargin(0.01);
- c->SetRightMargin(0.01);
-
- c->Divide(2, 2);
-
- // --- First pad: Data (multi, single, landau) + Fit -------------
- TVirtualPad* p = c->cd(1);
- p->SetLogy();
- p->SetBorderSize(0);
- p->SetBorderMode(0);
- p->SetTopMargin(0.01);
- p->SetRightMargin(0.01);
- fSummed->Draw();
- fSingle->Draw("same");
- fSimple->Draw("same");
-
- if (fF) {
- Double_t max = fF->GetMaximum();
- Double_t a[] = { 0 };
- TF1* fcopy = Function::MakeFunc(1,
- fF->GetParameter(Function::kC) /max,
- fF->GetParameter(Function::kDelta),
- fF->GetParameter(Function::kXi),
- fF->GetParameter(Function::kSigma),
- a, 0, 10);
- fcopy->SetLineColor(2);
- fcopy->SetLineStyle(2);
- fcopy->Draw("same");
- }
-
- Fitter* fitter = new Fitter(4);
- Function* fit = fitter->FitN(fSummed, 5);
- TF1* f = fit->GetF1();
- f->Draw("same");
- f->SetLineColor(kBlack);
- TF1* fst = static_cast<TF1*>(fitter->fFunctions.At(0));
- fst->SetLineWidth(3);
- fst->SetLineStyle(2);
- fst->SetLineColor(kMagenta);
- fst->Draw("same");
- PrintResult(f);
- DrawResult(.25, .95, f);
- TGraphErrors* gr = new TGraphErrors(fSummed->GetNbinsX());
- gr->SetName("fitError");
- gr->SetTitle("Fit w/errors");
- gr->SetFillColor(kYellow+1);
- gr->SetFillStyle(3001);
- for (Int_t j = 1; j <= fSummed->GetNbinsX(); j++) {
- Double_t x = fSummed->GetBinCenter(j);
- Double_t e = 0;
- Double_t y = fit->Evaluate(x, e);
- gr->SetPoint(j-1, x, y);
- gr->SetPointError(j-1, 0, e);
- }
- gr->Draw("same l e4");
-
- TLegend* l = new TLegend(.15, .11, .5, .4);
- l->SetBorderSize(0);
- l->SetFillColor(0);
- l->SetFillStyle(0);
- l->SetTextFont(132);
- l->AddEntry(fSimple, "L: Single-L", "lp");
- l->AddEntry(fSingle, "LG: Single-L#otimesG", "lp");
- l->AddEntry(fSummed,"NLG: Multi-L#otimesG", "lp");
- l->AddEntry(gr, "f_{N}: Fit to NLG");
- l->Draw();
-
- // --- Second pad: Data (1, 2, 3, ...) --------------------------
- p = c->cd(2);
- p->SetLogy();
- p->SetBorderSize(0);
- p->SetBorderMode(0);
- p->SetTopMargin(0.01);
- p->SetRightMargin(0.01);
- p->SetGridx();
- THStack* stack = new THStack(fNSum, "y");
- TIter next(stack->GetHists());
- TH1* h = 0;
- Int_t i = 2;
- Double_t max = 0;
- TLegend* l2 = new TLegend(.6, .5, .95, .95);
- l2->SetBorderSize(0);
- l2->SetFillStyle(0);
- l2->SetFillColor(0);
- l2->SetTextFont(132);
- while ((h = static_cast<TH1*>(next()))) {
- if (i == 2) max = h->GetMaximum();
- h->SetLineColor(i);
- h->SetFillColor(i);
- h->SetFillStyle(3001);
- h->Scale(1/max);
- l2->AddEntry(h, Form("%d particles", i - 1), "f");
- i++;
- }
- stack->Draw("hist nostack");
- stack->GetHistogram()->SetXTitle("Signal");
- stack->GetHistogram()->SetYTitle("Events");
- i = 2;
- for (Int_t j = 1; j <= 5; j++) {
- TF1* fi = Function::MakeIFunc(j, fRP[j-1], fDelta, fXi, fSigma, 0, 10);
- if (j == 1) max = fi->GetMaximum();
- fi->SetParameter(0, fRP[j-1]/max);
- fi->SetLineColor(i);
- fi->SetLineWidth(3);
- fi->SetLineStyle(2);
- fi->Draw("same");
-
- TF1* fj = Function::MakeIFunc(j, fit->GetC() * fit->GetA(j),
- fit->GetDelta(), fit->GetXi(),
- fit->GetSigma(), 0, 10);
- fj->SetLineColor(i);
- fj->SetLineWidth(2);
- fj->Draw("same");
-
- std::cout << "Component " << j << " scale=" << fi->GetParameter(0)
- << " (" << fRP[j-1] << ")" << " " << fj->GetParameter(0)
- << " (" << fit->GetC() << "*" << fit->GetA(j) << ")"
- << std::endl;
- i++;
- }
- TLegendEntry* ee = l2->AddEntry("", "Input f_{i}", "l");
- ee->SetLineStyle(2);
- ee->SetLineWidth(3);
- l2->AddEntry("", "Fit f_{i}", "l");
- l2->Draw();
- // fNSum->Draw("lego2");
-
- // --- Third pad: Mean multiplicity ------------------------------
- TH1* nhist1 = new TH1F("nhist1", "NHist", 5, 0.5, 5.5);
- nhist1->SetFillColor(kCyan+1);
- nhist1->SetMarkerColor(kCyan+1);
- nhist1->SetLineColor(kCyan+1);
- nhist1->SetFillStyle(3002);
- TH1* nhist2 = new TH1F("nhist2", "NHist", 5, 0.5, 5.5);
- nhist2->SetFillColor(kYellow+1);
- nhist2->SetMarkerColor(kYellow+1);
- nhist2->SetLineColor(kYellow+1);
- nhist2->SetFillStyle(3002);
-
- TH2* resp = new TH2F("resp", "Reponse", 100, 0, 10, 6, -.5, 5.5);
- resp->SetXTitle("Signal");
- resp->SetYTitle("# particles");
- for (Int_t j = 0; j < fNObs; j++) {
- if (((j+1) % (fNObs / 10)) == 0)
- std::cout << "Event # " << j+1 << std::endl;
- Double_t x = fSummed->GetRandom();
- Double_t y = fit->RandomEstimateNParticles(x);
- nhist1->Fill(y);
- resp->Fill(x, y);
- }
- TGraphErrors* graph = new TGraphErrors(fSummed->GetNbinsX());
- graph->SetName("evalWeighted");
- graph->SetTitle("Evaluated function");
- graph->SetLineColor(kBlack);
- graph->SetFillColor(kYellow+1);
- graph->SetFillStyle(3001);
- for (Int_t j = 1; j <= fSummed->GetNbinsX(); j++) {
- Double_t x = fSummed->GetBinCenter(j);
- Double_t e = 0;
- Double_t y = fit->EstimateNParticles(x, e);
- nhist2->Fill(y, fSummed->GetBinContent(j));
- graph->SetPoint(j-1, x, y);
- graph->SetPointError(j-1, 0, e);
- }
- nhist1->Scale(1. / nhist1->GetMaximum());
- nhist2->Scale(1. / nhist2->GetMaximum());
- fN->Scale(1. / fN->GetMaximum());
-
- p = c->cd(3);
- p->SetLogy();
- p->SetBorderSize(0);
- p->SetBorderMode(0);
- p->SetTopMargin(0.01);
- p->SetRightMargin(0.01);
- fN->Draw();
- nhist1->Draw("same");
- nhist2->Draw("same");
-
- TLegend* l3 = new TLegend(.3, .7, .95, .95);
- l3->SetBorderSize(0);
- l3->SetFillStyle(0);
- l3->SetFillColor(0);
- l3->SetTextFont(132);
- l3->AddEntry(fN,
- Form("Input n distribution: #bar{m}=%5.3f, s_{m}=%5.3f",
- fN->GetMean(), fN->GetRMS()), "f");
- l3->AddEntry(nhist1,
- Form("Random draw from fit: #bar{m}=%5.3f, s_{m}=%5.3f",
- nhist1->GetMean(), nhist1->GetRMS()), "f");
- l3->AddEntry(nhist2,
- Form("Weighted evaluation of fit: #bar{m}=%5.3f, s_{m}=%5.3f",
- nhist2->GetMean(), nhist2->GetRMS()), "f");
- l3->Draw();
-
- // --- Fourth pad: Reponse function ------------------------------
- p = c->cd(4);
- // p->SetLogy();
- p->SetBorderSize(0);
- p->SetBorderMode(0);
- p->SetTopMargin(0.01);
- p->SetRightMargin(0.01);
- p->SetGridx();
- TProfile* prof1 = fNSum->ProfileY();
- prof1->SetLineColor(kMagenta+1);
- prof1->Draw();
- TProfile* prof2 = resp->ProfileX();
- prof2->SetLineColor(kCyan+2);
- prof2->Draw("same");
- graph->SetLineWidth(2);
- graph->Draw("same LP E4");
-
- TLegend* l4 = new TLegend(.2, .11, .8, .4);
- l4->SetBorderSize(0);
- l4->SetFillStyle(0);
- l4->SetFillColor(0);
- l4->SetTextFont(132);
- l4->AddEntry(prof1, "Input distribution of N particles", "l");
- l4->AddEntry(prof2, "Random estimate of N particles", "l");
- l4->AddEntry(graph, "E: #sum_{i}^{N}i a_{i}f_{i}/"
- "#sum_{i}^{N}a_{i}f_{i} evaluated over NLG", "lf");
- l4->Draw();
-
- c->cd();
-
- if (type && type[0] != '\0')
- c->Print(Form("TestELossDist.%s", type));
- }
-};
-
-
-void
-TestELossDist(const char* type="png")
-{
- gRandom->SetSeed(12345);
- const Double_t a1[] = { 1, .05, .02, .005, .002 }; // 7TeV pp
- const Double_t a2[] = { 1, .1, .05, .01, .005 };
- const Double_t a3[] = { 1, .2, .10, .05, .01 };
-
- Generator* g = new Generator(0.55, 0.06, 0.06, 5, a1);
- g->Generate(100000);
- g->Draw(type);
-
-
-}
-
-
+++ /dev/null
-/**
- * Test script to fit the energy loss spectra
- *
- * @deprecated
- * This is a simple test script
- *
- * @ingroup pwg2_forward_analysis_scripts
- */
-#ifndef __CINT__
-# include "AliForwardUtil.h"
-# include <TFile.h>
-# include <TList.h>
-# include <TH1.h>
-# include <TF1.h>
-# include <TFitResult.h>
-# include <TMath.h>
-# include <TStyle.h>
-# include <TArrow.h>
-# include <TCanvas.h>
-#else
-class TCanvas;
-class TFile;
-class TH1;
-class TList;
-class TF1;
-#endif
-
-//____________________________________________________________________
-/**
- *
- *
- * @param ef
- * @param d
- * @param r
- * @param etabin
- *
- * @return
- *
- * @deprecated
- * This is a simple test script
- *
- * @ingroup pwg2_forward_analysis_scripts
- */
-TH1* GetEDist(TList* ef, UShort_t d, Char_t r, UShort_t etabin)
-{
- TList* dL = static_cast<TList*>(ef->FindObject(Form("FMD%d%c",d,r)));
- if (!dL) {
- Error("GetEDist", "Couldn't find list FMD%d%c",d,r);
- ef->ls();
- return 0;
- }
- if (etabin > 999) {
- TH1* hist = static_cast<TH1*>(dL->FindObject(Form("FMD%d%c_edist",d,r)));
- if (hist) {
- Error("GetEDist", "Couldn't find EDists histogram for FMD%d%c",d,r);
- return 0;
- }
- }
-
- TList* edL = static_cast<TList*>(dL->FindObject("EDists"));
- if (!edL) {
- Error("GetEDist", "Couldn't find list EDists for FMD%d%c",d,r);
- return 0;
- }
-
- TH1* hist = static_cast<TH1*>(edL->FindObject(Form("FMD%d%c_etabin%03d",
- d, r, etabin)));
- if (!hist) {
- Error("GetEDist", "Couldn't find histogra FMD%d%c_etabin%03d",
- d,r, etabin);
- return 0;
- }
-
- return hist;
-}
-
-//____________________________________________________________________
-/**
- *
- *
- * @param ef
- * @param d
- * @param r
- * @param eta
- *
- * @return
- *
- * @deprecated
- * This is a simple test script
- *
- * @ingroup pwg2_forward_analysis_scripts
-*/
-TH1* GetEDist(TList* ef, UShort_t d, Char_t r, Float_t eta)
-{
- if (!ef) {
- Error("GetEDist", "EF not set");
- return 0;
- }
- TAxis* etaAxis = static_cast<TAxis*>(ef->FindObject("etaAxis"));
- if (!etaAxis) {
- Error("GetEDist", "Couldn't find eta axis in list");
- return 0;
- }
-
- UShort_t bin = etaAxis->FindBin(eta);
- if (bin <= 0 || bin > etaAxis->GetNbins()) {
- Error("GetEDist", "eta=%f out of range [%f,%f] - getting ring histo",
- eta, etaAxis->GetXmin(), etaAxis->GetXmax());
- return GetEDist(ef, d, r, UShort_t(1000));
- }
-
- return GetEDist(ef, d, r, bin);
-}
-
-//____________________________________________________________________
-/**
- *
- *
- * @param file
- *
- * @return
- *
- * @deprecated
- * This is a simple test script
- *
- * @ingroup pwg2_forward_analysis_scripts
- */
-TList* GetEF(TFile* file)
-{
- TList* forward = static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
- if (!forward) {
- Error("GetEF", "Failed to get list PWG2forwardDnDeta/Forward from file");
- return 0;
- }
- TList* ef = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
- if (!ef) {
- Error("GetEF", "Failed to get energy fitter list");
- return 0;
- }
-
- return ef;
-}
-
-//____________________________________________________________________
-TList* ef = 0;
-
-//____________________________________________________________________
-/**
- *
- *
- *
- * @return
- *
- * @deprecated
- * This is a simple test script
- *
- * @ingroup pwg2_forward_analysis_scripts
- */
-TList* CheckEF()
-{
- if (ef) return ef;
-
- TFile* file = TFile::Open("AnalysisResults.root", "READ");
- if (!file) {
- Error("Fit1", "Failed to open file");
- return 0;
- }
- return GetEF(file);
-}
-
-//____________________________________________________________________
-TCanvas* c = 0;
-
-//____________________________________________________________________
-/**
- *
- *
- *
- * @return
- *
- * @deprecated
- * This is a simple test script
- *
- * @ingroup pwg2_forward_analysis_scripts
- */
-TCanvas* CheckC()
-{
- if (c) return c;
-
- gStyle->SetOptFit(1111);
- gStyle->SetCanvasColor(0);
- gStyle->SetCanvasBorderSize(0);
- gStyle->SetCanvasBorderMode(0);
- gStyle->SetCanvasDefW(800);
- gStyle->SetCanvasDefH(800);
- gStyle->SetPadTopMargin(0.05);
- gStyle->SetPadRightMargin(0.05);
- gStyle->SetTitleBorderSize(0);
- gStyle->SetTitleFillColor(0);
- gStyle->SetTitleStyle(0);
- gStyle->SetStatBorderSize(1);
- gStyle->SetStatColor(0);
- gStyle->SetStatStyle(0);
- gStyle->SetStatX(0.95);
- gStyle->SetStatY(0.95);
- gStyle->SetStatW(0.15);
- gStyle->SetStatH(0.15);
-
- c = new TCanvas("c", "c");
- c->SetLogy();
-
- return c;
-}
-
-//____________________________________________________________________
-/**
- *
- *
- * @param f
- *
- * @deprecated
- * This is a simple test script
- *
- * @ingroup pwg2_forward_analysis_scripts
- */
-void PrintFit(TF1* f)
-{
- Int_t nu = f->GetNDF();
- Double_t chi2 = f->GetChisquare();
- Double_t chi2nu = (nu > 0 ? chi2/nu : 0);
- Printf("%s: chi^2/nu=%f/%d=%f [%f,%f]",
- f->GetName(), chi2, nu, chi2nu,
- f->GetXmin(), f->GetXmax());
- for (Int_t i = 0; i < f->GetNpar(); i++) {
- Double_t v = f->GetParameter(i);
- Double_t e = f->GetParError(i);
- Double_t r = 100*(v == 0 ? 1 : e / v);
- Printf("%32s = %14.7f +/- %14.7f (%5.1f%%)",f->GetParName(i),v,e,r);
- }
-}
-
-//____________________________________________________________________
-/**
- *
- *
- * @param n
- * @param d
- * @param r
- * @param eta
- *
- * @deprecated
- * This is a simple test script
- *
- * @ingroup pwg2_forward_analysis_scripts
- */
-void TestFitELoss(Int_t n, UShort_t d, Char_t r, Float_t eta)
-{
- TList* ef1 = CheckEF();
- TCanvas* c1 = CheckC();
- if (!ef1) return;
- if (!c1) return;
-
- TH1* dist = GetEDist(ef1, d, r, eta);
- if (!dist) return;
-
- AliForwardUtil::ELossFitter f(0.4, 10, 4);
- TF1* landau1 = new TF1(*f.Fit1Particle(dist, 0));
- landau1->SetName("Landau1");
- PrintFit(landau1);
- TF1* landaun = new TF1(*f.FitNParticle(dist, n, 0));
- landau1->SetName(Form("Landau%d", n));
- PrintFit(landaun);
- landau1->SetRange(0,10);
- landaun->SetRange(0,10);
- landau1->SetLineWidth(4);
- landaun->SetLineWidth(4);
- landau1->SetLineColor(kBlack);
- landaun->SetLineColor(kBlack);
-
- dist->GetListOfFunctions()->Clear();
- dist->GetListOfFunctions()->Add(landau1);
- dist->GetListOfFunctions()->Add(landaun);
- dist->GetListOfFunctions()->ls();
- dist->Draw();
- landau1->Draw("same");
- landaun->Draw("same");
-
- Double_t mp = landaun->GetParameter(1);
- Double_t xi = landaun->GetParameter(2);
- for (Int_t i = 1; i <= n; i++) {
- Double_t x = i * (mp + xi * TMath::Log(i));
- Double_t y = landaun->Eval(x);
- Double_t y1 = y < 0.05 ? 1 : 0.01;
- TArrow* a = new TArrow(x,y1,x,y,0.03,"|>");
- Info("FitSteer", "Delta_{p,%d}=%f", i, x);
- a->SetLineWidth(2);
- a->SetAngle(30);
- a->Draw();
- }
-
- c1->cd();
-}
-//
-// EOF
-//
+++ /dev/null
-#ifndef __CINT__
-#include "AliForwardCorrectionManager.h"
-#include "AliFMDCorrELossFit.h"
-#include <TAxis.h>
-#include <TFile.h>
-#include <TList.h>
-#include <TError.h>
-#include <TH1.h>
-#include <TF1.h>
-#include <TClonesArray.h>
-#else
-class TAxis;
-class AliFMDCirrELossFit;
-class TH1;
-
-#endif
-
-/**
- * Class to make correction object and save to file
- *
- * Run like
- *
- * @verbatim
- * gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
- * gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/Compile.C");
- * Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/MakeELossFit.C");
- * MakeELossFit mef(sys, cms, field, mc, "AnalysisResults.root");
- * mef.Run();
- * @endverbatim
- * where @e sys, the collision system, is one of
- * - 1: pp
- * - 2: PbPb
- *
- * @e cms is the center of mass energy per nuclean in GeV (e.g., 2760
- * for first PbPb data), @e field is (signed) L3 magnetic in kG, and
- * @e mc is wether this correction applies to MC data or not.
- *
- * The class generates a file like
- * @verbatim
- * elossfits<sys>_<cms>GeV_<field>kG_<realmc>.root
- * @endverbatim
- * in the working directory. This file can be moved to
- * @verbatim
- * $(ALICE_ROOT)/PWG2/FORWARD/corrections/ELossFits
- * @endverbatim
- * and stored in SVN for later use.
- *
- * @depcrecated
- * The class AliFMDELossFitter automatically generates the
- * AliFMDCorrELossFit object.
- *
- * @ingroup pwg2_forward_analysis_scripts
- */
-
-class MakeELossFit
-{
-protected:
-public:
- TList* fFitter;
- TAxis* fAxis;
- TClonesArray fFits;
- UShort_t fSys;
- UShort_t fCMS;
- Short_t fField;
- Bool_t fMC;
-
- //__________________________________________________________________
- MakeELossFit(UShort_t sys,
- UShort_t cms,
- Short_t field,
- Bool_t mc,
- const char* filename="forward_eloss.root")
- : fFitter(0),
- fAxis(0),
- fFits("AliFMDCorrELossFit::ELossFit"),
- fSys(sys),
- fCMS(cms),
- fField(field),
- fMC(mc)
- {
- TFile* file = TFile::Open(filename, "READ");
- if (!file) {
- Error("MakeELossFit", "Failed to open file %s", filename);
- return;
- }
- TList* forward = static_cast<TList*>(file->Get("Forward"));
- // static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
- if (!forward) {
- Error("MakeELossFit", "Couldn't get forward list from %s", filename);
- return;
- }
-
- fFitter = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
- if (!fFitter) {
- Error("MakeELossFit", "Couldn't get fitter folder");
- return;
- }
-
- fAxis = static_cast<TAxis*>(fFitter->FindObject("etaAxis"));
- if (!fAxis) {
- Error("MakeELossFit", "Couldn't get eta axis");
- return;
- }
- file->Close();
-
-#if 0
- AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
- mgr.Init(sys, cms, field, mc, 0, true);
-#endif
- }
- //__________________________________________________________________
- TList* GetRing(UShort_t d, Char_t r) const
- {
- TList* rL = static_cast<TList*>(fFitter->FindObject(Form("FMD%d%c",d,r)));
- if (!rL) {
- Warning("DrawFits", "List FMD%d%c not found", d, r);
- return 0;
- }
- return rL;
- }
- //__________________________________________________________________
- TList* GetEDists(UShort_t d, Char_t r) const
- {
- TList* rList = GetRing(d, r);
- if (!rList)
- return 0;
- // rList->ls();
-
- TList* edists = static_cast<TList*>(rList->FindObject("EDists"));
- if (!edists) {
- Error("PrintFits", "Couldn't get EDists list for FMD%d%c", d,r);
- return 0;
- }
- return edists;
- }
- //__________________________________________________________________
- TH1* GetEDist(UShort_t d, Char_t r, UShort_t etabin)
- {
- TList* eList = GetEDists(d, r);
- if (!eList) {
- Warning("GetEDist", "No list for FMD%d%c", d, r);
- return 0;
- }
- // eList->ls();
-
- TString cmp(Form("FMD%d%c_etabin%03d", d, r, etabin));
-
- TIter next(eList);
- TObject* o = 0;
- while ((o = next())) {
- if (!cmp.CompareTo(o->GetName())) {
- return static_cast<TH1*>(o);
- }
- }
- return 0;
- }
-#ifndef __CINT__
- //__________________________________________________________________
- AliFMDCorrELossFit::ELossFit* FindBestFit(TH1* dist)
- {
- TList* funcs = dist->GetListOfFunctions();
- TIter next(funcs);
- TF1* func = 0;
- fFits.Clear();
- Int_t i = 0;
- Info("FindBestFit", "%s", dist->GetName());
- while ((func = static_cast<TF1*>(next()))) {
- AliFMDCorrELossFit::ELossFit* fit =
- new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func);
- fit->CalculateQuality(10, .20, 1e-7);
- }
- fFits.Sort(false);
- // fFits.Print();
- return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0));
- }
-#endif
- //__________________________________________________________________
- Bool_t Run()
- {
- if (!fFitter || !fAxis) {
- Error("Run", "Missing objects");
- return kFALSE;
- }
- AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
- obj->SetEtaAxis(*fAxis);
-
- for (UShort_t d = 1; d <= 3; d++) {
- Info("Run", "detector is FMD%d", d);
- UShort_t nQ = (d == 1 ? 1 : 2);
- for (UShort_t q = 0; q < nQ; q++) {
- Char_t r = (q == 0 ? 'I' : 'O');
- Int_t nBin = fAxis->GetNbins();
- Info("Run", " ring is FMD%d%c - %d bins", d, r, nBin);
-
- Bool_t oneSeen = kFALSE;
- for (UShort_t b = 1; b <= nBin; b++) {
- TH1* h = GetEDist(d, r, b);
- if (oneSeen && !h) break;
- if (!h) continue;
- if (!oneSeen) oneSeen = true;
-
- AliFMDCorrELossFit::ELossFit* best = FindBestFit(h);
- best->Print("");
- best->fDet = d;
- best->fRing = r;
- best->fBin = b;
- // Double_t eta = fAxis->GetBinCenter(b);
- Info("Run", " bin=%3d->%6.4f", b, eta);
- obj->SetFit(d, r, b, new AliFMDCorrELossFit::ELossFit(*best));
- }
- }
- }
- AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
- TString fname(mgr.GetFileName(AliForwardCorrectionManager::kELossFits,
- fSys, fCMS, fField, fMC));
- TFile* output = TFile::Open(fname.Data(), "RECREATE");
- if (!output) {
- Warning("Run", "Failed to open output file %s", fname.Data());
- return kFALSE;
- }
- obj->Write(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
- output->Write();
- output->Close();
- Info("Run", "File %s created. It should be copied to %s and stored in SVN",
- fname.Data(),mgr.GetFileDir(AliForwardCorrectionManager::kELossFits));
-
- return kTRUE;
- }
-};
-//
-// EOF
-//
+++ /dev/null
-namespace {
- enum {
- kSolid = 0x000,
- kHollow = 0x001,
- kCircle = 0x002,
- kSquare = 0x004,
- kUpTriangle = 0x006,
- kDownTriangle = 0x008,
- kDiamond = 0x00a,
- kCross = 0x00c,
- kStar = 0x00e
- };
- Int_t MarkerStyle(UInt_t bits)
- {
- Int_t base = bits & (0xFE);
- Bool_t hollow = bits & kHollow;
- switch (base) {
- case kCircle: return (hollow ? 24 : 20);
- case kSquare: return (hollow ? 25 : 21);
- case kUpTriangle: return (hollow ? 26 : 22);
- case kDownTriangle: return (hollow ? 32 : 23);
- case kDiamond: return (hollow ? 27 : 33);
- case kCross: return (hollow ? 28 : 34);
- case kStar: return (hollow ? 30 : 29);
- }
- return 1;
- }
- UShort_t MarkerBits(Int_t style)
- {
- UShort_t bits = 0;
- switch (style) {
- case 24: case 25: case 26: case 27: case 28: case 30: case 32:
- bits |= kHollow; break;
- }
- switch (style) {
- case 20: case 24: bits |= kCircle; break;
- case 21: case 25: bits |= kSquare; break;
- case 22: case 26: bits |= kUpTriangle; break;
- case 23: case 32: bits |= kDownTriangle; break;
- case 27: case 33: bits |= kDiamond; break;
- case 28: case 34: bits |= kCross; break;
- case 29: case 30: bits |= kStar; break;
- }
- return bits;
- }
- Int_t FlipHollow(Int_t style)
- {
- UShort_t bits = MarkerBits(style);
- Int_t ret = MarkerStyle(bits ^ kHollow);
- Info("FlipHollow", "style=%2d -> bits=0x%02x -> mask=0x%02x -> ret=%02d",
- style, bits, (bits ^ kHollow), ret);
- return ret;
- }
-}
-
-void DrawOne(const char* what, UShort_t base, Double_t y)
-{
- TLatex* l = new TLatex(.07, y, what);
- l->SetTextSize(0.03);
- l->Draw();
-
- Int_t filled = MarkerStyle(base);
- // Info("DrawOne", "%2d (%16s) -> %d", base, what, style);
- TMarker* p = new TMarker(.35, y, filled);
- p->SetMarkerSize(1.5);
- p->Draw();
-
- Int_t hollow = MarkerStyle(base|kHollow);
- p = new TMarker(.60, y, hollow);
- p->SetMarkerSize(1.5);
- p->Draw();
-
- p = new TMarker(.75, y, FlipHollow(filled));
- p->SetMarkerSize(1.5);
- p->Draw();
-
- p = new TMarker(.85, y, FlipHollow(hollow));
- p->SetMarkerSize(1.5);
- p->Draw();
-
-}
-
+++ /dev/null
-/**
- * Run the energy loss fit finder and generate corrections output file
- *
- * @param sys Collision system
- * @param cms Center of mass energy per nucleon in GeV
- * @param field Magnetic field
- * @param mc Whether this is for Monte-Carlo data
- * @param filename Input file name
- *
- * @ingroup pwg2_forward_analysis_scripts
- *
- * @depcrecated
- * The class AliFMDELossFitter automatically generates the
- * AliFMDCorrELossFit object.
- *
- */
-void
-TestRunMakeELossFit(UShort_t sys,
- UShort_t cms,
- Short_t field,
- Bool_t mc=false,
- const char* filename="forward_eloss.root")
-{
- std::cout << "Loading libraries ..." << std::endl;
- gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
-
- std::cout << "Loading compile script ..." << std::endl;
- gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/Compile.C");
-
- std::cout << "Compiling MakeELossFit.C script ..." << std::endl;
- Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/TestMakeELossFit.C");
-
- std::cout << "Making MakeELossFit object (sys=" << sys
- << ", cms=" << cms << ", field=" << field << ", mc=" << mc
- << ")" << std::endl;
- MakeELossFit mef(sys, cms, field, mc, filename);
-
- std::cout << "Running maker ..." << std::endl;
- mef.Run();
-}
-//
-// EOF
-//
--- /dev/null
+Bool_t
+BuildTrain(const char* script, const char* extra="", Bool_t useTmp=false)
+{
+ // --- Load needed libraries ---------------------------------------
+ gROOT->LoadClass("TVirtualMC", "libVMC");
+ gROOT->LoadClass("AliVEvent", "libSTEERBase");
+ gROOT->LoadClass("AliESDEvent", "libESD");
+ gROOT->LoadClass("AliAnalysisManager", "libANALYSIS");
+ gROOT->LoadClass("AliAnalysisTaskSE", "libANALYSISalice");
+ gROOT->LoadClass("AliAODForwardMult", "libPWG2forward2");
+
+ // --- Setup script path -------------------------------------------
+ const char* aliPath = gSystem->ExpandPathName("$ALICE_ROOT");
+ const char* fwd2Path =
+ gSystem->ExpandPathName("$ALICE_ROOT/PWG2/FORWARD/analysis2");
+ const char* macroPath = gROOT->GetMacroPath();
+ gROOT->SetMacroPath(".:%s:%s/trains:%s/scripts",macroPath,fwd2Path,fwd2Path);
+
+ // --- Setup include path ------------------------------------------
+ gSystem->AddIncludePath(Form("-I%s", fwd2Path));
+ gSystem->AddIncludePath(Form("-I%s", aliPath));
+ gSystem->AddIncludePath(Form("-I%s/include", aliPath));
+
+ // --- Check that we can find the script ---------------------------
+ TString path = gSystem->Which(gROOT->GetMacroPath(), script);
+ if (path.IsNull()) {
+ path = gSystem->Which(gROOT->GetMacroPath(), Form("%s.C", script));
+ if (path.IsNull()) {
+ Error("BuildTrain", "Cannot find script %s or %s.C in %s",
+ script, script, gROOT->GetMacroPath());
+ return false;
+ }
+ }
+
+ // --- Possibly make a temporary copy ------------------------------
+ if (useTmp) {
+ TString tmp("Train");
+ FILE* fp = gSystem->TempFileName(tmp, ".");
+ fclose(fp);
+ gSystem->Unlink(tmp);
+ tmp.Append(".C");
+
+ gSystem->CopyFile(path, tmp);
+ path = tmp;
+ }
+
+ // --- Now compile the thing ---------------------------------------
+ Int_t ret = gROOT->LoadMacro(Form("%s+%s", path.Data(), extra));
+
+ // --- If we did a temporary file, remove stuff --------------------
+ if (useTmp) {
+ gSystem->Unlink(path);
+ path.ReplaceAll(".C", "_C.d"); gSystem->Unlink(path);
+ path.ReplaceAll("_C.d", "_C.so"); gSystem->Unlink(path);
+ }
+
+ // --- Warning in case of problems ---------------------------------
+ if (ret != 0)
+ Warning("BuildSetup", "Failed to build %s (%s)", path.Data(), script);
+
+ return ret == 0;
+}
+//
+// EOF
+//
+
--- /dev/null
+//====================================================================
+/**
+ * Analysis train to make Forward and Central multiplicity
+ *
+ * To run, do
+ * @code
+ * gROOT->LoadMacro("TrainSetup.C");
+ * // Make train
+ * MakeAODTrain t("My Analysis");
+ * // Set variaous parameters on the train
+ * t.SetDataDir("/home/of/data");
+ * t.AddRun(118506)
+ * // Run it
+ * t.Run("LOCAL", "FULL", -1, false, false);
+ * @endcode
+ *
+ * @ingroup pwg2_forward_scripts_makers
+ * @ingroup pwg2_forward_aod
+ */
+class MakeAODTrain : public TrainSetup
+{
+public:
+ /**
+ * Constructor. Date and time must be specified when running this
+ * in Termiante mode on Grid
+ *
+ * @param name Name of train (free form)
+ * @param sys Collision system (1: pp, 2: PbPb)
+ * @param sNN Center of mass energy [GeV]
+ * @param field L3 magnetic field - one of {-5,0,+5} kG
+ * @param useCent Whether to use centrality
+ * @param dateTime Append date and time to name
+ * @param year Year - if not specified, current year
+ * @param month Month - if not specified, current month
+ * @param day Day - if not specified, current day
+ * @param hour Hour - if not specified, current hour
+ * @param min Minutes - if not specified, current minutes
+ */
+ MakeAODTrain(const char* name,
+ UShort_t sys = 0,
+ UShort_t sNN = 0,
+ Short_t field = 0,
+ Bool_t useCent = false,
+ Bool_t dateTime = false,
+ UShort_t year = 0,
+ UShort_t month = 0,
+ UShort_t day = 0,
+ UShort_t hour = 0,
+ UShort_t min = 0)
+ : TrainSetup(name, dateTime,
+ year, month, day, hour, min),
+ fSys(sys),
+ fSNN(sNN),
+ fField(field),
+ fUseCent(useCent)
+ {}
+ /**
+ * Run this analysis
+ *
+ * @param mode Mode - see TrainSetup::EMode
+ * @param oper Operation - see TrainSetup::EOperation
+ * @param nEvents Number of events (negative means all)
+ * @param mc If true, assume simulated events
+ * @param usePar If true, use PARs
+ */
+ void Run(const char* mode, const char* oper,
+ Int_t nEvents=-1, Bool_t mc=false,
+ Bool_t usePar=false)
+ {
+ Info("Run", "Running in mode=%s, oper=%s, events=%d, MC=%d, Par=%d",
+ mode, oper, nEvents, mc, usePar);
+ Exec("ESD", mode, oper, nEvents, mc, usePar);
+ }
+ /**
+ * Run this analysis
+ *
+ * @param mode Mode - see TrainSetup::EMode
+ * @param oper Operation - see TrainSetup::EOperation
+ * @param nEvents Number of events (negative means all)
+ * @param mc If true, assume simulated events
+ * @param usePar If true, use PARs
+ */
+ void Run(EMode mode, EOper oper, Int_t nEvents=-1, Bool_t mc=false,
+ Bool_t usePar = false)
+ {
+ Info("Run", "Running in mode=%d, oper=%d, events=%d, MC=%d, Par=%d",
+ mode, oper, nEvents, mc, usePar);
+ Exec(kESD, mode, oper, nEvents, mc, usePar);
+ }
+protected:
+ /**
+ * Create the tasks
+ *
+ * @param mode Processing mode
+ * @param par Whether to use par files
+ * @param mgr Analysis manager
+ */
+ void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager* mgr)
+ {
+ // --- Output file name ------------------------------------------
+ AliAnalysisManager::SetCommonFileName("forward.root");
+
+ // --- Load libraries/pars ---------------------------------------
+ LoadLibrary("PWG2forward2", mode, par, true);
+
+ // --- Set load path ---------------------------------------------
+ gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
+ gROOT->GetMacroPath()));
+
+ // --- Check if this is MC ---------------------------------------
+ Bool_t mc = mgr->GetMCtruthEventHandler() != 0;
+
+ // --- Add the task ----------------------------------------------
+ gROOT->Macro(Form("AddTaskForwardMult.C(%d,%d,%d,%d)",
+ mc, fSys, fSNN, fField));
+ AddExtraFile(gSystem->Which(gROOT->GetMacroPath(), "ForwardAODConfig.C"));
+
+ // --- Add the task ----------------------------------------------
+ gROOT->Macro(Form("AddTaskCentralMult.C(%d,%d,%d,%d)",
+ mc, fSys, fSNN, fField));
+ AddExtraFile(gSystem->Which(gROOT->GetMacroPath(), "CentralAODConfig.C"));
+ }
+ //__________________________________________________________________
+ /**
+ * Create physics selection , and add to manager
+ *
+ * @param mc Whether this is for MC
+ * @param mgr Manager
+ */
+ void CreatePhysicsSelection(Bool_t mc,
+ AliAnalysisManager* mgr)
+ {
+ TrainSetup::CreatePhysicsSelection(mc, mgr);
+
+ // --- Get input event handler -----------------------------------
+ AliInputEventHandler* ih =
+ dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
+ if (!ih)
+ Fatal("CreatePhysicsSelection", "Couldn't get input handler (%p)", ih);
+
+ // --- Get Physics selection -------------------------------------
+ AliPhysicsSelection* ps =
+ dynamic_cast<AliPhysicsSelection*>(ih->GetEventSelection());
+ if (!ps)
+ Fatal("CreatePhysicsSelection", "Couldn't get PhysicsSelection (%p)",ps);
+
+ // --- Ignore trigger class when selecting events. This means ---
+ // --- that we get offline+(A,C,E) events too --------------------
+ // ps->SetSkipTriggerClassSelection(true);
+ }
+ //__________________________________________________________________
+ /**
+ * Create the centrality selection only if requested
+ *
+ * @param mc Monte-Carlo truth flag
+ * @param mgr Manager
+ */
+ void CreateCentralitySelection(Bool_t mc, AliAnalysisManager* mgr)
+ {
+ if (!fUseCent) return;
+ TrainSetup::CreateCentralitySelection(mc, mgr);
+ }
+ UShort_t fSys;
+ UShort_t fSNN;
+ Short_t fField;
+ Bool_t fUseCent;
+};
+//
+// EOF
+//
--- /dev/null
+//====================================================================
+/**
+ * Analysis train to do energy loss fits
+ *
+ * @ingroup pwg2_forward_scripts_makers
+ */
+class MakeFMDELossTrain : public TrainSetup
+{
+public:
+ /**
+ * Constructor. Date and time must be specified when running this
+ * in Termiante mode on Grid
+ *
+ * @param name Name of train
+ * @param useCent Whether to use centrality or not
+ * @param dateTime Append date and time to name
+ * @param year Year
+ * @param month Month
+ * @param day Day
+ * @param hour Hour
+ * @param min Minutes
+ */
+ MakeFMDELossTrain(const char* name = "FMD Energy Loss",
+ Bool_t useCent = false,
+ Bool_t dateTime = false,
+ UShort_t year = 0,
+ UShort_t month = 0,
+ UShort_t day = 0,
+ UShort_t hour = 0,
+ UShort_t min = 0)
+ : TrainSetup(name, dateTime, year, month, day, hour, min),
+ fUseCent(useCent)
+ {}
+ //__________________________________________________________________
+ /**
+ * Run this analysis
+ *
+ * @param mode Mode
+ * @param oper Operation
+ * @param nEvents Number of events (negative means all)
+ * @param mc If true, assume simulated events
+ * @param par IF true, use par files
+ */
+ void Run(const char* mode, const char* oper,
+ Int_t nEvents=-1, Bool_t mc=false, Bool_t par=false)
+ {
+ Exec("ESD", mode, oper, nEvents, mc, par);
+ }
+ //__________________________________________________________________
+ /**
+ * Run this analysis
+ *
+ * @param mode Mode
+ * @param oper Operation
+ * @param nEvents Number of events (negative means all)
+ * @param mc If true, assume simulated events
+ * @param par IF true, use par files
+ */
+ void Run(EMode mode, EOper oper, Int_t nEvents=-1, Bool_t mc=false,
+ Bool_t par=false)
+ {
+ Exec(kESD, mode, oper, nEvents, mc, par);
+ }
+protected:
+ //__________________________________________________________________
+ /**
+ * Create the tasks
+ *
+ * @param mode Processing mode
+ * @param par Whether to use par files
+ * @param mgr Analysis manager
+ */
+ void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager* mgr)
+ {
+ // --- Output file name ------------------------------------------
+ AliAnalysisManager::SetCommonFileName("forward_eloss.root");
+
+ // --- Load libraries/pars ---------------------------------------
+ LoadLibrary("PWG2forward2", mode, par, true);
+
+ // --- Set load path ---------------------------------------------
+ gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
+ gROOT->GetMacroPath()));
+
+ // --- Check if this is MC ---------------------------------------
+ Bool_t mc = mgr->GetMCtruthEventHandler() != 0;
+
+ // --- Add the task ----------------------------------------------
+ gROOT->Macro(Form("AddTaskFMDELoss.C(%d,%d)", mc, fUseCent));
+ }
+ /**
+ * Create entrality selection if enabled
+ *
+ * @param mc Whether this is MC or not
+ * @param mgr Analysis manager
+ */
+ virtual void CreateCentralitySelection(Bool_t mc, AliAnalysisManager* mgr)
+ {
+ if (!fUseCent) return;
+
+ gROOT->Macro("AddTaskCentrality.C");
+ AliCentralitySelectionTask* ctask =
+ dynamic_cast<AliCentralitySelectionTask*>(mgr->GetTask("CentralitySelection"));
+ if (!ctask) return;
+ ctask->SetPass(fESDPass);
+ if (mc) ctask->SetMCInput();
+ }
+ /**
+ * Crete output handler - we don't want one here.
+ *
+ * @return 0
+ */
+ AliVEventHandler* CreateOutputHandler(EType) { return 0; }
+ Bool_t fUseCent; // Whether to use centrality or not
+};
+
+//
+// EOF
+//
--- /dev/null
+#include "TrainSetup.C"
+
+//====================================================================
+/**
+ * Analysis train to make Forward and Central MC corrections
+ *
+ * To run, do
+ * @code
+ * gROOT->LoadMacro("TrainSetup.C");
+ * // Make train
+ * MakeMCCorrTrain t("My Analysis");
+ * // Set variaous parameters on the train
+ * t.SetDataDir("/home/of/data");
+ * t.AddRun(118506)
+ * // Run it
+ * t.Run("LOCAL", "FULL", -1, false, false);
+ * @endcode
+ *
+ * @ingroup pwg2_forward_scripts_makers
+ * @ingroup pwg2_forward_mc
+ */
+class MakeMCCorrTrain : public TrainSetup
+{
+public:
+ /**
+ * Constructor. Date and time must be specified when running this
+ * in Termiante mode on Grid
+ *
+ * @param name Name of train (free form)
+ * @param dateTime Append date and time to name
+ * @param year Year - if not specified, current year
+ * @param month Month - if not specified, current month
+ * @param day Day - if not specified, current day
+ * @param hour Hour - if not specified, current hour
+ * @param min Minutes - if not specified, current minutes
+ */
+ MakeMCCorrTrain(const char* name,
+ Bool_t dateTime = false,
+ UShort_t year = 0,
+ UShort_t month = 0,
+ UShort_t day = 0,
+ UShort_t hour = 0,
+ UShort_t min = 0)
+ : TrainSetup(name, dateTime,
+ year, month, day, hour, min)
+ {}
+ /**
+ * Run this analysis
+ *
+ * @param mode Mode - see TrainSetup::EMode
+ * @param oper Operation - see TrainSetup::EOperation
+ * @param nEvents Number of events (negative means all)
+ * @param usePar If true, use PARs
+ */
+ void Run(const char* mode, const char* oper,
+ Int_t nEvents=-1, Bool_t usePar=false)
+ {
+ Info("Run", "Running in mode=%s, oper=%s, events=%d, Par=%d",
+ mode, oper, nEvents, usePar);
+ Exec("ESD", mode, oper, nEvents, true, usePar);
+ }
+ /**
+ * Run this analysis
+ *
+ * @param mode Mode - see TrainSetup::EMode
+ * @param oper Operation - see TrainSetup::EOperation
+ * @param nEvents Number of events (negative means all)
+ * @param usePar If true, use PARs
+ */
+ void Run(EMode mode, EOper oper, Int_t nEvents=-1,
+ Bool_t usePar = false)
+ {
+ Info("Run", "Running in mode=%d, oper=%d, events=%d, Par=%d",
+ mode, oper, nEvents, usePar);
+ Exec(kESD, mode, oper, nEvents, true, usePar);
+ }
+protected:
+ /**
+ * Create the tasks
+ *
+ * @param mode Processing mode
+ * @param par Whether to use par files
+ * @param mgr Analysis manager
+ */
+ void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager* mgr)
+ {
+ // --- Output file name ------------------------------------------
+ AliAnalysisManager::SetCommonFileName("forward_mccorr.root");
+
+ // --- Load libraries/pars ---------------------------------------
+ LoadLibrary("PWG2forward2", mode, par, true);
+
+ // --- Set load path ---------------------------------------------
+ gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
+ gROOT->GetMacroPath()));
+
+ // --- Check if this is MC ---------------------------------------
+ if (!mgr->GetMCtruthEventHandler()) return;
+
+ // --- Add the task ----------------------------------------------
+ gROOT->Macro("AddTaskForwardMCCorr.C");
+
+ // --- Add the task ----------------------------------------------
+ gROOT->Macro("AddTaskCentralMCCorr.C");
+ }
+ //__________________________________________________________________
+ /**
+ * Create physics selection , and add to manager
+ *
+ * @param mc Whether this is for MC
+ * @param mgr Manager
+ */
+ void CreatePhysicsSelection(Bool_t mc,
+ AliAnalysisManager* mgr)
+ {
+ TrainSetup::CreatePhysicsSelection(mc, mgr);
+
+ // --- Get input event handler -----------------------------------
+ AliInputEventHandler* ih =
+ dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
+ if (!ih)
+ Fatal("CreatePhysicsSelection", "Couldn't get input handler (%p)", ih);
+
+ // --- Get Physics selection -------------------------------------
+ AliPhysicsSelection* ps =
+ dynamic_cast<AliPhysicsSelection*>(ih->GetEventSelection());
+ if (!ps)
+ Fatal("CreatePhysicsSelection", "Couldn't get PhysicsSelection (%p)", ps);
+
+ // --- Ignore trigger class when selecting events. This means ---
+ // --- that we get offline+(A,C,E) events too --------------------
+ // ps->SetSkipTriggerClassSelection(true);
+ }
+ /**
+ * Do not the centrality selection
+ */
+ void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {}
+ Double_t fVzMin; // Least v_z
+ Double_t fVzMax; // Largest v_z
+};
+
+//
+// EOF
+//
--- /dev/null
+#include "TrainSetup.C"
+
+//====================================================================
+/**
+ * Analysis train to make @f$ dN/d\eta@f$
+ *
+ * To run, do
+ * @code
+ * gROOT->LoadMacro("TrainSetup.C");
+ * // Make train
+ * MakedNdetaTrain t("My Analysis");
+ * // Set variaous parameters on the train
+ * t.SetDataDir("/home/of/data");
+ * t.AddRun(118506)
+ * // Run it
+ * t.Run("LOCAL", "FULL", -1, false, false);
+ * @endcode
+ *
+ * @ingroup pwg2_forward_scripts_makers
+ * @ingroup pwg2_forward_dndeta
+ */
+class MakedNdetaTrain : public TrainSetup
+{
+public:
+ /**
+ * Constructor. Date and time must be specified when running this
+ * in Termiante mode on Grid
+ *
+ * @param name Name of train (free form)
+ * @param trig Trigger to use
+ * @param vzMin Least @f$ v_z@f$
+ * @param vzMax Largest @f$ v_z@f$
+ * @param scheme Normalisation scheme
+ * @param useCent Whether to use centrality
+ * @param dateTime Append date and time to name
+ * @param year Year - if not specified, current year
+ * @param month Month - if not specified, current month
+ * @param day Day - if not specified, current day
+ * @param hour Hour - if not specified, current hour
+ * @param min Minutes - if not specified, current minutes
+ */
+ MakedNdetaTrain(const char* name,
+ const char* trig="INEL",
+ Double_t vzMin=-10,
+ Double_t vzMax=10,
+ const char* scheme="FULL",
+ Bool_t useCent=false,
+ Bool_t dateTime=false,
+ UShort_t year = 0,
+ UShort_t month = 0,
+ UShort_t day = 0,
+ UShort_t hour = 0,
+ UShort_t min = 0)
+ : TrainSetup(name, dateTime, year, month, day, hour, min),
+ fTrig(trig),
+ fVzMin(vzMin),
+ fVzMax(vzMax),
+ fScheme(scheme),
+ fUseCent(useCent)
+ {}
+ /**
+ * Run this analysis
+ *
+ * @param mode Mode - see TrainSetup::EMode
+ * @param oper Operation - see TrainSetup::EOperation
+ * @param nEvents Number of events (negative means all)
+ * @param usePar If true, use PARs
+ */
+ void Run(const char* mode, const char* oper,
+ Int_t nEvents=-1, Bool_t usePar=false)
+ {
+ Exec("AOD", mode, oper, nEvents, false, usePar);
+ }
+ /**
+ * Run this analysis
+ *
+ * @param mode Mode - see TrainSetup::EMode
+ * @param oper Operation - see TrainSetup::EOperation
+ * @param nEvents Number of events (negative means all)
+ * @param usePar If true, use PARs
+ */
+ void Run(EMode mode, EOper oper, Int_t nEvents=-1,
+ Bool_t usePar=false)
+ {
+ Exec(kAOD, mode, oper, nEvents, false, usePar);
+ }
+ /**
+ * Set the trigger to use (INEL, INEL>0, NSD)
+ *
+ * @param trig Trigger to use
+ */
+ void SetTrigger(const char* trig) { fTrig = trig; }
+ /**
+ * Set the vertex range to accept
+ *
+ * @param min Minimum
+ * @param max Maximum
+ */
+ void SetVertexRange(Double_t min, Double_t max) { fVzMin=min; fVzMax=max; }
+ /**
+ * Set the normalisation scheme
+ *
+ * @param scheme Normalisation scheme options
+ */
+ void SetScheme(const char* scheme) { fScheme = scheme; }
+ /**
+ * Whether to use centrality or not
+ *
+ * @param use To use the centrality
+ */
+ void SetUseCentrality(Bool_t use) { fUseCent = use; }
+protected:
+ /**
+ * Create the tasks
+ *
+ * @param mode Processing mode
+ * @param par Whether to use par files
+ */
+ void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager*)
+ {
+ // --- Output file name ------------------------------------------
+ AliAnalysisManager::SetCommonFileName("forward_dndeta.root");
+
+ // --- Load libraries/pars ---------------------------------------
+ LoadLibrary("PWG2forward2", mode, par, true);
+
+ // --- Set load path ---------------------------------------------
+ gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
+ gROOT->GetMacroPath()));
+
+ // --- Add the task ----------------------------------------------
+ gROOT->Macro(Form("AddTaskForwarddNdeta.C(\"%s\",%f,%f,%d,\"%s\")",
+ fTrig.Data(), fVzMin, fVzMax, fUseCent, fScheme.Data()));
+
+ gROOT->Macro(Form("AddTaskCentraldNdeta.C(\"%s\",%f,%f,%d,\"%s\")",
+ fTrig.Data(), fVzMin, fVzMax, fUseCent, fScheme.Data()));
+
+ gROOT->Macro(Form("AddTaskMCTruthdNdeta.C(\"%s\",%f,%f,%d,\"%s\")",
+ fTrig.Data(), fVzMin, fVzMax, fUseCent, fScheme.Data()));
+ }
+ /**
+ * Do not the centrality selection
+ */
+ void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {}
+ /**
+ * Crete output handler - we don't want one here.
+ *
+ * @return 0
+ */
+ AliVEventHandler* CreateOutputHandler(EType) { return 0; }
+ TString fTrig; // Trigger to use
+ Double_t fVzMin; // Least v_z
+ Double_t fVzMax; // Largest v_z
+ TString fScheme; // Normalisation scheme
+ Bool_t fUseCent; // Use centrality
+};
+//
+// EOF
+//
--- /dev/null
+/**
+ * @file TrainSetup.C
+ * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
+ * @date Wed Mar 23 12:12:00 2011
+ *
+ * @brief
+ *
+ * @ingroup pwg2_forward_scripts_makers
+ *
+ */
+
+#ifndef __CINT__
+#include <fstream>
+#include <iostream>
+
+#include <TAlienCollection.h>
+#include <TArrayI.h>
+#include <TChain.h>
+#include <TDatime.h>
+#include <TEnv.h>
+#include <TFile.h>
+#include <TGrid.h>
+#include <TList.h>
+#include <TObjString.h>
+#include <TProof.h>
+#include <TString.h>
+#include <TSystem.h>
+#include <TSystemDirectory.h>
+#include <TSystemFile.h>
+#include <TROOT.h>
+
+#include <AliAODHandler.h>
+#include <AliAODInputHandler.h>
+#include <AliAnalysisDataContainer.h>
+#include <AliAnalysisManager.h>
+#include <AliAnalysisAlien.h>
+#include <AliESDInputHandler.h>
+#include <AliMCEventHandler.h>
+#include <AliVEventHandler.h>
+#include <AliPhysicsSelection.h>
+#include <AliCentralitySelectionTask.h>
+#else
+class TArrayI;
+class TChain;
+class AliAnalysisManager;
+#endif
+
+//====================================================================
+/**
+ * Generic set-up of an analysis train using the grid-handler (AliEn plugin).
+ *
+ * Users should define a class that derives from this. The class
+ * should implement the member function CreateTasks to add needed
+ * tasks to the train
+ *
+ * @code
+ * // MyTrain.C
+ * class MyTrain : public TrainSetup
+ * {
+ * public:
+ * MyTrain(Bool_t dateTime = false,
+ * UShort_t year = 0,
+ * UShort_t month = 0,
+ * UShort_t day = 0,
+ * UShort_t hour = 0,
+ * UShort_t min = 0)
+ * : TrainSetup("My train", dateTime, year, month, day, hour, min)
+ * {}
+ * void Run(const char* type, const char* mode, const char* oper,
+ * Int_t nEvents=-1, Bool_t mc=false,
+ * Bool_t usePar=false)
+ * {
+ * Exec(type, mode, oper, nEvents, mc, usePar);
+ * }
+ * protected:
+ * void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager* mgr)
+ * {
+ * AliAnalysisManager::SetCommonFileName("my_analysis.root");
+ * LoadLibrary("MyAnalysis", mode, par, true);
+ * Bool_t mc = mgr->GetMCtruthEventHandler() != 0;
+ * gROOT->Macro("MyAnalysis.C");
+ * }
+ * };
+ * @endcode
+ *
+ * This can then be run like
+ *
+ * @verbatim
+ * > aliroot
+ * Root> .L TrainSetup.C
+ * Root> .L MyTrain.C
+ * Root> MyTrain t;
+ * Root> t.Run();
+ * @endverbatim
+ *
+ * or as a script
+ *
+ * @code
+ * {
+ * gROOT->LoadMacro("TrainSetup.C");
+ * gROOT->LoadMacro("MyTrain.C");
+ * MyTrain t;
+ * t.Run();
+ * }
+ * @endcode
+ *
+ * To byte compile this, you need to
+ * - load the ROOT AliEn library
+ * - load the analysis libraries
+ * - add $ALICE_ROOT/include to header search
+ * first
+ *
+ * @verbatim
+ * > aliroot
+ * Root> gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
+ * Root> "$ALICE_ROOT/ANALYSIS/macros",
+ * Root> gROOT->GetMacroPath()));
+ * Root> gSystem->AddIncludePath("-I${ALICE_ROOT}/include");
+ * Root> gSystem->Load("libRAliEn");
+ * Root> gSystem->Load("libANALYSIS");
+ * Root> gSystem->Load("libANALYSISalice");
+ * Root> gROOT->LoadMacro("TrainSetup.C+");
+ * @endverbatim
+ *
+ * @ingroup pwg2_forward_scripts_makers
+ *
+ */
+struct TrainSetup
+{
+ /**
+ * Data type to process
+ */
+ enum EType {
+ /** Event Summary Data */
+ kESD,
+ /** Analysis Object Data */
+ kAOD
+ };
+ /**
+ * How to run the analysis
+ *
+ */
+ enum EMode {
+ /** Locally */
+ kLocal = 1,
+ /** In PROOF(-Lite) cluster */
+ kProof,
+ /** On the Grid */
+ kGrid
+ };
+ /**
+ * What stage of the analysis to run
+ *
+ */
+ enum EOper {
+ /** Testing. Local processing, a single copied from Grid */
+ kTest,
+ /** Off-line */
+ kOffline,
+ /** Submit to queue */
+ kSubmit,
+ /** Merge and terminate */
+ kTerminate,
+ /** Full run */
+ kFull
+ };
+
+ //__________________________________________________________________
+ /**
+ * Constructor
+ *
+ * @param name Name of analysis (free-form)
+ * @param useDateTime Whether to append date and time to the name
+ * @param year Year - if not specified, taken from current date
+ * @param month Month - if not specified, taken from current date
+ * @param day Day - if not specified, taken from current date
+ * @param hour Hour - if not specified, taken from current time
+ * @param min Minute - if not specified, taken from current time
+ */
+ TrainSetup(const char* name, Bool_t useDateTime=true,
+ UShort_t year=0, UShort_t month=0,
+ UShort_t day=0, UShort_t hour=0, UShort_t min=0)
+ : fName(name),
+ fRootVersion("v5-28-00a"),
+ fAliRootVersion("v4-21-18-AN"),
+ fAliEnAPIVersion("V1.1x"),
+ fProofServer("alicecaf.cern.ch"),
+ fDataDir("/alice/data/2010/LHC10c"),
+ fDataSet("/COMMON/COMMON/LHC09a4_run8100X#/esdTree"),
+ fXML(""),
+ fRunNumbers(0),
+ fListOfPARs(),
+ fListOfSources(),
+ fListOfLibraries(),
+ fListOfExtras(),
+ fNReplica(4),
+ fESDPass(3),
+ fPassPostfix(""),
+ fEscapedName(name),
+ fAllowOverwrite(kFALSE)
+ {
+ char c[] = { ' ', '/', '@', 0 };
+ char* p = c;
+ while (*p) {
+ fEscapedName.ReplaceAll(Form("%c", *p), "_");
+ p++;
+ }
+
+ if (useDateTime) {
+ if (year == 0 || month == 0 || day == 0) {
+ TDatime now;
+ year = now.GetYear();
+ month = now.GetMonth();
+ day = now.GetDay();
+ hour = now.GetHour();
+ min = now.GetMinute();
+ }
+ fEscapedName.Append(Form("_%04d%02d%02d_%02d%02d",
+ year, month, day, hour, min));
+ }
+
+ }
+
+ //__________________________________________________________________
+ /**
+ * Parse a string into a type enum
+ *
+ * @param type String to pass
+ *
+ * @return Enumaration value
+ */
+ static EType ParseType(const char* type, Bool_t& /*mc*/)
+ {
+ // mc = false;
+ TString sType(type);
+ sType.ToUpper();
+ EType eType = kESD;
+ // if (sType.Contains("MC")) mc = true;
+ if (sType.Contains("ESD")) eType = kESD;
+ else if (sType.Contains("AOD")) eType = kAOD;
+ else
+ Fatal("Run", "Unknown type '%s'", type);
+
+ return eType;
+ }
+ //__________________________________________________________________
+ /**
+ * Return a string that reflects the passed mode
+ *
+ * @param eMode Mode
+ *
+ * @return String representation of mode
+ */
+ static const char* ModeString(EMode eMode)
+ {
+ switch (eMode) {
+ case kLocal: return "LOCAL";
+ case kProof: return "PROOF";
+ case kGrid: return "GRID";
+ }
+ return 0;
+ }
+ //__________________________________________________________________
+ /**
+ * Parse a string for mode specifier
+ *
+ * @param mode Mode string
+ *
+ * @return EMode value
+ */
+ static EMode ParseMode(const char* mode)
+ {
+ TString sMode(mode);
+ sMode.ToUpper();
+ EMode eMode = kLocal;
+ if (sMode == "LOCAL") eMode = kLocal;
+ else if (sMode == "PROOF") eMode = kProof;
+ else if (sMode == "GRID") eMode = kGrid;
+ else
+ Fatal("Run", "Unknown mode '%s'", mode);
+ return eMode;
+ }
+
+ //__________________________________________________________________
+ /**
+ * Return a string that reflects the passed operation
+ *
+ * @param eOper Operation
+ *
+ * @return String representation of operation
+ */
+ static const char* OperString(EOper eOper)
+ {
+ switch (eOper) {
+ case kTest: return "TEST";
+ case kOffline: return "OFFLINE";
+ case kSubmit: return "SUBMIT";
+ case kTerminate: return "TERMINATE";
+ case kFull: return "FULL";
+ }
+ return 0;
+ }
+ //__________________________________________________________________
+ /**
+ * Parse an operation string
+ *
+ * @param oper Operation
+ *
+ * @return An EOper value
+ */
+ static EOper ParseOperation(const char* oper)
+ {
+ TString sOper(oper);
+ sOper.ToUpper();
+ EOper eOper = kFull;
+ if (sOper == "TEST") eOper = kTest;
+ else if (sOper == "OFFLINE") eOper = kOffline;
+ else if (sOper == "SUBMIT") eOper = kSubmit;
+ else if (sOper == "TERMINATE") eOper = kTerminate;
+ else if (sOper == "FULL") eOper = kFull;
+ else
+ Fatal("Run", "unknown operation '%s'", oper);
+ return eOper;
+ }
+
+ //__________________________________________________________________
+ /**
+ * Set ROOT version to use
+ *
+ * @param v Version string of ROOT
+ */
+ void SetROOTVersion(const char* v) { fRootVersion = v; }
+ //__________________________________________________________________
+ /**
+ * Set AliROOT version to use
+ *
+ * @param v Version string of AliROOT
+ */
+ void SetAliROOTVersion(const char* v) { fAliRootVersion = v; }
+ //__________________________________________________________________
+ /**
+ * Set the PROOF server URL
+ *
+ * @param s PROOF server URL
+ */
+ void SetProofServer(const char* s) { fProofServer = s; }
+ //__________________________________________________________________
+ /**
+ * Set the GRID/Local data dir
+ *
+ * @param d Directory with data
+ */
+ void SetDataDir(const char* d) { fDataDir = d; }
+ //__________________________________________________________________
+ /**
+ * Set the PROOF data set
+ *
+ * @param d PROOF registered data set
+ */
+ void SetDataSet(const char* d) { fDataSet = d; }
+ //__________________________________________________________________
+ /**
+ * Set the XML file to use
+ *
+ * @param x XML file
+ */
+ void SetXML(const char* x) { fXML = x; }
+ //__________________________________________________________________
+ /**
+ * Set how many replicas of the output we want
+ *
+ * @param n Number of replicas requested
+ */
+ void SetNReplica(Int_t n) { fNReplica = n; }
+ /**
+ * Set the ESD pass to use
+ *
+ * @param pass Pass number
+ */
+ void SetESDPass(Int_t pass) { fESDPass = pass; }
+ /**
+ * Set the ESD pass to use
+ *
+ * @param postfix Post fix to pass number
+ */
+ void SetPassPostfix(const char* postfix) { fPassPostfix = postfix; }
+ //__________________________________________________________________
+ /**
+ * Add a source file to be copied and byte compiled on slaves
+ *
+ * @param src Sources
+ * @param addToExtra If false, do not copy
+ */
+ void AddSource(const char* src, bool addToExtra=true)
+ {
+ fListOfSources.Add(new TObjString(src));
+ if (addToExtra) AddExtraFile(src); // Source code isn't copied!
+ }
+ //__________________________________________________________________
+ /**
+ * Add binary data to be uploaded to slaves
+ *
+ * @param lib Name of binary file
+ */
+ void AddLibrary(const char* lib) { fListOfLibraries.Add(new TObjString(lib));}
+ //__________________________________________________________________
+ /**
+ * Add a run to be analysed
+ *
+ * @param run Run number
+ */
+ void AddRun(Int_t run)
+ {
+ Int_t i = fRunNumbers.fN; fRunNumbers.Set(i+1); fRunNumbers[i] = run;
+ }
+ //__________________________________________________________________
+ /**
+ * Read run numbers from a file
+ *
+ * @param filename File name
+ */
+ void ReadRunNumbers(const char* filename)
+ {
+ std::ifstream file(filename);
+ if (!file)
+ Fatal("ReadRunNumbers", "Cannot read from %s", filename);
+
+ while (!file.eof()) {
+ Int_t run;
+ file >> run;
+ AddRun(run);
+ Char_t c;
+ file >> c;
+ if (file.bad()) break;
+ }
+ file.close();
+ }
+ //__________________________________________________________________
+ /**
+ * Add an extra file to be uploaded to slave
+ *
+ * @param file Extra file to be uploaded
+ */
+ void AddExtraFile(const char* file)
+ {
+ if (!file || file[0] == '\0') return;
+ fListOfExtras.Add(new TObjString(file));
+ }
+ //__________________________________________________________________
+ /**
+ * Set whether to allow overwritting existing files/directories
+ *
+ * @param allow If true, allow overwritting files/directories
+ */
+ void SetAllowOverwrite(Bool_t allow) { fAllowOverwrite = allow; }
+ //__________________________________________________________________
+ /**
+ * Print the setup
+ *
+ */
+ void Print() const
+ {
+ bool mc = AliAnalysisManager::GetAnalysisManager()
+ ->GetMCtruthEventHandler();
+ std::cout << fName << " train setup\n"
+ << std::boolalpha
+ << " ROOT version: " << fRootVersion << "\n"
+ << " AliROOT version: " << fAliRootVersion << "\n"
+ << " Name of proof server: " << fProofServer << "\n"
+ << " Grid Input directory: " << fDataDir << "\n"
+ << " Proof data set name: " << fDataSet << "\n"
+ << " XML collection: " << fXML << "\n"
+ << " Monte-Carlo input: " << mc << "\n"
+ << " Storage replication: " << fNReplica << "\n"
+ << " Run numbers: " << std::flush;
+ for (Int_t i = 0; i < fRunNumbers.GetSize(); i++)
+ std::cout << (i == 0 ? "" : ", ") << fRunNumbers.At(i);
+
+ std::cout << "\n"
+ << " PAR files: " << std::flush;
+ Bool_t first = true;
+ TObject* obj = 0;
+ TIter nextPar(&fListOfPARs);
+ while ((obj = nextPar())) {
+ std::cout << (first ? "" : ", ") << obj->GetName();
+ first = false;
+ }
+
+ std::cout << "\n"
+ << " Script sources: " << std::flush;
+ first = true;
+ TIter nextSrc(&fListOfSources);
+ while ((obj = nextSrc())) {
+ std::cout << (first ? "" : ", ") << obj->GetName();
+ first = false;
+ }
+
+ std::cout << "\n"
+ << " Libraries to load: " << std::flush;
+ first = true;
+ TIter nextLib(&fListOfLibraries);
+ while ((obj = nextLib())) {
+ std::cout << (first ? "" : ", ") << obj->GetName();
+ first = false;
+ }
+ std::cout << std::noboolalpha << std::endl;
+
+ AliAnalysisGrid* plugin =
+ AliAnalysisManager::GetAnalysisManager()->GetGridHandler();
+ if (!plugin) return;
+
+ }
+
+protected:
+ //__________________________________________________________________
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ TrainSetup(const TrainSetup& o)
+ : fName(o.fName),
+ fRootVersion(o.fRootVersion),
+ fAliRootVersion(o.fAliRootVersion),
+ fProofServer(o.fProofServer),
+ fDataDir(o.fDataDir),
+ fDataSet(o.fDataSet),
+ fXML(o.fXML),
+ fRunNumbers(o.fRunNumbers),
+ fListOfPARs(),
+ fListOfSources(),
+ fListOfLibraries(),
+ fListOfExtras(),
+ fNReplica(o.fNReplica),
+ fESDPass(o.fESDPass)
+ {
+ if (isdigit(fName[0])) {
+ Warning("TrainSetup", "Name starts with a digit, prepending 'a' to name");
+ fName = Form("a%s", fName.Data());
+ }
+ TObject* obj = 0;
+ TIter nextPar(&o.fListOfPARs);
+ while ((obj = nextPar())) fListOfPARs.Add(obj->Clone());
+ TIter nextSrc(&o.fListOfSources);
+ while ((obj = nextSrc())) fListOfSources.Add(obj->Clone());
+ TIter nextLib(&o.fListOfLibraries);
+ while ((obj = nextLib())) fListOfLibraries.Add(obj->Clone());
+ TIter nextExa(&o.fListOfExtras);
+ while ((obj = nextExa())) fListOfExtras.Add(obj->Clone());
+ }
+ //__________________________________________________________________
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object.
+ */
+ TrainSetup& operator=(const TrainSetup& o)
+ {
+ fName = o.fName;
+ fRootVersion = o.fRootVersion;
+ fAliRootVersion = o.fAliRootVersion;
+ fProofServer = o.fProofServer;
+ fDataDir = o.fDataDir;
+ fDataSet = o.fDataSet;
+ fXML = o.fXML;
+ fNReplica = o.fNReplica;
+ fESDPass = o.fESDPass;
+ fRunNumbers = o.fRunNumbers;
+ TObject* obj = 0;
+ TIter nextPar(&o.fListOfPARs);
+ while ((obj = nextPar())) fListOfPARs.Add(obj->Clone());
+ TIter nextSrc(&o.fListOfSources);
+ while ((obj = nextSrc())) fListOfSources.Add(obj->Clone());
+ TIter nextLib(&o.fListOfLibraries);
+ while ((obj = nextLib())) fListOfLibraries.Add(obj->Clone());
+ TIter nextExa(&o.fListOfExtras);
+ while ((obj = nextExa())) fListOfExtras.Add(obj->Clone());
+
+ return *this;
+ }
+
+ //__________________________________________________________________
+ /**
+ * Run this analysis
+ *
+ * @param type Type of input for analysis (kESD, kAOD)
+ * @param mode Mode of job (kLocal, kProof, kGrid)
+ * @param oper Operation
+ * @param nEvents Number of events to analyse (<0 means all)
+ * @param mc Whether to connect MC data
+ * @param usePar Whether to use PARs
+ * @param dbg Debug level
+ */
+ void Exec(const char* type,
+ const char* mode="GRID",
+ const char* oper="FULL",
+ Int_t nEvents=-1,
+ Bool_t mc=false,
+ Bool_t usePar=false,
+ Int_t dbg=0)
+ {
+ Info("Exec", "Doing exec with type=%s, mode=%s, oper=%s, events=%d "
+ "mc=%d, usePar=%d", type, mode, oper, nEvents, mc, usePar);
+ EType eType = ParseType(type, mc);
+ EMode eMode = ParseMode(mode);
+ EOper eOper = ParseOperation(oper);
+
+ Exec(eType, eMode, eOper, nEvents, mc, usePar, dbg);
+ }
+
+ //__________________________________________________________________
+ /**
+ * Run this analysis
+ *
+ * @param type Type of input for analysis (kESD, kAOD)
+ * @param mode Mode of job (kLocal, kProof, kGrid)
+ * @param oper Operation
+ * @param nEvents Number of events to analyse (<0 means all)
+ * @param mc Whether to connect MC data
+ * @param usePar Whether to use PARs
+ * @param dbg Debug level
+ */
+ void Exec(EType type,
+ EMode mode,
+ EOper oper,
+ Int_t nEvents,
+ Bool_t mc,
+ Bool_t usePar,
+ Int_t dbg=0)
+ {
+ Info("Exec", "Doing exec with type=%d, mode=%d, oper=%d, events=%d "
+ "mc=%d, usePar=%d", type, mode, oper, nEvents, mc, usePar);
+
+ if (mode == kProof) usePar = true;
+
+ if (!Connect(mode)) return;
+
+ TString cwd = gSystem->WorkingDirectory();
+ TString nam = EscapedName();
+ if (oper != kTerminate) {
+ if (!fAllowOverwrite && !gSystem->AccessPathName(nam.Data())) {
+ Error("Exec", "File/directory %s already exists", nam.Data());
+ return;
+ }
+ if (gSystem->AccessPathName(nam.Data())) {
+ if (gSystem->MakeDirectory(nam.Data())) {
+ Error("Exec", "Failed to make directory %s", nam.Data());
+ return;
+ }
+ }
+ }
+ else {
+ if (gSystem->AccessPathName(nam.Data())) {
+ Error("Exec", "File/directory %s does not exists", nam.Data());
+ return;
+ }
+ }
+
+ if (!gSystem->ChangeDirectory(nam.Data())) {
+ Error("Exec", "Failed to change directory to %s", nam.Data());
+ return;
+ }
+ Info("Exec", "Made subdirectory %s, and cd'ed there", nam.Data());
+
+ if (!LoadCommonLibraries(mode, usePar)) return;
+
+ // --- Create analysis manager -----------------------------------
+ AliAnalysisManager *mgr = new AliAnalysisManager(fName,"Analysis Train");
+
+ // In test mode, collect system information on every event
+ // if (oper == kTest) mgr->SetNSysInfo(1);
+ if (dbg > 0) mgr->SetDebugLevel(dbg);
+ if (mode == kLocal) mgr->SetUseProgressBar(kTRUE, 100);
+
+ // --- ESD input handler ------------------------------------------
+ AliVEventHandler* inputHandler = CreateInputHandler(type);
+ if (inputHandler) mgr->SetInputEventHandler(inputHandler);
+
+ // --- Monte-Carlo ------------------------------------------------
+ AliVEventHandler* mcHandler = CreateMCHandler(type,mc);
+ if (mcHandler) mgr->SetMCtruthEventHandler(mcHandler);
+
+ // --- AOD output handler -----------------------------------------
+ AliVEventHandler* outputHandler = CreateOutputHandler(type);
+ if (outputHandler) mgr->SetOutputEventHandler(outputHandler);
+
+ // --- Include analysis macro path in search path ----------------
+ gROOT->SetMacroPath(Form("%s:%s:$ALICE_ROOT/ANALYSIS/macros",
+ cwd.Data(), gROOT->GetMacroPath()));
+
+ // --- Physics selction - only for ESD ---------------------------
+ if (type == kESD) CreatePhysicsSelection(mc, mgr);
+
+ // --- Create centrality task ------------------------------------
+ CreateCentralitySelection(mc, mgr);
+
+ // --- Create tasks ----------------------------------------------
+ CreateTasks(mode, usePar, mgr);
+
+ // --- Create Grid handler ----------------------------------------
+ // _must_ be done after all tasks has been added
+ AliAnalysisAlien* gridHandler = CreateGridHandler(type, mode, oper);
+ if (gridHandler) mgr->SetGridHandler(gridHandler);
+
+ // --- Create the chain ------------------------------------------
+ TChain* chain = CreateChain(type, mode, oper, mc);
+ if (mode == kLocal && !chain) {
+ Error("Exec", "No chain defined in local mode!");
+ return;
+ }
+
+ // --- Print setup -----------------------------------------------
+ Print();
+
+ // --- Initialise the train --------------------------------------
+ if (!mgr->InitAnalysis()) {
+ gSystem->ChangeDirectory(cwd.Data());
+ Error("Run","Failed to initialise train");
+ return;
+ }
+
+ // --- Show status -----------------------------------------------
+ mgr->PrintStatus();
+
+ Long64_t ret = StartAnalysis(mgr, mode, chain, nEvents);
+
+ // Make sure we go back
+ gSystem->ChangeDirectory(cwd.Data());
+
+ if (ret < 0) Error("Exec", "Analysis failed");
+ }
+ //__________________________________________________________________
+ /**
+ * Start the analysis
+ *
+ * @param mgr Analysis manager
+ * @param mode Run mode
+ * @param chain Input data (local and proof only)
+ * @param nEvents Number of events to analyse
+ */
+ Long64_t StartAnalysis(AliAnalysisManager* mgr,
+ EMode mode,
+ TChain* chain,
+ Int_t nEvents)
+ {
+ // --- Run the analysis ------------------------------------------
+ switch (mode) {
+ case kLocal:
+ if (!chain) {
+ Error("StartAnalysis", "No chain defined");
+ return -1;
+ }
+ if (nEvents < 0) nEvents = chain->GetEntries();
+ return mgr->StartAnalysis(ModeString(mode), chain, nEvents);
+ case kProof:
+ if (fDataSet.IsNull()) {
+ if (!chain) {
+ Error("StartAnalysis", "No chain defined");
+ return -1;
+ }
+ if (nEvents < 0) nEvents = chain->GetEntries();
+ return mgr->StartAnalysis(ModeString(mode), chain, nEvents);
+ }
+ return mgr->StartAnalysis(ModeString(mode), fDataSet);
+ case kGrid:
+ if (nEvents < 0)
+ return mgr->StartAnalysis(ModeString(mode));
+ return mgr->StartAnalysis(ModeString(mode), nEvents);
+ }
+ // We should never get here
+ return -1;
+ }
+ //__________________________________________________________________
+ /**
+ * Return the escaped name
+ *
+ *
+ * @return Escaped name
+ */
+ const TString& EscapedName() const
+ {
+ return fEscapedName;
+ }
+ //__________________________________________________________________
+ /**
+ * Create a grid handler
+ *
+ * @param type Data type
+ * @param mode Run mode
+ * @param oper Operation
+ *
+ * @return Grid handler
+ */
+ virtual AliAnalysisAlien*
+ CreateGridHandler(EType type, EMode mode, EOper oper)
+ {
+ if (mode != kGrid) return 0;
+
+ TString name = EscapedName();
+
+ // Create the plug-in object, and set run mode
+ AliAnalysisAlien* plugin = new AliAnalysisAlien();
+ plugin->SetRunMode(OperString(oper));
+
+ // Production mode - not used here
+ // plugin->SetProductionMode();
+
+ // Set output to be per run
+ plugin->SetOutputToRunNo();
+
+ // Set the job tag
+ plugin->SetJobTag(fName);
+
+ // Set number of test files - used in test mode only
+ plugin->SetNtestFiles(1);
+
+ // Set required version of software
+ plugin->SetAPIVersion(fAliEnAPIVersion);
+ plugin->SetROOTVersion(fRootVersion);
+ plugin->SetAliROOTVersion(fAliRootVersion);
+
+ // Keep log files
+ plugin->SetKeepLogs();
+
+ // Declare root of input data directory
+ plugin->SetGridDataDir(fDataDir);
+
+ // Data search patterns
+ TString pat;
+ if (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) {
+ pat = "*/";
+ plugin->SetRunPrefix("");
+ }
+ else {
+ pat = Form("*ESDs/pass%d%s/*/", fESDPass, fPassPostfix.Data());
+ plugin->SetRunPrefix("000");
+ }
+ pat.Append(Form("*%s.root", type == kESD ? "ESDs" : "AOD"));
+ plugin->SetDataPattern(pat);
+
+ // Add the run numbers
+ for (Int_t i = 0; i < fRunNumbers.fN; i++) {
+ if (fRunNumbers[i] < 0) continue;
+ plugin->AddRunNumber(fRunNumbers[i]);
+ }
+
+ // Set the working directory to be the trains name (with special
+ // characters replaced by '_' and the date appended), and also set
+ // the output directory (relative to working directory)
+ plugin->SetGridWorkingDir(name.Data());
+ plugin->SetGridOutputDir("output");
+
+ // Enable configured PARs
+ TIter nextPar(&fListOfPARs);
+ TObject* parName;
+ while ((parName = nextPar()))
+ plugin->EnablePackage(parName->GetName());
+
+ // Add sources that need to be compiled on the workers using
+ // AcLIC.
+ TString addSources = SetupSources();
+ if (!addSources.IsNull()) plugin->SetAnalysisSource(addSources.Data());
+
+ // Add binary libraries that should be uploaded to the workers
+ TString addLibs = SetupLibraries();
+ if (!addLibs.IsNull()) plugin->SetAdditionalLibs(addLibs.Data());
+
+ // Disable default outputs
+ plugin->SetDefaultOutputs(true);
+
+ // Merge parameters
+ plugin->SetMaxMergeFiles(20);
+ plugin->SetMergeExcludes("AliAOD.root "
+ "*EventStat*.root "
+ "*event_stat*.root");
+
+ // Set number of runs per master - set to one to per run
+ plugin->SetNrunsPerMaster(1);
+
+ // Loop over defined containers in the analysis manager,
+ // and declare these as outputs
+ TString listOfAODs = "";
+ TString listOfHists = "";
+ AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
+ AliAnalysisDataContainer* cont = 0;
+ TIter nextCont(mgr->GetOutputs());
+ while ((cont = static_cast<AliAnalysisDataContainer*>(nextCont()))) {
+ TString outName(cont->GetFileName());
+ TString& list = (outName == "default" ? listOfAODs : listOfHists);
+ if (outName == "default") {
+ if (!mgr->GetOutputEventHandler()) continue;
+
+ outName = mgr->GetOutputEventHandler()->GetOutputFileName();
+ }
+ if (list.Contains(outName)) continue;
+ if (!list.IsNull()) list.Append(",");
+ list.Append(outName);
+ }
+ if (!mgr->GetExtraFiles().IsNull()) {
+ if (!listOfAODs.IsNull()) listOfAODs.Append("+");
+ TString extra = mgr->GetExtraFiles();
+ extra.ReplaceAll(" ", ",");
+ listOfAODs.Append(extra);
+ }
+ TString outArchive = Form("stderr, stdout@disk=%d", fNReplica);
+ if (!listOfHists.IsNull())
+ outArchive.Append(Form(" hist_archive.zip:%s@disk=%d",
+ listOfHists.Data(), fNReplica));
+ if (!listOfAODs.IsNull())
+ outArchive.Append(Form(" aod_archive.zip:%s@disk=%d",
+ listOfAODs.Data(), fNReplica));
+ if (listOfAODs.IsNull() && listOfHists.IsNull())
+ Fatal("CreateGridHandler", "No outputs defined");
+ // Disabled for now
+ // plugin->SetOutputArchive(outArchive);
+
+ // Set name of generated analysis macro
+ plugin->SetAnalysisMacro(Form("%s.C", name.Data()));
+
+ // Maximum number of sub-jobs
+ // plugin->SetSplitMaxInputFileNumber(25);
+
+ // Set the Time-To-Live
+ plugin->SetTTL(70000);
+
+ // Re-submit failed jobs as long as the ratio of failed jobs is
+ // below this percentage.
+ plugin->SetMasterResubmitThreshold(95);
+
+ // Set the input format
+ plugin->SetInputFormat("xml-single");
+
+ // Set the name of the generated jdl
+ plugin->SetJDLName(Form("%s.jdl", name.Data()));
+
+ // Set the name of the generated executable
+ plugin->SetExecutable(Form("%s.sh", name.Data()));
+
+ // Set the job price !?
+ plugin->SetPrice(1);
+
+ // Set whether to merge via JDL
+ plugin->SetMergeViaJDL(true);
+
+ // Fast read otion
+ plugin->SetFastReadOption(false);
+
+ // Whether to overwrite existing output
+ plugin->SetOverwriteMode(true);
+
+ // Set the executable binary name and options
+ plugin->SetExecutableCommand("aliroot -b -q -x");
+
+ // Split by storage element - must be lower case!
+ plugin->SetSplitMode("se");
+
+ return plugin;
+ }
+ //__________________________________________________________________
+ /**
+ * Create input handler
+ *
+ * @param type
+ *
+ * @return
+ */
+ virtual AliVEventHandler* CreateInputHandler(EType type)
+ {
+ switch (type) {
+ case kESD: return new AliESDInputHandler();
+ case kAOD: return new AliAODInputHandler();
+ }
+ return 0;
+ }
+ //__________________________________________________________________
+ /**
+ * Create input handler
+ *
+ * @param type Run type (ESD or AOD)
+ * @param mc Assume monte-carlo input
+ *
+ * @return
+ */
+ virtual AliVEventHandler* CreateMCHandler(EType type, bool mc)
+ {
+ if (!mc) return 0;
+ if (type != kESD) return 0;
+ Info("CreateMCHandler", "Making MC handler");
+ AliMCEventHandler* mcHandler = new AliMCEventHandler();
+ mcHandler->SetReadTR(true);
+ return mcHandler;
+ }
+ //__________________________________________________________________
+ /**
+ * Create output event handler
+ *
+ * @param type
+ *
+ * @return
+ */
+ virtual AliVEventHandler* CreateOutputHandler(EType type)
+ {
+ AliAODHandler* ret = new AliAODHandler();
+ switch (type) {
+ case kESD:
+ ret->SetOutputFileName("AliAOD.root");
+ break;
+ case kAOD:
+ ret->SetOutputFileName("AliAOD.pass2.root");
+ break;
+ }
+ return ret;
+ }
+ //__________________________________________________________________
+ /**
+ * Create physics selection , and add to manager
+ *
+ * @param mc Whether this is for MC
+ * @param mgr Manager
+ */
+ virtual void CreatePhysicsSelection(Bool_t mc,
+ AliAnalysisManager* mgr)
+ {
+ gROOT->Macro(Form("AddTaskPhysicsSelection.C(%d)", mc));
+ mgr->RegisterExtraFile("event_stat.root");
+ }
+ //__________________________________________________________________
+ /**
+ * Create physics selection , and add to manager
+ *
+ * @param mc Whether this is for MC
+ * @param mgr Manager
+ */
+ virtual void CreateCentralitySelection(Bool_t mc, AliAnalysisManager* mgr)
+ {
+ gROOT->Macro("AddTaskCentrality.C");
+ AliCentralitySelectionTask* ctask =
+ dynamic_cast<AliCentralitySelectionTask*>(mgr->GetTask("CentralitySelection"));
+ if (!ctask) return;
+ ctask->SetPass(fESDPass);
+ if (mc) ctask->SetMCInput();
+ }
+ //__________________________________________________________________
+ /**
+ * Create analysis tasks
+ *
+ * @param mode Run mode
+ * @param mgr Manager
+ * @param par Whether to use pars
+ */
+ virtual void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager* mgr)=0;
+ //__________________________________________________________________
+ /**
+ * Connect to external services (Proof and/or grid)
+ *
+ * @param mode Running mode
+ *
+ * @return true on success
+ */
+ virtual Bool_t Connect(EMode mode)
+ {
+ if (mode == kLocal) return true;
+
+ // --- Set-up connections to Proof cluster and alien -------------
+ if (mode == kProof) {
+ // --- Find user name ------------------------------------------
+ TString userName(gSystem->Getenv("alien_API_USER"));
+ if (userName.IsNull()) {
+ userName = gSystem->GetUserInfo()->fUser;
+ Warning("Connect",
+ "environment variable 'alien_API_USER' not set, using %s",
+ userName.Data());
+ }
+
+ // --- Set prefered GSI method ---------------------------------
+ gEnv->SetValue("XSec.GSI.DelegProxy", "2");
+
+ // --- Now open connection to PROOF cluster --------------------
+ TString serv = "";
+ Bool_t lite = false;
+ if (fProofServer.BeginsWith("workers=") || fProofServer.IsNull()) {
+ lite = true;
+ serv = fProofServer;
+ }
+ else
+ serv = Form("%s@%s", userName.Data(), fProofServer.Data());
+ TProof::Open(serv);
+ if (!gProof) {
+ Error("Connect", "Failed to connect to Proof cluster %s as %s",
+ fProofServer.Data(), userName.Data());
+ return false;
+ }
+ if (lite) return true;
+ }
+
+ // --- Open a connection to the grid -----------------------------
+ TGrid::Connect("alien://");
+ if (!gGrid || !gGrid->IsConnected()) {
+ // This is only fatal in grid mode
+ Error("Connect", "Failed to connect to AliEN");
+ if (mode == kGrid) return false;
+ return true;
+ }
+ if (mode == kGrid) return true;
+
+
+ // --- Set and make output directory -----------------------------
+ TString name = EscapedName();
+ TString homeDir(gGrid->GetHomeDirectory());
+ TString workDir(homeDir);
+ workDir.Append("/");
+ workDir.Append(name);
+
+ // Make working directory
+ if (!gGrid->Cd(workDir)) {
+ gGrid->Cd(homeDir);
+ if (gGrid->Mkdir(workDir)) {
+ gGrid->Cd(name);
+ Info("Connect", "Directory %s created", workDir.Data());
+ }
+ }
+ // Make output directory
+ gGrid->Mkdir("proof_output");
+ gGrid->Cd("proof_output");
+
+ return true;
+ }
+ //__________________________________________________________________
+ /**
+ * Load common libraries
+ *
+ * @param mode Running mode
+ * @param par If true, load as PARs
+ *
+ * @return true on success
+ */
+ Bool_t LoadCommonLibraries(EMode mode, Bool_t par)
+ {
+ if (!gSystem->Getenv("ALICE_ROOT")) {
+ Error("LoadCommonLibraries", "Local AliROOT not available");
+ return false;
+ }
+ gSystem->Load("libTree.so");
+ gSystem->Load("libGeom.so");
+ gSystem->Load("libVMC.so");
+ gSystem->Load("libPhysics.so");
+ gSystem->Load("libMinuit.so");
+
+ Bool_t ret = true;
+ Bool_t basic = mode == kGrid ? false : par;
+
+ ret &= LoadLibrary("STEERBase", mode, basic, false);
+ ret &= LoadLibrary("ESD", mode, basic, false);
+ ret &= LoadLibrary("AOD", mode, basic, false);
+ ret &= LoadLibrary("ANALYSIS", mode, basic, true);
+ ret &= LoadLibrary("ANALYSISalice", mode, basic, true);
+
+ return ret;
+ }
+ //__________________________________________________________________
+ /**
+ * Load a library
+ *
+ * @param what What library to load
+ * @param mode Mode (local, proof, grid)
+ * @param par If true, load as PAR
+ * @param rec If true, also load on slaves
+ *
+ * @return true on success
+ */
+ Bool_t LoadLibrary(const char* what, EMode mode, Bool_t par, Bool_t rec=false)
+ {
+ if (!what || what[0] == '\0') return true;
+
+ TString module(what);
+ TString libName(what);
+ if (!libName.BeginsWith("lib")) libName = Form("lib%s", libName.Data());
+ if (!libName.EndsWith(".so")) libName.Append(".so");
+
+ Int_t ret = 0;
+
+ switch (mode) {
+ case kLocal: // Just load and exit
+ gSystem->Load(libName.Data());
+ break;
+ case kGrid:
+ if (par) {
+ ret = SetupPAR(what) ? 0 : -1;
+ if (rec) fListOfPARs.Add(new TObjString(what));
+ } else {
+ ret = gSystem->Load(libName.Data());
+ if (rec) fListOfLibraries.Add(new TObjString(libName));
+ }
+ break;
+ case kProof:
+ ret = gProof->UploadPackage(what);
+ if (ret < 0) {
+ ret = gProof->UploadPackage(gSystem->ExpandPathName(Form("../%s.par",
+ what)));
+ if (ret < 0) {
+ ret =
+ gProof->UploadPackage(gSystem
+ ->ExpandPathName(Form("$ALICE_ROOT/%s.par",
+ what)));
+ if (ret < 0) {
+ Error("LoadLibrary",
+ "Could not find module %s.par in current directory nor "
+ "in $ALICE_ROOT", module.Data());
+ return false;
+ }
+ }
+ }
+ ret = gProof->EnablePackage(what);
+ break;
+ }
+ if (ret < 0) {
+ Error("LoadLibrary", "Couldn't load %s", what);
+ return false;
+ }
+ return true;
+ }
+
+ //__________________________________________________________________
+ Bool_t SetupPAR(const char* what)
+ {
+ if (!what || what[0] == '\0') return -1;
+
+ TString parFile(Form("%s.par", what));
+ if (gSystem->AccessPathName(parFile.Data())) {
+ if (gSystem->AccessPathName(Form("../%s.par", what))) {
+ // If not found
+ TString aliParFile =
+ gSystem->ExpandPathName(Form("$(ALICE_ROOT)/%s.par", what));
+ if (gSystem->AccessPathName(aliParFile.Data())) {
+ Error("SetupPAR", "PAR file %s not found in current directory or "
+ "$(ALICE_ROOT)", what);
+ return false;
+ }
+ // Copy to current directory
+ TFile::Cp(aliParFile, parFile);
+ }
+ else
+ gSystem->Exec(Form("ln -s ../%s.par .", what));
+ }
+
+ // Extract archive
+ gSystem->Exec(Form("tar xvzf %s", parFile.Data()));
+
+ // Change directory into par archive
+ TString cwd = gSystem->WorkingDirectory();
+
+ if (!gSystem->ChangeDirectory(what)) {
+ Error("SetupPAR", "Failed to change directory to %s", what);
+ return false;
+ }
+
+ // Test the build
+ if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
+ Info("SetupPar", "Building in PAR archive %s", what);
+ if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
+ Error("SetupPar", "Failed to build in PAR directory %s", what);
+ gSystem->ChangeDirectory(cwd.Data());
+ return false;
+ }
+ }
+
+ // Check for setup script
+ if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
+ Info("SetupPAR", "Setting up for PAR %s", what);
+ gROOT->Macro("PROOF-INF/SETUP.C");
+ }
+ if (!gSystem->ChangeDirectory(cwd.Data())) return false;
+
+ return true;
+ }
+ //__________________________________________________________________
+ TString SetupExtras()
+ {
+ TString ret;
+ TIter next(&fListOfExtras);
+ TObjString* obj = 0;
+ while ((obj = static_cast<TObjString*>(next()))) {
+ TString path = gSystem->ExpandPathName(obj->GetName());
+ if (!path.BeginsWith("/"))
+ // If not an absolute path, prepend to up-one
+ path = Form("../%s", path.Data());
+ if (gSystem->AccessPathName(path.Data())) {
+ // File not accessible
+ Warning("SetupExtras", "File %s not accessible", path.Data());
+ continue;
+ }
+ ret.Append(Form("%s ", gSystem->BaseName(obj->GetName())));
+ gSystem->Exec(Form("ln -s %s .", path.Data()));
+ }
+ ret = ret.Strip();
+ return ret;
+ }
+ //__________________________________________________________________
+ TString SetupSources()
+ {
+ TString nam = EscapedName();
+ TString ret;
+ TIter next(&fListOfSources);
+ TObject* src;
+ while ((src = next())) {
+ TString path = gSystem->ExpandPathName(src->GetName());
+ if (!path.BeginsWith("/"))
+ // If not an absolute path, prepend to up-one
+ path = Form("../%s", path.Data());
+ if (gSystem->AccessPathName(path.Data())) {
+ // File not accessible
+ Warning("SetupSources", "File %s not accessible", path.Data());
+ continue;
+ }
+ ret.Append(Form("%s ", gSystem->BaseName(src->GetName())));
+ gSystem->Exec(Form("ln -s %s .", path.Data()));
+ }
+ ret = ret.Strip();
+ return ret;
+ }
+ //__________________________________________________________________
+ TString SetupLibraries()
+ {
+ TString ret;
+ TIter next(&fListOfLibraries);
+ TObject* lib;
+ while ((lib = next())) {
+ ret.Append(lib->GetName());
+ ret.Append(" ");
+ }
+ // Also add extra files to this variable
+ ret.Append(SetupExtras());
+ ret = ret.Strip();
+ return ret;
+ }
+ //__________________________________________________________________
+ /**
+ * Scan directory @a dir (possibly recursive) for tree files to add
+ * to the chain. This does not follow sym-links
+ *
+ * @param dir Directory to scan
+ * @param chain Chain to add to
+ * @param type Type of tree (ESD or AOD)
+ * @param recursive Whether to scan recursively
+ * @param mc Look also for MC files if true
+ *
+ * @return true if any files where added
+ */
+ Bool_t ScanDirectory(TSystemDirectory* dir, TChain* chain,
+ EType type, bool recursive, bool mc)
+ {
+ TString fnPattern;
+ switch (type) {
+ case kESD: fnPattern = "AliESD"; break;
+ case kAOD: fnPattern = "AliAOD"; break;
+ }
+
+ // Assume failure
+ Bool_t ret = false;
+
+ // Get list of files, and go back to old working directory
+ TString oldDir(gSystem->WorkingDirectory());
+ TList* files = dir->GetListOfFiles();
+ if (!gSystem->ChangeDirectory(oldDir)) {
+ Error("ScanDirectory", "Failed to go back to %s", oldDir.Data());
+ return false;
+ }
+ if (!files) return false;
+
+ TList toAdd;
+ toAdd.SetOwner();
+ Bool_t hasGAlice = (!mc ? true : false);
+ Bool_t hasKine = (!mc ? true : false);
+ Bool_t hasTrRef = (!mc ? true : false);
+
+ // Sort list of files and check if we should add it
+ files->Sort();
+ TIter next(files);
+ TSystemFile* file = 0;
+ while ((file = static_cast<TSystemFile*>(next()))) {
+ TString name(file->GetName());
+ TString title(file->GetTitle());
+ TString full(gSystem->ConcatFileName(file->GetTitle(), name.Data()));
+ if (dynamic_cast<TSystemDirectory*>(file)) full = title;
+ // Ignore special links
+ if (name == "." || name == "..") {
+ // Info("ScanDirectory", "Ignoring %s", name.Data());
+ continue;
+ }
+
+ FileStat_t fs;
+ if (gSystem->GetPathInfo(full.Data(), fs)) {
+ Warning("ScanDirectory", "Cannot stat %s (%s)", full.Data(),
+ gSystem->WorkingDirectory());
+ continue;
+ }
+ // Check if this is a directory
+ if (file->IsDirectory(full)) {
+ if (recursive) {
+ // if (title[0] == '/')
+ TSystemDirectory* d = new TSystemDirectory(file->GetName(),
+ full.Data());
+ if (ScanDirectory(d,chain,type,recursive,mc))
+ ret = true;
+ delete d;
+ }
+ continue;
+ }
+
+ // If this is not a root file, ignore
+ if (!name.EndsWith(".root")) continue;
+
+ // If this file does not contain AliESDs, ignore
+ if (!name.Contains(fnPattern)) {
+ // Info("ScanDirectory", "%s does not match pattern %s",
+ // name.Data(), fnPattern.Data());
+ if (mc) {
+ if (name.CompareTo("galice.root") == 0) hasGAlice = true;
+ if (name.CompareTo("Kinematics.root") == 0) hasKine = true;
+ if (name.CompareTo("TrackRefs.root") == 0) hasTrRef = true;
+ }
+ continue;
+ }
+
+ // Add
+ // Info("ScanDirectory", "Adding %s", full.Data());
+ toAdd.Add(new TObjString(full));
+ }
+
+ if (mc && toAdd.GetEntries() > 0 &&
+ (!hasGAlice || !hasKine || !hasTrRef)) {
+ Warning("ScanDirectory",
+ "one or more of {galice,Kinematics,TrackRefs}.root missing from "
+ "%s, not adding anything from this directory",
+ dir->GetTitle());
+ toAdd.Delete();
+ }
+
+ TIter nextAdd(&toAdd);
+ TObjString* s = 0;
+ while ((s = static_cast<TObjString*>(nextAdd()))) {
+ // Info("ScanDirectory", "Adding %s", s->GetString().Data());
+ chain->Add(s->GetString());
+ }
+ if (toAdd.GetEntries() > 0) ret = true;
+
+ gSystem->ChangeDirectory(oldDir);
+ return ret;
+ }
+ //__________________________________________________________________
+ /**
+ * Create a chain from an XML containing an collection
+ *
+ * @param treeName Name of tree's
+ * @param xmlFile XML collection
+ *
+ * @return Newly allocated chain or null
+ */
+ TChain* CreateChainFromXML(const char* treeName,
+ const char* xmlFile)
+ {
+ TGridCollection* collection = TAlienCollection::Open(xmlFile);
+ if (!collection) {
+ Error("CreateChainFromXML", "Cannot create AliEn collection from "
+ "XML file %s", xmlFile);
+ return 0;
+ }
+
+ TChain* chain = new TChain(treeName);
+ collection->Reset();
+ while (collection->Next()) chain->Add(collection->GetTURL(""));
+
+ return chain;
+ }
+ //__________________________________________________________________
+ /**
+ * Create a chain of data
+ *
+ * @param type Type of data
+ * @param mode Operation mode
+ * @param mc Assume MC input if true
+ *
+ * @return TChain of data
+ */
+ TChain* CreateChain(EType type, EMode mode, EOper /* oper */, Bool_t mc)
+ {
+ TString treeName;
+ switch (type) {
+ case kESD: treeName = "esdTree"; break;
+ case kAOD: treeName = "aodTree"; break;
+ }
+
+ TChain* chain = 0;
+ switch (mode) {
+ case kProof:
+ if (!fDataSet.IsNull()) break;
+ // Otherwise fall through
+ case kLocal:
+ if (fXML.IsNull()) {
+ chain = new TChain(treeName.Data());
+ TString dir(fDataDir);
+ if (dir == ".") dir = "";
+ if (!dir.BeginsWith("/")) dir = Form("../%s", dir.Data());
+ TString savdir(gSystem->WorkingDirectory());
+ TSystemDirectory d(gSystem->BaseName(dir.Data()), dir.Data());
+ if (!ScanDirectory(&d, chain, type, true, mc)) {
+ delete chain;
+ chain = 0;
+ }
+ gSystem->ChangeDirectory(savdir);
+ }
+ else
+ chain = CreateChainFromXML(treeName.Data(), fXML.Data());
+ break;
+ case kGrid: break; // Do nothing - we use plugin
+ }
+
+ if (chain && chain->GetNtrees() <= 0) {
+ delete chain;
+ return 0;
+ }
+ return chain;
+ }
+ //__________________________________________________________________
+ TString fName; // Name of analysis
+ TString fRootVersion; // ROOT version to use
+ TString fAliRootVersion; // AliROOT version to use
+ TString fAliEnAPIVersion; // AliEn API version to use
+ TString fProofServer; // Name of proof server
+ TString fDataDir; // Grid Input directory
+ TString fDataSet; // Proof data set name
+ TString fXML; // XML collection for local/proof mode
+ TArrayI fRunNumbers; // List of run number
+ TList fListOfPARs; // List of PAR files to use
+ TList fListOfSources; // List of sources to upload and AcLIC
+ TList fListOfLibraries; // List of libraries to load
+ TList fListOfExtras; // List of extra files to upload
+ Int_t fNReplica; // Storage replication
+ Int_t fESDPass;
+ TString fPassPostfix; // Possible pass postfix
+ TString fEscapedName;
+ Bool_t fAllowOverwrite;
+};
+
+
+void
+BuildTrainSetup()
+{
+ gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
+ "$ALICE_ROOT/ANALYSIS/macros",
+ gROOT->GetMacroPath()));
+ gSystem->AddIncludePath("-I${ALICE_ROOT}/include");
+ gSystem->Load("libRAliEn");
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+ TString path = gSystem->Which(gROOT->GetMacroPath(), "TrainSetup.C");
+ Info("BuildTrainSetup", "Path=%s", path.Data());
+ TString tmp("TrainSetup");
+ FILE* fp = gSystem->TempFileName(tmp, ".");
+ fclose(fp);
+ gSystem->Unlink(tmp);
+ tmp.Append(".C");
+ Info("BuildTrainSetup", "Copy %s -> %s", path.Data(), tmp.Data());
+ gSystem->CopyFile(path, tmp);
+ gROOT->LoadMacro(Form("%s+g", tmp.Data()));
+ gSystem->Unlink(tmp);
+ tmp.ReplaceAll(".C", "_C.so");
+ gSystem->Unlink(tmp);
+ tmp.ReplaceAll("_C.so", "_C.d");
+ gSystem->Unlink(tmp);
+}
+
+
+//____________________________________________________________________
+//
+// EOF
+//
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff