9777b85957865892ebd2a856438f9acf298d2ab7
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRawUtils.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 /* History of cvs commits:
18  *
19  * $Log$ */
20
21 //*-- Author: Marco van Leeuwen (LBL)
22 #include "AliEMCALRawUtils.h"
23
24 #include "TF1.h"
25 #include "TGraph.h"
26 #include "TSystem.h"
27
28 #include "AliLog.h"
29 #include "AliRunLoader.h"
30 #include "AliCaloAltroMapping.h"
31 #include "AliAltroBuffer.h"
32 #include "AliRawReader.h"
33 #include "AliCaloRawStream.h"
34 #include "AliDAQ.h"
35
36 #include "AliEMCALLoader.h"
37 #include "AliEMCALGeometry.h"
38 #include "AliEMCALDigitizer.h"
39 #include "AliEMCALDigit.h"
40
41
42 ClassImp(AliEMCALRawUtils)
43
44 // Signal shape parameters
45 Int_t    AliEMCALRawUtils::fgOrder       = 2 ;      // Order of gamma function 
46 Double_t AliEMCALRawUtils::fgTimeMax     = 2.56E-5 ;  // each sample is over 100 ns fTimeMax/fTimeBins
47 Double_t AliEMCALRawUtils::fgTau         = 165E-9 ;   // 165 ns (from testbeam; not very accurate)
48 Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ;   // 15 time bins ~ 1.5 musec
49
50 // some digitization constants
51 Int_t    AliEMCALRawUtils::fgThreshold = 1;
52 Int_t    AliEMCALRawUtils::fgDDLPerSuperModule = 2;  // 2 ddls per SuperModule
53
54 AliEMCALRawUtils::AliEMCALRawUtils(): fHighLowGainFactor(0.) {
55   fHighLowGainFactor = 16. ;          // adjusted for a low gain range of 82 GeV (10 bits) 
56 }
57 //____________________________________________________________________________
58 AliEMCALRawUtils::~AliEMCALRawUtils() {
59 }
60 //____________________________________________________________________________
61 void AliEMCALRawUtils::Digits2Raw()
62 {
63   // convert digits of the current event to raw data
64   
65   AliRunLoader *rl = AliRunLoader::GetRunLoader();
66   AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
67
68   // get the digits
69   loader->LoadDigits("EMCAL");
70   loader->GetEvent();
71   TClonesArray* digits = loader->Digits() ;
72   
73   if (!digits) {
74     Warning("Digits2Raw", "no digits found !");
75     return;
76   }
77     
78   // get the geometry
79   AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
80   if (!geom) {
81     AliError(Form("No geometry found !"));
82     return;
83   }
84   
85   static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here
86   AliAltroBuffer* buffers[nDDL];
87   for (Int_t i=0; i < nDDL; i++)
88     buffers[i] = 0;
89
90   Int_t adcValuesLow[fgkTimeBins];
91   Int_t adcValuesHigh[fgkTimeBins];
92
93   //Load Mapping RCU files once
94   TString path = gSystem->Getenv("ALICE_ROOT");
95   path += "/EMCAL/mapping/RCU";
96   TString path0 = path+"0.data";//This file will change in future
97   TString path1 = path+"1.data";//This file will change in future
98   AliAltroMapping * mapping[2] ; // For the moment only 2
99   mapping[0] = new AliCaloAltroMapping(path0.Data());
100   mapping[1] = new AliCaloAltroMapping(path1.Data());
101
102   // loop over digits (assume ordered digits)
103   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
104     AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
105     if (digit->GetAmp() < fgThreshold) 
106       continue;
107
108     //get cell indices
109     Int_t nSM = 0;
110     Int_t nIphi = 0;
111     Int_t nIeta = 0;
112     Int_t iphi = 0;
113     Int_t ieta = 0;
114     Int_t nModule = 0;
115     geom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
116     geom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
117     
118     //Check which is the RCU of the cell.
119     Int_t iRCU = -111;
120     //RCU0
121     if (0<=iphi&&iphi<8) iRCU=0; // first cable row
122     else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; 
123     //second cable row
124     //RCU1
125     else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; 
126     //second cable row
127     else if(16<=iphi&&iphi<24) iRCU=1; // third cable row
128     
129     //Which DDL?
130     Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU;
131     if (iDDL >= nDDL)
132       Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL);
133
134     if (buffers[iDDL] == 0) {      
135       // open new file and write dummy header
136       TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
137       buffers[iDDL] = new AliAltroBuffer(fileName.Data(),mapping[iRCU]);
138       buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
139     }
140     
141     // out of time range signal (?)
142     if (digit->GetTimeR() > GetRawFormatTimeMax() ) {
143       AliInfo("Signal is out of time range.\n");
144       buffers[iDDL]->FillBuffer((Int_t)digit->GetAmp());
145       buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() );  // time bin
146       buffers[iDDL]->FillBuffer(3);          // bunch length      
147       buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM);  // trailer
148       // calculate the time response function
149     } else {
150       Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh, adcValuesLow) ; 
151       if (lowgain) 
152         buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow, fgThreshold);
153       else 
154         buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh, fgThreshold);
155     }
156   }
157   
158   // write headers and close files
159   for (Int_t i=0; i < nDDL; i++) {
160     if (buffers[i]) {
161       buffers[i]->Flush();
162       buffers[i]->WriteDataHeader(kFALSE, kFALSE);
163       delete buffers[i];
164     }
165   }
166   mapping[0]->Delete();
167   mapping[1]->Delete();
168   loader->UnloadDigits();
169 }
170
171 //____________________________________________________________________________
172 void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader)
173 {
174   // convert raw data of the current event to digits
175   AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
176   if (!geom) {
177     AliError(Form("No geometry found !"));
178     return;
179   }
180
181   AliRunLoader *rl = AliRunLoader::GetRunLoader();
182   AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
183
184   // get the digits
185   TClonesArray* digits = loader->Digits() ;
186   digits->Clear(); 
187
188   if (!digits) {
189     Error("Raw2Digits", "no digits found !");
190     return;
191   }
192   if (!reader) {
193     Error("Raw2Digits", "no raw reader found !");
194     return;
195   }
196
197   // Use AliAltroRawStream to read the ALTRO format.  No need to
198   // reinvent the wheel :-) 
199   AliCaloRawStream in(reader,"EMCAL");
200   // Select EMCAL DDL's;
201   reader->Select("EMCAL");
202
203   // reading is from previously existing AliEMCALGetter.cxx
204   // ReadRaw method
205   TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4);
206   signalF->SetParNames("Charge", "Gain", "Amplitude", "TimeZero"); 
207   
208   Int_t id =  -1;
209   Int_t idigit = 0 ; 
210   Double_t time = 0. ; 
211   Double_t amp = 0. ; 
212
213   TGraph * gSig = new TGraph(GetRawFormatTimeBins()) ; 
214
215   Int_t eofReached = 0;
216   Int_t lowGain = 0;
217
218   in.Next(); // Go to first digit
219   do {
220     id =  geom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
221     lowGain = in.IsLowGain();
222     gSig->SetPoint(in.GetTime(), 
223                    in.GetTime()* GetRawFormatTimeMax() / GetRawFormatTimeBins(), 
224                    in.GetSignal()) ;
225
226     Int_t iTime = 1;
227     do {
228       if (!in.Next())
229         eofReached = 1;
230       else {
231         gSig->SetPoint(in.GetTime(), 
232                        in.GetTime()* GetRawFormatTimeMax() / GetRawFormatTimeBins(), 
233                        in.GetSignal()) ;
234       }
235       iTime++;
236     } while (!eofReached && !in.IsNewRow() && !in.IsNewColumn() && !in.IsNewModule());
237
238     FitRaw(gSig, signalF, amp, time) ; 
239     if (lowGain) 
240       amp *= fHighLowGainFactor; 
241     
242     if (amp > 0) {
243       AliDebug(2,Form("id %d amp %g", id, amp));
244       new((*digits)[idigit]) AliEMCALDigit( -1, -1, id, (Int_t)amp, time, idigit) ;     
245           idigit++ ; 
246     }
247     Int_t index ; 
248         
249     // Reset graph
250     for (index = 0; index < GetRawFormatTimeBins(); index++) {
251       gSig->SetPoint(index, index * GetRawFormatTimeMax() / GetRawFormatTimeBins(), 0) ;  
252     } 
253   } while (!eofReached); // EMCAL entries loop
254   
255   delete signalF ; 
256   delete gSig;
257   
258   return ; 
259 }
260
261 //____________________________________________________________________________ 
262 void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Double_t & amp, Double_t & time)
263 {
264   // Fits the raw signal time distribution; from AliEMCALGetter 
265
266   const Int_t kNoiseThreshold = 0 ;
267   Double_t timezero1 = 0., timezero2 = 0., timemax = 0. ;
268   Double_t signal = 0., signalmax = 0. ;       
269   amp = time = 0. ; 
270
271   timezero1 = signalmax = timemax = 0. ;
272   timezero2 = GetRawFormatTimeMax();
273   Int_t index ; 
274   for (index = 0; index < GetRawFormatTimeBins(); index++) {
275     gSig->GetPoint(index, time, signal) ; 
276     if (signal > kNoiseThreshold && timezero1 == 0.) 
277       timezero1 = time ;
278     if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
279       timezero2 = time ; 
280     if (signal > signalmax) {
281       signalmax = signal ; 
282       timemax   = time ; 
283     }
284   }
285
286   if ( timemax < GetRawFormatTimeMax() * 0.4 ) { // else its noise 
287     signalF->SetParameter(0, signalmax) ; 
288     signalF->SetParameter(1, timemax) ; 
289     gSig->Fit(signalF, "QRON", "", 0., timezero2); //, "QRON") ; 
290     amp = signalF->GetParameter(0) ; 
291     time = signalF->GetParameter(1) - fgTimeTrigger;
292   }
293   return;
294 }
295 //__________________________________________________________________
296 Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
297 {
298   // Shape of the electronics raw reponse:
299   // It is a semi-gaussian, 2nd order Gamma function of the general form
300   //
301   // t' = (t - t0 + tau) / tau
302   // F = A * t**N * exp( N * ( 1 - t) )   for t >= 0
303   // F = 0                                for t < 0 
304   //
305   // parameters:
306   // A:   par[0]   // Amplitude = peak value
307   // t0:  par[1]
308   // tau: fgTau
309   // N:   fgOrder
310   //
311   Double_t signal ;
312   Double_t xx = ( x[0] - par[1] + fgTau ) / fgTau ; 
313
314   if (xx < 0) 
315     signal = 0. ;  
316   else {  
317     signal = par[0] * TMath::Power(xx , fgOrder) * TMath::Exp(fgOrder * (1 - xx )) ; 
318   }
319   return signal ;  
320 }
321
322 //__________________________________________________________________
323 Bool_t AliEMCALRawUtils::RawSampledResponse(
324 const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const 
325 {
326   // for a start time dtime and an amplitude damp given by digit, 
327   // calculates the raw sampled response AliEMCAL::RawResponseFunction
328
329   const Int_t kRawSignalOverflow = 0x3FF ; 
330   Bool_t lowGain = kFALSE ; 
331
332   TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4);
333   signalF.SetParameter(0, damp) ; 
334   signalF.SetParameter(1, dtime + fgTimeTrigger) ; 
335
336   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
337     Double_t time = iTime * GetRawFormatTimeMax() / GetRawFormatTimeBins() ;
338     Double_t signal = signalF.Eval(time) ;     
339     adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
340     if ( adcH[iTime] > kRawSignalOverflow ){  // larger than 10 bits 
341       adcH[iTime] = kRawSignalOverflow ;
342       lowGain = kTRUE ; 
343     }
344
345     signal /= fHighLowGainFactor;
346
347     adcL[iTime] =  static_cast<Int_t>(signal + 0.5) ;
348     if ( adcL[iTime] > kRawSignalOverflow)  // larger than 10 bits 
349       adcL[iTime] = kRawSignalOverflow ;
350   }
351   return lowGain ; 
352 }