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