bug fixed
[u/mrichter/AliRoot.git] / MUON / AliMUONPedestal.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 #include "AliMUONPedestal.h"
19 #include "AliMUONErrorCounter.h"
20 #include "AliMUONVStore.h"
21 #include "AliMUON2DMap.h"
22 #include "AliMUONCalibParamND.h"
23 #include "AliMpConstants.h"
24 #include <TString.h>
25 #include <TTimeStamp.h>
26 #include <TMath.h>
27 #include <TTree.h>
28 #include <TFile.h>
29 #include <TH1F.h>
30 #include <THashList.h>
31 #include <Riostream.h>
32
33 #include <sstream>
34
35 //-----------------------------------------------------------------------------
36 /// \class AliMUONPedestal
37 ///
38 /// Implementation of the pedestal computing
39 ///
40 /// add
41 /// 
42 ///
43 /// \author Alberto Baldisseri, JL Charvet (05/05/2009)
44 //-----------------------------------------------------------------------------
45
46 /// \cond CLASSIMP
47 ClassImp(AliMUONPedestal)
48 /// \endcond
49
50 //______________________________________________________________________________
51 AliMUONPedestal::AliMUONPedestal()
52 : TObject(),
53 //fN(0),
54 fNCurrentEvents(0),
55 fNEvthreshold(0),
56 fSorting(0),
57 fNEvents(0),
58 fRunNumber(0),
59 fNChannel(0),
60 fNManu(0),
61 fNManuConfig(0),
62 fConfig(1),
63 fErrorBuspatchTable(new AliMUON2DMap(kFALSE)),
64 fManuBuspatchTable(new AliMUON2DMap(kFALSE)),
65 fManuBPoutofconfigTable(new AliMUON2DMap(kFALSE)),
66 fDate(new TTimeStamp()),
67 fFilcout(0),
68 fHistoFileName(),
69 fPedestalStore(new AliMUON2DMap(kTRUE)),
70 fIndex(-1),
71 fPrefixDA()
72 {
73 /// Default constructor
74 }
75
76 //______________________________________________________________________________
77 AliMUONPedestal::~AliMUONPedestal()
78 {
79 /// Destructor
80   delete fErrorBuspatchTable;
81   delete fManuBuspatchTable;
82   delete fPedestalStore;
83   delete fManuBPoutofconfigTable;
84 }
85
86 //______________________________________________________________________________
87 const char* 
88 AliMUONPedestal::GetHistoFileName() const
89 {
90   /// Return the name of file we use to store histograms
91   return fHistoFileName.Data();
92 }
93
94 //______________________________________________________________________________
95 void AliMUONPedestal::LoadConfig(const char* dbfile)
96 {
97   /// Load MUONTRK configuration from ascii file "dbfile" (in DetDB)
98
99   Int_t manuId;
100   Int_t busPatchId;
101
102   ifstream filein(dbfile,ios::in);
103   
104   while (!filein.eof())
105     { 
106       filein >> busPatchId >> manuId;
107
108       AliMUONErrorCounter* manuCounter;
109       AliMUONVCalibParam* ped = 
110         static_cast<AliMUONVCalibParam*>(fPedestalStore ->FindObject(busPatchId, manuId));
111
112       if (!ped) {
113         fNManuConfig++;
114         fNChannel+=64;
115   ped = new AliMUONCalibParamND(2, AliMpConstants::ManuNofChannels(),busPatchId, manuId, -1.); // put default wise -1, not connected channel
116         fPedestalStore ->Add(ped);  
117
118         if (!(manuCounter = static_cast<AliMUONErrorCounter*>(fManuBuspatchTable->FindObject(busPatchId,manuId))))
119           {
120             // New (buspatch,manu)
121             manuCounter = new AliMUONErrorCounter(busPatchId,manuId);
122             fManuBuspatchTable->Add(manuCounter);
123           }
124       }
125     }
126
127 //______________________________________________________________________________
128 void AliMUONPedestal::MakePed(Int_t busPatchId, Int_t manuId, Int_t channelId, Int_t charge)
129 {
130   static Int_t warn=0;
131   /// Compute pedestals values
132   AliMUONVCalibParam* ped = 
133     static_cast<AliMUONVCalibParam*>(fPedestalStore ->FindObject(busPatchId, manuId));
134
135   if (!ped)   
136     {
137       if(fConfig) 
138         {  // Fill out_of_config (buspatch,manu) table
139           if (!(static_cast<AliMUONErrorCounter*>(fManuBPoutofconfigTable->FindObject(busPatchId,manuId))))
140             fManuBPoutofconfigTable->Add(new AliMUONErrorCounter(busPatchId,manuId));
141           if(warn<10) cout << " !!! WARNING  : busPatchId = " << busPatchId << " manuId = " << manuId << " not in the Detector configuration " << endl;
142           else if(warn==10) cout << " !!! see .log file for an exhaustive list of (busPatchId, manuId) out of Detector configuration \n" << endl; 
143            warn++;
144            (*fFilcout) << " !!! WARNING  : busPatchId = " << busPatchId << " manuId = " << manuId << " not in the Detector configuration " << endl; 
145         }
146       else {fNManu++;}
147       fNChannel+=64;
148       // put default wise -1, not connected channel
149       ped = new AliMUONCalibParamND(2, AliMpConstants::ManuNofChannels(),busPatchId, manuId, -1.); 
150       fPedestalStore ->Add(ped);  
151     }
152
153   // Initialization for the first value
154   if (ped->ValueAsDouble(channelId, 0) == -1)  
155     { 
156       if(fConfig && channelId == 0){fNManu++;}
157       ped->SetValueAsDouble(channelId, 0, 0.);
158     }
159   if (ped->ValueAsDouble(channelId, 1) == -1) ped->SetValueAsDouble(channelId, 1, 0.);
160
161   Double_t pedMean  = ped->ValueAsDouble(channelId, 0) + (Double_t) charge;
162   Double_t pedSigma = ped->ValueAsDouble(channelId, 1) + (Double_t) charge*charge;
163
164   ped->SetValueAsDouble(channelId, 0, pedMean);
165   ped->SetValueAsDouble(channelId, 1, pedSigma);
166
167   AliMUONErrorCounter* manuCounter;
168   if (!(manuCounter = static_cast<AliMUONErrorCounter*>(fManuBuspatchTable->FindObject(busPatchId,manuId))))
169     {
170       // New (buspatch,manu)
171       manuCounter = new AliMUONErrorCounter(busPatchId,manuId);
172       fManuBuspatchTable->Add(manuCounter);
173     }
174   else
175     {
176       // Existing buspatch
177       manuCounter->Increment();
178     }   
179 }
180 //______________________________________________________________________________
181 void AliMUONPedestal::Finalize()
182 {
183   /// final polishing of the store
184   
185   Double_t pedMean;
186   Double_t pedSigma;
187   Double_t pedSigmalimit=0.5;
188   Int_t busPatchId;
189   Int_t manuId;
190   Int_t channelId;
191
192   // print in logfile
193   if (fErrorBuspatchTable->GetSize())
194     {
195       cout<<"\nWarning: Buspatches with less statistics (due to parity errors)"<<endl;
196       (*fFilcout)<<"\nWarning: Buspatches with less statistics (due to parity errors)"<<endl;
197       TIter nextParityError(fErrorBuspatchTable->CreateIterator());
198       AliMUONErrorCounter* parityerror;
199       while((parityerror = static_cast<AliMUONErrorCounter*>(nextParityError())))
200         {
201           cout<<"  bp "<<parityerror->BusPatch()<<": events used = "<<fNEvents-parityerror->Events()<<endl;
202           (*fFilcout)<<"  bp "<<parityerror->BusPatch()<<": events used = "<<fNEvents-parityerror->Events()<<endl;
203         }
204     }
205
206   Int_t nADC4090=0;
207   Int_t nADCmax=0;
208   // iterator over pedestal
209   TIter next(fPedestalStore ->CreateIterator());
210   AliMUONVCalibParam* ped;
211
212   while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
213     {
214       busPatchId              = ped->ID0();
215       manuId                  = ped->ID1();
216       if(manuId==0)
217         {
218           cout << "Warning: ManuId = " << manuId << " !!! in  BP = " << busPatchId << endl;
219           (*fFilcout) << "Warning: ManuId = " << manuId << " !!! in  BP = " << busPatchId << endl;
220         }
221       Int_t eventCounter;
222       // Correct the number of events for buspatch with errors
223       AliMUONErrorCounter* errorCounter;
224       if ((errorCounter = (AliMUONErrorCounter*)fErrorBuspatchTable->FindObject(busPatchId)))
225         {
226           eventCounter = fNEvents - errorCounter->Events();
227         }
228       else
229         {
230           eventCounter = fNEvents;
231         }
232       Int_t occupancy=0; // channel missing in raw data or read but rejected (case of parity error)
233       // value of (buspatch, manu) occupancy
234       AliMUONErrorCounter* manuCounter;
235       manuCounter = static_cast<AliMUONErrorCounter*>(fManuBuspatchTable->FindObject(busPatchId,manuId));
236       if(eventCounter>0)occupancy = manuCounter->Events()/64/eventCounter;
237       if(occupancy>1)
238         {
239           cout << "Warning: ManuId = " << manuId << " !!! in  BP = " << busPatchId << " occupancy (>1) = " << occupancy << endl;
240           (*fFilcout) << "Warning: ManuId = " << manuId << " !!! in  BP = " << busPatchId << " occupancy (>1) = " << occupancy <<endl;
241         }
242
243       for (channelId = 0; channelId < ped->Size() ; ++channelId) 
244         {
245           pedMean  = ped->ValueAsDouble(channelId, 0);
246
247           if (pedMean >= 0) // connected channels
248             {
249               ped->SetValueAsDouble(channelId, 0, pedMean/(Double_t)eventCounter);
250               pedMean  = ped->ValueAsDouble(channelId, 0);
251               pedSigma = ped->ValueAsDouble(channelId, 1);
252               ped->SetValueAsDouble(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Double_t)eventCounter - pedMean*pedMean)));
253
254               if(eventCounter < fNEvthreshold )
255                 { nADCmax++; ped->SetValueAsDouble(channelId, 0, ADCMax());
256                   ped->SetValueAsDouble(channelId, 1, ADCMax());}
257               if( ped->ValueAsDouble(channelId, 1) < pedSigmalimit )
258                 { nADC4090++; ped->SetValueAsDouble(channelId, 0, ADCMax()-5);
259                   ped->SetValueAsDouble(channelId, 1, ADCMax()-5);}
260               if(manuId == 0 || occupancy>1)
261                 { nADCmax++; ped->SetValueAsDouble(channelId, 0, ADCMax());
262                   ped->SetValueAsDouble(channelId, 1, ADCMax());
263                   if(occupancy>1 && channelId==0)ped->SetValueAsDouble(channelId, 0, ADCMax()+occupancy);}
264             }
265           else
266             { nADCmax++; ped->SetValueAsDouble(channelId, 0, ADCMax());
267               ped->SetValueAsDouble(channelId, 1, ADCMax());}
268         }
269     }
270   if(nADCmax>0)
271     { char* detail=Form("Warning: Nb of Channels with bad Pedestal (Ped=4095) = %d over %d",nADCmax,fNChannel);
272       printf("%s\n",detail);
273       (*fFilcout) <<  detail << endl; }
274   if(nADC4090>0)
275     { char* detail=Form("Warning: Nb of Channels with PedSigma<0.5 (Ped=4090) = %d over %d",nADC4090,fNChannel);
276       printf("%s\n",detail);
277       (*fFilcout) <<  detail << endl; }
278 }
279 //______________________________________________________________________________
280 void AliMUONPedestal::MakeASCIIoutput(ostream& out) const
281 {
282   /// put pedestal store in the output stream
283
284   out<<"//===========================================================================" << endl;
285   out<<"//                 Pedestal file calculated by "<< fPrefixDA.Data() << endl;
286   out<<"//===========================================================================" << endl;
287   out<<"//     * Run           : " << fRunNumber << endl; 
288   out<<"//     * Date          : " << fDate->AsString("l") <<endl;
289   out<<"//     * Statictics    : " << fNEvents << endl;
290   if(fConfig)
291     out<<"//     * Nb of MANUS   : " << fNManuConfig << " read in the Det. config. " << endl;
292   out<<"//     * Nb of MANUS   : " << fNManu << " read in raw data " << endl;
293   out<<"//     * Nb of MANUS   : " << fNChannel/64 << " written in pedestal file " << endl;
294   out<<"//     * Nb of channels: " << fNChannel << endl;
295   out<<"//"<<endl;
296   out<<"//     * Below " << fNEvthreshold << " events=> Ped.&sig.=4095" << endl;
297   out<<"//     * PedSigma < 0.5 => Ped.&sig.=4090" << endl;
298
299   if (fErrorBuspatchTable->GetSize())
300     {
301       out<<"//"<<endl;
302       out<<"//    * Buspatches with less statistics (due to parity errors)"<<endl;
303       TIter next(fErrorBuspatchTable->CreateIterator());
304       AliMUONErrorCounter* parityerror;
305       while((parityerror = static_cast<AliMUONErrorCounter*>(next())))
306         {
307           if(fNEvents-parityerror->Events()>fNEvthreshold)
308             { out<<"//      BusPatch = "<<parityerror->BusPatch()<<"\t Nevents used = "<<fNEvents-parityerror->Events()<<endl; }
309           else
310             { out<<"//      BusPatch = "<<parityerror->BusPatch()<<"\t Nevents used = "<<fNEvents-parityerror->Events()<< " (Ped.&sig.=4095)" << endl; }
311         }
312     }  
313   Int_t writitle=0;
314   Int_t occupancy=1;
315   if(occupancy)
316     {
317       TIter next(fPedestalStore ->CreateIterator());
318       AliMUONVCalibParam* ped;
319       while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
320         {
321           Int_t busPatchId = ped->ID0();
322           Int_t manuId = ped->ID1();
323           Double_t pedMean  = ped->ValueAsDouble(0, 0); // check pedestal value for channelId=0
324
325           if(pedMean>ADCMax()) 
326             {
327               writitle++;
328               if(writitle==1){ 
329                 out<<"//"<<endl;
330                 out<<"//    * Puzzling (Buspatch,Manu) read in raw data ? (Ped.&sig.=4095)"<<endl;}
331               occupancy=TMath::Nint(pedMean-ADCMax());
332               ped->SetValueAsDouble(0, 0, ADCMax());
333               out<<"//      BusPatch = "<< busPatchId <<"\t ManuId =  "<< manuId << "\t occupancy = " << occupancy  << endl;
334             }
335
336           if (manuId==0 || (fConfig && static_cast<AliMUONErrorCounter*>(fManuBPoutofconfigTable->FindObject(busPatchId,manuId))))
337             {
338               writitle++;
339               if(writitle==1){ 
340                 out<<"//"<<endl;
341                 out<<"//    * Puzzling (Buspatch,Manu) read in raw data ? (Ped.&sig.=4095)"<<endl;}
342               out<<"//      BusPatch = "<< busPatchId <<"\t ManuId =  "<< manuId << "\t missing in the mapping" << endl;
343             }
344         }
345     }
346   out<<"//"<<endl;
347   out<<"//---------------------------------------------------------------------------" << endl;
348   out<<"//---------------------------------------------------------------------------" << endl;
349   out<<"//      BP     MANU     CH.      MEAN    SIGMA"<<endl;
350   out<<"//---------------------------------------------------------------------------" << endl;
351
352   TIter next(fPedestalStore->CreateIterator());
353   AliMUONVCalibParam* ped;  
354
355   // Sorting 
356   if  (fSorting)
357     {
358       cout << " ..... sorting pedestal values ....."  << endl;
359       THashList pedtable(100,2);
360       while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
361       {
362         pedtable.Add(ped);
363       }
364       pedtable.Sort();
365       //      iterator over sorted pedestal
366       TIter nextSorted(&pedtable);
367       while ( (ped = (AliMUONVCalibParam*)(nextSorted()) ) )
368       {
369         Int_t busPatchId = ped->ID0();
370         Int_t manuId = ped->ID1();
371         for ( Int_t channelId = 0; channelId < ped->Size(); ++channelId ) 
372         {
373           Double_t pedMean  = ped->ValueAsDouble(channelId, 0);
374           Double_t pedSigma = ped->ValueAsDouble(channelId, 1);
375           out << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t" << pedMean <<"\t"<< pedSigma << endl;
376         }
377       }
378     }
379   else
380     {
381       while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
382         {
383           Int_t busPatchId = ped->ID0();
384           Int_t manuId = ped->ID1();
385           for ( Int_t channelId = 0; channelId < ped->Size(); ++channelId ) 
386             {
387               Double_t pedMean  = ped->ValueAsDouble(channelId, 0);
388               Double_t pedSigma = ped->ValueAsDouble(channelId, 1);
389               out << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t" << pedMean <<"\t"<< pedSigma << endl;
390             }
391         }
392     }
393 }
394
395 //______________________________________________________________________________
396 void AliMUONPedestal::MakeControlHistos()
397 {
398   /// Create control histograms
399   if (fIndex>=0) return; // Pedestal run (fIndex=-1)
400
401   Double_t pedMean;
402   Double_t pedSigma;
403   Double_t evt;
404   Int_t busPatchId;
405   Int_t manuId;
406   Int_t channelId;
407
408 // histo
409   TFile*  histoFile = 0;
410   TTree* tree = 0;
411   TH1F* pedMeanHisto = 0;
412   TH1F* pedSigmaHisto = 0;
413     
414   fHistoFileName=Form("%s.root",fPrefixDA.Data());
415   histoFile = new TFile(fHistoFileName,"RECREATE","MUON Tracking pedestals");
416
417   Char_t name[255];
418   Char_t title[255];
419   sprintf(name,"pedmean_allch");
420   sprintf(title,"Pedestal mean all channels");
421   Int_t nx = ADCMax()+1;
422   Int_t xmin = 0;
423   Int_t xmax = ADCMax(); 
424   pedMeanHisto = new TH1F(name,title,nx,xmin,xmax);
425   pedMeanHisto->SetDirectory(histoFile);
426
427   sprintf(name,"pedsigma_allch");
428   sprintf(title,"Pedestal sigma all channels");
429   nx = 201;
430   xmin = 0;
431   xmax = 200; 
432   pedSigmaHisto = new TH1F(name,title,nx,xmin,xmax);
433   pedSigmaHisto->SetDirectory(histoFile);
434
435   tree = new TTree("t","Pedestal tree");
436   tree->Branch("bp",&busPatchId,"bp/I");
437   tree->Branch("manu",&manuId,",manu/I");
438   tree->Branch("channel",&channelId,",channel/I");
439   tree->Branch("pedMean",&pedMean,",pedMean/D");
440   tree->Branch("pedSigma",&pedSigma,",pedSigma/D");
441   tree->Branch("nevt",&evt,",evt/D");
442
443   // iterator over pedestal
444   TIter next(fPedestalStore ->CreateIterator());
445   AliMUONVCalibParam* ped;
446   AliMUONErrorCounter* manuCounter;
447   
448   while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
449   {
450     busPatchId = ped->ID0();
451     manuId = ped->ID1();
452     
453     for ( channelId = 0; channelId < ped->Size(); ++channelId ) 
454     {
455       pedMean  = ped->ValueAsDouble(channelId, 0);
456       pedSigma = ped->ValueAsDouble(channelId, 1);
457       manuCounter = static_cast<AliMUONErrorCounter*>(fManuBuspatchTable->FindObject(busPatchId,manuId));
458       evt = manuCounter->Events()/64;
459           
460       pedMeanHisto->Fill(pedMean);
461       pedSigmaHisto->Fill(pedSigma);
462       tree->Fill();  
463     }
464   }
465     
466   histoFile->Write();  
467   histoFile->Close(); 
468
469 }