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