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