AliPhysics  b76e98e (b76e98e)
RerunELossFits.C
Go to the documentation of this file.
1 
10 TCollection* GetCollection(TDirectory* dir, const TString& name,
11  Bool_t verbose=true)
12 {
13  if (!dir) {
14  Error("GetCollection", "No parent directory for %s", name.Data());
15  return 0;
16  }
17  TCollection* ret = 0;
18  dir->GetObject(name, ret);
19  if (!ret) {
20  if (verbose)
21  Error("GetCollection", "Couldn't find %s in %s",
22  name.Data(), dir->GetName());
23  return 0;
24  }
25  return ret;
26 }
39 TObject* GetObject(const TCollection* parent, const TString& name,
40  const TClass* cls=0, Bool_t verbose=true)
41 {
42  if (!parent) {
43  if (verbose)
44  Error("GetObject", "No parent collection for %s", name.Data());
45  return 0;
46  }
47  TObject* ret = parent->FindObject(name);
48  if (!ret) {
49  if (verbose)
50  Error("GetObject", "Couldn't find %s in %s",
51  name.Data(), parent->GetName());
52  return 0;
53  }
54  if (cls && !ret->IsA()->InheritsFrom(cls)) {
55  if (verbose)
56  Error("GetObject", "%s in %s is a %s, not a %s", name.Data(),
57  parent->GetName(), ret->ClassName(), cls->GetName());
58  return 0;
59  }
60  return ret;
61 }
62 void RemoveObject(TCollection* parent, const TString& name)
63 {
64  while (parent->GetEntries() > 0) {
65  TObject* o = GetObject(parent, name, 0, false);
66  if (!o) return;
67  Info("RemoveObject", "Removing %s from %s", name.Data(), parent->GetName());
68  TObject* ret = parent->Remove(o);
69  if (!ret || ret != o)
70  Warning("RemoveObject", "Failed to remove %s from %s (%p %p)",
71  name.Data(), parent->GetName(), o, ret);
72  // delete o;
73  }
74 }
75 
85 TCollection* GetCollection(const TCollection* parent, const TString& name,
86  Bool_t verbose=true)
87 {
88  TObject* o = GetObject(parent, name, TCollection::Class(), verbose);
89  if (!o) return 0;
90  return static_cast<TCollection*>(o);
91 }
92 
94 {
95  const char* other[] = { "fmdESDFixer",
96  "fmdSharingFilter",
97  "fmdDensityCalculator",
98  0 };
99  const char** ptr = other;
100  while (*ptr) { RemoveObject(c, *ptr); ptr++; }
101 
102  TCollection* ef = GetCollection(c, "fmdEnergyFitter");
103  const char* stacks[] = { "chi2", "c", "delta", "xi", "sigma", "sigman",
104  "a2", "a3", "a4", "a5", 0 };
105  ptr = stacks;
106  while (*ptr) { RemoveObject(ef, *ptr); ptr++; }
107 
108  const char* dets[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
109  ptr = dets;
110  while (*ptr) {
111  TCollection* det = GetCollection(ef, *ptr);
112  if (!det) {
113  ptr++;
114  continue;
115  }
116  const char* subs[] = { "elossDists", "elossResults", "elossResiduals", 0 };
117  const char** sub = subs;
118  while (*sub) { RemoveObject(det, *sub); sub++; }
119 
120  // TH1* eloss = static_cast<TH1*>(GetObject(det,"eloss",TH1::Class()));
121  // if (eloss) Printf("%s %d entries", *ptr, eloss->GetEntries());
122  ptr++;
123  }
124  // c->ls();
125 }
126 
141 void RerunELossFits(Bool_t forceSet=false,
142  const TString& input="forward_eloss.root",
143  Bool_t shift=true,
144  const TString& output="",
145  UShort_t flags=0x1)
146 {
147  const char* fwd = "$ALICE_PHYSICS/PWGLF/FORWARD/analysis2";
148  gROOT->Macro(Form("%s/scripts/LoadLibs.C", fwd));
149 
150  TFile* inFile = 0;
151  TFile* outFile = 0;
152  TString outName(output);
153  if (outName.IsNull()) {
154  outName = input;
155  outName.ReplaceAll(".root", "_rerun.root");
156  }
157  Bool_t allOk = false;
158  try {
159  // --- Open input file ---------------------------------------------
160  inFile = TFile::Open(input, "READ");
161  if (!inFile)
162  throw TString::Format("Failed to open %s", input.Data());
163 
164  // --- InFiled input collections --------------------------------------
165  TCollection* inFwdSum = GetCollection(inFile, "ForwardELossSums");
166  if (!inFwdSum) inFwdSum = GetCollection(inFile, "forwardQAResults");
167  if (!inFwdSum) throw new TString("Cannot proceed without sums");
168  inFwdSum->SetName("ForwardELossSums");
169 
170  TCollection* inFwdRes = GetCollection(inFile, "ForwardELossResults");
171  if (!inFwdRes) inFwdRes = GetCollection(inFile, "forwardQAResults");
172  if (!inFwdRes) {
173  inFwdRes =
174  static_cast<TCollection*>(inFwdSum->Clone("ForwardELossResults"));
175  // throw new TString("Cannot proceed with merged list");
176  }
177  inFwdRes->SetName("ForwardELossResults");
178 
179  TCollection* inEFSum = GetCollection(inFwdRes, "fmdEnergyFitter");
180  if (!inEFSum) throw new TString("Cannot proceed without accumulated data");
181 
182  TCollection* inEFRes = GetCollection(inFwdRes, "fmdEnergyFitter");
183  if (!inEFRes) throw new TString("Cannot proceed without previous results");
184 
185  TCollection* inESSum = GetCollection(inFwdRes, "fmdEventInspector");
186  Int_t sys = 0;
187  Long_t nEvAcc = 0;
188  if (inESSum) {
189  TObject* oSys = GetObject(inESSum, "sys", TParameter<int>::Class());
190  if (oSys) sys = (static_cast<TParameter<int>*>(oSys))->GetVal();
191  TH1* oEvAcc = static_cast<TH1*>(GetObject(inESSum, "nEventsAccepted",
192  TH1::Class()));
193  if (oEvAcc) nEvAcc = oEvAcc->GetEntries();
194  }
195  if (nEvAcc > 0 && sys > 0) {
196  Long_t minEvents = 0;
197  switch (sys) {
198  //case 1: minEvents = 1000000; break; // At least 1M
199  case 1: minEvents = 500000; break; // At least 1M
200  case 2: minEvents = 10000; break; // At least 10k
201  case 3: // Fall-through
202  case 4: minEvents = 100000; break; // At least 100k
203  }
204  if (nEvAcc < minEvents) {
205  Error("RerunELossFits", "Too few events (%ld<%ld) for %d",
206  nEvAcc, minEvents, sys);
207  return;
208  }
209  }
210  // --- Open output file --------------------------------------------
211  outFile = TFile::Open(outName, "RECREATE");
212  if (!outFile)
213  throw TString::Format("Failed to open %s", outName.Data());
214 
215  // --- Make our fitter object ------------------------------------
217  fitter->SetDoFits(true);
218  fitter->SetEnableDeltaShift(shift);
219  fitter->Init();
220  if (forceSet || !fitter->ReadParameters(inEFSum)) {
221  Printf("Forced settings");
222 
223  const TAxis* etaAxis = static_cast<TAxis*>(GetObject(inEFSum,"etaAxis"));
224  if (!etaAxis) throw new TString("Cannot proceed without eta axis");
225  fitter->SetEtaAxis(*etaAxis);
226 
227  // Set maximum energy loss to consider
228  fitter->SetMaxE(15);
229  // Set number of energy loss bins
230  fitter->SetNEbins(500);
231  // Set whether to use increasing bin sizes
232  // fitter->SetUseIncreasingBins(true);
233  // Set whether to do fit the energy distributions
234  fitter->SetDoFits(kTRUE);
235  // Set whether to make the correction object
236  fitter->SetDoMakeObject(kTRUE);
237  // Set the low cut used for energy
238  fitter->SetLowCut(0.4);
239  // Set the number of bins to subtract from maximum of distributions
240  // to get the lower bound of the fit range
241  fitter->SetFitRangeBinWidth(4);
242  // Set the maximum number of landaus to try to fit (max 5)
243  fitter->SetNParticles(5);
244  // Set the minimum number of entries in the distribution before
245  // trying to fit to the data - 10k seems the least we can do
246  fitter->SetMinEntries(10000);
247  // fitter->SetMaxChi2PerNDF(10);
248  // Enable debug
249  }
250  // Set the maximum number of landaus to try to fit (max 5)
251  Int_t maxPart = sys == 1 ? 3 : 5;
252  fitter->SetNParticles(maxPart);
253  fitter->SetDoMakeObject(true);
254  fitter->SetMinEntries(10000);
255  if (fitter->GetLowCut() < .5) fitter->SetLowCut(0.55);
256  if (flags & 0x2) fitter->SetDebug(3);
257  if (flags & 0x1)
258  fitter->SetStoreResiduals(AliFMDEnergyFitter::kResidualSquareDifference);
259  // fitter->SetRegularizationCut(1e8); // Lower by factor 3
260  // Set the number of bins to subtract from maximum of distributions
261  // to get the lower bound of the fit range
262  // fitter->SetFitRangeBinWidth(2);
263  // Skip all of FMD2 and 3
264  // fitter->SetSkips(AliFMDEnergyFitter::kFMD2|AliFMDEnergyFitter::kFMD3);
265 
266  // --- Write copy of sum collection to output --------------------
267  TCollection* outFwdSum = static_cast<TCollection*>(inFwdSum->Clone());
268  outFile->cd();
269  CleanCollection(outFwdSum);
270  TCollection* outEFSum = GetCollection(outFwdSum, "fmdEnergyFitter");
271  if (outEFSum) {
272  TObject* np = GetObject(outEFSum,"nParticles",TParameter<int>::Class());
273  if (np) (static_cast<TParameter<int>*>(np))->SetVal(maxPart);
274  }
275  outFwdSum->Write(inFwdSum->GetName(), TObject::kSingleKey);
276  outFwdSum->ls();
277 
278  // --- Now do the fits -------------------------------------------
279  fitter->Print();
280  TStopwatch timer;
281  timer.Start();
282  fitter->Fit(static_cast<TList*>(outFwdSum));
283  timer.Print();
284 
285  // --- Copy full result folder -----------------------------------
286  TCollection* outFwdRes = static_cast<TCollection*>(inFwdRes->Clone());
287  CleanCollection(outFwdRes);
288  // Remove old fits
289  while (true) {
290  TCollection* outEFRes = GetCollection(outFwdRes,"fmdEnergyFitter",false);
291  if (!outEFRes) break;
292  outFwdRes->Remove(outEFRes);
293  }
294  // outFwdRes->ls();
295  // Make our new fit results folder, and add it to results folder
296  TCollection* tmp = GetCollection(outFwdSum, "fmdEnergyFitter");
297  outEFRes = static_cast<TCollection*>(tmp->Clone());
298  outEFRes->Add(new TNamed("refitted", "Refit of the data"));
299  outFwdRes->Add(outEFRes);
300  // outFwdRes->ls();
301 
302  // --- Write out new results folder ------------------------------
303  outFile->cd();
304  outFwdRes->Write(inFwdRes->GetName(), TObject::kSingleKey);
305  Printf("Wrote results to \"%s\" (%s)", outName.Data(), outFile->GetName());
306  allOk = true;
307  }
308  catch (const TString* e) {
309  Error("RerunELossFits", e->Data());
310  }
311  catch (const TString& e) {
312  Error("RerunELossFits", e.Data());
313  }
314  if (inFile) inFile->Close();
315  if (outFile) {
316  Printf("Wrote new output to \"%s\"", outName.Data());
317  outFile->Close();
318  }
319 
320  if (allOk)
321  gROOT->Macro(Form("%s/corrs/DrawCorrELoss.C(false,false,\"%s\")",
322  fwd, outName.Data()));
323 }
324 
325 
TObject * GetObject(const TCollection *parent, const TString &name, const TClass *cls=0, Bool_t verbose=true)
void RerunELossFits(Bool_t forceSet=false, const TString &input="forward_eloss.root", Bool_t shift=true, const TString &output="", UShort_t flags=0x1)
TCanvas * c
Definition: TestFitELoss.C:172
int Int_t
Definition: External.C:63
TList * fitter
Definition: DrawAnaELoss.C:26
void CleanCollection(TCollection *c)
const char * fwd
const char * dets[]
Definition: PlotSysInfo.C:138
unsigned short UShort_t
Definition: External.C:28
TList * ef
Definition: TestFitELoss.C:145
bool Bool_t
Definition: External.C:53
TCollection * GetCollection(TDirectory *dir, const TString &name, Bool_t verbose=true)
Definition: External.C:196
void RemoveObject(TCollection *parent, const TString &name)
TDirectoryFile * dir