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.10 2005/05/28 14:19:04 schutz
22 * Compilation warnings fixed by T.P.
26 //_________________________________________________________________________
27 // Class to calculate calibration parameters from beam tests etc.
28 // First pass - one should calcuate pedestals, in the second pass - gains.
30 // To calculate pedestals we scan pedestals events and fill histos for each
31 // channel. Then in each histo we find maximum and fit it with Gaussian in the
32 // visinity of the maximum. Extracted mean and width of distribution are put into
33 // resulting histogram which cheks results for 'reasonability'.
35 // To evaluate gains: scans beam events and calculate gain in the cristall with
36 // maximal energy deposition, assuming, that 90% of energy deposited in it.
37 // For each channel separate histogramm is filled. When scan is finished,
38 // these histograms are fitted with Gaussian and finds mean and width, which
39 // are put in final histogram. Finaly gains are checked for deviation from mean.
41 // Finally fills database.
44 // AliPHOSCalibrator * c = new AliPHOSCalibrator("path/galice.root") ;
45 // c->AddRun("path2/galice.root") ;
46 // c->ScanPedestals();
47 // c->CalculatePedestals();
48 // c->WritePedestals();
50 // c->CalculateGains() ;
53 //*-- Author : D.Peressounko (RRC KI)
54 //////////////////////////////////////////////////////////////////////////////
56 // --- ROOT system ---
59 #include "TObjString.h"
61 #include "TClonesArray.h"
63 // --- Standard library ---
65 // --- AliRoot header files ---
67 #include "AliPHOSCalibrManager.h"
68 #include "AliPHOSCalibrationData.h"
69 #include "AliPHOSCalibrator.h"
70 #include "AliPHOSConTableDB.h"
71 #include "AliRawReaderDate.h"
72 #include "AliPHOSRawStream2004.h"
73 #include "AliPHOSDigit.h"
78 #define PHYSICS_EVENT 7
81 ClassImp(AliPHOSCalibrator)
84 //____________________________________________________________________________
85 AliPHOSCalibrator::AliPHOSCalibrator():TTask("AliPHOSCalibrator","Default")
87 //Default constuctor for root. Normally should not be used
96 fConTableDB = "Beamtest2002" ;
97 fConTableDBFile = "ConTableDB.root" ;
99 //____________________________________________________________________________
100 AliPHOSCalibrator::AliPHOSCalibrator(const char* file, const char* title):
101 TTask("AliPHOSCalibrator",title)
103 //Constructor which should normally be used.
104 //file: path/galice.root - header file
105 //title: branch name of PHOS reconstruction (e.g. "Default")
108 fRunList = new TList() ;
109 fRunList->SetOwner() ;
110 fRunList->Add(new TObjString(file)) ;
117 fAcceptCorr = 10 ; //Maximal deviation from mean, considered as normal
119 fGainAcceptCorr = 5 ; //Factor for gain deviation
125 fConTableDB = "Beamtest2002" ;
126 fConTableDBFile = "ConTableDB.root" ;
129 //____________________________________________________________________________
130 AliPHOSCalibrator::~AliPHOSCalibrator()
140 delete fhPedestalsWid ;
146 //____________________________________________________________________________
147 void AliPHOSCalibrator::AddRun(const char * filename)
149 //Adds one more run to list of runs, which will be scanned in ScanXXX methods
151 TObjString * fn = new TObjString(filename) ;
153 fRunList=new TList() ;
154 fRunList->SetOwner() ;
159 TIter next(fRunList) ;
161 while((r=(TObjString *)(next()))){
162 if(fn->String().CompareTo(r->String())==0){
163 AliError(Form("Run already in list: %s",filename)) ;
171 //____________________________________________________________________________
172 void AliPHOSCalibrator::Exec(Option_t * option)
174 // reads parameters and does the calibration
175 ScanPedestals(option);
176 CalculatePedestals();
182 //____________________________________________________________________________
183 void AliPHOSCalibrator::Init(void)
185 // intializes everything
187 //check if ConTableDB already read
189 SetConTableDB(fConTableDBFile) ;
192 fNch = fctdb->GetNchanels() ;
193 fhPedestals = new TH1F("hPedestals","Pedestals mean",fNch,0.,fNch) ;
194 fhPedestalsWid= new TH1F("hPedestalsWid","Pedestals width",fNch,0.,fNch) ;
195 fhGains = new TH1F("hGains","Gains ",fNch,0.,fNch) ;
196 fhGainsWid = new TH1F("hGainsWid","Gains width",fNch,0.,fNch) ;
198 //____________________________________________________________________________
199 void AliPHOSCalibrator::SetConTableDB(const char * file,const char * name)
201 //Reads Connection Table database with name "name" from file "file"
203 if(file==0 || name == 0){
204 AliError(Form("Please, specify file with database and its title")) ;
207 if(fctdb && strcmp(fctdb->GetTitle(),name)==0) //already read
216 TFile * v = gROOT->GetFile(fConTableDBFile) ;
218 v = TFile::Open(fConTableDBFile) ;
220 AliError(Form("Can not open file with Connection Table DB: %s",fConTableDBFile.Data())) ;
223 fctdb = new AliPHOSConTableDB(*(dynamic_cast<AliPHOSConTableDB *>(v->Get("AliPHOSConTableDB")))) ;
227 //____________________________________________________________________________
228 void AliPHOSCalibrator::PlotPedestal(Int_t chanel)
230 //Plot histogram for a given channel, filled in Scan method
231 if(fPedHistos && fPedHistos->GetEntriesFast()){
232 static_cast<TH1F*>(fPedHistos->At(chanel))->Draw() ;
235 AliInfo(Form("Histograms not created yet! \n")) ;
238 //____________________________________________________________________________
239 void AliPHOSCalibrator::PlotPedestals(void)
241 // draws pedestals distribution
242 fhPedestals->Draw() ;
244 //____________________________________________________________________________
245 void AliPHOSCalibrator::PlotGain(Int_t chanel)
247 //Plot histogram for a given channel, filled in Scan method
248 if(fGainHistos && fGainHistos->GetEntriesFast()){
249 static_cast<TH1F*>(fGainHistos->At(chanel))->Draw() ;
252 AliInfo(Form("Histograms not created yet! \n")) ;
255 //____________________________________________________________________________
256 void AliPHOSCalibrator::PlotGains(void)
258 // draws gains distribution
261 //____________________________________________________________________________
262 void AliPHOSCalibrator::ScanPedestals(Option_t * option )
264 //scan all files in list fRunList and fill pedestal hisgrams
265 //option: "clear" - clear pedestal histograms filled up to now
266 // "deb" - plot file name currently processed
271 if(fPedHistos && strstr(option,"clear"))
272 fPedHistos->Delete() ;
274 fPedHistos = new TObjArray(fNch) ;
276 //Create histos for each channel, fills them and extracts mean values.
277 //First - prepare histos
279 for(ich=0;ich<fNch ;ich++){
280 TH1F * h = static_cast<TH1F *>(fPedHistos->At(ich)) ;
284 TString name("Pedestal for channel ") ;
286 fPedHistos->AddAt(new TH1F(n,name,fNChan,0,fNChan),ich) ;
290 TIter next(fRunList) ;
292 while((file = static_cast<TObjString *>(next()))){
293 if(strstr(option,"deb"))
294 printf("Processing file %s \n ",file->String().Data()) ;
297 AliRawReaderDate *rawReader = new AliRawReaderDate(file->String().Data()) ;
298 AliPHOSRawStream2004 *rawStream = new AliPHOSRawStream2004(rawReader) ;
299 rawStream->SetConTableDB(fctdb) ;
300 TClonesArray * digits = new TClonesArray("AliPHOSDigit",300) ;
302 //Scan all event in file
303 while(rawReader->NextEvent()){
304 //Is it PHYSICAL event
305 if(rawReader->GetType() == PHYSICS_EVENT){
307 if(rawStream->ReadDigits(digits)){
308 if(rawStream->IsPEDevent()){
309 for(Int_t idigit = 0; idigit<digits->GetEntriesFast() ; idigit++){
310 AliPHOSDigit * digit = static_cast<AliPHOSDigit *>(digits->At(idigit) ) ;
311 ich = fctdb->AbsId2Raw(digit->GetId());
313 Float_t amp = digit->GetAmp() ;
314 TH1F * hh = dynamic_cast<TH1F*>(fPedHistos->At(ich)) ;
322 if(strstr(option,"deb"))
323 AliInfo(Form(" found %d events \n ",nevents)) ;
329 //____________________________________________________________________________
330 void AliPHOSCalibrator::CalculatePedestals()
332 //Fit histograms, filled in ScanPedestals method with Gaussian
333 //find mean and width, check deviation from mean for each channel.
335 if(!fPedHistos || !fPedHistos->At(0)){
336 AliError(Form("You should run ScanPedestals first!")) ;
340 //Now fit results with Gauss
341 TF1 * gs = new TF1("gs","gaus",0.,10000.) ;
343 for(ich=0;ich<fNch ;ich++){
344 TH1F * h = static_cast<TH1F *>(fPedHistos->At(ich)) ;
345 Int_t max = h->GetMaximumBin() ;
346 Axis_t xmin = max/2. ;
347 Axis_t xmax = max*3/2 ;
348 gs->SetRange(xmin,xmax) ;
350 par[0] = h->GetBinContent(max) ;
353 gs->SetParameters(par[0],par[1],par[2]) ;
355 gs->GetParameters(par) ;
356 fhPedestals->SetBinContent(ich,par[1]) ;
357 fhPedestals->SetBinError(ich,par[2]) ;
358 fhPedestalsWid->Fill(ich,par[2]) ;
362 //now check reasonability of results
363 TF1 * p0 = new TF1("p0","pol0",0.,fNch) ;
364 fhPedestals->Fit("p0","Q") ;
366 p0->GetParameters(&meanPed);
367 for(ich=0;ich<fNch ;ich++){
368 Float_t ped = fhPedestals->GetBinContent(ich) ;
369 if(ped < 0 || ped > meanPed+fAcceptCorr){
370 TString out("Pedestal of channel ") ;
374 out+= "it is too far from mean " ;
376 AliError(Form("PHOSCalibrator %s",out.Data())) ;
382 //____________________________________________________________________________
383 void AliPHOSCalibrator::ScanGains(Option_t * option)
385 //Scan all runs, listed in fRunList and fill histograms for all channels
386 //options: "clear" - clean histograms, filled up to now
387 // "deb" - print current file name
388 // "narrow" - scan only narrow beam events
392 if(fGainHistos && strstr(option,"clear"))
393 fGainHistos->Delete() ;
395 if(strstr(option,"deball"))
396 AliInfo(Form("creating array for %d channels \n",fNch)) ;
397 fGainHistos = new TObjArray(fNch) ;
400 //Create histos for each channel, fills them and extracts mean values.
401 //First - prepare histos
403 if(!fGainHistos->GetEntriesFast()){
405 for(ich=0;ich<fNch ;ich++){
408 TString name("Gains for channel ") ;
410 fGainHistos->AddAt(new TH1F(n,name,fNGainBins,0,fGainMax),ich) ;
411 // static_cast<TH1F*>(fGainHistos->At(ich))->Sumw2() ;
415 TIter next(fRunList) ;
417 while((file = static_cast<TObjString *>(next()))){
419 AliRawReaderDate *rawReader = new AliRawReaderDate(file->String().Data()) ;
420 AliPHOSRawStream2004 *rawStream = new AliPHOSRawStream2004(rawReader) ;
421 rawStream->SetConTableDB(fctdb) ;
423 TClonesArray * digits = new TClonesArray("AliPHOSDigit",300) ;
425 //Scan all event in file
426 while(rawReader->NextEvent()){
427 //Is it PHYSICAL event
428 if(rawReader->GetType() == PHYSICS_EVENT){
429 if(rawStream->ReadDigits(digits)){
431 if(rawStream->IsNELevent() || rawStream->IsWELevent()){
433 AliPHOSDigit * digit ;
436 for(Int_t idigit = 0; idigit<digits->GetEntriesFast() ; idigit++){
437 digit = static_cast<AliPHOSDigit *>(digits->At(idigit) ) ;
438 if(digit->GetAmp() > max){
440 max = digit->GetAmp() ;
443 digit = static_cast<AliPHOSDigit *>(digits->At(imax) ) ;
444 Int_t ich = fctdb->AbsId2Raw(digit->GetId());
446 Float_t pedestal = fhPedestals->GetBinContent(ich) ;
447 const Float_t kshowerInCrystall = 0.9 ;
448 Float_t gain = fBeamEnergy*kshowerInCrystall/
449 (digit->GetAmp() - pedestal) ;
450 static_cast<TH1F*>(fGainHistos->At(ich))->Fill(gain) ;
459 if(strstr(option,"deb"))
460 AliInfo(Form(" found %d events \n",nevents)) ;
463 //____________________________________________________________________________
464 void AliPHOSCalibrator::CalculateGains(void)
468 if(!fGainHistos || !fGainHistos->GetEntriesFast()){
469 AliError(Form("You should run ScanGains first!")) ;
473 //Fit results with Landau
474 TF1 * gs = new TF1("gs","landau",0.,10000.) ;
476 for(ich=0;ich<fNch ;ich++){
477 TH1F * h = static_cast<TH1F *>(fGainHistos->At(ich)) ;
478 Int_t bmax = h->GetMaximumBin() ;
479 Axis_t center = h->GetBinCenter(bmax) ;
480 Axis_t xmin = center - 0.01 ;
481 Axis_t xmax = center + 0.02 ;
482 gs->SetRange(xmin,xmax) ;
484 par[0] = h->GetBinContent(bmax) ;
487 gs->SetParameters(par[0],par[1],par[2]) ;
489 gs->GetParameters(par) ;
490 fhGains->SetBinContent(ich,par[1]) ;
491 fhGains->SetBinError(ich,par[2]) ;
492 fhGainsWid->Fill(ich,par[2]) ;
496 //now check reasonability of results
497 TF1 * p0 = new TF1("p0","pol0",0.,fNch) ;
498 fhGains->Fit("p0","Q") ;
500 p0->GetParameters(&meanGain);
501 for(ich=0;ich<fNch ;ich++){
502 Float_t gain = fhGains->GetBinContent(ich) ;
503 if(gain < meanGain/fGainAcceptCorr || gain > meanGain*fGainAcceptCorr){
504 TString out("Gain of channel ") ;
508 out+= "it is too far from mean " ;
510 AliError(Form("PHOSCalibrator %s",out.Data())) ;
516 //____________________________________________________________________________
517 void AliPHOSCalibrator::ReadFromASCII(const char * filename){
518 // We read pedestals and gains from *.dat file with following format:
519 // 0 0 0 0 37.09 1972. // next nmodrows*nmodcols*ncryrows*ncrycols lines
520 // 0 0 0 1 28.53 2072. // contains <RR CC r c ped peak>
521 // 0 0 0 2 30.93 1938. //
522 // where module is an array of 8*8 crystals and RR and CC are module raw and column position
523 FILE * file = fopen(filename, "r");
525 Error("ReadFromASCII", "could not open file %s", filename);
528 if(!fctdb || !fhPedestals || !fhGains){
536 Int_t modRaw,modCol,raw,col;
539 while(fscanf(file,"%d %d %d %d %f %f",&modRaw,&modCol,&raw,&col,&ped,&pik)==6){
540 //Calculate plain crystal position:
541 Int_t rawPosition = (modRaw*8+raw)*fctdb->GetNColumns()+modCol*8+col ;
542 fhPedestals->SetBinContent(rawPosition,ped) ;
544 fhGains->SetBinContent(rawPosition,1./pik);
546 fhGains->SetBinContent(rawPosition,0.);
549 if(nread != fctdb->GetNColumns()*fctdb->GetNRaws()){
550 Error("ReadFromASCII","Read %d parameters instead of %d\n",nread,fctdb->GetNColumns()*fctdb->GetNRaws()) ;
554 //_____________________________________________________________________________
555 void AliPHOSCalibrator::WritePedestals(const char * version)
557 //Write calculated data to file using AliPHOSCalibrManager
558 //version and validitirange (begin-end) will be used to identify data
561 AliError(Form("\n Please, supply Connection Table DB (use SetConTableDB()) \n" )) ;
565 AliPHOSCalibrationData ped("Pedestals",version);
566 for(Int_t i=0; i<fNch;i++){
567 Int_t absid=fctdb->Raw2AbsId(i) ;
568 ped.SetData(absid,fhPedestals->GetBinContent(i)) ;
569 ped.SetDataCheck(absid,fhPedestalsWid->GetBinContent(i)) ;
572 // //evaluate validity range
574 // TIter next(fRunList) ;
575 // Int_t ibegin=99999;
577 // TObjString * file ;
578 // while((file=((TObjString*)next()))){
579 // TString s = file->GetString() ;
580 // TString ss = s(s.Last('_'),s.Last('.'));
582 // if(sscanf(ss.Data(),"%d",&tmp)){
589 // ped.SetValidityRange(ibegin,iend) ;
592 // ped.SetValidityRange(begin,end) ;
594 //check, may be Manager instance already configured?
595 AliPHOSCalibrManager * cmngr = AliPHOSCalibrManager::GetInstance() ;
597 AliWarning(Form("Using database file 'PHOSBTCalibration.root'")) ;
598 cmngr = AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
600 cmngr->WriteData(ped) ;
602 //_____________________________________________________________________________
603 void AliPHOSCalibrator::ReadPedestals(const char * version)
605 //Read data from file using AliPHOSCalibrManager
606 //version and range will be used to choose proper data
608 AliPHOSCalibrationData ped("Pedestals",version);
609 AliPHOSCalibrManager * cmngr = AliPHOSCalibrManager::GetInstance() ;
611 AliWarning(Form("Using database file 'PHOSBTCalibration.root'")) ;
612 cmngr = AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
614 cmngr->GetParameters(ped) ;
615 Int_t npeds=ped.NChannels() ;
616 fNch = fctdb->GetNchanels() ;
619 fhPedestals = new TH1F("hPedestals","Pedestals mean",fNch,0.,fNch) ;
620 for(Int_t i=0;i<npeds;i++){
621 Int_t raw =fctdb->AbsId2Raw(i) ;
623 fhPedestals->SetBinContent(raw-1,ped.Data(i)) ;
624 fhPedestals->SetBinError(raw-1,ped.DataCheck(i)) ;
628 //_____________________________________________________________________________
629 void AliPHOSCalibrator::ReadGains(const char * version)
631 //Read data from file using AliPHOSCalibrManager
632 //version and range will be used to choose proper data
634 AliPHOSCalibrationData gains("Gains",version);
635 AliPHOSCalibrManager * cmngr = AliPHOSCalibrManager::GetInstance() ;
637 AliWarning(Form("Using database file 'PHOSBTCalibration.root'")) ;
638 cmngr = AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
640 cmngr->GetParameters(gains) ;
641 Int_t npeds=gains.NChannels() ;
642 fNch = fctdb->GetNchanels() ;
645 fhGains = new TH1F("hGainss","Gains mean",fNch,0.,fNch) ;
646 for(Int_t i=0;i<npeds;i++){
647 Int_t raw =fctdb->AbsId2Raw(i) ;
649 fhGains->SetBinContent(raw-1,gains.Data(i)) ;
650 fhGains->SetBinError(raw-1,gains.DataCheck(i)) ;
654 //_____________________________________________________________________________
655 void AliPHOSCalibrator::WriteGains(const char * version)
657 //Write gains through AliPHOSCalibrManager
658 //version and validity range(begin-end) are used to identify data
661 AliError(Form("\n Please, supply Connection Table DB (use SetConTableDB()) \n" )) ;
665 AliPHOSCalibrationData gains("Gains",version);
666 for(Int_t i=0; i<fNch;i++){
667 Int_t absid=fctdb->Raw2AbsId(i) ;
668 gains.SetData(absid,fhGains->GetBinContent(i)) ;
669 gains.SetDataCheck(absid,fhGainsWid->GetBinContent(i)) ;
672 // TIter next(fRunList) ;
673 // Int_t ibegin=99999;
675 // TObjString * file ;
676 // while((file=((TObjString*)next()))){
677 // TString s = file->GetString() ;
678 // TSubString ss = s(s.Last('_'),s.Last('.'));
680 // if(sscanf(ss.Data(),"%d",&tmp)){
687 // gains.SetValidityRange(ibegin,iend) ;
690 // gains.SetValidityRange(begin,end) ;
691 AliPHOSCalibrManager * cmngr = AliPHOSCalibrManager::GetInstance() ;
693 AliWarning(Form("Using database file 'PHOSBTCalibration.root'")) ;
694 cmngr = AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
696 cmngr->WriteData(gains) ;
698 //_____________________________________________________________________________
699 void AliPHOSCalibrator::Print(const Option_t *)const
702 AliInfo(Form("--------------PHOS Calibrator-----------------\n")) ;
703 printf("Files to handle:\n") ;
704 TIter next(fRunList) ;
706 while((r=(TObjString *)(next())))
707 printf(" %s\n",r->GetName()) ;
709 printf("Name of ConTableDB:.....................%s\n",fConTableDB.Data()) ;
710 printf("File of ConTableDB:.....................%s\n",fConTableDBFile.Data() ) ;
711 printf("Maximal deviation from mean Gain (factor):.%f\n",fGainAcceptCorr) ;
712 printf("Maximal deviation of Pedestal from mean:...%f\n",fAcceptCorr) ;
713 printf("Range used in Gain histos:..............%f\n",fGainMax) ;
714 printf("Number of bins in Gain histos:..........%d\n",fNGainBins) ;
715 printf("Number of channels to calibrate:........%d\n",fNch) ;
716 printf("Number of bins in pedestal histos:......%d\n",fNChan) ;
717 printf("--------------------------------------------------\n") ;