]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx
corrected obvious typo. which value to chose between 0.045 and 0.050 is a bit matter...
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDEnergyFitter.cxx
index caedd41a879d72ed2c939818218c445b5acf8255..8c7eb80e004ca4aa33c8106c21492eecc8d33776 100644 (file)
@@ -16,7 +16,6 @@
 #include <TROOT.h>
 #include <iostream>
 #include <iomanip>
-#include <TFile.h>
 
 ClassImp(AliFMDEnergyFitter)
 #if 0
@@ -31,13 +30,14 @@ namespace {
 AliFMDEnergyFitter::AliFMDEnergyFitter()
   : TNamed(), 
     fRingHistos(),
-    fLowCut(0.3),
-    fNParticles(3),
-    fMinEntries(100),
+    fLowCut(0.4),
+    fNParticles(5),
+    fMinEntries(1000),
     fFitRangeBinWidth(4),
     fDoFits(false),
     fDoMakeObject(false),
     fEtaAxis(),
+    fCentralityAxis(),
     fMaxE(10),
     fNEbins(300), 
     fUseIncreasingBins(true),
@@ -55,13 +55,14 @@ AliFMDEnergyFitter::AliFMDEnergyFitter()
 AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
   : TNamed("fmdEnergyFitter", title), 
     fRingHistos(), 
-    fLowCut(0.3),
-    fNParticles(3),
-    fMinEntries(100),
+    fLowCut(0.4),
+    fNParticles(5),
+    fMinEntries(1000),
     fFitRangeBinWidth(4),
     fDoFits(false),
     fDoMakeObject(false),
     fEtaAxis(0,0,0),
+    fCentralityAxis(0,0,0),
     fMaxE(10),
     fNEbins(300), 
     fUseIncreasingBins(true),
@@ -78,6 +79,8 @@ AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
   //
   fEtaAxis.SetName("etaAxis");
   fEtaAxis.SetTitle("#eta");
+  fCentralityAxis.SetName("centralityAxis");
+  fCentralityAxis.SetTitle("Centrality [%]");
   fRingHistos.Add(new RingHistos(1, 'I'));
   fRingHistos.Add(new RingHistos(2, 'I'));
   fRingHistos.Add(new RingHistos(2, 'O'));
@@ -96,6 +99,7 @@ AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
     fDoFits(o.fDoFits),
     fDoMakeObject(o.fDoMakeObject),
     fEtaAxis(o.fEtaAxis),
+    fCentralityAxis(o.fCentralityAxis),
     fMaxE(o.fMaxE),
     fNEbins(o.fNEbins), 
     fUseIncreasingBins(o.fUseIncreasingBins),
@@ -145,7 +149,17 @@ AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
   fFitRangeBinWidth= o.fFitRangeBinWidth;
   fDoFits        = o.fDoFits;
   fDoMakeObject  = o.fDoMakeObject;
-  fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax());
+  fEtaAxis.Set(o.fEtaAxis.GetNbins(),
+              o.fEtaAxis.GetXmin(),
+              o.fEtaAxis.GetXmax());
+  if (o.fCentralityAxis.GetXbins()) {
+    const TArrayD* bins = o.fCentralityAxis.GetXbins();
+    fCentralityAxis.Set(bins->GetSize()-1,bins->GetArray());
+  }
+  else 
+    fCentralityAxis.Set(o.fCentralityAxis.GetNbins(),
+                       o.fCentralityAxis.GetXmin(),
+                       o.fCentralityAxis.GetXmax());
   fDebug         = o.fDebug;
   fMaxRelParError= o.fMaxRelParError;
   fMaxChi2PerNDF = o.fMaxChi2PerNDF;
@@ -199,10 +213,16 @@ AliFMDEnergyFitter::Init(const TAxis& eAxis)
   if (fEtaAxis.GetNbins() == 0 || 
       TMath::Abs(fEtaAxis.GetXmax() - fEtaAxis.GetXmin()) < 1e-7) 
     SetEtaAxis(eAxis);
+  if (fCentralityAxis.GetNbins() == 0) {
+    UShort_t n = 12;
+    Double_t bins[] = {  0.,  5., 10., 15., 20., 30., 
+                        40., 50., 60., 70., 80., 100. };
+    SetCentralityAxis(n, bins);
+  }
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
   while ((o = static_cast<RingHistos*>(next())))
-    o->Init(fEtaAxis, fMaxE, fNEbins, fUseIncreasingBins);
+    o->Init(fEtaAxis, fCentralityAxis, fMaxE, fNEbins, fUseIncreasingBins);
 }  
 //____________________________________________________________________
 void
@@ -238,10 +258,18 @@ AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax)
   //
   fEtaAxis.Set(nBins,etaMin,etaMax);
 }
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::SetCentralityAxis(UShort_t n, Double_t* bins)
+{
+  fCentralityAxis.Set(n-1, bins);
+}
+
 
 //____________________________________________________________________
 Bool_t
-AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, 
+AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
+                              Double_t         cent,
                               Bool_t           empty)
 {
   // 
@@ -249,11 +277,15 @@ AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
   // 
   // Parameters:
   //    input     Input 
+  //    cent      Centrality 
   //    empty     Whether the event is 'empty'
   // 
   // Return:
   //    True on success, false otherwise 
   //
+  Int_t icent = fCentralityAxis.FindBin(cent);
+  if (icent < 1 || icent > fCentralityAxis.GetNbins()) icent = 0;
+
   for(UShort_t d = 1; d <= 3; d++) {
     Int_t nRings = (d == 1 ? 1 : 2);
     for (UShort_t q = 0; q < nRings; q++) {
@@ -276,7 +308,7 @@ AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
          Int_t ieta = fEtaAxis.FindBin(eta1);
          if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue; 
 
-         histos->Fill(empty, ieta-1, mult);
+         histos->Fill(empty, ieta-1, icent, mult);
 
        } // for strip
       } // for sector
@@ -334,13 +366,6 @@ AliFMDEnergyFitter::Fit(const TList* dir)
   if (!fDoMakeObject) return;
 
   MakeCorrectionsObject(d);
-  d->ls();
-  TDirectory* savdir = gDirectory;
-  TFile* tmp = TFile::Open("elossfits.root", "RECREATE");
-  d->Write();
-  tmp->Write();
-  tmp->Close();
-  savdir->cd();
 }
 
 //____________________________________________________________________
@@ -370,6 +395,7 @@ AliFMDEnergyFitter::MakeCorrectionsObject(TList* d)
   TString fileName(mgr.GetFilePath(AliForwardCorrectionManager::kELossFits));
   AliInfo(Form("Object %s created in output - should be extracted and copied "
               "to %s", oName.Data(), fileName.Data()));
+  d->Add(new TNamed("filename", fileName.Data()));
   d->Add(obj, oName.Data());
 }
 
@@ -425,7 +451,7 @@ AliFMDEnergyFitter::Print(Option_t*) const
   char ind[gROOT->GetDirLevel()+1];
   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
   ind[gROOT->GetDirLevel()] = '\0';
-  std::cout << ind << "AliFMDEnergyFitter: " << GetName() << '\n'
+  std::cout << ind << ClassName() << ": " << GetName() << '\n'
            << ind << " Low cut:                " << fLowCut << " E/E_mip\n"
            << ind << " Max(particles):         " << fNParticles << '\n'
            << ind << " Min(entries):           " << fMinEntries << '\n'
@@ -555,14 +581,16 @@ AliFMDEnergyFitter::RingHistos::~RingHistos()
 
 //____________________________________________________________________
 void
-AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, Double_t mult)
+AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, 
+                                    Int_t /* icent */,  Double_t mult)
 {
   // 
   // Fill histogram 
   // 
   // Parameters:
   //    empty  True if event is empty
-  //    ieta   Eta bin
+  //    ieta   Eta bin (0-based)
+  //    icent  Centrality bin (1-based)
   //    mult   Signal 
   //
   if (empty) { 
@@ -571,6 +599,9 @@ AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, Double_t mult)
   }
   fEDist->Fill(mult);
   
+  // Here, we should first look up in a centrality array, and then in
+  // that array look up the eta bin - or perhaps we should do it the
+  // other way around?
   if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return;
   
   TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta));
@@ -729,6 +760,7 @@ AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name,
 //____________________________________________________________________
 void
 AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis, 
+                                    const TAxis& /* cAxis */,
                                     Double_t maxDE, Int_t nDEbins, 
                                     Bool_t useIncrBin)
 {
@@ -749,13 +781,67 @@ AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis,
                         fName.Data()), 200, 0, 6);
   fList->Add(fEDist);
   fList->Add(fEmpty);
+  // Here, we should make an array of centrality bins ...
+  // In that case, the structure will be 
+  // 
+  //    -+- Ring1 -+- Centrality1 -+- Eta1
+  //     |         |               +- Eta2
+  //     ...       ...             ...
+  //     |         |               +- EtaM
+  //     |         +- Centrality2 -+- Eta1
+  //     ...       ...             ...
+  //     |         |               +- EtaM
+  //     ...       ...
+  //     |         +- CentralityN -+- Eta1
+  //     ...       ...             ...
+  //     |                         +- EtaM
+  //    -+- Ring2 -+- Centrality1 -+- Eta1
+  //     |         |               +- Eta2
+  //     ...       ...             ...
+  //     |         |               +- EtaM
+  //     |         +- Centrality2 -+- Eta1
+  //     ...       ...             ...
+  //     |         |               +- EtaM
+  //     ...       ...
+  //     |         +- CentralityN -+- Eta1
+  //     ...       ...             ...
+  //     |                         +- EtaM
+  //     ...       ...             ...
+  //    -+- RingO -+- Centrality1 -+- Eta1
+  //               |               +- Eta2
+  //               ...             ...
+  //               |               +- EtaM
+  //               +- Centrality2 -+- Eta1
+  //               ...             ...
+  //               |               +- EtaM
+  //               ...
+  //               +- CentralityN -+- Eta1
+  //               ...             ...
+  //                               +- EtaM
+  // 
   // fEtaEDists.Expand(eAxis.GetNbins());
   for (Int_t i = 1; i <= eAxis.GetNbins(); i++) { 
     Double_t min = eAxis.GetBinLowEdge(i);
     Double_t max = eAxis.GetBinUpEdge(i);
+    // Or may we should do it here?  In that case we'd have the
+    // following structure:
+    //     Ring1 -+- Eta1 -+- Centrality1 
+    //            |        +- Centrality2
+    //            ...      ...
+    //            |        +- CentralityN
+    //            +- Eta2 -+- Centrality1 
+    //            |        +- Centrality2
+    //            ...      ...
+    //            |        +- CentralityN
+    //            ...
+    //            +- EtaM -+- Centrality1 
+    //                     +- Centrality2
+    //                     ...
+    //                     +- CentralityN
     Make(i, min, max, maxDE, nDEbins, useIncrBin);
   }
   fList->Add(&fEtaEDists);
+  // fEtaEDists.ls();
 }
 //____________________________________________________________________
 TObjArray*
@@ -1005,7 +1091,9 @@ AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
   if (nParticles == 1) {
     TF1* r = f.Fit1Particle(dist, 0);
     if (!r) return 0;
-    return new TF1(*r);
+    TF1* ret = new TF1(*r);
+    dist->GetListOfFunctions()->Add(ret);
+    return ret;
   }
 
   // Fit from 2 upto n particles  
@@ -1013,15 +1101,15 @@ AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
 
 
   // Now, we need to select the best fit 
-  Int_t nFits = f.fFitResults.GetEntriesFast();
+  Int_t nFits = f.GetFitResults().GetEntriesFast();
   TF1*  good[nFits];
   for (Int_t i = nFits-1; i >= 0; i--) { 
     good[i] = 0;
-    TF1* ff = static_cast<TF1*>(f.fFunctions.At(i));
+    TF1* ff = static_cast<TF1*>(f.GetFunctions().At(i));
     // ff->SetLineColor(Color());
     ff->SetRange(0, maxRange);
     dist->GetListOfFunctions()->Add(new TF1(*ff));
-    if (CheckResult(static_cast<TFitResult*>(f.fFitResults.At(i)),
+    if (CheckResult(static_cast<TFitResult*>(f.GetFitResults().At(i)),
                    relErrorCut, chi2nuCut)) {
       good[i] = ff;
       ff->SetLineWidth(2);
@@ -1029,20 +1117,20 @@ AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
     }
   }
   // If all else fails, use the 1 particle fit 
-  TF1* ret = static_cast<TF1*>(f.fFunctions.At(0));
+  TF1* ret = static_cast<TF1*>(f.GetFunctions().At(0));
 
   // Find the fit with the most valid particles 
   for (Int_t i = nFits-1; i > 0; i--) {
     if (!good[i]) continue;
     if (fDebug > 1) {
       AliInfo(Form("Choosing fit with n=%d", i+1));
-      f.fFitResults.At(i)->Print();
+      f.GetFitResults().At(i)->Print();
     }
     ret = good[i];
     break;
   }
   // Give a warning if we're using fall-back 
-  if (ret == f.fFunctions.At(0)) {
+  if (ret == f.GetFunctions().At(0)) {
     AliWarning("Choosing fall-back 1 particle fit");
   }
   // Copy our result and return (the functions are owned by the fitter)
@@ -1133,7 +1221,7 @@ AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r,
 
 //__________________________________________________________________
 void
-AliFMDEnergyFitter::RingHistos::FindBestFits(TList*              d, 
+AliFMDEnergyFitter::RingHistos::FindBestFits(const TList*        d, 
                                             AliFMDCorrELossFit& obj,
                                             const TAxis&        eta,     
                                             Double_t            relErrorCut, 
@@ -1189,7 +1277,7 @@ AliFMDEnergyFitter::RingHistos::FindBestFits(TList*              d,
 
 //__________________________________________________________________
 AliFMDCorrELossFit::ELossFit* 
-AliFMDEnergyFitter::RingHistos::FindBestFit(TH1* dist,
+AliFMDEnergyFitter::RingHistos::FindBestFit(const TH1* dist,
                                            Double_t relErrorCut, 
                                            Double_t chi2nuCut,
                                            Double_t minWeightCut)