Added some more scripts
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 16 Dec 2010 11:46:56 +0000 (11:46 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 16 Dec 2010 11:46:56 +0000 (11:46 +0000)
16 files changed:
PWG2/FORWARD/analysis2/scripts/Compile.C
PWG2/FORWARD/analysis2/scripts/DrawFits.C
PWG2/FORWARD/analysis2/scripts/ExtractELoss.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/FitELoss.C
PWG2/FORWARD/analysis2/scripts/LoadLibs.C
PWG2/FORWARD/analysis2/scripts/LoadPars.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/MakeCorrRepository.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/MakeELossFit.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/MakeESDChain.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/MoveCorrections.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/RunCopyELossFit.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/RunCopyMergeEff.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/RunCopySecMap.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/RunCopyVtxBias.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/RunMakeELossFit.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/scripts/TestAcc.C [new file with mode: 0644]

index dbfc35e38b69df3e89998268b604b896fdee6a9f..26f618cbae3348bb9e9fde90cea8e77f20ed2c12 100644 (file)
@@ -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} " 
index d3276ac7dc6157018f8ef4d8998830ed791726f6..8aaf0adc2a98840ba8f3c6e7173635255d49810d 100644 (file)
@@ -1,3 +1,8 @@
+/**
+ * Script to draw the energy loss fits 
+ * 
+ * @ingroup pwg2_forward_analysis_scripts
+ */
 #include <TFile.h>
 #include <THStack.h>
 #include <TList.h>
@@ -24,8 +29,8 @@ TList* OpenFile(const char* fname)
     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;
diff --git a/PWG2/FORWARD/analysis2/scripts/ExtractELoss.C b/PWG2/FORWARD/analysis2/scripts/ExtractELoss.C
new file mode 100644 (file)
index 0000000..2ff53e2
--- /dev/null
@@ -0,0 +1,96 @@
+/**
+ * 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
+//
index c287933545d8166505a1e263786eefc349150e53..8bb4a9e4f1e9ea4a62eafd91661fd1fd3ac8f7ec 100644 (file)
@@ -1,3 +1,8 @@
+/**
+ * Test script to fit the energy loss spectra 
+ * 
+ * @ingroup pwg2_forward_analysis_scripts
+ */
 #ifndef __CINT__
 # include "AliForwardUtil.h"
 # include <TFile.h>
index 1799bcfdb0aed20eb7f7bc913408a690e363b609..1532541e539a9b1f74bf2e3973d752a1010f1025 100644 (file)
@@ -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 (file)
index 0000000..590d712
--- /dev/null
@@ -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 (file)
index 0000000..34d94f8
--- /dev/null
@@ -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 (file)
index 0000000..e8509ab
--- /dev/null
@@ -0,0 +1,228 @@
+#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
+//
diff --git a/PWG2/FORWARD/analysis2/scripts/MakeESDChain.C b/PWG2/FORWARD/analysis2/scripts/MakeESDChain.C
new file mode 100644 (file)
index 0000000..39de328
--- /dev/null
@@ -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<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;
+}
diff --git a/PWG2/FORWARD/analysis2/scripts/MoveCorrections.C b/PWG2/FORWARD/analysis2/scripts/MoveCorrections.C
new file mode 100644 (file)
index 0000000..e449501
--- /dev/null
@@ -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<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);       
+
+}
+  
diff --git a/PWG2/FORWARD/analysis2/scripts/RunCopyELossFit.C b/PWG2/FORWARD/analysis2/scripts/RunCopyELossFit.C
new file mode 100644 (file)
index 0000000..09aded1
--- /dev/null
@@ -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 (file)
index 0000000..641f6a2
--- /dev/null
@@ -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 (file)
index 0000000..6bc8ec4
--- /dev/null
@@ -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 (file)
index 0000000..1da2b01
--- /dev/null
@@ -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 (file)
index 0000000..6ffb77a
--- /dev/null
@@ -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 (file)
index 0000000..b012f93
--- /dev/null
@@ -0,0 +1,244 @@
+#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();
+}