PWGLF/FORWARD/analysis2/corrs/RerunELossFits.C

   1 /**
   2  * Get a collection from a file directory
   3  *
   4  * @param dir   Parent directory
   5  * @param name  Name of collection
   6  *
   7  * @return collection or null
   8  */
   9 TCollection* GetCollection(TDirectory* dir, const TString& name)
  10 {
  11   if (!dir) {
  12     Error("GetCollection", "No parent directory for %s", name.Data());
  13     return 0;
  14   }
  15   TCollection* ret = 0;
  16   dir->GetObject(name, ret);
  17   if (!ret) {
  18     Error("GetCollection", "Couldn't find %s in %s",
  19           name.Data(), dir->GetName());
  20     return 0;
  21   }
  22   return ret;
  23 }
  24 /**
  25  * Get an object from a collection.  Optionally, we check that the
  26  * type of the possibly found object matches the request.
  27  *
  28  * @param parent Parent collection
  29  * @param name   Name of object
  30  * @param cls If specified, check that the found object (if any) is of
  31  * this class.
  32  *
  33  * @return Found object (possibly type-checked) or null
  34  */
  35 TObject* GetObject(const TCollection* parent, const TString& name,
  36                    const TClass* cls=0)
  37 {
  38   if (!parent) {
  39     Error("GetObject", "No parent collection for %s", name.Data());
  40     return 0;
  41   }
  42   TObject* ret = parent->FindObject(name);
  43   if (!ret) {
  44     Error("GetObject", "Couldn't find %s in %s",
  45           name.Data(), parent->GetName());
  46     return 0;
  47   }
  48   if (cls && !ret->IsA()->InheritsFrom(cls)) {
  49     Error("GetObject", "%s in %s is a %s, not a %s", name.Data(),
  50           parent->GetName(), ret->ClassName(), cls->GetName());
  51     return 0;
  52   }
  53   return ret;
  54 }
  55 /**
  56  * Get a collection contained in another collection
  57  *
  58  * @param parent Parent collection
  59  * @param name   Name of collection to find
  60  *
  61  * @return Found collection or null
  62  */
  63 TCollection* GetCollection(const TCollection* parent, const TString& name)
  64 {
  65   TObject* o = GetObject(parent, name, TCollection::Class());
  66   if (!o) return 0;
  67   return static_cast<TCollection*>(o);
  68 }
  69 /**
  70  * Re-run the energy loss fitter on a merged output file
  71  *
  72  * @param input  File name of merged output file
  73  * @param output If specified, the file the new results are written
  74  * to.  If this is not specified, it defaults to the name of the input
  75  * file with "_rerun" attached to the base name
  76  *
  77  * @param forceSet  Forcibly set things
  78  * @param input     Input file
  79  * @param output    Output file
  80  */
  81 void RerunELossFits(Bool_t forceSet=false,
  82                     const TString& input="forward_eloss.root",
  83                     Bool_t shift=true,
  84                     const TString& output="")
  85 {
  86   const char* fwd = "$ALICE_ROOT/PWGLF/FORWARD/analysis2";
  87   gROOT->Macro(Form("%s/scripts/LoadLibs.C", fwd));
  88
  89   TFile*  inFile  = 0;
  90   TFile*  outFile = 0;
  91   TString outName(output);
  92   if (outName.IsNull()) {
  93     outName = input;
  94     outName.ReplaceAll(".root", "_rerun.root");
  95   }
  96   Bool_t  allOk = false;
  97   try {
  98     // --- Open input file ---------------------------------------------
  99     inFile = TFile::Open(input, "READ");
 100     if (!inFile)
 101       throw TString::Format("Failed to open %s", input.Data());
 102
 103     // --- InFiled input collections --------------------------------------
 104     TCollection* inFwdSum = GetCollection(inFile, "ForwardELossSums");
 105     if (!inFwdSum) throw TString("Cannot proceed without sums");
 106
 107     TCollection* inFwdRes = GetCollection(inFile, "ForwardELossResults");
 108     if (!inFwdRes) throw TString("Cannot proceed with merged list");
 109
 110     TCollection* inEFSum = GetCollection(inFwdRes, "fmdEnergyFitter");
 111     if (!inEFSum) throw TString("Cannot proceed without accumulated data");
 112
 113     TCollection* inEFRes = GetCollection(inFwdRes, "fmdEnergyFitter");
 114     if (!inEFRes) throw TString("Cannot proceed without previous results");
 115
 116     // --- Open output file --------------------------------------------
 117     outFile = TFile::Open(outName, "RECREATE");
 118     if (!outFile)
 119       throw TString::Format("Failed to open %s", outName.Data());
 120
 121     // --- Write copy of sum collection to output --------------------
 122     TCollection* outFwdSum = static_cast<TCollection*>(inFwdSum->Clone());
 123     outFile->cd();
 124     outFwdSum->Write(inFwdSum->GetName(), TObject::kSingleKey);
 125
 126     // --- Make our fitter object ------------------------------------
 127     AliFMDEnergyFitter* fitter = new AliFMDEnergyFitter("energy");
 128     fitter->SetDoFits(true);
 129     fitter->SetEnableDeltaShift(shift);
 130     fitter->Init();
 131     if (forceSet || !fitter->ReadParameters(inEFSum)) {
 132       Printf("Forced settings");
 133
 134       const TAxis* etaAxis = static_cast<TAxis*>(GetObject(inEFSum,"etaAxis"));
 135       if (!etaAxis) throw TString("Cannot proceed without eta axis");
 136       fitter->SetEtaAxis(*etaAxis);
 137
 138       // Set maximum energy loss to consider
 139       fitter->SetMaxE(15);
 140       // Set number of energy loss bins
 141       fitter->SetNEbins(500);
 142       // Set whether to use increasing bin sizes
 143       // fitter->SetUseIncreasingBins(true);
 144       // Set whether to do fit the energy distributions
 145       fitter->SetDoFits(kTRUE);
 146       // Set whether to make the correction object
 147       fitter->SetDoMakeObject(kTRUE);
 148       // Set the low cut used for energy
 149       fitter->SetLowCut(0.4);
 150       // Set the number of bins to subtract from maximum of distributions
 151       // to get the lower bound of the fit range
 152       fitter->SetFitRangeBinWidth(4);
 153       // Set the maximum number of landaus to try to fit (max 5)
 154       fitter->SetNParticles(5);
 155       // Set the minimum number of entries in the distribution before
 156       // trying to fit to the data - 10k seems the least we can do
 157       fitter->SetMinEntries(10000);
 158       // fitter->SetMaxChi2PerNDF(10);
 159       // Enable debug
 160     }
 161     fitter->SetDebug(3);
 162     fitter->SetStoreResiduals(AliFMDEnergyFitter::kResidualSquareDifference);
 163     fitter->SetRegularizationCut(1e6); // Lower by factor 3
 164     // Set the number of bins to subtract from maximum of distributions
 165     // to get the lower bound of the fit range
 166     // fitter->SetFitRangeBinWidth(2);
 167     // Skip all of FMD2 and 3
 168     // fitter->SetSkips(AliFMDEnergyFitter::kFMD2|AliFMDEnergyFitter::kFMD3);
 169
 170     // --- Now do the fits -------------------------------------------
 171     fitter->Print();
 172     outFwdSum->ls("R");
 173     fitter->Fit(static_cast<TList*>(outFwdSum));
 174
 175     // --- Copy full result folder -----------------------------------
 176     TCollection* outFwdRes = static_cast<TCollection*>(inFwdRes->Clone());
 177     // Remove old fits
 178     TCollection* outEFRes = GetCollection(outFwdRes, "fmdEnergyFitter");
 179     outFwdRes->Remove(outEFRes);
 180     // Make our new fit results folder, and add it to results folder
 181     TCollection* tmp = GetCollection(outFwdSum, "fmdEnergyFitter");
 182     outEFRes = static_cast<TCollection*>(tmp->Clone());
 183     outEFRes->Add(new TNamed("refitted", "Refit of the data"));
 184     outFwdRes->Add(outEFRes);
 185
 186     // --- Write out new results folder ------------------------------
 187     outFile->cd();
 188     outFwdRes->Write(inFwdRes->GetName(), TObject::kSingleKey);
 189     Printf("Wrote results to \"%s\" (%s)", outName.Data(), outFile->GetName());
 190     allOk = true;
 191   }
 192   catch (const TString& e) {
 193     Error("RerunELossFits", e);
 194   }
 195   if (inFile)  inFile->Close();
 196   if (outFile) {
 197     Printf("Wrote new output to \"%s\"", outName.Data());
 198     outFile->Close();
 199   }
 200
 201   if (allOk)
 202     gROOT->Macro(Form("%s/corrs/DrawCorrELoss.C(false,false,\"%s\")",
 203                       fwd, outName.Data()));
 204 }
 205
 206