AliPhysics  32e057f (32e057f)
ExtractAcceptance.C
Go to the documentation of this file.
1 //_____________________________________________________________________
16  Char_t r,
17  Double_t vz,
18  Int_t& nDead,
19  std::ostream* deadScript)
20 {
21  AliFMDGeometry* geom = AliFMDGeometry::Instance();
22  AliFMDParameters* pars = AliFMDParameters::Instance();
23  Float_t konst = pars->GetDACPerMIP();
24 
25  UShort_t nS = (r == 'I' || r == 'i' ? 20 : 40);
26  UShort_t nT = (r == 'I' || r == 'i' ? 512 : 256);
27 
28  // Make our two histograms
29  TH2D* hAll = new TH2D("all","All",200,-4,6,nS,0,2*TMath::Pi());
30  hAll->SetXTitle("#eta");
31  hAll->SetYTitle("#phi");
32  hAll->Sumw2();
33  hAll->SetDirectory(0);
34  TH2D* hOK = static_cast<TH2D*>(hAll->Clone());
35  hOK->SetDirectory(0);
36 
37  if (deadScript)
38  *deadScript << "\n // === FMD" << d << r << " === " << std::endl;
39  // Loop over all sectors and strips in this ring
40  Int_t nOK = 0;
41  Int_t nAll = 0;
42  Int_t nPhi = hAll->GetNbinsY();
43  for (UShort_t s = 0; s < nS; s++) {
44  for (UShort_t t = 0; t < nT; t++) {
45  // Get eta,phi by quering the geometry (first for (x,y,z), then
46  // correcting for the vertex position, and then calculating
47  // (eta, phi))
48  Double_t x, y, z;
49  geom->Detector2XYZ(d, r, s, t, x, y, z);
50  z -= vz;
51  Double_t q, eta, phi, theta;
52  AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z, q, eta, phi, theta);
53  if (phi < 0) phi += 2*TMath::Pi();
54  TString reason;
55 
56  // Check if this is a dead channel or not
57  Bool_t isDead = pars->IsDead(d, r, s, t);
58  if (isDead) {
59  if (deadScript)
60  Warning("", "FMD%d%c[%2d,%3d] is dead because OCDB says so");
61  reason = "OCDB";
62  }
63 
64  // Pedestal, noise, gain
65  Double_t ped = pars->GetPedestal (d, r, s, t);
66  Double_t noise = pars->GetPedestalWidth(d, r, s, t);
67  Double_t gain = pars->GetPulseGain (d, r, s, t);
68  if (ped < 10) {
69  if (!isDead && deadScript) {
70  reason = "Low pedestal";
71  Warning("", "FMD%d%c[%2d,%3d] is dead because pedestal %f < 10",
72  d, r, s, t, ped);
73  }
74  isDead = true;
75  }
76  Float_t corr = 0;
77  if (noise > .5 && gain > .5) corr = noise / (gain * konst);
78  if (corr > 0.05 || corr <= 0) {
79  if (!isDead && deadScript) {
80  reason = "high noise/low gain";
81  Warning("", "FMD%d%c[%2d,%3d] is dead because %f/(%f*%f)=%f > 0.05 or negative",
82  d, r, s, t, noise, gain, konst, corr);
83  }
84  isDead = true;
85  }
86 
87 
88 
89 
90  // Special check for FMD2i - upper part of sectors 16/17 have
91  // have anomalous gains and/or noise - common sources are
92  // power regulartors for bias currents and the like
93 #if 0
94  Int_t VA = t/128;
95  if(d==2 && r=='I' && VA>1 && (s==16 || s==17)) {
96  if (!isDead && deadScript) {
97  reason = "I say so";
98  Warning("", "FMD%d%c[%2d,%3d] is dead because I say so", d, r, s, t);
99  }
100  isDead =true;
101  }
102 #endif
103 
104  // Find the eta bin number and corresponding overflow bin
105  Int_t etaBin = hAll->GetXaxis()->FindBin(eta);
106  Int_t ovrBin = hAll->GetBin(etaBin, nPhi+1);
107 
108  // Increment all histogram
109  hAll->Fill(eta, phi);
110  hAll->AddBinContent(ovrBin);
111  nAll++;
112 
113  // If not dead, increment OK histogram
114  if (!isDead) {
115  hOK->Fill(eta, phi);
116  hOK->AddBinContent(ovrBin);
117  nOK++;
118  }
119  else {
120  nDead++;
121  if (deadScript)
122  *deadScript << " filter->AddDead(" << d << ",'" << r << "',"
123  << std::setw(2) << s << ','
124  << std::setw(3) << t << "); // "
125  << reason << std::endl;
126  }
127  }
128  }
129  // Divide out the efficiency.
130  // Note, that the overflow bins along eta now contains the ratio
131  // nOK/nAll Strips for a given eta bin.
132  hOK->Divide(hOK,hAll,1,1,"B");
133 
134  // Invert overflow bin
135  for (Int_t etaBin = 1; etaBin <= hOK->GetNbinsX(); etaBin++) {
136  Double_t ovr = hOK->GetBinContent(etaBin, nPhi+1);
137  Double_t novr = (ovr < 1e-12 ? 0 : 1./ovr);
138  hOK->SetBinContent(etaBin, nPhi+1, novr);
139 #if 0
140  if (ovr > 0 && ovr != 1)
141  Info("", "Setting overflow bin (%3d,%3d) to 1/%f=%f", etaBin, nPhi+1,
142  ovr, hOK->GetBinContent(etaBin, nPhi+1));
143 #endif
144  }
145  // Clean up
146  delete hAll;
147 
148  Printf("=== FMD%d%c at vz=%+5.1f - %d/%d=%3d%% OK (w/overflow)",
149  d, r, vz, nOK, nAll, (100*nOK)/nAll);
150 
151  // Return result
152  return hOK;
153 }
154 
155 //_____________________________________________________________________
166 void ExtractAcceptance(Int_t runNo=121526,
167  Int_t nVtxBins=10,
168  Float_t vtxLow=-10,
169  Float_t vtxHigh=10)
170 {
171  const char* fwd = "$ALICE_PHYSICS/PWGLF/FORWARD/analysis2";
172  gSystem->AddIncludePath(Form("-I%s", fwd));
173  gROOT->Macro(Form("%s/scripts/LoadLibs.C", fwd));
174 
175  // gSystem->Load("libANALYSIS");
176  // gSystem->Load("libANALYSISalice");
177  // gSystem->Load("libPWGLFforward2");
178  // Float_t delta = (vtxHigh - vtxLow) / (Float_t)nVtxBins;
179 
180  Bool_t gridOnline = kTRUE;
181  if(!(TGrid::Connect("alien://",0,0,"t")))
182  gridOnline = kFALSE;
183 
184  // --- Figure out the year --------------------------------------
185  UShort_t year = 0;
186  if (runNo <= 99999) year = 2009;
187  else if (runNo <= 139667) year = 2010;
188  else if (runNo <= 170718) year = 2011;
189  else if (runNo <= 194306) year = 2012;
190  else if (runNo <= 197709) year = 2013;
191  else if (runNo <= 208364) year = 2014;
192  else year = 2015;
193  if (year <= 0) {
194  Warning("", "Couldn't deduce the year from the run number");
195  // return;
196  }
197 
198  // --- Initialisations ------------------------------------------
199  //Set up CDB manager
200  Printf("=== Setting up OCDB");
201 
202  AliCDBManager* cdb = AliCDBManager::Instance();
203  if(gridOnline) {
204  cdb->SetDefaultStorageFromRun(runNo);
205  // cdb->SetDefaultStorage(Form("alien://Folder=/alice/data/%4d/OCDB", year));
206  }
207  else
208  cdb->SetDefaultStorage("local://$(ALICE_PHYSICS)/OCDB");
209  cdb->SetRun(runNo);
210 
211  // Get the geometry
212  Printf("=== Loading geometry");
213  AliGeomManager::LoadGeometry();
214 
215  // Get an initialize parameters
216  Printf("=== Intialising parameters");
217  AliFMDParameters* pars = AliFMDParameters::Instance();
218  pars->Init();
219 
220  // Get an initialise geometry
221  Printf("=== Initialising geomtry");
222  AliFMDGeometry* geom = AliFMDGeometry::Instance();
223  geom->Init();
224  geom->InitTransformations();
225 
226  // --- Get the general run parameters ------------------------------
227  AliCDBEntry* grpE = cdb->Get("GRP/GRP/Data");
228  if (!grpE) {
229  AliWarningF("No GRP entry found for run %d", runNo);
230  return;
231  }
232  AliGRPObject* grp = static_cast<AliGRPObject*>(grpE->GetObject());
233  if (!grp) {
234  AliWarningF("No GRP object found for run %d", runNo);
235  return;
236  }
237  Float_t beamE = grp->GetBeamEnergy();
238  TString beamT = grp->GetBeamType();
239 # if 0
240  // This isn't really needed as the acceptance map is indifferent to
241  // the field settings.
242  Float_t l3cur = grp->GetL3Current(AliGRPObject::kMean);
243  Char_t l3pol = grp->GetL3Polarity();
244  Bool_t l3lhc = grp->IsPolarityConventionLHC();
245  Bool_t l3uni = grp->IsUniformBMap();
246  AliMagF* fldM =
247  AliMagF::CreateFieldMap(TMath::Abs(l3cur) * (l3pol ? -1:1), 0,
248  (l3lhc ? 0 : 1), l3uni, beamE, beamT.Data());
249  Float_t l3fld = fldM->SolenoidField();
250 #endif
251 
254  Short_t fld = 0; // AliForwardUtil::ParseMagneticField(l3fld);
255  Printf("=== Run=%d, year=%d, sys=%d, sNN=%d, fld=%d",
256  runNo, year, sys, sNN, fld);
257 
258  // --- Output object -----------------------------------------------
259  // Make our correction object
261  corr->SetVertexAxis(nVtxBins, vtxLow, vtxHigh);
262 
263  // --- Output script -----------------------------------------------
264  std::ofstream deadScript("deadstrips.C");
265  deadScript << "// Automatically generaeted by ExtractAcceptance.C\n"
266  << "// Add additional dead strips to sharing filter\n"
267  << "// Information taken from OCDB entry for\n"
268  << "//\n"
269  << "// run = " << runNo << "\n"
270  << "// year = " << year << "\n"
271  << "// system = " << sys << "\n"
272  << "// sqrt{sNN} = " << sNN << "GeV\n"
273  << "// L3 field = " << fld << "kG\n"
274  << "void deadstrips(AliFMDESDFixer* filter)\n"
275  << "{" << std::endl;
276 
277  // --- Loop over verticies and rings -------------------------------
278  Int_t nDead = 0;
279  Bool_t gotDead = false;
280  Float_t dV = (vtxHigh - vtxLow) / nVtxBins;
281  Printf("=== Looping over vertices: %d bins from %+6.2f to %+6.2f "
282  "in steps of %+6.2f",
283  nVtxBins, vtxLow, vtxHigh, dV);
284  for (Double_t v = vtxLow+dV/2; v < vtxHigh; v += dV) {
285  for(UShort_t d = 1; d <= 3;d++) {
286  UShort_t nR = (d == 1 ? 1 : 2);
287  for (UShort_t q = 0; q < nR; q++) {
288  Char_t r = (q == 0 ? 'I' : 'O');
289 
290  // Delegate to other function
291  Int_t nLocal = 0;
292  TH2D* ratio = MakeOneRing(d, r, v, nLocal,
293  !gotDead ? &deadScript : 0);
294  nDead += nLocal;
295  if (!ratio) {
296  Warning("ExtractAcceptance", "Didn't get correction from "
297  "FMD%d%c @ vz=%+6.2fcm", d, r, v);
298  continue;
299  }
300  // Printf("v=%+6.2f FMD%d%c, got %d dead strips", v, d, r, nLocal);
301 
302  // Set the correction
303  corr->SetCorrection(d, r, v, ratio);
304  }
305  }
306  gotDead = true;
307  }
308  corr->SetHasOverflow();
309  corr->Print();
310  // corr->ls();
311 
312  deadScript << "}\n"
313  << "// EOF" << std::endl;
314  deadScript.close();
315 
316  // Write to a file
317  Printf("=== Writing to disk");
319  if (!cm.Store(corr, runNo, sys, sNN, fld, false, false,
320  "fmd_corrections.root")) {
321  Error("", "Failed to store acceptance correction in local file");
322  return;
323  }
324 
325  std::ofstream f("Upload.C");
326  f << "// Generated by ExtractELoss.C\n"
327  << "TString MakeDest(const TString& dest, const TString& fname)\n"
328  << "{\n"
329  << " TString tmp(dest);\n"
330  << " if (!tmp.IsNull()) {\n"
331  << " if (!tmp.EndsWith(\"/\")) tmp.Append(\"/\");\n"
332  << " tmp.Append(fname);\n"
333  << " }\n"
334  << " return tmp;\n"
335  << "}\n\n"
336  << "void Upload(const TString& dest=\"\")\n"
337  << "{\n"
338  << " gROOT->Macro(\"" << fwd << "/scripts/LoadLibs.C\");\n"
339  << " \n"
340  << " const char* fmdFile = \"fmd_corrections.root\";\n"
341  << " TString fdest = MakeDest(dest, fmdFile);\n"
342  << " \n"
343  << " AliForwardCorrectionManager::Instance().Append(fmdFile, fdest);\n"
344  << "}\n"
345  << "// EOF\n"
346  << std::endl;
347  f.close();
348 
349  Info("ExtracAcceptance",
350  "Run generated Upload.C(DEST) script to copy files in place");
351 }
352 
353 //
354 // EOF
355 //
356 
357 
double Double_t
Definition: External.C:58
TH2D * MakeOneRing(UShort_t d, Char_t r, Double_t vz, Int_t &nDead, std::ostream *deadScript)
static UShort_t ParseCenterOfMassEnergy(UShort_t sys, Float_t cms)
TSystem * gSystem
char Char_t
Definition: External.C:18
void SetVertexAxis(const TAxis &axis)
int Int_t
Definition: External.C:63
virtual Bool_t Store(TObject *o, ULong_t runNo, UShort_t sys, UShort_t sNN, Short_t field, Bool_t mc, Bool_t sat, const char *file, const char *meth="NEAR") const
float Float_t
Definition: External.C:68
Definition: External.C:228
void SetHasOverflow(Bool_t present=true)
void ExtractAcceptance(Int_t runNo=121526, Int_t nVtxBins=10, Float_t vtxLow=-10, Float_t vtxHigh=10)
short Short_t
Definition: External.C:23
const char * fwd
static UShort_t ParseCollisionSystem(const char *sys)
GRPData * grp
Definition: GRP.C:361
void Print(Option_t *option="R") const
unsigned short UShort_t
Definition: External.C:28
bool Bool_t
Definition: External.C:53
Bool_t SetCorrection(UShort_t d, Char_t r, Double_t v, TH2D *h)
static AliForwardCorrectionManager & Instance()