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