gSystem->Load("libANALYSISalice.so");
gSystem->Load("libPWG2forward2.so");
TString macroPath(gROOT->GetMacroPath());
- macroPath.Append(":${ALICE_ROOT}/FMD/scripts");
+ macroPath.Append(":${ALICE_ROOT}/PWG2/FORWARD/analysis2");
+ macroPath.Append(":${ALICE_ROOT}/PWG2/FORWARD/analysis2/scripts");
gROOT->SetMacroPath(macroPath.Data());
gSystem->SetIncludePath("-I`root-config --incdir` "
"-I${ALICE_ROOT} "
+/**
+ * Script to draw the energy loss fits
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
#include <TFile.h>
#include <THStack.h>
#include <TList.h>
return 0;
}
- TList* forward =
- static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
+ TList* forward = static_cast<TList*>(file->Get("Forward"));
+ // static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
if (!forward) {
Error("DrawFits", "Couldn't get forward list from %s", fname);
return 0;
--- /dev/null
+/**
+ * Script to draw the energy loss fits
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+#ifndef __CINT__
+#include <TFile.h>
+#include <TList.h>
+#include <TError.h>
+#include "AliFMDCorrELossFit.h"
+#include "AliForwardCorrectionManager.h"
+#endif
+
+//____________________________________________________________________
+void
+ExtractELoss(const char* fname="energyFits.root",
+ UShort_t sys=1, UShort_t sNN=900, Short_t field=5, Bool_t mc=false)
+{
+#ifdef __CINT__
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+#endif
+
+ TFile* file = TFile::Open(fname, "READ");
+ if (!file) {
+ Error("ExtractELoss", "Couldn't open %s", fname);
+ return;
+ }
+
+ TList* forward = static_cast<TList*>(file->Get("Forward"));
+ // static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
+ if (!forward) {
+ Error("ExtractELoss", "Couldn't get forward list from %s", fname);
+ return;
+ }
+
+ TList* fitter = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
+ if (!fitter) {
+ Error("ExtractELoss", "Couldn't get fitter folder");
+ return;
+ }
+
+ TString cName(AliFMDCorrELossFit::Class()->GetName());
+
+ AliFMDCorrELossFit* obj =
+ static_cast<AliFMDCorrELossFit*>(fitter->FindObject(cName));
+ if (!obj) {
+ Error("ExtractELoss", "Couldn't get %s correction object", cName.Data());
+ return;
+ }
+
+ AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+
+ TString ofName(mgr.GetFileName(AliForwardCorrectionManager::kELossFits,
+ sys, sNN, field, mc));
+ TFile* output = TFile::Open(ofName.Data(), "RECREATE");
+ if (!output) {
+ Error("ExtractELoss", "Failed to open file %s", ofName.Data());
+ return;
+ }
+
+ TString oName(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
+ obj->Write(oName);
+
+ output->Write();
+ output->ls();
+ output->Close();
+
+ TString dName(mgr.GetFileDir(AliForwardCorrectionManager::kELossFits));
+ Info("ExtractELoss", "Wrote %s object %s to %s",cName.Data(),oName.Data(),
+ ofName.Data());
+ Info("ExtractELoss", "%s should be copied to %s",ofName.Data(),dName.Data());
+ Info("ExtractELoss", "Do for example\n\t"
+ "aliroot $ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MoveCorrections.C\(0,0,0,0,1\)\n\t"
+ "cp %s %s/", ofName.Data(), gSystem->ExpandPathName(dName.Data()));
+}
+
+
+
+//____________________________________________________________________
+void
+ExtractELoss(const char* fname="energyFits.root",
+ const char* sys="p-p",
+ Float_t sNN=900,
+ Float_t field=5)
+{
+ UShort_t uSys = AliForwardUtil::ParseCollisionSystem(sys);
+ UShort_t usNN = AliForwardUtil::ParseCenterOfMassEnergy(uSys,sNN);
+ Short_t sField = AliForwardUtil::ParseMagneticField(field);
+
+ ExtractELoss(fname, uSys, usNN, sField);
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+/**
+ * Test script to fit the energy loss spectra
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
#ifndef __CINT__
# include "AliForwardUtil.h"
# include <TFile.h>
+/**
+ * Load the libraries of PWG2/FORWARD/analsysis2
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
void
LoadLibs()
{
gSystem->Load("libANALYSIS");
gSystem->Load("libANALYSISalice");
gSystem->Load("libPWG0base");
- gSystem->Load("libPWG2forward");
gSystem->Load("libPWG2forward2");
}
-
+//
+// EOF
+//
--- /dev/null
+/**
+ * Set-up for a PROOF analysis job. Make TProof object and load pars.
+ *
+ */
+void
+LoadPars()
+{
+ TProof::Open("workers=2");
+ const char* pkgs[] = { "STEERBase", "ESD", "AOD", "ANALYSIS",
+ "ANALYSISalice", "PWG2forward2", 0};
+ const char** pkg = pkgs;
+ while (*pkg) {
+ gProof->UploadPackage(Form("${ALICE_ROOT}/%s.par",*pkg));
+ gProof->EnablePackage(*pkg);
+ pkg++;
+ }
+}
+//
+// EOF
+//
--- /dev/null
+void
+MakeCorrRepository()
+{
+ AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+
+ UInt_t what[] = {
+ AliForwardCorrectionManager::kSecondaryMap,
+ AliForwardCorrectionManager::kELossFits,
+ AliForwardCorrectionManager::kVertexBias,
+ AliForwardCorrectionManager::kMergingEfficiency,
+ AliForwardCorrectionManager::kDoubleHit,
+ 0
+ };
+ UInt_t* ptr = what;
+ while (*ptr) {
+ TString dir(gSystem->ExpandPathName(mgr.GetFileDir(*ptr)));
+ if (dir.IsNull()) {
+ ptr++;
+ continue;
+ }
+
+ Info("MakeCorrRepository", "Making directory %s", dir.Data());
+ gSystem->MakeDirectory(dir.Data());
+ ptr++;
+ }
+}
+
--- /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.
+ *
+ * @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="AnalysisResults.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
+TChain*
+MakeESDChain(const char* esddir)
+{
+ // --- Our data chain ----------------------------------------------
+ TChain* chain = new TChain("esdTree");
+
+ // --- Get list of ESDs --------------------------------------------
+ // Open source directory, and make sure we go back to were we were
+ TString oldDir(gSystem->WorkingDirectory());
+ TSystemDirectory d(esddir, esddir);
+ TList* files = d.GetListOfFiles();
+ gSystem->ChangeDirectory(oldDir);
+
+ // 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()))) {
+ if (file->IsDirectory()) continue;
+ TString name(file->GetName());
+ if (!name.EndsWith(".root")) continue;
+ if (!name.Contains("AliESDs")) continue;
+ TString esd(Form("%s/%s", file->GetTitle(), name.Data()));
+ Info("RunManager", "Adding %s to chain", esd.Data());
+ chain->Add(esd);
+ }
+ return chain;
+}
--- /dev/null
+Bool_t
+Backup(const TString& fname, Int_t n=10)
+{
+ TString fn(fname); fn.Append(Form(".%d", n));
+ if (!gSystem->AccessPathName(fn.Data())) {
+ Error("Backup", "Last possible backup slot (%d) filled for %s",
+ n, fname.Data());
+ return false;
+ }
+
+ for (Int_t i=n-1; i >= 0; i--) {
+ TString fi(fname);
+ if (i > 0) fi.Append(Form(".%d", i));
+ if (gSystem->AccessPathName(fi.Data())) continue;
+
+ TString fb(fname); fb.Append(Form(".%d", i+1));
+
+ if (gSystem->Rename(fi.Data(), fb.Data())) {
+ Error("Backup", "Failed to backup %s to %s", fi.Data(), fb.Data());
+ return false;
+ }
+ }
+ return true;
+}
+
+
+Bool_t
+MoveWhat(UInt_t what)
+{
+ TString nWhat;
+ switch (what) {
+ case AliForwardCorrectionManager::kSecondaryMap:
+ nWhat = "secondary map"; break;
+ case AliForwardCorrectionManager::kDoubleHit:
+ nWhat = "double hit"; break;
+ case AliForwardCorrectionManager::kVertexBias:
+ nWhat = "vertex bias"; break;
+ case AliForwardCorrectionManager::kMergingEfficiency:
+ nWhat = "merging efficiency";break;
+ case AliForwardCorrectionManager::kELossFits:
+ nWhat = "energy loss fits"; break;
+ }
+ Info("MakeWhat", " Copying %s corrections", nWhat.Data());
+
+ AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+ TString dir = gSystem->ExpandPathName(mgr.GetFileDir(what));
+
+ // Make the directory if it doesn't exist
+ if (gSystem->AccessPathName(dir.Data())) {
+ Info("MakeWhat", " Making directory %s ... ", dir.Data());
+ if (gSystem->mkdir(dir.Data(), true)) {
+ Error("MakeWhat", "couldn't make directory %s", dir.Data());
+ return false;
+ }
+ }
+
+ TString pattern = mgr.GetFileName(what, 1, 900, 5, false);
+ pattern.ReplaceAll("0900GeV", "[0-9]+GeV");
+ pattern.ReplaceAll("pp", "[a-zA-Z]+");
+ pattern.ReplaceAll("p5kG", "[pm][0-9]+kG");
+ pattern.ReplaceAll("real", "[a-zA-Z]+");
+
+ TPRegexp regex(pattern);
+
+ TSystemDirectory sysdir(".", ".");
+ TIter next(sysdir.GetListOfFiles());
+ TSystemFile* file = 0;
+ while ((file = static_cast<TSystemFile*>(next()))) {
+ if (file->IsDirectory()) continue;
+ TString fname(file->GetName());
+
+ if (!regex.Match(fname)) continue;
+
+ // Info("MakeWhat", " match: %s", fname.Data());
+ TString to(gSystem->ConcatFileName(dir.Data(), fname.Data()));
+
+ if (!Backup(to)) return false;
+
+ Info("MakeWhat", " copying %s\n to %s",
+ fname.Data(), to.Data());
+ if (gSystem->CopyFile(fname.Data(), to.Data(), false)) {
+ Error("MakeWhat", "Failed to copy %s to %s", fname.Data(), to.Data());
+ return false;
+ }
+ }
+ return true;
+}
+
+void
+MoveCorrections(bool sec=true,
+ bool dbl=true,
+ bool vtx=true,
+ bool merge=true,
+ bool eloss=true)
+{
+ Info("MoveCorrections", "Loadingl libraries ...");
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+
+ Info("MoveCorrections", "Moving selected corrections ...");
+ if (sec) MoveWhat(AliForwardCorrectionManager::kSecondaryMap);
+ if (dbl) MoveWhat(AliForwardCorrectionManager::kDoubleHit);
+ if (vtx) MoveWhat(AliForwardCorrectionManager::kVertexBias);
+ if (merge) MoveWhat(AliForwardCorrectionManager::kMergingEfficiency);
+ if (eloss) MoveWhat(AliForwardCorrectionManager::kELossFits);
+
+}
+
--- /dev/null
+Color_t Color(UShort_t d, Char_t r ) const
+{
+ return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
+ + ((r == 'I' || r == 'i') ? 2 : -2));
+}
+
+/**
+ *
+ * @param sys Collision system
+ * @param cms Center of mass energy per nucleon in GeV
+ * @param field Magnetic field
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+void
+RunCopyELossFit(UShort_t sys, UShort_t cms, Short_t field, bool mc=false)
+{
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+ gSystem->Load("libPWG2forward.so");
+
+ AliFMDAnaParameters* p = AliFMDAnaParameters::Instance();
+ p->SetRealData(!mc);
+ p->SetEnergy(Float_t(cms));
+ p->SetMagField(Float_t(field));
+ p->SetCollisionSystem(sys == 1 ? "pp" : "PbPb");
+ p->Init(true, AliFMDAnaParameters::kBackgroundCorrection|
+ AliFMDAnaParameters::kEnergyDistributions);
+
+ Int_t nVtx = p->GetNvtxBins();
+ Double_t minVtx = -p->GetVtxCutZ();
+ Double_t maxVtx = -p->GetVtxCutZ();
+ Int_t nEta = p->GetNetaBins();
+ Double_t minEta = p->GetEtaMin();
+ Double_t maxEta = p->GetEtaMax();
+
+ TAxis vtxAxis(nVtx, minVtx, maxVtx);
+ TAxis etaAxis(nEta, minEta, maxEta);
+ AliFMDCorrELossFit* m = new AliFMDCorrELossFit;
+ m->SetEtaAxis(etaAxis);
+
+ for (UShort_t d = 1; d <= 3; d++) {
+ UShort_t nQ = (d == 1 ? 1 : 2);
+ for (UShort_t q = 0; q < nQ; q++) {
+ Char_t r = (q == 0 ? 'I' : 'O');
+
+ for (UShort_t b = 1; b < nEta; b++) {
+ Double_t eta = etaAxis.GetBinCenter(b);
+ Info("RunCopyELossFit", "FMD%d%c, bin %3d, eta %+8.4f", d, r, b, eta);
+ if (eta < minEta || eta > maxEta) continue;
+ TH1F* oldhist = p->GetEnergyDistribution(d,r,eta);
+ if (!oldhist) {
+ Warning("RunCopyELossFit",
+ "Didn't find energy distribution for FMD%d%c, eta bin %3d",
+ d, r, b);
+ continue;
+ }
+ TF1* oldfunc = oldhist->GetFunction("FMDfitFunc");
+ if (!oldfunc) {
+ // Warning("RunCopyELossFit",
+ // "Didn't find energy fit for FMD%d%c, eta bin %3d",
+ // d, r, b);
+ continue;
+ }
+ Double_t chi2 = oldfunc->GetChisquare();
+ UShort_t nu = oldfunc->GetNDF();
+ Double_t c = oldfunc->GetParameter(0);
+ Double_t delta = oldfunc->GetParameter(1);
+ Double_t xi = oldfunc->GetParameter(2);
+ Double_t a2 = 0, a3 = 0, ea2 = 0, ea3 = 0;
+ if (oldfunc->GetNpar() >= 4) {
+ a2 = oldfunc->GetParameter(3);
+ ea2 = oldfunc->GetParError(3);
+ }
+ if (oldfunc->GetNpar() >= 5) {
+ a3 = oldfunc->GetParameter(4);
+ ea3 = oldfunc->GetParError(4);
+ }
+ Int_t n = (a3 < 1e-5 ? (a2 < 1e-5 ? 1 : 2) : 3);
+ Double_t ec = oldfunc->GetParError(0);
+ Double_t edelta = oldfunc->GetParError(1);
+ Double_t exi = oldfunc->GetParError(2);
+ delta = delta - xi * -0.22278298;
+ edelta = TMath::Sqrt(TMath::Power(edelta,2) -
+ TMath::Power(0.22278298*exi,2));
+ Info("RunCopyELossFit", "FMD%d%c etabin=%3d: n=%d\n"
+ " C=%f+/-%f, Delta=%f+/-%f, xi=%f+/-%f, a2=%f+/-%f, a3=%f+/-%f",
+ d, r, b, n, c, ec, delta, edelta, xi, exi, a2, ea2, a3, ea3);
+ Double_t a[] = { a2, a3, -9999 };
+ Double_t ea[] = { ea2, ea3, -9999 };
+
+ AliFMDCorrELossFit::ELossFit* fit = new
+ AliFMDCorrELossFit::ELossFit(0, n,
+ chi2, nu,
+ c, ec,
+ delta, edelta,
+ xi, exi,
+ 0, 0,
+ 0, 0,
+ a, ea);
+ fit->CalculateQuality();
+ fit->fBin = b;
+ fit->fDet = d;
+ fit->fRing = r;
+ fit->Print();
+
+ Info("RunCopyELossFit", "Setting fit FMD%d%c etabin=%3d: %p",
+ d, r, b, fit);
+ m->SetFit(d, r, b, fit);
+ }
+ }
+ }
+ AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+ TString fname(mgr.GetFileName(AliForwardCorrectionManager::kELossFits,
+ sys, cms, field, mc));
+ TFile* output = TFile::Open(fname.Data(), "RECREATE");
+ if (!output) {
+ Warning("RunCopyELossFit", "Failed to open output file %s", fname.Data());
+ return kFALSE;
+ }
+ m->Write(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
+ output->Write();
+ output->Close();
+ Info("RunCopyELossFit",
+ "File %s created. It should be copied to %s and stored in SVN",
+ fname.Data(), mgr.GetFileDir(AliForwardCorrectionManager::kELossFits));
+
+}
+//
+// EOF
+//
--- /dev/null
+Color_t Color(UShort_t d, Char_t r ) const
+{
+ return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
+ + ((r == 'I' || r == 'i') ? 2 : -2));
+}
+
+/**
+ *
+ * @param sys Collision system
+ * @param cms Center of mass energy per nucleon in GeV
+ * @param field Magnetic field
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+void
+RunCopyMergeEff(UShort_t sys, UShort_t cms, Short_t field)
+{
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+ gSystem->Load("libPWG2forward.so");
+
+ AliFMDAnaParameters* p = AliFMDAnaParameters::Instance();
+ p->SetEnergy(Float_t(cms));
+ p->SetMagField(Float_t(field));
+ p->SetCollisionSystem(sys == 1 ? "pp" : "PbPb");
+ p->Init(true, AliFMDAnaParameters::kBackgroundCorrection|
+ AliFMDAnaParameters::kSharingEfficiency);
+
+ Int_t nVtx = p->GetNvtxBins();
+ Double_t minVtx = -p->GetVtxCutZ();
+ Double_t maxVtx = -p->GetVtxCutZ();
+ Int_t nEta = p->GetNetaBins();
+ Double_t minEta = p->GetEtaMin();
+ Double_t maxEta = p->GetEtaMax();
+
+ TAxis vtxAxis(nVtx, minVtx, maxVtx);
+ AliFMDCorrMergingEfficiency* m = new AliFMDCorrMergingEfficiency;
+ m->SetVertexAxis(nVtx, minVtx, maxVtx);
+
+ for (UShort_t d = 1; d <= 3; d++) {
+ UShort_t nQ = (d == 1 ? 1 : 2);
+ for (UShort_t q = 0; q < nQ; q++) {
+ Char_t r = (q == 0 ? 'I' : 'O');
+
+ for (UShort_t b = 1; b <= nVtx; b++) {
+ TH1F* oldmap = p->GetSharingEfficiency(d, r, b-1);
+ if (!oldmap) {
+ Warning("RunCopyMergeEff",
+ "Didn't find secondary map correction "
+ "for FMD%d%c, vertex bin %3d", d, r, b);
+ continue;
+ }
+
+ TH1D* newmap = new TH1D(Form("FMD%d%c_vtxbin%03d", d, r, b),
+ Form("Merging efficiency for FMD%d%c "
+ "in vertex bin %d [%+8.4f,%+8.4f]",
+ d, r, b, vtxAxis.GetBinLowEdge(b),
+ vtxAxis.GetBinUpEdge(b)),
+ nEta, minEta, maxEta);
+ newmap->SetXTitle("#eta");
+ newmap->SetYTitle("dN_{ch,i,incl}/d#eta / #sum_{i} N_{ch,i,FMD}");
+ newmap->SetDirectory(0);
+ newmap->SetStats(0);
+ newmap->Sumw2();
+ newmap->Add(oldmap);
+
+ Info("RunCopyMergeEff",
+ "Copying %s to %s", oldmap->GetName(), newmap->GetName());
+
+ m->SetCorrection(d, r, b, newmap);
+ }
+ }
+ }
+ AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+ TString fname(mgr.GetFileName(AliForwardCorrectionManager::kMergingEfficiency,
+ sys, cms, field, false));
+ TFile* output = TFile::Open(fname.Data(), "RECREATE");
+ if (!output) {
+ Warning("Run", "Failed to open output file %s", fname.Data());
+ return kFALSE;
+ }
+ m->Write(mgr.GetObjectName(AliForwardCorrectionManager::kMergingEfficiency));
+ output->Write();
+ output->Close();
+ Info("Run", "File %s created. It should be copied to %s and stored in SVN",
+ fname.Data(),
+ mgr.GetFileDir(AliForwardCorrectionManager::kMergingEfficiency));
+
+}
+//
+// EOF
+//
--- /dev/null
+Color_t Color(UShort_t d, Char_t r ) const
+{
+ return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
+ + ((r == 'I' || r == 'i') ? 2 : -2));
+}
+
+/**
+ *
+ * @param sys Collision system
+ * @param cms Center of mass energy per nucleon in GeV
+ * @param field Magnetic field
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+void
+RunCopySecMap(UShort_t sys, UShort_t cms, Short_t field)
+{
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+ gSystem->Load("libPWG2forward.so");
+
+ AliFMDAnaParameters* p = AliFMDAnaParameters::Instance();
+ p->SetEnergy(Float_t(cms));
+ p->SetMagField(Float_t(field));
+ p->SetCollisionSystem(sys == 1 ? "pp" : "PbPb");
+ p->Init(true, AliFMDAnaParameters::kBackgroundCorrection);
+
+ Int_t nVtx = p->GetNvtxBins();
+ Double_t minVtx = -p->GetVtxCutZ();
+ Double_t maxVtx = p->GetVtxCutZ();
+ Int_t nEta = p->GetNetaBins();
+ Double_t minEta = p->GetEtaMin();
+ Double_t maxEta = p->GetEtaMax();
+
+ TAxis vtxAxis(nVtx, minVtx, maxVtx);
+ AliFMDCorrSecondaryMap* m = new AliFMDCorrSecondaryMap;
+ m->SetVertexAxis(nVtx, minVtx, maxVtx);
+ m->SetEtaAxis(nEta,minEta,maxEta);
+ AliFMDCorrDoubleHit* dh = new AliFMDCorrDoubleHit;
+
+ for (UShort_t d = 1; d <= 3; d++) {
+ UShort_t nQ = (d == 1 ? 1 : 2);
+ for (UShort_t q = 0; q < nQ; q++) {
+ Char_t r = (q == 0 ? 'I' : 'O');
+
+ TH1F* old2hit = p->GetDoubleHitCorrection(d, r);
+ if (!old2hit)
+ Warning("RunCopySecMap",
+ "Didn't find double hit correction for FMD%d%c", d, r);
+ else {
+ TH1D* new2hit = new TH1D(Form("FMD%d%c", d, r),
+ Form("Double hit correction for FMD%d%c",d,r),
+ nEta, minEta, maxEta);
+ new2hit->SetXTitle("#eta");
+ new2hit->SetYTitle("#sum_{i} N_{i,strips hit}(#eta)/"
+ "#sum_{i} N_{i,total hits}(#eta)");
+ new2hit->SetFillColor(Color(d,r));
+ new2hit->SetFillStyle(3001);
+ new2hit->SetDirectory(0);
+ new2hit->SetStats(0);
+ new2hit->Sumw2();
+ new2hit->Add(old2hit);
+
+ Info("RunCopySecMap",
+ "Copying %s to %s", old2hit->GetName(), new2hit->GetName());
+
+ dh->SetCorrection(d, r, new2hit);
+ }
+
+ for (UShort_t b = 1; b <= nVtx; b++) {
+ TH2F* oldmap = p->GetBackgroundCorrection(d, r, b-1);
+ if (!oldmap) {
+ Warning("RunCopySecMap",
+ "Didn't find secondary map correction "
+ "for FMD%d%c, vertex bin %3d", d, r, b);
+ continue;
+ }
+
+ TH2D* newmap = new TH2D(Form("FMD%d%c_vtxbin%03d", d, r, b),
+ Form("Secondary map correction for FMD%d%c "
+ "in vertex bin %d [%+8.4f,%+8.4f]",
+ d, r, b, vtxAxis.GetBinLowEdge(b),
+ vtxAxis.GetBinUpEdge(b)),
+ nEta, minEta, maxEta,
+ oldmap->GetYaxis()->GetNbins(),
+ oldmap->GetYaxis()->GetXmin(),
+ oldmap->GetYaxis()->GetXmax());
+ newmap->SetXTitle("#eta");
+ newmap->SetYTitle("#phi [radians]");
+ newmap->SetZTitle("#sum_{i} N_{ch,i,primary} / #sum_{i} N_{ch,i,FMD}");
+ newmap->SetDirectory(0);
+ newmap->SetStats(0);
+ newmap->Sumw2();
+ newmap->Add(oldmap);
+
+ Info("RunCopySecMap",
+ "Copying %s to %s", oldmap->GetName(), newmap->GetName());
+
+ m->SetCorrection(d, r, b, newmap);
+ }
+ }
+ }
+ AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+ TString fname(mgr.GetFileName(AliForwardCorrectionManager::kSecondaryMap,
+ sys, cms, field, false));
+ TFile* output = TFile::Open(fname.Data(), "RECREATE");
+ if (!output) {
+ Warning("Run", "Failed to open output file %s", fname.Data());
+ return kFALSE;
+ }
+ m->Write(mgr.GetObjectName(AliForwardCorrectionManager::kSecondaryMap));
+ output->Write();
+ output->Close();
+ Info("Run", "File %s created. It should be copied to %s and stored in SVN",
+ fname.Data(),mgr.GetFileDir(AliForwardCorrectionManager::kSecondaryMap));
+
+ fname = mgr.GetFileName(AliForwardCorrectionManager::kDoubleHit,
+ sys, cms, field, false);
+ output = TFile::Open(fname.Data(), "RECREATE");
+ if (!output) {
+ Warning("Run", "Failed to open output file %s", fname.Data());
+ return kFALSE;
+ }
+ dh->Write(mgr.GetObjectName(AliForwardCorrectionManager::kDoubleHit));
+ output->Write();
+ output->Close();
+ Info("Run", "File %s created. It should be copied to %s and stored in SVN",
+ fname.Data(),mgr.GetFileDir(AliForwardCorrectionManager::kDoubleHit));
+}
+//
+// EOF
+//
--- /dev/null
+Color_t Color(UShort_t d, Char_t r ) const
+{
+ return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
+ + ((r == 'I' || r == 'i') ? 2 : -2));
+}
+
+/**
+ *
+ * @param sys Collision system
+ * @param cms Center of mass energy per nucleon in GeV
+ * @param field Magnetic field
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+void
+RunCopyVtxBias(UShort_t sys, UShort_t cms, Short_t field)
+{
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+ gSystem->Load("libPWG2forward.so");
+
+ AliFMDAnaParameters* p = AliFMDAnaParameters::Instance();
+ p->SetEnergy(Float_t(cms));
+ p->SetMagField(Float_t(field));
+ p->SetCollisionSystem(sys == 1 ? "pp" : "PbPb");
+ p->Init(true, AliFMDAnaParameters::kBackgroundCorrection|
+ AliFMDAnaParameters::kEventSelectionEfficiency);
+
+ Int_t nVtx = p->GetNvtxBins();
+ Double_t minVtx = -p->GetVtxCutZ();
+ Double_t maxVtx = -p->GetVtxCutZ();
+ Int_t nEta = p->GetNetaBins();
+ Double_t minEta = p->GetEtaMin();
+ Double_t maxEta = p->GetEtaMax();
+
+ TAxis vtxAxis(nVtx, minVtx, maxVtx);
+ AliFMDCorrVertexBias* obj = new AliFMDCorrVertexBias;
+ obj->SetVertexAxis(vtxAxis);
+
+ for (UShort_t b = 1; b <= nVtx; b++) {
+ for (UShort_t q = 0; q < 2; q++) {
+ Char_t r = q == 0 ? 'I' : 'O';
+
+ TH2F* oldcorr = p->GetEventSelectionEfficiency("INEL",b-1,r);
+ if (!oldcorr) {
+ Warning("RunCopyVtxBias",
+ "Didn't find secondary map correction "
+ "for ring type %c, vertex bin %3d", r, b);
+ continue;
+ }
+
+ TH2D* newcorr = new TH2D(Form("FMDX%c", r),
+ Form("Vertex bias correction for %c rings "
+ "in vertex bin %d [%+8.4f,%+8.4f]",
+ r, b, vtxAxis.GetBinLowEdge(b),
+ vtxAxis.GetBinUpEdge(b)),
+ nEta, minEta, maxEta,
+ oldcorr->GetYaxis()->GetNbins(),
+ oldcorr->GetYaxis()->GetXmin(),
+ oldcorr->GetYaxis()->GetXmax());
+ newcorr->SetXTitle("#eta");
+ newcorr->SetYTitle("#phi [radians]");
+ newcorr->SetZTitle("1/N_{t}#sum_{i}^{N_{tv}} N_{ch,i,primary} / "
+ "1/N_{v}#sum_{i}^{N_{v}} N_{ch,i,primary}");
+ newcorr->SetDirectory(0);
+ newcorr->SetStats(0);
+ newcorr->Sumw2();
+ newcorr->Add(oldcorr);
+
+ Info("RunCopyVtxBias",
+ "Copying %s to %s", oldcorr->GetName(), newcorr->GetName());
+
+ obj->SetCorrection(r, b, newcorr);
+ }
+ }
+ AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+ TString fname(mgr.GetFileName(AliForwardCorrectionManager::kVertexBias,
+ sys, cms, field, false));
+ 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::kVertexBias));
+ output->Write();
+ output->Close();
+ Info("Run", "File %s created. It should be copied to %s and stored in SVN",
+ fname.Data(),mgr.GetFileDir(AliForwardCorrectionManager::kVertexBias));
+
+}
+//
+// EOF
+//
--- /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
+ */
+void
+RunMakeELossFit(UShort_t sys,
+ UShort_t cms,
+ Short_t field,
+ Bool_t mc=false,
+ const char* filename="AnalysisResults.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/MakeELossFit.C");
+
+ std::cout << "Making MakeELossFit object (sys=" << sys
+ << ", cms=" << cms << ", field=" << field << ", mc=" << mc
+ << ")" << std::endl;
+ MakeELossFit mef(sys, cms, field, mc, "AnalysisResults.root");
+
+ std::cout << "Runing maker ..." << std::endl;
+ mef.Run();
+}
+//
+// 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();
+}