minor coverity defect: added protection for self-assignment
[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   status.SetLimits(*recoParam);
281   
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 //______________________________________________________________________________
406 bool 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   
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 //______________________________________________________________________________
473 TChain* 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