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