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