Write only detector coefficients from HLT (Raphaelle)
[u/mrichter/AliRoot.git] / MUON / AliMUONPadStatusMapMaker.cxx
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 /// \class AliMUONPadStatusMapMaker
20 /// 
21 /// Convert a pad statuses into pad status maps.
22 /// 
23 /// A pad status is one 32-bits word describing whether this pad pedestal, gains
24 /// hv is correct or not.
25 ///
26 /// A pad status *map* is one 32-bits (of which 24 only are used)
27 /// word describing whether this pad neighbours are ok or not
28 /// (whether a pad is ok or not is determined by applying a given
29 /// bitmask to the pad status word). Each bit in this word is related to one
30 /// neighbour, assuming the pad itself is at bit 0
31 ///
32 /// ----------------
33 /// |  3 |  5 |  8 |
34 /// ----------------
35 /// |  2 |  0 |  7 |
36 /// ----------------
37 /// |  1 |  4 |  6 |
38 /// ----------------
39 ///
40 /// Note that for instance in NonBending plane of slats, at the boundaries
41 /// between two pad densities, the pictures is a bit different, e.g.
42 /// (bits in () are always zero)
43 ///
44 /// so some care must be taken when designing a mask to be tested ;-) if you
45 /// want to go farther than immediate neighbours...
46 ///
47 /// If a pad is at a physical boundary, is will for sure have some bits at 1
48 /// (i.e. a non-existing neighbour is considered = bad).
49 ///
50 ///
51 /// add something about the reject list/probabilities here... (LA)
52 ///
53 /// \author Laurent Aphecetche
54 //-----------------------------------------------------------------------------
55
56 #include "AliMUONPadStatusMapMaker.h"
57
58 #include "AliCodeTimer.h"
59 #include "AliLog.h"
60 #include "AliMpDDLStore.h"
61 #include "AliMpDetElement.h"
62 #include "AliMpManuIterator.h"
63 #include "AliMUON2DMap.h"
64 #include "AliMUONCalibParamNF.h"
65 #include "AliMUONCalibParamNI.h"
66 #include "AliMUONCalibrationData.h"
67 #include "AliMUONPadStatusMaker.h"
68 #include "AliMUONRejectList.h"
69 #include "AliMUONVCalibParam.h"
70 #include "AliMUONVStore.h"
71 #include "AliMpConstants.h"
72 #include <Riostream.h>
73 #include <TList.h>
74 #include "TRandom.h"
75 #include <cassert>
76
77 /// \cond CLASSIMP
78 ClassImp(AliMUONPadStatusMapMaker)
79 /// \endcond
80
81 Int_t AliMUONPadStatusMapMaker::fgkSelfDead = 1;
82
83 //_____________________________________________________________________________
84 AliMUONPadStatusMapMaker::AliMUONPadStatusMapMaker(const AliMUONPadStatusMaker& padStatusMaker,
85                                                    Int_t mask,
86                                                    Bool_t deferredInitialization) 
87 : TObject(),
88 fkStatusMaker(padStatusMaker),
89 fMask(mask),
90 fStatusMap(new AliMUON2DMap(true)),
91 fRejectProbabilities(new AliMUON2DMap(true)),
92 fRejectList(0x0),
93 fComputeOnDemand(deferredInitialization)
94 {
95   /// ctor
96   if (!deferredInitialization)
97   {
98     AliCodeTimerAuto("Computing complete status map at once",0);
99     AliMUONVStore* neighboursStore = padStatusMaker.NeighboursStore();
100     AliMUONVCalibParam* param;
101     TIter next(neighboursStore->CreateIterator());
102     while ( ( param = static_cast<AliMUONVCalibParam*>(next()) ) )
103     {
104       Int_t detElemId = param->ID0();
105       Int_t manuId = param->ID1();
106       ComputeStatusMap(detElemId,manuId);
107     }
108   }
109   
110   /// Whatever the deferred flag is, we *have* to compute the reject 
111   /// probabilities here and now, for *all* channels.
112   
113   AliMUONRejectList* rl = padStatusMaker.CalibrationData().RejectList();
114   
115   if (rl)
116   {
117     AliMpManuIterator it;
118     Int_t detElemId;
119     Int_t manuId;
120     
121     while ( it.Next(detElemId,manuId) )
122     {
123       AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(detElemId);
124       Int_t busPatchId = AliMpDDLStore::Instance()->GetBusPatchId(detElemId,manuId);
125       
126       AliMUONVCalibParam* param = new AliMUONCalibParamNF(1,AliMpConstants::ManuNofChannels(),detElemId,manuId,0);
127       
128       Int_t n(0);
129       
130       for ( Int_t i = 0; i < AliMpConstants::ManuNofChannels(); ++i ) 
131       {
132         Float_t proba(0.0);
133         
134         if ( de->IsConnectedChannel(manuId,i) )
135         {
136           proba = TMath::Max(rl->DetectionElementProbability(detElemId),rl->BusPatchProbability(busPatchId));
137           
138           proba = TMath::Max(proba,rl->ManuProbability(detElemId,manuId));
139           
140           proba = TMath::Max(proba,rl->ChannelProbability(detElemId,manuId,i));
141           
142           if ( proba > 0 ) 
143           {
144             ++n;
145             param->SetValueAsFloat(i,0,proba);
146           }
147         }
148       }
149       
150       if ( n > 0 ) 
151       {
152         fRejectProbabilities->Add(param);
153       }
154       else
155       {
156         // no need to add empty stuff...
157         delete param;
158       }
159     }
160   
161     if ( rl->IsBinary())
162     {
163       fRejectList = fRejectProbabilities;
164       fRejectProbabilities = 0x0;
165       AliDebug(1,"RejectList = RejectProbabilities");
166       StdoutToAliDebug(1,fRejectList->Print("","MEAN"));
167     }
168     else
169     {
170       AliWarning("Will run with non trivial survival probabilities for channels, manus, etc... Better check this is a simulation and not real data !");
171       fRejectList = new AliMUON2DMap(true);
172     }
173   }
174   else
175   {
176     fRejectList = fRejectProbabilities;
177     fRejectProbabilities = 0x0;
178     AliInfo("No RejectList found, so no RejectList will be used.");
179   }
180 }
181
182 //_____________________________________________________________________________
183 AliMUONPadStatusMapMaker::~AliMUONPadStatusMapMaker()
184 {
185   /// dtor
186   delete fStatusMap;
187   delete fRejectProbabilities;
188   delete fRejectList;
189 }
190
191 //_____________________________________________________________________________
192 AliMUONVCalibParam*
193 AliMUONPadStatusMapMaker::ComputeStatusMap(Int_t detElemId, Int_t manuId) const
194 {
195   /// Compute the status map for a given manu, and add it to our internal
196   /// fStatusMap internal storage
197   
198   AliCodeTimerAuto("(Int_t,Int_t)",0)
199     
200   AliMUONVCalibParam* param = new AliMUONCalibParamNI(1,AliMpConstants::ManuNofChannels(),
201                                                       detElemId,manuId,-1);    
202                                     
203   Bool_t ok = fStatusMap->Add(param);
204   if (!ok)
205   {
206     AliFatal(Form("Could not add manu %d of de %d",manuId,detElemId));
207   }
208                                   
209   AliMUONVCalibParam* neighbours = fkStatusMaker.Neighbours(detElemId,manuId);
210   
211   AliMUONVCalibParam* statusParam = fkStatusMaker.PadStatus(detElemId,manuId);
212   
213   Int_t n = neighbours->Dimension();
214   
215   for ( Int_t manuChannel = 0; manuChannel < param->Size(); ++manuChannel )
216   {
217     Int_t statusMap(0);
218     
219     Int_t x = neighbours->ValueAsIntFast(manuChannel,0);
220     if ( x < 0 ) 
221     {
222       // channel is not a valid one (i.e. (manuId,manuChannel) is not an existing pad)
223       statusMap = -1;//fgkSelfDead;
224       continue;
225     }
226         
227     for ( Int_t i = 0; i < n; ++i )
228     {
229       // Compute the statusmap related to the status of neighbouring
230       // pads. An invalid pad means "outside of edges".
231             
232       Int_t y = neighbours->ValueAsIntFast(manuChannel,i);      
233       Int_t m,c;
234       neighbours->UnpackValue(y,m,c);
235       if ( c < 0 ) continue;
236       Int_t status = 0;
237       if ( !m )
238       {
239         status = -1;
240       }
241       else
242       {
243         status = statusParam->ValueAsIntFast(c); //fkStatusMaker.PadStatus(detElemId,m,c);
244       }
245       if ( ( fMask != 0 ) && ( (status & fMask) != 0 ) )
246       {
247         statusMap |= (1<<i);
248       }
249     }    
250     param->SetValueAsIntFast(manuChannel,0,statusMap);
251   }
252   return param;
253 }
254
255 //_____________________________________________________________________________
256 void 
257 AliMUONPadStatusMapMaker::RefreshRejectProbabilities()
258 {
259   /// From the (fixed) fRejectProbabilities, compute
260   /// a fRejectList that will be valid for one event
261   /// If fRejectProbabilities=0x0 it means we're dealing with
262   /// trivial probabilities (0 or 1) and those are assumed to be already
263   /// in fRejectList then.
264   
265   if ( !fRejectProbabilities ) return;
266   
267   AliCodeTimerAuto("",0);
268   
269   fRejectList->Clear();
270   
271   TIter next(fRejectProbabilities->CreateIterator());
272   AliMUONVCalibParam* paramProba;
273   AliMUONVCalibParam* paramReject;
274   
275   while ( ( paramProba = static_cast<AliMUONVCalibParam*>(next()) ) )
276   {
277     paramReject = new AliMUONCalibParamNF(1,paramProba->Size(),paramProba->ID0(),paramProba->ID1(),0.0);
278     
279     Int_t n(0);
280     
281     for ( Int_t i = 0; i < paramProba->Size(); ++i ) 
282     {
283       Float_t proba = paramProba->ValueAsFloat(i);
284       Float_t x(proba);
285       
286       if ( proba > 0.0 && proba < 1.0 ) 
287       {
288         x = gRandom->Rndm();
289         proba = ( x < proba ) ? 1.0 : 0.0;
290       }
291       
292       if (proba>0.0)
293       {
294         ++n;
295         paramReject->SetValueAsFloat(i,0,proba);
296       }
297     }
298     if (n) fRejectList->Add(paramReject);
299   }
300 }
301
302 //_____________________________________________________________________________
303 Int_t
304 AliMUONPadStatusMapMaker::StatusMap(Int_t detElemId, Int_t manuId, 
305                                     Int_t manuChannel) const
306                                       
307 {
308   /// Get the pad status map
309   
310   AliMUONVCalibParam* param = static_cast<AliMUONVCalibParam*>(fStatusMap->FindObject(detElemId,manuId));
311   if (!param)
312   {
313     if ( fComputeOnDemand ) 
314     {
315       // not yet computed, so do it now
316       param = ComputeStatusMap(detElemId,manuId);
317     }
318     else
319     {
320       // we're locked. probably a bad manuId ?
321       return fgkSelfDead;
322     }
323   }
324   
325   Int_t statusMap = param->ValueAsInt(manuChannel);
326   
327   AliMUONVCalibParam* r = static_cast<AliMUONVCalibParam*>(fRejectList->FindObject(detElemId,manuId));
328   
329   if (r)
330   {
331     Float_t v= r->ValueAsFloat(manuChannel);
332     
333     assert (v==0.0 || v==1.0 ); 
334
335     if ( v > 0 ) 
336     {
337       statusMap |= fgkSelfDead;
338     }
339   }
340   
341   return statusMap;
342   
343 }