From ab0f914c4be5ce381fa8e7eb383a0ae8363975a0 Mon Sep 17 00:00:00 2001 From: cholm Date: Thu, 16 Dec 2010 11:46:56 +0000 Subject: [PATCH] Added some more scripts --- PWG2/FORWARD/analysis2/scripts/Compile.C | 3 +- PWG2/FORWARD/analysis2/scripts/DrawFits.C | 9 +- PWG2/FORWARD/analysis2/scripts/ExtractELoss.C | 96 +++++++ PWG2/FORWARD/analysis2/scripts/FitELoss.C | 5 + PWG2/FORWARD/analysis2/scripts/LoadLibs.C | 10 +- PWG2/FORWARD/analysis2/scripts/LoadPars.C | 20 ++ .../analysis2/scripts/MakeCorrRepository.C | 27 ++ PWG2/FORWARD/analysis2/scripts/MakeELossFit.C | 228 ++++++++++++++++ PWG2/FORWARD/analysis2/scripts/MakeESDChain.C | 28 ++ .../analysis2/scripts/MoveCorrections.C | 107 ++++++++ .../analysis2/scripts/RunCopyELossFit.C | 130 ++++++++++ .../analysis2/scripts/RunCopyMergeEff.C | 91 +++++++ .../FORWARD/analysis2/scripts/RunCopySecMap.C | 131 ++++++++++ .../analysis2/scripts/RunCopyVtxBias.C | 92 +++++++ .../analysis2/scripts/RunMakeELossFit.C | 38 +++ PWG2/FORWARD/analysis2/scripts/TestAcc.C | 244 ++++++++++++++++++ 16 files changed, 1254 insertions(+), 5 deletions(-) create mode 100644 PWG2/FORWARD/analysis2/scripts/ExtractELoss.C create mode 100644 PWG2/FORWARD/analysis2/scripts/LoadPars.C create mode 100644 PWG2/FORWARD/analysis2/scripts/MakeCorrRepository.C create mode 100644 PWG2/FORWARD/analysis2/scripts/MakeELossFit.C create mode 100644 PWG2/FORWARD/analysis2/scripts/MakeESDChain.C create mode 100644 PWG2/FORWARD/analysis2/scripts/MoveCorrections.C create mode 100644 PWG2/FORWARD/analysis2/scripts/RunCopyELossFit.C create mode 100644 PWG2/FORWARD/analysis2/scripts/RunCopyMergeEff.C create mode 100644 PWG2/FORWARD/analysis2/scripts/RunCopySecMap.C create mode 100644 PWG2/FORWARD/analysis2/scripts/RunCopyVtxBias.C create mode 100644 PWG2/FORWARD/analysis2/scripts/RunMakeELossFit.C create mode 100644 PWG2/FORWARD/analysis2/scripts/TestAcc.C diff --git a/PWG2/FORWARD/analysis2/scripts/Compile.C b/PWG2/FORWARD/analysis2/scripts/Compile.C index dbfc35e38b6..26f618cbae3 100644 --- a/PWG2/FORWARD/analysis2/scripts/Compile.C +++ b/PWG2/FORWARD/analysis2/scripts/Compile.C @@ -24,7 +24,8 @@ Compile(const char* script, Option_t* option="g") 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} " diff --git a/PWG2/FORWARD/analysis2/scripts/DrawFits.C b/PWG2/FORWARD/analysis2/scripts/DrawFits.C index d3276ac7dc6..8aaf0adc2a9 100644 --- a/PWG2/FORWARD/analysis2/scripts/DrawFits.C +++ b/PWG2/FORWARD/analysis2/scripts/DrawFits.C @@ -1,3 +1,8 @@ +/** + * Script to draw the energy loss fits + * + * @ingroup pwg2_forward_analysis_scripts + */ #include #include #include @@ -24,8 +29,8 @@ TList* OpenFile(const char* fname) return 0; } - TList* forward = - static_cast(file->Get("PWG2forwardDnDeta/Forward")); + TList* forward = static_cast(file->Get("Forward")); + // static_cast(file->Get("PWG2forwardDnDeta/Forward")); if (!forward) { Error("DrawFits", "Couldn't get forward list from %s", fname); return 0; diff --git a/PWG2/FORWARD/analysis2/scripts/ExtractELoss.C b/PWG2/FORWARD/analysis2/scripts/ExtractELoss.C new file mode 100644 index 00000000000..2ff53e2d551 --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/ExtractELoss.C @@ -0,0 +1,96 @@ +/** + * Script to draw the energy loss fits + * + * @ingroup pwg2_forward_analysis_scripts + */ +#ifndef __CINT__ +#include +#include +#include +#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(file->Get("Forward")); + // static_cast(file->Get("PWG2forwardDnDeta/Forward")); + if (!forward) { + Error("ExtractELoss", "Couldn't get forward list from %s", fname); + return; + } + + TList* fitter = static_cast(forward->FindObject("fmdEnergyFitter")); + if (!fitter) { + Error("ExtractELoss", "Couldn't get fitter folder"); + return; + } + + TString cName(AliFMDCorrELossFit::Class()->GetName()); + + AliFMDCorrELossFit* obj = + static_cast(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 +// diff --git a/PWG2/FORWARD/analysis2/scripts/FitELoss.C b/PWG2/FORWARD/analysis2/scripts/FitELoss.C index c287933545d..8bb4a9e4f1e 100644 --- a/PWG2/FORWARD/analysis2/scripts/FitELoss.C +++ b/PWG2/FORWARD/analysis2/scripts/FitELoss.C @@ -1,3 +1,8 @@ +/** + * Test script to fit the energy loss spectra + * + * @ingroup pwg2_forward_analysis_scripts + */ #ifndef __CINT__ # include "AliForwardUtil.h" # include diff --git a/PWG2/FORWARD/analysis2/scripts/LoadLibs.C b/PWG2/FORWARD/analysis2/scripts/LoadLibs.C index 1799bcfdb0a..1532541e539 100644 --- a/PWG2/FORWARD/analysis2/scripts/LoadLibs.C +++ b/PWG2/FORWARD/analysis2/scripts/LoadLibs.C @@ -1,3 +1,8 @@ +/** + * Load the libraries of PWG2/FORWARD/analsysis2 + * + * @ingroup pwg2_forward_analysis_scripts + */ void LoadLibs() { @@ -8,7 +13,8 @@ LoadLibs() gSystem->Load("libANALYSIS"); gSystem->Load("libANALYSISalice"); gSystem->Load("libPWG0base"); - gSystem->Load("libPWG2forward"); gSystem->Load("libPWG2forward2"); } - +// +// EOF +// diff --git a/PWG2/FORWARD/analysis2/scripts/LoadPars.C b/PWG2/FORWARD/analysis2/scripts/LoadPars.C new file mode 100644 index 00000000000..590d712932d --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/LoadPars.C @@ -0,0 +1,20 @@ +/** + * 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 +// diff --git a/PWG2/FORWARD/analysis2/scripts/MakeCorrRepository.C b/PWG2/FORWARD/analysis2/scripts/MakeCorrRepository.C new file mode 100644 index 00000000000..34d94f84733 --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/MakeCorrRepository.C @@ -0,0 +1,27 @@ +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++; + } +} + diff --git a/PWG2/FORWARD/analysis2/scripts/MakeELossFit.C b/PWG2/FORWARD/analysis2/scripts/MakeELossFit.C new file mode 100644 index 00000000000..e8509ab6624 --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/MakeELossFit.C @@ -0,0 +1,228 @@ +#ifndef __CINT__ +#include "AliForwardCorrectionManager.h" +#include "AliFMDCorrELossFit.h" +#include +#include +#include +#include +#include +#include +#include +#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_GeV_kG_.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(file->Get("Forward")); + // static_cast(file->Get("PWG2forwardDnDeta/Forward")); + if (!forward) { + Error("MakeELossFit", "Couldn't get forward list from %s", filename); + return; + } + + fFitter = static_cast(forward->FindObject("fmdEnergyFitter")); + if (!fFitter) { + Error("MakeELossFit", "Couldn't get fitter folder"); + return; + } + + fAxis = static_cast(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(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(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(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(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(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 +// diff --git a/PWG2/FORWARD/analysis2/scripts/MakeESDChain.C b/PWG2/FORWARD/analysis2/scripts/MakeESDChain.C new file mode 100644 index 00000000000..39de328e8de --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/MakeESDChain.C @@ -0,0 +1,28 @@ +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(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; +} diff --git a/PWG2/FORWARD/analysis2/scripts/MoveCorrections.C b/PWG2/FORWARD/analysis2/scripts/MoveCorrections.C new file mode 100644 index 00000000000..e4495015a3a --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/MoveCorrections.C @@ -0,0 +1,107 @@ +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(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); + +} + diff --git a/PWG2/FORWARD/analysis2/scripts/RunCopyELossFit.C b/PWG2/FORWARD/analysis2/scripts/RunCopyELossFit.C new file mode 100644 index 00000000000..09aded15d31 --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/RunCopyELossFit.C @@ -0,0 +1,130 @@ +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 +// diff --git a/PWG2/FORWARD/analysis2/scripts/RunCopyMergeEff.C b/PWG2/FORWARD/analysis2/scripts/RunCopyMergeEff.C new file mode 100644 index 00000000000..641f6a2afab --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/RunCopyMergeEff.C @@ -0,0 +1,91 @@ +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 +// diff --git a/PWG2/FORWARD/analysis2/scripts/RunCopySecMap.C b/PWG2/FORWARD/analysis2/scripts/RunCopySecMap.C new file mode 100644 index 00000000000..6bc8ec4da48 --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/RunCopySecMap.C @@ -0,0 +1,131 @@ +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 +// diff --git a/PWG2/FORWARD/analysis2/scripts/RunCopyVtxBias.C b/PWG2/FORWARD/analysis2/scripts/RunCopyVtxBias.C new file mode 100644 index 00000000000..1da2b01b1b9 --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/RunCopyVtxBias.C @@ -0,0 +1,92 @@ +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 +// diff --git a/PWG2/FORWARD/analysis2/scripts/RunMakeELossFit.C b/PWG2/FORWARD/analysis2/scripts/RunMakeELossFit.C new file mode 100644 index 00000000000..6ffb77aca3c --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/RunMakeELossFit.C @@ -0,0 +1,38 @@ +/** + * 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 +// diff --git a/PWG2/FORWARD/analysis2/scripts/TestAcc.C b/PWG2/FORWARD/analysis2/scripts/TestAcc.C new file mode 100644 index 00000000000..b012f93701d --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/TestAcc.C @@ -0,0 +1,244 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//_____________________________________________________________________ +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(); +} -- 2.43.0