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