- Adding check and flagging for HG present
[u/mrichter/AliRoot.git] / MUON / CreateWeightedRejectList.C
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 ///
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 ///
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
82 class RunInfo
83 {
84 public:
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
93 private:
94   Int_t fNumber;
95   Long64_t fNevents;
96 };
97
98 std::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
104 AliMUONRejectList* CheckDE_BP_ManuPedestals(Int_t rejectMask, AliMUONPadStatusMaker& status);
105
106 Int_t AddEventsSingleRun(Int_t index, 
107                          Int_t run_number, 
108                          AliMUONRejectList& rejectedEvents);
109
110 std::vector<RunInfo> runs;
111
112 //______________________________________________________________________________
113 bool 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 //____________________________________________________________________________________________
228 AliMUONRejectList* 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 //____________________________________________________________________________________________
268 Int_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   
280   Long64_t nEvents = runs[index].Nevents();
281   
282   Int_t rejectMask = recoParam->PadGoodnessMask();
283   
284   // Functions to compute the number of bad channels at the DE, BP and Manu level
285   AliMUONRejectList* nBadChannels = CheckDE_BP_ManuPedestals(rejectMask, status);
286   
287   AliMpDEIterator DEiter;
288   for (DEiter.First(); !DEiter.IsDone(); DEiter.Next())
289     if (DEiter.CurrentDEId() < 1100)
290     {
291       Int_t DEid = DEiter.CurrentDEId();
292       Int_t nBusPatches = DEiter.CurrentDE()->GetNofBusPatches();
293       
294       Double_t nChannelsDE = static_cast<Double_t>(AliMpDDLStore::Instance(kFALSE)->GetDetElement(DEid)->NofChannels());
295       Double_t nBadChannelsDE = nBadChannels->DetectionElementProbability(DEid);
296       Double_t ratioBadChannelsDE = 0.0;
297       if (nChannelsDE != 0.0)
298         ratioBadChannelsDE = 100.0 * nBadChannelsDE / nChannelsDE;
299       
300       
301       Bool_t goodDE = kTRUE;
302       if (ratioBadChannelsDE >= maxBadChannels)
303       {
304         Double_t newRejectedEventsDE = rejectedEvents.DetectionElementProbability(DEid) + static_cast<Double_t>(nEvents);
305         rejectedEvents.SetDetectionElementProbability(DEid, newRejectedEventsDE);
306         goodDE = kFALSE;
307       }
308                 
309       if (goodDE)      
310         for (Int_t ii = 0; ii < nBusPatches; ii++)
311         {
312           Int_t BPid = DEiter.CurrentDE()->GetBusPatchId(ii);
313           Int_t nManus = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetNofManus();
314           
315           Double_t nChannelsBP = static_cast<Double_t>(nManus * AliMpConstants::ManuNofChannels());
316           Double_t nBadChannelsBP = nBadChannels->BusPatchProbability(BPid);
317           Double_t ratioBadChannelsBP = 100.0 * nBadChannelsBP / nChannelsBP;
318           
319           Bool_t goodBP = kTRUE;
320           
321           if (goodBP)
322             if (ratioBadChannelsBP >= maxBadChannels)
323             {
324               Double_t newRejectedEventsBP = rejectedEvents.BusPatchProbability(BPid) + static_cast<Double_t>(nEvents);
325               rejectedEvents.SetBusPatchProbability(BPid, newRejectedEventsBP);
326               goodBP = kFALSE;
327             }
328           
329           if (goodBP)
330             for (Int_t jj = 0; jj < nManus; jj++)
331             {
332               Int_t manuId = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetManuId(jj);
333               
334               Double_t nChannelsManu = DEiter.CurrentDE()->NofChannelsInManu(manuId);;
335               Double_t nBadChannelsManu = nBadChannels->ManuProbability(DEid, manuId);
336               Double_t ratioBadChannelsManu = 100.0 * nBadChannelsManu / nChannelsManu;
337               
338               Bool_t goodManu = kTRUE;
339               
340               if (ratioBadChannelsManu >= maxBadChannels)
341               {
342                 Double_t newRejectedEventsManu = rejectedEvents.ManuProbability(DEid, manuId) + static_cast<Double_t>(nEvents);
343                 rejectedEvents.SetManuProbability(DEid, manuId, newRejectedEventsManu);
344                 goodManu = kFALSE;
345               }
346               
347               
348               Int_t nChannels = DEiter.CurrentDE()->NofChannelsInManu(manuId);
349               
350               if (goodManu)
351                 for (Int_t channel = 0; channel < nChannels; channel++)
352                   if (nBadChannels->ChannelProbability(DEid, manuId, channel) > 0.0)
353                   {
354                     Double_t newRejectedEventsChannel = rejectedEvents.ChannelProbability(DEid, manuId, channel) + static_cast<Double_t>(nEvents);
355                     rejectedEvents.SetChannelProbability(DEid, manuId, channel, newRejectedEventsChannel);
356                   } // end of Channels loop
357             } // end of MANUs loop          
358         } // end of BPs loop
359     } // end of DEs loop
360   
361   delete nBadChannels;
362   
363   return nEvents;
364 }
365
366 //______________________________________________________________________________
367 //
368 //
369 // The code below might be used someday to retrieve the number of events
370 // with at least one muon directly from a runlist, instead of requiring
371 // as input both the run numbers and the number of events.
372 //
373 // This approach though is not without its problems :
374 //
375 // - even if we are using tags, as the tags are chunk-based (for the moment at
376 //   least), it's pretty slow to loop on all the chunks (it is becoming usual
377 //   to get runs with some 600 chunks or so...)
378 // - if part of a run is not reconstructed (for whatever reason) with the pass
379 //   we use to compute the number of events with at least one muon, we might
380 //   underestimate its weight in the final rejectlist for the next passes (
381 //   if for the next passes the full run is correctly reconstructed for instance)
382 // - a muon in the tag = either a trigger track or a tracker track, i.e. we don't
383 //   know for sure it's a good track...
384 //
385 // Given all those reasons, we for the moment use a dumb approach : 
386 //
387 // - go to alice-logbook.cern.ch . Select the runs you want and the trigger class
388 //   you want (most probably CMUS1B), and export the list as (run,nevents)
389 // - run this macro with this (run,nevents) list
390 //
391 //______________________________________________________________________________
392
393 #if 0
394
395 #include "TChain.h"
396 #include "AliDetectorTagCuts.h"
397 #include "AliEventTagCuts.h"
398 #include "AliLHCTagCuts.h"
399 #include "AliRunTagCuts.h"
400 #include "AliTagAnalysis.h"
401 #include "TGridResult.h"
402
403 //______________________________________________________________________________
404 bool CreateXMLCollectionFromRunList(const char* collectionName, 
405                                     const char* runlist,
406                                     const char* type,
407                                     int passNumber)
408 {
409   /// From a list of run, create a collection with only events containing at least one muon
410   /// \param collectionName : output name of the collection (without extension, which will be .xml)
411   /// \param runlist : text file containing one integer per line = one run number per line
412   /// \param type : files to consider, either ESD or AD
413   /// \param passNumber : 1 or 2 most probably (to distinguish which reco pass is used)
414   
415   if (!gGrid) TGrid::Connect("alien://");
416   if (!gGrid) 
417   {
418     return 0x0;
419   }
420   
421   TString stype(type);
422   stype.ToUpper();
423   
424   if ( stype != "ESD" && stype != "AOD" )
425   {
426     cout << "Only ESD or AOD type supported" << endl;
427     return false;
428   }
429   
430   ifstream in(gSystem->ExpandPathName(runlist));
431   Int_t runNumber;
432   Int_t ntagfiles(0);
433   
434   AliTagAnalysis tagAnalysis("ESD");
435   
436   while ( in >> runNumber ) 
437   {
438     TGridResult *res = gGrid->Query("/alice/data",
439                                     Form("%09d/%ss/pass%d/*%d*/Run%d.Event*.ESD.tag.root",
440                                          runNumber,stype.Data(),passNumber,runNumber,runNumber));
441     Int_t nFiles = res->GetEntries();
442     if (!nFiles) 
443     {
444       continue;
445     }
446     
447     for (Int_t i = 0; i < nFiles; ++i) 
448     {
449       TString filename = res->GetKey(i, "turl");
450       if(filename == "") continue;
451       tagAnalysis.AddTagsFile(filename.Data(),kFALSE);
452       ++ntagfiles;
453     }
454     delete res;
455   }
456   
457   cout << ntagfiles << " tag files added" << endl;
458   
459   AliRunTagCuts runCuts;
460   AliEventTagCuts eventCuts;
461   AliLHCTagCuts lhcCuts;
462   AliDetectorTagCuts detCuts;
463   
464   eventCuts.SetNMuonRange(1,99999);
465   
466   return tagAnalysis.CreateXMLCollection(collectionName,&runCuts,&lhcCuts,&detCuts,&eventCuts);  
467   
468 }
469
470 //______________________________________________________________________________
471 TChain* CreateChainFromXMLCollection(const char* collectionName, const char* type)
472 {
473   /// Create a chain from an XML collection file.
474   
475   if (!gGrid) TGrid::Connect("alien://");
476   if (!gGrid) 
477   {
478     return 0x0;
479   }
480   
481   TString stype(type);
482   stype.ToUpper();
483   
484   if ( stype != "ESD" && stype != "AOD" )
485   {
486     cout << "Only ESD or AOD type supported" << endl;
487     return 0x0;
488   }
489   
490   AliTagAnalysis tagAnalysis(stype.Data());
491   
492   return tagAnalysis.CreateChainFromCollection(collectionName,stype=="ESD" ? "esdTree":"aodTree");
493 }
494
495 #endif