PWG2/FORWARD/analysis2/scripts/RunCopyELossFit.C

   1 Color_t Color(UShort_t d, Char_t r ) const
   2 {
   3   return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
   4           + ((r == 'I' || r == 'i') ? 2 : -2));
   5 }
   6
   7 void
   8 RunCopyELossFit(UShort_t sys, UShort_t cms, Short_t field, Bool_t mc=false,
   9                 const Char_t* path=0)
  10 {
  11   RunCopyELossFit(sys == 1 ? "pp" : "PbPb",
  12                   cms,
  13                   field,
  14                   mc,
  15                   path);
  16 }
  17
  18 /**
  19  *
  20  * @param sys       Collision system
  21  * @param cms       Center of mass energy per nucleon in GeV
  22  * @param field     Magnetic field
  23  *
  24  * @ingroup pwg2_forward_analysis_scripts
  25  */
  26 void
  27 RunCopyELossFit(const char* sys, UShort_t cms, Short_t field, bool mc=false,
  28                 const char* path=0)
  29 {
  30   gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
  31   gSystem->Load("libPWG2forward.so");
  32
  33   AliFMDAnaParameters* p = AliFMDAnaParameters::Instance();
  34   p->SetRealData(!mc);
  35   p->SetEnergy(Float_t(cms));
  36   p->SetMagField(Float_t(field));
  37   p->SetCollisionSystem(sys);
  38   if (path) {
  39     p->SetBackgroundPath(path);
  40     p->SetEnergyPath(path);
  41     p->SetEventSelectionPath(path);
  42     p->SetSharingEfficiencyPath(path);
  43   }
  44   p->Init(true, AliFMDAnaParameters::kBackgroundCorrection|
  45           AliFMDAnaParameters::kEnergyDistributions);
  46
  47   Int_t    nVtx   = p->GetNvtxBins();
  48   Double_t minVtx = -p->GetVtxCutZ();
  49   Double_t maxVtx = -p->GetVtxCutZ();
  50   Int_t    nEta   = p->GetNetaBins();
  51   Double_t minEta = p->GetEtaMin();
  52   Double_t maxEta = p->GetEtaMax();
  53
  54   TAxis vtxAxis(nVtx, minVtx, maxVtx);
  55   TAxis etaAxis(nEta, minEta, maxEta);
  56   AliFMDCorrELossFit* m = new AliFMDCorrELossFit;
  57   m->SetEtaAxis(etaAxis);
  58
  59   for (UShort_t d = 1; d <= 3; d++) {
  60     UShort_t nQ = (d == 1 ? 1 : 2);
  61     for (UShort_t q = 0; q < nQ; q++) {
  62       Char_t r = (q == 0 ? 'I' : 'O');
  63
  64       for (UShort_t b = 1; b < nEta; b++) {
  65         Double_t eta = etaAxis.GetBinCenter(b);
  66         Info("RunCopyELossFit", "FMD%d%c, bin %3d, eta %+8.4f", d, r, b, eta);
  67         if (eta < minEta || eta > maxEta) continue;
  68         TH1F* oldhist = p->GetEnergyDistribution(d,r,eta);
  69         if (!oldhist) {
  70           Warning("RunCopyELossFit",
  71                   "Didn't find energy distribution for FMD%d%c, eta bin %3d",
  72                   d, r, b);
  73           continue;
  74         }
  75         TF1*  oldfunc = oldhist->GetFunction("FMDfitFunc");
  76         if (!oldfunc) {
  77           // Warning("RunCopyELossFit",
  78           // "Didn't find energy fit for FMD%d%c, eta bin %3d",
  79           // d, r, b);
  80           continue;
  81         }
  82         Double_t chi2   = oldfunc->GetChisquare();
  83         UShort_t nu     = oldfunc->GetNDF();
  84         Double_t c      = oldfunc->GetParameter(0);
  85         Double_t delta  = oldfunc->GetParameter(1);
  86         Double_t xi     = oldfunc->GetParameter(2);
  87         Double_t a2 = 0, a3 = 0, ea2 = 0, ea3 = 0;
  88         if (oldfunc->GetNpar() >= 4) {
  89           a2     = oldfunc->GetParameter(3);
  90           ea2    = oldfunc->GetParError(3);
  91         }
  92         if (oldfunc->GetNpar() >= 5) {
  93           a3     = oldfunc->GetParameter(4);
  94           ea3    = oldfunc->GetParError(4);
  95         }
  96         Int_t    n      = (a3 < 1e-5 ? (a2 < 1e-5 ? 1 : 2) : 3);
  97         Double_t ec     = oldfunc->GetParError(0);
  98         Double_t edelta = oldfunc->GetParError(1);
  99         Double_t exi    = oldfunc->GetParError(2);
 100         delta           = delta - xi * -0.22278298;
 101         edelta          = TMath::Sqrt(TMath::Power(edelta,2) -
 102                                       TMath::Power(0.22278298*exi,2));
 103         Info("RunCopyELossFit", "FMD%d%c etabin=%3d: n=%d\n"
 104              "  C=%f+/-%f, Delta=%f+/-%f, xi=%f+/-%f, a2=%f+/-%f, a3=%f+/-%f",
 105              d, r, b, n, c, ec, delta, edelta, xi, exi, a2, ea2, a3, ea3);
 106         Double_t a[]  = { a2, a3, -9999 };
 107         Double_t ea[] = { ea2, ea3, -9999 };
 108
 109         AliFMDCorrELossFit::ELossFit* fit = new
 110           AliFMDCorrELossFit::ELossFit(0, n,
 111                                        chi2,   nu,
 112                                        c,      ec,
 113                                        delta,  edelta,
 114                                        xi,     exi,
 115                                        0,      0,
 116                                        0,      0,
 117                                        a,      ea);
 118         fit->CalculateQuality();
 119         fit->fBin = b;
 120         fit->fDet = d;
 121         fit->fRing = r;
 122         fit->Print();
 123
 124         Info("RunCopyELossFit", "Setting fit FMD%d%c etabin=%3d: %p",
 125              d, r, b, fit);
 126         m->SetFit(d, r, b, fit);
 127       }
 128     }
 129   }
 130   UShort_t isys = AliForwardUtil::ParseCollisionSystem(sys);
 131   AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
 132   TString fname(mgr.GetFileName(AliForwardCorrectionManager::kELossFits,
 133                                 isys, cms, field, mc));
 134   TFile* output = TFile::Open(fname.Data(), "RECREATE");
 135   if (!output) {
 136     Warning("RunCopyELossFit", "Failed to open output file %s", fname.Data());
 137     return kFALSE;
 138   }
 139   m->Write(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
 140   output->Write();
 141   output->Close();
 142   Info("RunCopyELossFit",
 143        "File %s created.  It should be copied to %s and stored in SVN",
 144        fname.Data(), mgr.GetFileDir(AliForwardCorrectionManager::kELossFits));
 145
 146 }
 147 //
 148 // EOF
 149 //