Fixes for #86059: Install data when ALICE_ROOT does not point to source (Christian)
[u/mrichter/AliRoot.git] / MUON / CreateWeightedRejectList.C
CommitLineData
f0938a5b 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16// $Id$
17
18///
19/// \ingroup macros
20///
21/// \file CreateWeightedRejectList.C
22///
23/// \brief Macro test for the creation of a RejectList object (for simulations) in the OCDB, for the MUON Tracke
24///
25/// Usage:
26///
27/// root[0] .L CreateWeightedRejectList.C+
28/// root[1] CreateWeightedRejectList("runlist.txt","local://$HOME/myLocalOCDB");
29///
30/// where runlist.txt has 2 integers per line = "run nevents"
31/// where nevents is the number of events where there's (for instance) a CMUS1B trigger
32///
951dd287 33///
34/// Assuming the file coming from the Export feature of the logbook is runlist.logbook.txt =>
35///
36/// 115521;PHYSICS_1;1357;0.9
37/// 115516;PHYSICS_1;944;0.9
38///
39/// The awk command below will output the needed format for this macro :
40///
41/// awk '{split ($0,a,";"); print a[1] " " a[3];}' runlist.logbook.txt =>
42///
43/// 115521 1357
44/// 115516 944
45///
f0938a5b 46/// \author Matthieu Lenhardt, Subatech
47///
48
49#include "AliCDBEntry.h"
50#include "AliCDBManager.h"
51#include "AliDAQ.h"
52#include "AliDCSValue.h"
53#include "AliMUONCDB.h"
54#include "AliMUONCalibParamNI.h"
55#include "AliMUONCalibrationData.h"
56#include "AliMUONPadStatusMaker.h"
57#include "AliMUONPadStatusMapMaker.h"
58#include "AliMUONRecoParam.h"
59#include "AliMUONRejectList.h"
60#include "AliMUONRejectList.h"
61#include "AliMUONTrackerData.h"
62#include "AliMUONVCalibParam.h"
63#include "AliMpBusPatch.h"
64#include "AliMpCDB.h"
65#include "AliMpConstants.h"
66#include "AliMpDCSNamer.h"
67#include "AliMpDDLStore.h"
68#include "AliMpDEIterator.h"
69#include "AliMpDEManager.h"
70#include "AliMpDetElement.h"
71#include "AliSysInfo.h"
72#include "TEnv.h"
73#include <Riostream.h>
74#include <TFile.h>
75#include <TGrid.h>
76#include <TMap.h>
77#include <TObjArray.h>
78#include <TObject.h>
79#include <TTree.h>
80#include <vector>
81
82class RunInfo
83{
84public:
85 RunInfo(Int_t run, Long64_t nevents) : fNumber(run), fNevents(nevents) { }
86
87 Int_t Number() const { return fNumber; }
88
89 Long64_t Nevents() const { return fNevents; }
90
91 void SetNevents(Long64_t n) { fNevents = n; }
92
93private:
94 Int_t fNumber;
95 Long64_t fNevents;
96};
97
98std::ostream& operator<<(std::ostream& out, const RunInfo& run)
99{
100 out << Form("RUN %09d NEVENTS %10d",run.Number(),run.Nevents());
101 return out;
102}
103
104AliMUONRejectList* CheckDE_BP_ManuPedestals(Int_t rejectMask, AliMUONPadStatusMaker& status);
105
106Int_t AddEventsSingleRun(Int_t index,
107 Int_t run_number,
108 AliMUONRejectList& rejectedEvents);
109
110std::vector<RunInfo> runs;
111
112//______________________________________________________________________________
113bool CreateWeightedRejectList(const char* runlistfile="runlist.txt",
114 const char* rejectListPath="local://$HOME/OCDB")
115{
116 /// Create a weighted RejectList for the runs included in the run list
117 /// The cuts are the same that the one applied in the RecoParam used to create the ESDs.
118 ///
119 /// \param runlistfile : an ASCII file where each line is a pair "run nevents"
120 /// \param rejectListPath : set this to the folder where the RejectList will be created (must have right access to it, of course)
121
122 TGrid::Connect("alien://");
123
124 ifstream in(runlistfile);
125 int run, nevents;
126
127 while ( in >> run >> nevents )
128 {
129 runs.push_back(RunInfo(run,nevents));
130 };
131
132 for ( std::vector<RunInfo>::size_type i = 0; i < runs.size(); ++i )
133 {
134 cout << runs[i] << endl;
135 }
136
137 if (runs.empty()) return false;
138
139 Int_t firstRun = runs[0].Number();
140 Int_t lastRun = runs[runs.size()-1].Number();
141
142 const char* rawOCDB = "raw://";
143
144 AliCDBManager* manager = AliCDBManager::Instance();
145
146 manager->SetCacheFlag(kFALSE);
147
148 manager->SetDefaultStorage(rawOCDB);
149
150 manager->SetSpecificStorage("MUON/Calib/RejectList", rejectListPath);
151
152 AliMUONRejectList weightedRejectList;
153 AliMUONRejectList rejectedEvents;
154
155 Long64_t nEventsTotal = 0;
156
157 for (std::vector<RunInfo>::size_type ii = 0; ii < runs.size(); ++ii)
158 {
159 Int_t runNumber = runs[ii].Number();
160 AliSysInfo::Instance()->AddStamp(Form("RUN%d",runNumber));
161 nEventsTotal += AddEventsSingleRun(ii, runNumber, rejectedEvents);
162 }
163
164 AliMpDEIterator DEiter;
165
166 for (DEiter.First(); !DEiter.IsDone(); DEiter.Next())
167 if (DEiter.CurrentDEId() < 1100)
168 {
169 Int_t DEid = DEiter.CurrentDEId();
170
171 Double_t nBadEventsDE = rejectedEvents.DetectionElementProbability(DEid);
172 Double_t nEventsTotalDE = static_cast<Double_t>(nEventsTotal);
173 Float_t probaDE = 0.0;
174 if (nEventsTotalDE != 0.0)
175 probaDE = nBadEventsDE / nEventsTotalDE;
176 weightedRejectList.SetDetectionElementProbability(DEid, probaDE);
177
178 Int_t nBusPatches = DEiter.CurrentDE()->GetNofBusPatches();
179 for (Int_t ii = 0; ii < nBusPatches; ii++)
180 {
181 Int_t BPid = DEiter.CurrentDE()->GetBusPatchId(ii);
182
183 Double_t nBadEventsBP = rejectedEvents.BusPatchProbability(BPid);
184 Double_t nEventsTotalBP = nEventsTotalDE - nBadEventsDE;
185 Float_t probaBP = 0.0;
186 if (nEventsTotalBP != 0.0)
187 probaBP = nBadEventsBP / nEventsTotalBP;
188 weightedRejectList.SetBusPatchProbability(BPid, probaBP);
189
190 Int_t nManus = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetNofManus();
191 for (Int_t jj = 0; jj < nManus; jj++)
192 {
193 Int_t manuId = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetManuId(jj);
194
195 Double_t nBadEventsManu = rejectedEvents.ManuProbability(DEid, manuId);
196 Double_t nEventsTotalManu = nEventsTotalBP - nBadEventsBP;
197 Float_t probaManu = 0.0;
198 if (nEventsTotalManu != 0.0)
199 probaManu = nBadEventsManu / nEventsTotalManu;
200 weightedRejectList.SetManuProbability(DEid, manuId, probaManu);
201
202 for (Int_t channel = 0; channel < AliMpConstants::ManuNofChannels(); channel++)
203 {
204 Double_t nBadEventsChannel = rejectedEvents.ChannelProbability(DEid, manuId, channel);
205 Double_t nEventsTotalChannel = nEventsTotalManu - nBadEventsManu;
206 Float_t probaChannel = 0.0;
207 if (nEventsTotalChannel != 0.0)
208 {
209 probaChannel = nBadEventsChannel / nEventsTotalChannel;
210 }
211
212 weightedRejectList.SetChannelProbability(DEid, manuId, channel, probaChannel);
213 }
214 }
215 }
216 }
217
218
219 AliMUONCDB::WriteToCDB(&weightedRejectList, "MUON/Calib/RejectList",
220 firstRun , lastRun,
221 "Weighted reject List for MCH, for simulations only",
222 "Matthieu Lenhardt");
223
224 return true;
225}
226
227//____________________________________________________________________________________________
228AliMUONRejectList* CheckDE_BP_ManuPedestals(Int_t rejectMask, AliMUONPadStatusMaker& status)
229{
230 AliMUONRejectList* nBadChannels = new AliMUONRejectList;
231
232 AliMpDEIterator DEiter;
233 for (DEiter.First(); !DEiter.IsDone(); DEiter.Next())
234 if (DEiter.CurrentDEId() < 1100)
235 {
236 Int_t DEid = DEiter.CurrentDEId();
237 Int_t nBusPatches = DEiter.CurrentDE()->GetNofBusPatches();
238 for (Int_t ii = 0; ii < nBusPatches; ii++)
239 {
240 Int_t BPid = DEiter.CurrentDE()->GetBusPatchId(ii);
241
242 Int_t nManus = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetNofManus();
243 for (Int_t jj = 0; jj < nManus; jj++)
244 {
245 Int_t manuId = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetManuId(jj);
246
247 Int_t nChannels = DEiter.CurrentDE()->NofChannelsInManu(manuId);
248 for (Int_t channel = 0; channel < nChannels; channel++)
249 {
250 Int_t padStatus = status.PadStatus(DEid, manuId, channel);
251 if ((rejectMask == 0 && padStatus != 0) || (padStatus & rejectMask) != 0)
252 {
253 nBadChannels->SetDetectionElementProbability(DEid, nBadChannels->DetectionElementProbability(DEid) + 1.0);
254 nBadChannels->SetBusPatchProbability(BPid, nBadChannels->BusPatchProbability(BPid) + 1.0);
255 nBadChannels->SetManuProbability(DEid, manuId, nBadChannels->ManuProbability(DEid, manuId) + 1.0);
256 nBadChannels->SetChannelProbability(DEid, manuId, channel, 1.0);
257 }
258 }
259 }
260 }
261 }
262
263 return nBadChannels;
264}
265
266
267//____________________________________________________________________________________________
268Int_t AddEventsSingleRun(Int_t index, Int_t runNumber, AliMUONRejectList& rejectedEvents)
269{
270 AliCDBManager::Instance()->SetRun(runNumber);
271
272 AliMpCDB::LoadAll();
273
274 Double_t maxBadChannels = 80.0;
275
276 AliMUONRecoParam *recoParam = AliMUONCDB::LoadRecoParam();
277 AliMUONCalibrationData calibrationData(runNumber);
278 AliMUONPadStatusMaker status(calibrationData);
279
beffc4d5 280 status.SetLimits(*recoParam);
281
f0938a5b 282 Long64_t nEvents = runs[index].Nevents();
283
284 Int_t rejectMask = recoParam->PadGoodnessMask();
285
286 // Functions to compute the number of bad channels at the DE, BP and Manu level
287 AliMUONRejectList* nBadChannels = CheckDE_BP_ManuPedestals(rejectMask, status);
288
289 AliMpDEIterator DEiter;
290 for (DEiter.First(); !DEiter.IsDone(); DEiter.Next())
291 if (DEiter.CurrentDEId() < 1100)
292 {
293 Int_t DEid = DEiter.CurrentDEId();
294 Int_t nBusPatches = DEiter.CurrentDE()->GetNofBusPatches();
295
296 Double_t nChannelsDE = static_cast<Double_t>(AliMpDDLStore::Instance(kFALSE)->GetDetElement(DEid)->NofChannels());
297 Double_t nBadChannelsDE = nBadChannels->DetectionElementProbability(DEid);
298 Double_t ratioBadChannelsDE = 0.0;
299 if (nChannelsDE != 0.0)
300 ratioBadChannelsDE = 100.0 * nBadChannelsDE / nChannelsDE;
301
302
303 Bool_t goodDE = kTRUE;
304 if (ratioBadChannelsDE >= maxBadChannels)
305 {
306 Double_t newRejectedEventsDE = rejectedEvents.DetectionElementProbability(DEid) + static_cast<Double_t>(nEvents);
307 rejectedEvents.SetDetectionElementProbability(DEid, newRejectedEventsDE);
308 goodDE = kFALSE;
309 }
310
311 if (goodDE)
312 for (Int_t ii = 0; ii < nBusPatches; ii++)
313 {
314 Int_t BPid = DEiter.CurrentDE()->GetBusPatchId(ii);
315 Int_t nManus = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetNofManus();
316
317 Double_t nChannelsBP = static_cast<Double_t>(nManus * AliMpConstants::ManuNofChannels());
318 Double_t nBadChannelsBP = nBadChannels->BusPatchProbability(BPid);
319 Double_t ratioBadChannelsBP = 100.0 * nBadChannelsBP / nChannelsBP;
320
321 Bool_t goodBP = kTRUE;
322
323 if (goodBP)
324 if (ratioBadChannelsBP >= maxBadChannels)
325 {
326 Double_t newRejectedEventsBP = rejectedEvents.BusPatchProbability(BPid) + static_cast<Double_t>(nEvents);
327 rejectedEvents.SetBusPatchProbability(BPid, newRejectedEventsBP);
328 goodBP = kFALSE;
329 }
330
331 if (goodBP)
332 for (Int_t jj = 0; jj < nManus; jj++)
333 {
334 Int_t manuId = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetManuId(jj);
335
336 Double_t nChannelsManu = DEiter.CurrentDE()->NofChannelsInManu(manuId);;
337 Double_t nBadChannelsManu = nBadChannels->ManuProbability(DEid, manuId);
338 Double_t ratioBadChannelsManu = 100.0 * nBadChannelsManu / nChannelsManu;
339
340 Bool_t goodManu = kTRUE;
341
342 if (ratioBadChannelsManu >= maxBadChannels)
343 {
344 Double_t newRejectedEventsManu = rejectedEvents.ManuProbability(DEid, manuId) + static_cast<Double_t>(nEvents);
345 rejectedEvents.SetManuProbability(DEid, manuId, newRejectedEventsManu);
346 goodManu = kFALSE;
347 }
348
349
350 Int_t nChannels = DEiter.CurrentDE()->NofChannelsInManu(manuId);
351
352 if (goodManu)
353 for (Int_t channel = 0; channel < nChannels; channel++)
354 if (nBadChannels->ChannelProbability(DEid, manuId, channel) > 0.0)
355 {
356 Double_t newRejectedEventsChannel = rejectedEvents.ChannelProbability(DEid, manuId, channel) + static_cast<Double_t>(nEvents);
357 rejectedEvents.SetChannelProbability(DEid, manuId, channel, newRejectedEventsChannel);
358 } // end of Channels loop
359 } // end of MANUs loop
360 } // end of BPs loop
361 } // end of DEs loop
362
363 delete nBadChannels;
364
365 return nEvents;
366}
367
368//______________________________________________________________________________
369//
370//
371// The code below might be used someday to retrieve the number of events
372// with at least one muon directly from a runlist, instead of requiring
373// as input both the run numbers and the number of events.
374//
375// This approach though is not without its problems :
376//
377// - even if we are using tags, as the tags are chunk-based (for the moment at
378// least), it's pretty slow to loop on all the chunks (it is becoming usual
379// to get runs with some 600 chunks or so...)
380// - if part of a run is not reconstructed (for whatever reason) with the pass
381// we use to compute the number of events with at least one muon, we might
382// underestimate its weight in the final rejectlist for the next passes (
383// if for the next passes the full run is correctly reconstructed for instance)
384// - a muon in the tag = either a trigger track or a tracker track, i.e. we don't
385// know for sure it's a good track...
386//
387// Given all those reasons, we for the moment use a dumb approach :
388//
389// - go to alice-logbook.cern.ch . Select the runs you want and the trigger class
390// you want (most probably CMUS1B), and export the list as (run,nevents)
391// - run this macro with this (run,nevents) list
392//
393//______________________________________________________________________________
394
395#if 0
396
397#include "TChain.h"
398#include "AliDetectorTagCuts.h"
399#include "AliEventTagCuts.h"
400#include "AliLHCTagCuts.h"
401#include "AliRunTagCuts.h"
402#include "AliTagAnalysis.h"
403#include "TGridResult.h"
404
405//______________________________________________________________________________
406bool CreateXMLCollectionFromRunList(const char* collectionName,
407 const char* runlist,
408 const char* type,
409 int passNumber)
410{
411 /// From a list of run, create a collection with only events containing at least one muon
412 /// \param collectionName : output name of the collection (without extension, which will be .xml)
413 /// \param runlist : text file containing one integer per line = one run number per line
414 /// \param type : files to consider, either ESD or AD
415 /// \param passNumber : 1 or 2 most probably (to distinguish which reco pass is used)
416
417 if (!gGrid) TGrid::Connect("alien://");
418 if (!gGrid)
419 {
420 return 0x0;
421 }
422
f0938a5b 423 TString stype(type);
424 stype.ToUpper();
425
426 if ( stype != "ESD" && stype != "AOD" )
427 {
428 cout << "Only ESD or AOD type supported" << endl;
429 return false;
430 }
431
432 ifstream in(gSystem->ExpandPathName(runlist));
433 Int_t runNumber;
434 Int_t ntagfiles(0);
435
436 AliTagAnalysis tagAnalysis("ESD");
437
438 while ( in >> runNumber )
439 {
440 TGridResult *res = gGrid->Query("/alice/data",
441 Form("%09d/%ss/pass%d/*%d*/Run%d.Event*.ESD.tag.root",
442 runNumber,stype.Data(),passNumber,runNumber,runNumber));
443 Int_t nFiles = res->GetEntries();
444 if (!nFiles)
445 {
446 continue;
447 }
448
449 for (Int_t i = 0; i < nFiles; ++i)
450 {
451 TString filename = res->GetKey(i, "turl");
452 if(filename == "") continue;
453 tagAnalysis.AddTagsFile(filename.Data(),kFALSE);
454 ++ntagfiles;
455 }
456 delete res;
457 }
458
459 cout << ntagfiles << " tag files added" << endl;
460
461 AliRunTagCuts runCuts;
462 AliEventTagCuts eventCuts;
463 AliLHCTagCuts lhcCuts;
464 AliDetectorTagCuts detCuts;
465
466 eventCuts.SetNMuonRange(1,99999);
467
468 return tagAnalysis.CreateXMLCollection(collectionName,&runCuts,&lhcCuts,&detCuts,&eventCuts);
469
470}
471
472//______________________________________________________________________________
473TChain* CreateChainFromXMLCollection(const char* collectionName, const char* type)
474{
475 /// Create a chain from an XML collection file.
476
477 if (!gGrid) TGrid::Connect("alien://");
478 if (!gGrid)
479 {
480 return 0x0;
481 }
482
483 TString stype(type);
484 stype.ToUpper();
485
486 if ( stype != "ESD" && stype != "AOD" )
487 {
488 cout << "Only ESD or AOD type supported" << endl;
489 return 0x0;
490 }
491
492 AliTagAnalysis tagAnalysis(stype.Data());
493
494 return tagAnalysis.CreateChainFromCollection(collectionName,stype=="ESD" ? "esdTree":"aodTree");
495}
496
497#endif