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