1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /* History of cvs commits:
21 * Revision 1.12 2006/09/07 18:31:08 kharlov
22 * Effective c++ corrections (T.Pocheptsov)
24 * Revision 1.11 2006/01/13 16:59:39 kharlov
25 * Rename classes to read raw data 2004
27 * Revision 1.10 2005/05/28 14:19:04 schutz
28 * Compilation warnings fixed by T.P.
32 //_________________________________________________________________________
33 // Class to calculate calibration parameters from beam tests etc.
34 // First pass - one should calcuate pedestals, in the second pass - gains.
36 // To calculate pedestals we scan pedestals events and fill histos for each
37 // channel. Then in each histo we find maximum and fit it with Gaussian in the
38 // visinity of the maximum. Extracted mean and width of distribution are put into
39 // resulting histogram which cheks results for 'reasonability'.
41 // To evaluate gains: scans beam events and calculate gain in the cristall with
42 // maximal energy deposition, assuming, that 90% of energy deposited in it.
43 // For each channel separate histogramm is filled. When scan is finished,
44 // these histograms are fitted with Gaussian and finds mean and width, which
45 // are put in final histogram. Finaly gains are checked for deviation from mean.
47 // Finally fills database.
50 // AliPHOSCalibrator * c = new AliPHOSCalibrator("path/galice.root") ;
51 // c->AddRun("path2/galice.root") ;
52 // c->ScanPedestals();
53 // c->CalculatePedestals();
54 // c->WritePedestals();
56 // c->CalculateGains() ;
59 //*-- Author : D.Peressounko (RRC KI)
60 //////////////////////////////////////////////////////////////////////////////
62 // --- ROOT system ---
65 #include "TObjString.h"
67 #include "TClonesArray.h"
69 // --- Standard library ---
71 // --- AliRoot header files ---
73 #include "AliPHOSCalibrManager.h"
74 #include "AliPHOSCalibrationData.h"
75 #include "AliPHOSCalibrator.h"
76 #include "AliPHOSConTableDB.h"
77 #include "AliRawReaderDate.h"
78 #include "AliPHOSRawStream2004.h"
79 #include "AliPHOSDigit.h"
81 #include "AliRawEventHeaderBase.h"
83 ClassImp(AliPHOSCalibrator)
86 //____________________________________________________________________________
87 AliPHOSCalibrator::AliPHOSCalibrator() :
88 TTask("AliPHOSCalibrator","Default"),
97 fConTableDB("Beamtest2002"),
98 fConTableDBFile("ConTableDB.root"),
107 //Default constuctor for root. Normally should not be used
110 //____________________________________________________________________________
111 AliPHOSCalibrator::AliPHOSCalibrator(const char* file, const char* title):
112 TTask("AliPHOSCalibrator",title),
121 fConTableDB("Beamtest2002"),
122 fConTableDBFile("ConTableDB.root"),
132 //Constructor which should normally be used.
133 //file: path/galice.root - header file
134 //title: branch name of PHOS reconstruction (e.g. "Default")
135 fRunList->SetOwner();
136 fRunList->Add(new TObjString(file));
139 //____________________________________________________________________________
140 AliPHOSCalibrator::AliPHOSCalibrator(const AliPHOSCalibrator & ctor) :
150 fConTableDB("Beamtest2002"),
151 fConTableDBFile("ConTableDB.root"),
160 // cpy ctor: no implementation yet
161 // requested by the Coding Convention
162 Fatal("cpy ctor", "not implemented") ;
165 //____________________________________________________________________________
166 AliPHOSCalibrator::~AliPHOSCalibrator()
176 delete fhPedestalsWid ;
182 //____________________________________________________________________________
183 void AliPHOSCalibrator::AddRun(const char * filename)
185 //Adds one more run to list of runs, which will be scanned in ScanXXX methods
187 TObjString * fn = new TObjString(filename) ;
189 fRunList=new TList() ;
190 fRunList->SetOwner() ;
195 TIter next(fRunList) ;
197 while((r=(TObjString *)(next()))){
198 if(fn->String().CompareTo(r->String())==0){
199 AliError(Form("Run already in list: %s",filename)) ;
207 //____________________________________________________________________________
208 void AliPHOSCalibrator::Exec(Option_t * option)
210 // reads parameters and does the calibration
211 ScanPedestals(option);
212 CalculatePedestals();
218 //____________________________________________________________________________
219 void AliPHOSCalibrator::Init(void)
221 // intializes everything
223 //check if ConTableDB already read
225 SetConTableDB(fConTableDBFile) ;
228 fNch = fctdb->GetNchanels() ;
229 fhPedestals = new TH1F("hPedestals","Pedestals mean",fNch,0.,fNch) ;
230 fhPedestalsWid= new TH1F("hPedestalsWid","Pedestals width",fNch,0.,fNch) ;
231 fhGains = new TH1F("hGains","Gains ",fNch,0.,fNch) ;
232 fhGainsWid = new TH1F("hGainsWid","Gains width",fNch,0.,fNch) ;
234 //____________________________________________________________________________
235 void AliPHOSCalibrator::SetConTableDB(const char * file,const char * name)
237 //Reads Connection Table database with name "name" from file "file"
239 if(file==0 || name == 0){
240 AliError(Form("Please, specify file with database and its title")) ;
243 if(fctdb && strcmp(fctdb->GetTitle(),name)==0) //already read
252 TFile * v = gROOT->GetFile(fConTableDBFile) ;
254 v = TFile::Open(fConTableDBFile) ;
256 AliError(Form("Can not open file with Connection Table DB: %s",fConTableDBFile.Data())) ;
259 fctdb = new AliPHOSConTableDB(*(dynamic_cast<AliPHOSConTableDB *>(v->Get("AliPHOSConTableDB")))) ;
263 //____________________________________________________________________________
264 void AliPHOSCalibrator::PlotPedestal(Int_t chanel)
266 //Plot histogram for a given channel, filled in Scan method
267 if(fPedHistos && fPedHistos->GetEntriesFast()){
268 static_cast<TH1F*>(fPedHistos->At(chanel))->Draw() ;
271 AliInfo(Form("Histograms not created yet! \n")) ;
274 //____________________________________________________________________________
275 void AliPHOSCalibrator::PlotPedestals(void)
277 // draws pedestals distribution
278 fhPedestals->Draw() ;
280 //____________________________________________________________________________
281 void AliPHOSCalibrator::PlotGain(Int_t chanel)
283 //Plot histogram for a given channel, filled in Scan method
284 if(fGainHistos && fGainHistos->GetEntriesFast()){
285 static_cast<TH1F*>(fGainHistos->At(chanel))->Draw() ;
288 AliInfo(Form("Histograms not created yet! \n")) ;
291 //____________________________________________________________________________
292 void AliPHOSCalibrator::PlotGains(void)
294 // draws gains distribution
297 //____________________________________________________________________________
298 void AliPHOSCalibrator::ScanPedestals(Option_t * option )
300 //scan all files in list fRunList and fill pedestal hisgrams
301 //option: "clear" - clear pedestal histograms filled up to now
302 // "deb" - plot file name currently processed
307 if(fPedHistos && strstr(option,"clear"))
308 fPedHistos->Delete() ;
310 fPedHistos = new TObjArray(fNch) ;
312 //Create histos for each channel, fills them and extracts mean values.
313 //First - prepare histos
315 for(ich=0;ich<fNch ;ich++){
316 TH1F * h = static_cast<TH1F *>(fPedHistos->At(ich)) ;
320 TString name("Pedestal for channel ") ;
322 fPedHistos->AddAt(new TH1F(n,name,fNChan,0,fNChan),ich) ;
326 TIter next(fRunList) ;
328 while((file = static_cast<TObjString *>(next()))){
329 if(strstr(option,"deb"))
330 printf("Processing file %s \n ",file->String().Data()) ;
333 AliRawReaderDate *rawReader = new AliRawReaderDate(file->String().Data()) ;
334 AliPHOSRawStream2004 *rawStream = new AliPHOSRawStream2004(rawReader) ;
335 rawStream->SetConTableDB(fctdb) ;
336 TClonesArray * digits = new TClonesArray("AliPHOSDigit",300) ;
338 //Scan all event in file
339 while(rawReader->NextEvent()){
340 //Is it PHYSICAL event
341 if(rawReader->GetType() == AliRawEventHeaderBase::kPhysicsEvent){
343 if(rawStream->ReadDigits(digits)){
344 if(rawStream->IsPEDevent()){
345 for(Int_t idigit = 0; idigit<digits->GetEntriesFast() ; idigit++){
346 AliPHOSDigit * digit = static_cast<AliPHOSDigit *>(digits->At(idigit) ) ;
347 ich = fctdb->AbsId2Raw(digit->GetId());
349 Float_t amp = digit->GetAmp() ;
350 TH1F * hh = dynamic_cast<TH1F*>(fPedHistos->At(ich)) ;
358 if(strstr(option,"deb"))
359 AliInfo(Form(" found %d events \n ",nevents)) ;
365 //____________________________________________________________________________
366 void AliPHOSCalibrator::CalculatePedestals()
368 //Fit histograms, filled in ScanPedestals method with Gaussian
369 //find mean and width, check deviation from mean for each channel.
371 if(!fPedHistos || !fPedHistos->At(0)){
372 AliError(Form("You should run ScanPedestals first!")) ;
376 //Now fit results with Gauss
377 TF1 * gs = new TF1("gs","gaus",0.,10000.) ;
379 for(ich=0;ich<fNch ;ich++){
380 TH1F * h = static_cast<TH1F *>(fPedHistos->At(ich)) ;
381 Int_t max = h->GetMaximumBin() ;
382 Axis_t xmin = max/2. ;
383 Axis_t xmax = max*3/2 ;
384 gs->SetRange(xmin,xmax) ;
386 par[0] = h->GetBinContent(max) ;
389 gs->SetParameters(par[0],par[1],par[2]) ;
391 gs->GetParameters(par) ;
392 fhPedestals->SetBinContent(ich,par[1]) ;
393 fhPedestals->SetBinError(ich,par[2]) ;
394 fhPedestalsWid->Fill(ich,par[2]) ;
398 //now check reasonability of results
399 TF1 * p0 = new TF1("p0","pol0",0.,fNch) ;
400 fhPedestals->Fit("p0","Q") ;
402 p0->GetParameters(&meanPed);
403 for(ich=0;ich<fNch ;ich++){
404 Float_t ped = fhPedestals->GetBinContent(ich) ;
405 if(ped < 0 || ped > meanPed+fAcceptCorr){
406 TString out("Pedestal of channel ") ;
410 out+= "it is too far from mean " ;
412 AliError(Form("PHOSCalibrator %s",out.Data())) ;
418 //____________________________________________________________________________
419 void AliPHOSCalibrator::ScanGains(Option_t * option)
421 //Scan all runs, listed in fRunList and fill histograms for all channels
422 //options: "clear" - clean histograms, filled up to now
423 // "deb" - print current file name
424 // "narrow" - scan only narrow beam events
428 if(fGainHistos && strstr(option,"clear"))
429 fGainHistos->Delete() ;
431 if(strstr(option,"deball"))
432 AliInfo(Form("creating array for %d channels \n",fNch)) ;
433 fGainHistos = new TObjArray(fNch) ;
436 //Create histos for each channel, fills them and extracts mean values.
437 //First - prepare histos
439 if(!fGainHistos->GetEntriesFast()){
441 for(ich=0;ich<fNch ;ich++){
444 TString name("Gains for channel ") ;
446 fGainHistos->AddAt(new TH1F(n,name,fNGainBins,0,fGainMax),ich) ;
447 // static_cast<TH1F*>(fGainHistos->At(ich))->Sumw2() ;
451 TIter next(fRunList) ;
453 while((file = static_cast<TObjString *>(next()))){
455 AliRawReaderDate *rawReader = new AliRawReaderDate(file->String().Data()) ;
456 AliPHOSRawStream2004 *rawStream = new AliPHOSRawStream2004(rawReader) ;
457 rawStream->SetConTableDB(fctdb) ;
459 TClonesArray * digits = new TClonesArray("AliPHOSDigit",300) ;
461 //Scan all event in file
462 while(rawReader->NextEvent()){
463 //Is it PHYSICAL event
464 if(rawReader->GetType() == AliRawEventHeaderBase::kPhysicsEvent){
465 if(rawStream->ReadDigits(digits)){
467 if(rawStream->IsNELevent() || rawStream->IsWELevent()){
469 AliPHOSDigit * digit ;
472 for(Int_t idigit = 0; idigit<digits->GetEntriesFast() ; idigit++){
473 digit = static_cast<AliPHOSDigit *>(digits->At(idigit) ) ;
474 if(digit->GetAmp() > max){
476 max = digit->GetAmp() ;
479 digit = static_cast<AliPHOSDigit *>(digits->At(imax) ) ;
480 Int_t ich = fctdb->AbsId2Raw(digit->GetId());
482 Float_t pedestal = fhPedestals->GetBinContent(ich) ;
483 const Float_t kshowerInCrystall = 0.9 ;
484 Float_t gain = fBeamEnergy*kshowerInCrystall/
485 (digit->GetAmp() - pedestal) ;
486 static_cast<TH1F*>(fGainHistos->At(ich))->Fill(gain) ;
495 if(strstr(option,"deb"))
496 AliInfo(Form(" found %d events \n",nevents)) ;
499 //____________________________________________________________________________
500 void AliPHOSCalibrator::CalculateGains(void)
504 if(!fGainHistos || !fGainHistos->GetEntriesFast()){
505 AliError(Form("You should run ScanGains first!")) ;
509 //Fit results with Landau
510 TF1 * gs = new TF1("gs","landau",0.,10000.) ;
512 for(ich=0;ich<fNch ;ich++){
513 TH1F * h = static_cast<TH1F *>(fGainHistos->At(ich)) ;
514 Int_t bmax = h->GetMaximumBin() ;
515 Axis_t center = h->GetBinCenter(bmax) ;
516 Axis_t xmin = center - 0.01 ;
517 Axis_t xmax = center + 0.02 ;
518 gs->SetRange(xmin,xmax) ;
520 par[0] = h->GetBinContent(bmax) ;
523 gs->SetParameters(par[0],par[1],par[2]) ;
525 gs->GetParameters(par) ;
526 fhGains->SetBinContent(ich,par[1]) ;
527 fhGains->SetBinError(ich,par[2]) ;
528 fhGainsWid->Fill(ich,par[2]) ;
532 //now check reasonability of results
533 TF1 * p0 = new TF1("p0","pol0",0.,fNch) ;
534 fhGains->Fit("p0","Q") ;
536 p0->GetParameters(&meanGain);
537 for(ich=0;ich<fNch ;ich++){
538 Float_t gain = fhGains->GetBinContent(ich) ;
539 if(gain < meanGain/fGainAcceptCorr || gain > meanGain*fGainAcceptCorr){
540 TString out("Gain of channel ") ;
544 out+= "it is too far from mean " ;
546 AliError(Form("PHOSCalibrator %s",out.Data())) ;
552 //____________________________________________________________________________
553 void AliPHOSCalibrator::ReadFromASCII(const char * filename){
554 // We read pedestals and gains from *.dat file with following format:
555 // 0 0 0 0 37.09 1972. // next nmodrows*nmodcols*ncryrows*ncrycols lines
556 // 0 0 0 1 28.53 2072. // contains <RR CC r c ped peak>
557 // 0 0 0 2 30.93 1938. //
558 // where module is an array of 8*8 crystals and RR and CC are module raw and column position
559 FILE * file = fopen(filename, "r");
561 Error("ReadFromASCII", "could not open file %s", filename);
564 if(!fctdb || !fhPedestals || !fhGains){
572 Int_t modRaw,modCol,raw,col;
575 while(fscanf(file,"%d %d %d %d %f %f",&modRaw,&modCol,&raw,&col,&ped,&pik)==6){
576 //Calculate plain crystal position:
577 Int_t rawPosition = (modRaw*8+raw)*fctdb->GetNColumns()+modCol*8+col ;
578 fhPedestals->SetBinContent(rawPosition,ped) ;
580 fhGains->SetBinContent(rawPosition,1./pik);
582 fhGains->SetBinContent(rawPosition,0.);
585 if(nread != fctdb->GetNColumns()*fctdb->GetNRaws()){
586 Error("ReadFromASCII","Read %d parameters instead of %d\n",nread,fctdb->GetNColumns()*fctdb->GetNRaws()) ;
590 //_____________________________________________________________________________
591 void AliPHOSCalibrator::WritePedestals(const char * version)
593 //Write calculated data to file using AliPHOSCalibrManager
594 //version and validitirange (begin-end) will be used to identify data
597 AliError(Form("\n Please, supply Connection Table DB (use SetConTableDB()) \n" )) ;
601 AliPHOSCalibrationData ped("Pedestals",version);
602 for(Int_t i=0; i<fNch;i++){
603 Int_t absid=fctdb->Raw2AbsId(i) ;
604 ped.SetData(absid,fhPedestals->GetBinContent(i)) ;
605 ped.SetDataCheck(absid,fhPedestalsWid->GetBinContent(i)) ;
608 // //evaluate validity range
610 // TIter next(fRunList) ;
611 // Int_t ibegin=99999;
613 // TObjString * file ;
614 // while((file=((TObjString*)next()))){
615 // TString s = file->GetString() ;
616 // TString ss = s(s.Last('_'),s.Last('.'));
618 // if(sscanf(ss.Data(),"%d",&tmp)){
625 // ped.SetValidityRange(ibegin,iend) ;
628 // ped.SetValidityRange(begin,end) ;
630 //check, may be Manager instance already configured?
631 AliPHOSCalibrManager * cmngr = AliPHOSCalibrManager::GetInstance() ;
633 AliWarning(Form("Using database file 'PHOSBTCalibration.root'")) ;
634 cmngr = AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
636 cmngr->WriteData(ped) ;
638 //_____________________________________________________________________________
639 void AliPHOSCalibrator::ReadPedestals(const char * version)
641 //Read data from file using AliPHOSCalibrManager
642 //version and range will be used to choose proper data
644 AliPHOSCalibrationData ped("Pedestals",version);
645 AliPHOSCalibrManager * cmngr = AliPHOSCalibrManager::GetInstance() ;
647 AliWarning(Form("Using database file 'PHOSBTCalibration.root'")) ;
648 cmngr = AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
650 cmngr->GetParameters(ped) ;
651 Int_t npeds=ped.NChannels() ;
652 fNch = fctdb->GetNchanels() ;
655 fhPedestals = new TH1F("hPedestals","Pedestals mean",fNch,0.,fNch) ;
656 for(Int_t i=0;i<npeds;i++){
657 Int_t raw =fctdb->AbsId2Raw(i) ;
659 fhPedestals->SetBinContent(raw-1,ped.Data(i)) ;
660 fhPedestals->SetBinError(raw-1,ped.DataCheck(i)) ;
664 //_____________________________________________________________________________
665 void AliPHOSCalibrator::ReadGains(const char * version)
667 //Read data from file using AliPHOSCalibrManager
668 //version and range will be used to choose proper data
670 AliPHOSCalibrationData gains("Gains",version);
671 AliPHOSCalibrManager * cmngr = AliPHOSCalibrManager::GetInstance() ;
673 AliWarning(Form("Using database file 'PHOSBTCalibration.root'")) ;
674 cmngr = AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
676 cmngr->GetParameters(gains) ;
677 Int_t npeds=gains.NChannels() ;
678 fNch = fctdb->GetNchanels() ;
681 fhGains = new TH1F("hGainss","Gains mean",fNch,0.,fNch) ;
682 for(Int_t i=0;i<npeds;i++){
683 Int_t raw =fctdb->AbsId2Raw(i) ;
685 fhGains->SetBinContent(raw-1,gains.Data(i)) ;
686 fhGains->SetBinError(raw-1,gains.DataCheck(i)) ;
690 //_____________________________________________________________________________
691 void AliPHOSCalibrator::WriteGains(const char * version)
693 //Write gains through AliPHOSCalibrManager
694 //version and validity range(begin-end) are used to identify data
697 AliError(Form("\n Please, supply Connection Table DB (use SetConTableDB()) \n" )) ;
701 AliPHOSCalibrationData gains("Gains",version);
702 for(Int_t i=0; i<fNch;i++){
703 Int_t absid=fctdb->Raw2AbsId(i) ;
704 gains.SetData(absid,fhGains->GetBinContent(i)) ;
705 gains.SetDataCheck(absid,fhGainsWid->GetBinContent(i)) ;
708 // TIter next(fRunList) ;
709 // Int_t ibegin=99999;
711 // TObjString * file ;
712 // while((file=((TObjString*)next()))){
713 // TString s = file->GetString() ;
714 // TSubString ss = s(s.Last('_'),s.Last('.'));
716 // if(sscanf(ss.Data(),"%d",&tmp)){
723 // gains.SetValidityRange(ibegin,iend) ;
726 // gains.SetValidityRange(begin,end) ;
727 AliPHOSCalibrManager * cmngr = AliPHOSCalibrManager::GetInstance() ;
729 AliWarning(Form("Using database file 'PHOSBTCalibration.root'")) ;
730 cmngr = AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
732 cmngr->WriteData(gains) ;
734 //_____________________________________________________________________________
735 void AliPHOSCalibrator::Print(const Option_t *)const
738 AliInfo(Form("--------------PHOS Calibrator-----------------\n")) ;
739 printf("Files to handle:\n") ;
740 TIter next(fRunList) ;
742 while((r=(TObjString *)(next())))
743 printf(" %s\n",r->GetName()) ;
745 printf("Name of ConTableDB:.....................%s\n",fConTableDB.Data()) ;
746 printf("File of ConTableDB:.....................%s\n",fConTableDBFile.Data() ) ;
747 printf("Maximal deviation from mean Gain (factor):.%f\n",fGainAcceptCorr) ;
748 printf("Maximal deviation of Pedestal from mean:...%f\n",fAcceptCorr) ;
749 printf("Range used in Gain histos:..............%f\n",fGainMax) ;
750 printf("Number of bins in Gain histos:..........%d\n",fNGainBins) ;
751 printf("Number of channels to calibrate:........%d\n",fNch) ;
752 printf("Number of bins in pedestal histos:......%d\n",fNChan) ;
753 printf("--------------------------------------------------\n") ;