[u/mrichter/AliRoot.git] / FMD / scripts / DrawCalibRaw.C

//____________________________________________________________________
//
// $Id: DrawDigits.C 30718 2009-01-22 16:07:40Z cholm $
//
// Script that contains a class to draw eloss from hits, versus ADC
// counts from digits, using the AliFMDInputHits class in the util library. 
//
// It draws the energy loss versus the p/(mq^2).  It can be overlayed
// with the Bethe-Bloc curve to show how the simulation behaves
// relative to the expected. 
//
// Use the script `Compile.C' to compile this class using ACLic. 
//
#include <TH1D.h>
#include <AliCDBManager.h>
#include <AliFMDHit.h>
#include <AliFMDDigit.h>
#include <AliFMDInput.h>
#include <AliFMDEdepMap.h>
#include <AliFMDParameters.h>
#include <AliFMDCalibStripRange.h>
#include <AliFMDCalibSampleRate.h>
#include <AliFMDCalibGain.h>
#include <AliFMDCalibPedestal.h>
#include <AliFMDRawReader.h>
#include <iostream>
#include <fstream>
#include <TStyle.h>
#include <TArrayF.h>
#include <AliLog.h>
#include <TSystem.h>
#include <TFile.h>

/** @class DrawDigits
    @brief Draw hit energy loss versus digit ADC
    @code 
    Root> .L Compile.C
    Root> Compile("DrawCalibRaw.C")
    Root> DrawCalibRaw cr;
    Root> cr.SetRawFile("file");
    Root> cr.Run();
    @endcode
    @ingroup FMD_script
 */
class DrawCalibRaw : public AliFMDInput
{
private:
  Double_t   fFactor;
  TString    fCalibDir;
  TH1D*      fELoss; // Histogram 
  TObjArray* fDets;
  TFile*     fOut;
  Bool_t     fHasData;
  Int_t      fGotNEvents;

public:
  //__________________________________________________________________
  DrawCalibRaw(const       char*  file, 
	       const char* calibDir    = 0,
	       Double_t    noiseFactor = 5, 
	       Bool_t      save        = kTRUE,
	       Int_t       m           = 420, 
	       Double_t    mmin        = -0.5, 
	       Double_t    mmax        = 20.5) 
    : AliFMDInput(), 
      fFactor(noiseFactor),
      fCalibDir(""), 
      fELoss(0), 
      fDets(0), 
      fOut(0), 
      fHasData(kFALSE),
      fGotNEvents(0)
  { 
    if (calibDir) fCalibDir = calibDir;
    AddLoad(kRawCalib);
    fELoss = new TH1D("eLoss", "Scaled Energy loss", m, mmin, mmax);
    fELoss->SetXTitle("#Delta E/#Delta E_{mip}");
    fELoss->Sumw2();
    fELoss->SetDirectory(0);
    
    SetRawFile(file);

    if (!save) return;
    fDets = new TObjArray(3);

  }
  //__________________________________________________________________
  Bool_t CheckFile(const char* prefix, int number, TString& f)
  {
    f = (Form("%s%d.csv", prefix, number));
    f = gSystem->Which(fCalibDir.Data(), f.Data());
    return !f.IsNull();
  }

  //__________________________________________________________________
  Bool_t Init()
  {
    AliCDBManager* cdb = AliCDBManager::Instance();
    cdb->SetRun(0);
    cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");

    AliFMDCalibStripRange* range = 0;
    AliFMDCalibSampleRate* rate  = 0;
    AliFMDCalibPedestal*   peds  = 0;
    AliFMDCalibGain*       gains = 0;
    for (Int_t i = 1; i <= 3; i++) { 
      TString f;
      if (CheckFile("conditions", i, f)) {
	Info("Init", "Reading conditions for FMD%d from %s", i, f.Data());
	if (!range) range = new AliFMDCalibStripRange;
	if (!rate)  rate  = new AliFMDCalibSampleRate;
	std::ifstream in(f.Data());
	range->ReadFromFile(in);
	rate->ReadFromFile(in);
      }
      if (CheckFile("peds", i, f)) {
	Info("Init", "Reading pedestals for FMD%d from %s", i, f.Data());
	if (!peds) peds = new AliFMDCalibPedestal;
	std::ifstream in(f.Data());
	peds->ReadFromFile(in);
      }
      if (CheckFile("gains", i, f)) {
	Info("Init", "Reading gains for FMD%d from %s", i, f.Data());
	if (!gains) gains = new AliFMDCalibGain;
	std::ifstream in(f.Data());
	gains->ReadFromFile(in);
      }
    }

    Int_t mask = (AliFMDParameters::kDeadMap|
		  AliFMDParameters::kZeroSuppression|
		  AliFMDParameters::kAltroMap);

    if (!range) mask |= AliFMDParameters::kStripRange;
    if (!rate)  mask |= AliFMDParameters::kSampleRate;
    if (!peds)  mask |= AliFMDParameters::kPedestal;
    if (!gains) mask |= AliFMDParameters::kPulseGain;

    AliFMDParameters* pars = AliFMDParameters::Instance();
    pars->Init(kFALSE, mask);

    if (range)  pars->SetStripRange(range);
    if (rate)   pars->SetSampleRate(rate);
    if (peds)   pars->SetPedestal(peds);
    if (gains)  pars->SetGain(gains);
    
    Bool_t ret = AliFMDInput::Init();
    
    if (!fDets) return ret;

    fOut  = TFile::Open(Form("histo_%s", fRawFile.Data()), "RECREATE");

    Int_t m    = fELoss->GetXaxis()->GetNbins();
    Int_t mmin = fELoss->GetXaxis()->GetXmin();
    Int_t mmax = fELoss->GetXaxis()->GetXmax();
    for (Int_t d = 1; d <= 3; d++) {
      Int_t      nRng = (d == 1 ? 1 : 2);
      TObjArray* det  = 0;
      fDets->AddAt(det = new TObjArray(nRng), d-1);
      det->SetName(Form("FMD%d", d));
      TDirectory* detD = fOut->mkdir(det->GetName());
      for (Int_t q = 0; q < nRng; q++) { 
	Char_t r       = q == 0 ? 'I' : 'O';
	Int_t  nSec    = q == 0 ?  20 :  40;
	Int_t  nStr    = q == 0 ? 512 : 256;
	TObjArray* rng = 0;
	det->AddAt(rng = new TObjArray(nSec), q);
	rng->SetName(Form("FMD%d%c", d, r));
	TDirectory* rngD = detD->mkdir(rng->GetName());
	for (Int_t s = 0; s < nSec; s++) { 
	  TObjArray* sec = 0;
	  rng->AddAt(sec = new TObjArray(nStr), s);
	  sec->SetName(Form("FMD%d%c_%02d", d, r, s));
	  TDirectory* secD = rngD->mkdir(sec->GetName());
	  for (Int_t t = 0; t < nStr; t++) { 
	    secD->cd();
	    TH1* str = new TH1D(Form("FMD%d%c_%02d_%03d", d, r, s, t), 
				 Form("Scaled energy loss in FMD%d%c[%2d,%3d]",
				      d, r, s, t), m, mmin, mmax);
	    str->SetXTitle("#Delta E/#Delta E_{mip}");
	    // str->SetDirectory(secD);
	    sec->AddAt(str, t);
	  }
	}
      }
    }

    return ret;
  }

  //__________________________________________________________________
  Bool_t Begin(Int_t e)
  {
    fHasData = kFALSE;
    return AliFMDInput::Begin(e);
  }

  //__________________________________________________________________
  Bool_t ProcessRawCalibDigit(AliFMDDigit* digit)
  {
    if (!digit) return kTRUE;

    AliFMDParameters* parm = AliFMDParameters::Instance();
    UShort_t d             =  digit->Detector();
    Char_t   r             =  digit->Ring();
    UShort_t s             =  digit->Sector();
    UShort_t t             =  digit->Strip();
    Double_t gain          =  parm->GetPulseGain(d, r, s, t);
    Double_t ped           =  parm->GetPedestal(d, r, s, t);
    Double_t pedW          =  parm->GetPedestalWidth(d, r, s, t);
    Double_t adc           =  digit->Counts();
    Double_t threshold     =  pedW * fFactor;
    if (gain < 0.1 || gain > 10) return kTRUE;
    if (pedW > 10) { 
      Warning("ProcessRawCalibDigit", "FMD%d%c[%2d,%3d] is noisy: %f",
	      d, r, s, t, pedW);
      return kTRUE;
    }

    if (fFMDReader && fFMDReader->IsZeroSuppressed(d-1))
      adc += fFMDReader->NoiseFactor(d-1) * pedW;
    else 
      threshold            += ped;

    if (adc < threshold) return kTRUE;
    
    Double_t mult = (adc-ped) / (gain * parm->GetDACPerMIP());

    fHasData = kTRUE;
    fELoss->Fill(mult);

    // if (t >= 10) return kTRUE;
    TObjArray* det = static_cast<TObjArray*>(fDets->At(d-1));
    TObjArray* rng = static_cast<TObjArray*>(det->At(r == 'I' ? 0 : 1));
    TObjArray* sec = static_cast<TObjArray*>(rng->At(s));
    TH1*       str = static_cast<TH1*>(sec->At(t));
    str->Fill(mult);

    return kTRUE;
  }
  //__________________________________________________________________
  Bool_t End()
  {
    if (fHasData) fGotNEvents++;
    return AliFMDInput::End();
  }
  //__________________________________________________________________
  Bool_t Finish()
  {
    std::cout << "A total of " << fGotNEvents << " with data" << std::endl;
    gStyle->SetPalette(1);
    gStyle->SetOptTitle(0);
    gStyle->SetCanvasColor(0);
    gStyle->SetCanvasBorderSize(0);
    gStyle->SetPadColor(0);
    gStyle->SetPadBorderSize(0);
    fELoss->SetStats(kFALSE);
    fELoss->SetFillColor(kRed);
    fELoss->SetFillStyle(3001);
    fELoss->Scale(1. / fELoss->GetEntries());
    fELoss->DrawCopy("e1 bar");

    if (fDets && fOut) { 
      std::cout << "Flusing to disk ... " << std::flush;
      fOut->cd();
      fELoss->Write();
      fOut->Write();
      fOut->Close();
      std::cout << "done" << std::endl;
    }
    return kTRUE;
  }

  ClassDef(DrawCalibRaw,0);
};

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
818fff8d	1	//____________________________________________________________________
	2	//
	3	// $Id: DrawDigits.C 30718 2009-01-22 16:07:40Z cholm $
	4	//
	5	// Script that contains a class to draw eloss from hits, versus ADC
	6	// counts from digits, using the AliFMDInputHits class in the util library.
	7	//
	8	// It draws the energy loss versus the p/(mq^2). It can be overlayed
	9	// with the Bethe-Bloc curve to show how the simulation behaves
	10	// relative to the expected.
	11	//
	12	// Use the script `Compile.C' to compile this class using ACLic.
	13	//
	14	#include <TH1D.h>
	15	#include <AliCDBManager.h>
	16	#include <AliFMDHit.h>
	17	#include <AliFMDDigit.h>
	18	#include <AliFMDInput.h>
	19	#include <AliFMDEdepMap.h>
	20	#include <AliFMDParameters.h>
	21	#include <AliFMDCalibStripRange.h>
	22	#include <AliFMDCalibSampleRate.h>
	23	#include <AliFMDCalibGain.h>
	24	#include <AliFMDCalibPedestal.h>
	25	#include <AliFMDRawReader.h>
	26	#include <iostream>
	27	#include <fstream>
	28	#include <TStyle.h>
	29	#include <TArrayF.h>
	30	#include <AliLog.h>
	31	#include <TSystem.h>
	32	#include <TFile.h>
	33
	34	/** @class DrawDigits
	35	@brief Draw hit energy loss versus digit ADC
	36	@code
	37	Root> .L Compile.C
	38	Root> Compile("DrawCalibRaw.C")
	39	Root> DrawCalibRaw cr;
	40	Root> cr.SetRawFile("file");
	41	Root> cr.Run();
	42	@endcode
	43	@ingroup FMD_script
	44	*/
	45	class DrawCalibRaw : public AliFMDInput
	46	{
	47	private:
	48	Double_t fFactor;
	49	TString fCalibDir;
	50	TH1D* fELoss; // Histogram
	51	TObjArray* fDets;
	52	TFile* fOut;
	53	Bool_t fHasData;
	54	Int_t fGotNEvents;
	55
	56	public:
	57	//__________________________________________________________________
	58	DrawCalibRaw(const char* file,
	59	const char* calibDir = 0,
	60	Double_t noiseFactor = 5,
	61	Bool_t save = kTRUE,
	62	Int_t m = 420,
	63	Double_t mmin = -0.5,
	64	Double_t mmax = 20.5)
65	: AliFMDInput(),
66	fFactor(noiseFactor),
67	fCalibDir(""),
68	fELoss(0),
69	fDets(0),
70	fOut(0),
71	fHasData(kFALSE),
72	fGotNEvents(0)
73	{
74	if (calibDir) fCalibDir = calibDir;
75	AddLoad(kRawCalib);
76	fELoss = new TH1D("eLoss", "Scaled Energy loss", m, mmin, mmax);
77	fELoss->SetXTitle("#Delta E/#Delta E_{mip}");
78	fELoss->Sumw2();
79	fELoss->SetDirectory(0);
80
81	SetRawFile(file);
82
83	if (!save) return;
84	fDets = new TObjArray(3);
85
86	}
87	//__________________________________________________________________
88	Bool_t CheckFile(const char* prefix, int number, TString& f)
89	{
90	f = (Form("%s%d.csv", prefix, number));
91	f = gSystem->Which(fCalibDir.Data(), f.Data());
92	return !f.IsNull();
93	}
94
95	//__________________________________________________________________
96	Bool_t Init()
97	{
98	AliCDBManager* cdb = AliCDBManager::Instance();
99	cdb->SetRun(0);
100	cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
101
102	AliFMDCalibStripRange* range = 0;
103	AliFMDCalibSampleRate* rate = 0;
104	AliFMDCalibPedestal* peds = 0;
105	AliFMDCalibGain* gains = 0;
106	for (Int_t i = 1; i <= 3; i++) {
107	TString f;
108	if (CheckFile("conditions", i, f)) {
109	Info("Init", "Reading conditions for FMD%d from %s", i, f.Data());
110	if (!range) range = new AliFMDCalibStripRange;
111	if (!rate) rate = new AliFMDCalibSampleRate;
112	std::ifstream in(f.Data());
113	range->ReadFromFile(in);
114	rate->ReadFromFile(in);
115	}
116	if (CheckFile("peds", i, f)) {
117	Info("Init", "Reading pedestals for FMD%d from %s", i, f.Data());
118	if (!peds) peds = new AliFMDCalibPedestal;
119	std::ifstream in(f.Data());
120	peds->ReadFromFile(in);
121	}
122	if (CheckFile("gains", i, f)) {
123	Info("Init", "Reading gains for FMD%d from %s", i, f.Data());
124	if (!gains) gains = new AliFMDCalibGain;
125	std::ifstream in(f.Data());
126	gains->ReadFromFile(in);
127	}
128	}
129
130	Int_t mask = (AliFMDParameters::kDeadMap\|
131	AliFMDParameters::kZeroSuppression\|
132	AliFMDParameters::kAltroMap);
133
134	if (!range) mask \|= AliFMDParameters::kStripRange;
135	if (!rate) mask \|= AliFMDParameters::kSampleRate;
136	if (!peds) mask \|= AliFMDParameters::kPedestal;
137	if (!gains) mask \|= AliFMDParameters::kPulseGain;
138
139	AliFMDParameters* pars = AliFMDParameters::Instance();
140	pars->Init(kFALSE, mask);
141
142	if (range) pars->SetStripRange(range);
143	if (rate) pars->SetSampleRate(rate);
144	if (peds) pars->SetPedestal(peds);
145	if (gains) pars->SetGain(gains);
146
147	Bool_t ret = AliFMDInput::Init();
148
149	if (!fDets) return ret;
150
151	fOut = TFile::Open(Form("histo_%s", fRawFile.Data()), "RECREATE");
152
153	Int_t m = fELoss->GetXaxis()->GetNbins();
154	Int_t mmin = fELoss->GetXaxis()->GetXmin();
155	Int_t mmax = fELoss->GetXaxis()->GetXmax();
156	for (Int_t d = 1; d <= 3; d++) {
157	Int_t nRng = (d == 1 ? 1 : 2);
158	TObjArray* det = 0;
159	fDets->AddAt(det = new TObjArray(nRng), d-1);
160	det->SetName(Form("FMD%d", d));
161	TDirectory* detD = fOut->mkdir(det->GetName());
162	for (Int_t q = 0; q < nRng; q++) {
163	Char_t r = q == 0 ? 'I' : 'O';
164	Int_t nSec = q == 0 ? 20 : 40;
165	Int_t nStr = q == 0 ? 512 : 256;
166	TObjArray* rng = 0;
167	det->AddAt(rng = new TObjArray(nSec), q);
168	rng->SetName(Form("FMD%d%c", d, r));
169	TDirectory* rngD = detD->mkdir(rng->GetName());
170	for (Int_t s = 0; s < nSec; s++) {
171	TObjArray* sec = 0;
172	rng->AddAt(sec = new TObjArray(nStr), s);
173	sec->SetName(Form("FMD%d%c_%02d", d, r, s));
174	TDirectory* secD = rngD->mkdir(sec->GetName());
175	for (Int_t t = 0; t < nStr; t++) {
176	secD->cd();
177	TH1* str = new TH1D(Form("FMD%d%c_%02d_%03d", d, r, s, t),
178	Form("Scaled energy loss in FMD%d%c[%2d,%3d]",
179	d, r, s, t), m, mmin, mmax);
180	str->SetXTitle("#Delta E/#Delta E_{mip}");
181	// str->SetDirectory(secD);
182	sec->AddAt(str, t);
183	}
184	}
185	}
186	}
187
188	return ret;
189	}
190
191	//__________________________________________________________________
192	Bool_t Begin(Int_t e)
193	{
194	fHasData = kFALSE;
195	return AliFMDInput::Begin(e);
196	}
197
198	//__________________________________________________________________
199	Bool_t ProcessRawCalibDigit(AliFMDDigit* digit)
200	{
201	if (!digit) return kTRUE;
202
203	AliFMDParameters* parm = AliFMDParameters::Instance();
204	UShort_t d = digit->Detector();
205	Char_t r = digit->Ring();
206	UShort_t s = digit->Sector();
207	UShort_t t = digit->Strip();
208	Double_t gain = parm->GetPulseGain(d, r, s, t);
209	Double_t ped = parm->GetPedestal(d, r, s, t);
210	Double_t pedW = parm->GetPedestalWidth(d, r, s, t);
211	Double_t adc = digit->Counts();
212	Double_t threshold = pedW * fFactor;
213	if (gain < 0.1 \|\| gain > 10) return kTRUE;
214	if (pedW > 10) {
215	Warning("ProcessRawCalibDigit", "FMD%d%c[%2d,%3d] is noisy: %f",
216	d, r, s, t, pedW);
217	return kTRUE;
218	}
219
220	if (fFMDReader && fFMDReader->IsZeroSuppressed(d-1))
221	adc += fFMDReader->NoiseFactor(d-1) * pedW;
222	else
223	threshold += ped;
224
225	if (adc < threshold) return kTRUE;
226
227	Double_t mult = (adc-ped) / (gain * parm->GetDACPerMIP());
228
229	fHasData = kTRUE;
230	fELoss->Fill(mult);
231
232	// if (t >= 10) return kTRUE;
233	TObjArray* det = static_cast<TObjArray*>(fDets->At(d-1));
234	TObjArray* rng = static_cast<TObjArray*>(det->At(r == 'I' ? 0 : 1));
235	TObjArray* sec = static_cast<TObjArray*>(rng->At(s));
236	TH1* str = static_cast<TH1*>(sec->At(t));
237	str->Fill(mult);
238
239	return kTRUE;
240	}
241	//__________________________________________________________________
242	Bool_t End()
243	{
244	if (fHasData) fGotNEvents++;
245	return AliFMDInput::End();
246	}
247	//__________________________________________________________________
248	Bool_t Finish()
249	{
250	std::cout << "A total of " << fGotNEvents << " with data" << std::endl;
251	gStyle->SetPalette(1);
252	gStyle->SetOptTitle(0);
253	gStyle->SetCanvasColor(0);
254	gStyle->SetCanvasBorderSize(0);
255	gStyle->SetPadColor(0);
256	gStyle->SetPadBorderSize(0);
257	fELoss->SetStats(kFALSE);
258	fELoss->SetFillColor(kRed);
259	fELoss->SetFillStyle(3001);
260	fELoss->Scale(1. / fELoss->GetEntries());
261	fELoss->DrawCopy("e1 bar");
262
263	if (fDets && fOut) {
264	std::cout << "Flusing to disk ... " << std::flush;
265	fOut->cd();
266	fELoss->Write();
267	fOut->Write();
268	fOut->Close();
269	std::cout << "done" << std::endl;
270	}
271	return kTRUE;
272	}
273
274	ClassDef(DrawCalibRaw,0);
275	};
276
277	//____________________________________________________________________
278	//
279	// EOF
280	//