Increased error checking and possibility to extract some diagnostics
[u/mrichter/AliRoot.git] / FMD / AliFMDDigitizer.cxx
1 //************************************************************************
2 // Copyright(c) 2004, 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 // $Id$ */
16 /**
17  * @file    AliFMDDigitizer.cxx
18  * 
19  * @author  Christian Holm Christensen <cholm@nbi.dk>
20  * @date    Mon Mar 27 12:38:26 2006
21  * @brief   FMD Digitizers implementation
22  *
23  * @ingroup FMD_sim
24  */
25 //////////////////////////////////////////////////////////////////////////////
26 //
27 //  This class contains the procedures simulation ADC  signal for the
28 //  Forward Multiplicity detector  : SDigits->Digits
29 // 
30 //  Digits consists of
31 //   - Detector #
32 //   - Ring ID                                             
33 //   - Sector #     
34 //   - Strip #
35 //   - ADC count in this channel                                  
36 //
37 //  Digits consists of
38 //   - Detector #
39 //   - Ring ID                                             
40 //   - Sector #     
41 //   - Strip #
42 //   - Total energy deposited in the strip
43 //   - ADC count in this channel                                  
44 //
45 // As the Digits and SDigits have so much in common, the classes
46 // AliFMDDigitizer and AliFMDSDigitizer are implemented via a base
47 // class AliFMDBaseDigitizer.
48 //
49 //                 +---------------------+
50 //                 | AliFMDBaseDigitizer |
51 //                 +---------------------+
52 //                           ^
53 //                           |
54 //                +----------+---------+
55 //                |                    |
56 //      +-----------------+     +------------------+
57 //      | AliFMDDigitizer |     | AliFMDSDigitizer |
58 //      +-----------------+     +------------------+
59 //                |
60 //     +-------------------+
61 //     | AliFMDSSDigitizer |
62 //     +-------------------+
63 //
64 // These classes has several paramters: 
65 //
66 //     fPedestal
67 //     fPedestalWidth
68 //         (Only AliFMDDigitizer)
69 //         Mean and width of the pedestal.  The pedestal is simulated
70 //         by a Guassian, but derived classes my override MakePedestal
71 //         to simulate it differently (or pick it up from a database).
72 //
73 //     fVA1MipRange
74 //         The dymamic MIP range of the VA1_ALICE pre-amplifier chip 
75 //
76 //     fAltroChannelSize
77 //         The largest number plus one that can be stored in one
78 //         channel in one time step in the ALTRO ADC chip. 
79 //
80 //     fSampleRate
81 //         How many times the ALTRO ADC chip samples the VA1_ALICE
82 //         pre-amplifier signal.   The VA1_ALICE chip is read-out at
83 //         10MHz, while it's possible to drive the ALTRO chip at
84 //         25MHz.  That means, that the ALTRO chip can have time to
85 //         sample each VA1_ALICE signal up to 2 times.  Although it's
86 //         not certain this feature will be used in the production,
87 //         we'd like have the option, and so it should be reflected in
88 //         the code.
89 //
90 //
91 // The shaping function of the VA1_ALICE is generally given by 
92 //
93 //      f(x) = A(1 - exp(-Bx))
94 //
95 // where A is the total charge collected in the pre-amp., and B is a
96 // paramter that depends on the shaping time of the VA1_ALICE circut.
97 // 
98 // When simulating the shaping function of the VA1_ALICe
99 // pre-amp. chip, we have to take into account, that the shaping
100 // function depends on the previous value of read from the pre-amp. 
101 //
102 // That results in the following algorithm:
103 //
104 //    last = 0;
105 //    FOR charge IN pre-amp. charge train DO 
106 //      IF last < charge THEN 
107 //        f(t) = (charge - last) * (1 - exp(-B * t)) + last
108 //      ELSE
109 //        f(t) = (last - charge) * exp(-B * t) + charge)
110 //      ENDIF
111 //      FOR i IN # samples DO 
112 //        adc_i = f(i / (# samples))
113 //      DONE
114 //      last = charge
115 //   DONE
116 //
117 // Here, 
118 //
119 //   pre-amp. charge train 
120 //       is a series of 128 charges read from the VA1_ALICE chip
121 //
122 //   # samples
123 //       is the number of times the ALTRO ADC samples each of the 128
124 //       charges from the pre-amp. 
125 //
126 // Where Q is the total charge collected by the VA1_ALICE
127 // pre-amplifier.   Q is then given by 
128 //
129 //           E S 
130 //      Q =  - -
131 //           e R
132 //
133 // where E is the total energy deposited in a silicon strip, R is the
134 // dynamic range of the VA1_ALICE pre-amp (fVA1MipRange), e is the
135 // energy deposited by a single MIP, and S ALTRO channel size in each
136 // time step (fAltroChannelSize).  
137 //
138 // The energy deposited per MIP is given by 
139 //
140 //      e = M * rho * w 
141 //
142 // where M is the universal number 1.664, rho is the density of
143 // silicon, and w is the depth of the silicon sensor. 
144 //
145 // The final ADC count is given by 
146 //
147 //      C' = C + P
148 //
149 // where P is the (randomized) pedestal (see MakePedestal)
150 //
151 // This class uses the class template AliFMDMap<Type> to make an
152 // internal cache of the energy deposted of the hits.  The class
153 // template is instantasized as 
154 //
155 //  typedef AliFMDMap<std::pair<Float_t, UShort_t> > AliFMDEdepMap;
156 //
157 // The first member of the values is the summed energy deposition in a
158 // given strip, while the second member of the values is the number of
159 // hits in a given strip.  Using the second member, it's possible to
160 // do some checks on just how many times a strip got hit, and what
161 // kind of error we get in our reconstructed hits.  Note, that this
162 // information is currently not written to the digits tree.  I think a
163 // QA (Quality Assurance) digit tree is better suited for that task.
164 // However, the information is there to be used in the future. 
165 //
166 //
167 // Latest changes by Christian Holm Christensen
168 //
169 //////////////////////////////////////////////////////////////////////////////
170
171 //      /1
172 //      |           A(-1 + B + exp(-B))
173 //      | f(x) dx = ------------------- = 1
174 //      |                    B
175 //      / 0
176 //
177 // and B is the a parameter defined by the shaping time (fShapingTime).  
178 //
179 // Solving the above equation, for A gives
180 //
181 //                 B
182 //      A = ----------------
183 //          -1 + B + exp(-B)
184 //
185 // So, if we define the function g: [0,1] -> [0:1] by 
186 //
187 //               / v
188 //               |              Bu + exp(-Bu) - Bv - exp(-Bv) 
189 //      g(u,v) = | f(x) dx = -A -----------------------------
190 //               |                            B
191 //               / u
192 //
193 // we can evaluate the ALTRO sample of the VA1_ALICE pre-amp between
194 // any two times (u, v), by 
195 //       
196 //
197 //                                B         Bu + exp(-Bu) - Bv - exp(-Bv)
198 //      C = Q g(u,v) = - Q ---------------- -----------------------------
199 //                         -1 + B + exp(-B)              B                  
200 //
201 //               Bu + exp(-Bu) - Bv - exp(-Bv) 
202 //        = -  Q -----------------------------
203 //                    -1 + B + exp(-B)
204 //
205
206 #include <TTree.h>              // ROOT_TTree
207 #include <TFile.h>
208 #include "AliFMDDebug.h"        // Better debug macros
209 #include "AliFMDDigitizer.h"    // ALIFMDSSDIGITIZER_H
210 #include "AliFMD.h"             // ALIFMD_H
211 #include "AliFMDSDigit.h"       // ALIFMDDIGIT_H
212 #include "AliFMDDigit.h"        // ALIFMDDIGIT_H
213 #include "AliFMDParameters.h"   // ALIFMDPARAMETERS_H
214 #include <AliRunDigitizer.h>    // ALIRUNDIGITIZER_H
215 #include <AliRun.h>             // ALIRUN_H
216 #include <AliLoader.h>          // ALILOADER_H
217 #include <AliRunLoader.h>       // ALIRUNLOADER_H
218     
219 //====================================================================
220 ClassImp(AliFMDDigitizer)
221 #if 0
222 ;
223 #endif
224
225 //____________________________________________________________________
226 Bool_t
227 AliFMDDigitizer::Init()
228 {
229   // 
230   // Initialisation
231   // 
232   if (!AliFMDBaseDigitizer::Init()) return kFALSE;
233   
234 #if 0
235   // Get the AliRun object 
236   AliRun* run = fRunLoader->GetAliRun();
237   if (!run) { 
238     AliWarning("Loading gAlice");
239     fRunLoader->LoadgAlice();
240     if (!run) { 
241       AliError("Can not get Run from Run Loader");
242       return kFALSE;
243     }
244   }
245   
246   // Get the AliFMD object 
247   fFMD = static_cast<AliFMD*>(run->GetDetector("FMD"));
248   if (!fFMD) {
249     AliError("Can not get FMD from gAlice");
250     return kFALSE;
251   }  
252 #endif
253   return kTRUE;
254 }
255
256
257 //____________________________________________________________________
258 void
259 AliFMDDigitizer::Exec(Option_t*)
260 {
261   // 
262   // Execute this digitizer.  
263   // This member function will be called once per event by the passed
264   // AliRunDigitizer manager object. 
265   // 
266   // Parameters:
267   //    options Not used 
268   //
269   if (!fManager) { 
270     AliError("No digitisation manager defined");
271     return;
272   }
273
274   // Clear array of deposited energies 
275   fEdep.Reset();
276
277   AliRunLoader* runLoader = 0;
278   if (!gAlice) { 
279     TString folderName(fManager->GetInputFolderName(0));
280     runLoader = AliRunLoader::GetRunLoader(folderName.Data());
281     if (!runLoader) { 
282       AliError(Form("Failed at getting run loader from %s",
283                     folderName.Data()));
284       return;
285     }
286     if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
287     runLoader->GetAliRun();
288   }
289   if (!gAlice) { 
290     AliError("Can not get Run from Run Loader");
291     return;
292   }
293   
294   // Get the AliFMD object 
295   fFMD = static_cast<AliFMD*>(gAlice->GetDetector("FMD"));
296   if (!fFMD) {
297     AliError("Can not get FMD from gAlice");
298     return;
299   }  
300
301
302   // Loop over input files
303   Int_t nFiles= fManager->GetNinputs();
304   AliFMDDebug(1, (" Digitizing event number %d, got %d inputs",
305                   fManager->GetOutputEventNr(), nFiles));
306   for (Int_t inputFile = 0; inputFile < nFiles; inputFile++) {
307     AliFMDDebug(5, ("Now reading input # %d", inputFile));
308     // Get the current loader 
309     AliRunLoader* currentLoader = 
310       AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
311     if (!currentLoader) { 
312       Error("Exec", "no run loader for input file # %d", inputFile);
313       continue;
314     }
315
316     // Cache contriutions 
317     AliFMDDebug(5, ("Now summing the contributions from input # %d",inputFile));
318
319     // Get the FMD loader 
320     AliLoader* inFMD = currentLoader->GetLoader("FMDLoader");
321     // And load the summable digits
322     inFMD->LoadSDigits("READ");
323   
324     // Get the tree of summable digits
325     TTree* sdigitsTree = inFMD->TreeS();
326     if (!sdigitsTree)  {
327       AliError("No sdigit tree from manager");
328       continue;
329     }
330     if (AliLog::GetDebugLevel("FMD","") >= 10) {
331       TFile* file = sdigitsTree->GetCurrentFile();
332       if (!file) {
333         AliWarning("Input tree has no file!");
334       }
335       else { 
336         AliFMDDebug(10, ("Input tree file %s content:", file->GetName()));
337         file->ls();
338       }
339       // AliFMDDebug(5, ("Input tree %s file structure:", 
340       //                 sdigitsTree->GetName()));
341       // sdigitsTree->Print();
342     }
343
344     // Get the FMD branch 
345     TBranch* sdigitsBranch = sdigitsTree->GetBranch("FMD");
346     if (!sdigitsBranch) {
347       AliError("Failed to get sdigit branch");
348       return;
349     }
350
351     // Set the branch addresses 
352     fFMD->SetTreeAddress();
353
354     // Sum contributions from the sdigits
355     AliFMDDebug(3, ("Will now sum contributions from SDigits"));
356     SumContributions(sdigitsBranch);
357
358     // Unload the sdigits
359     inFMD->UnloadSDigits();
360   }  
361
362   TString       outFolder(fManager->GetOutputFolderName());
363   AliRunLoader* out    = AliRunLoader::GetRunLoader(outFolder.Data());
364   AliLoader*    outFMD = out->GetLoader("FMDLoader");
365   if (!outFMD) { 
366     AliError("Cannot get the FMDLoader output folder");
367     return;
368   }
369   TTree* outTree = MakeOutputTree(outFMD);
370   if (!outTree) { 
371     AliError("Failed to get output tree");
372     return;
373   }
374   // Set the branch address 
375   fFMD->SetTreeAddress();
376   
377   // And digitize the cached data 
378   DigitizeHits();
379   
380   // Fill the tree
381   Int_t write = outTree->Fill();
382   AliFMDDebug(5, ("Wrote %d bytes to digit tree", write));
383   
384   // Store the digits
385   StoreDigits(outFMD);
386 }
387
388 //____________________________________________________________________
389 void
390 AliFMDDigitizer::SumContributions(TBranch* sdigitsBranch) 
391 {
392   // 
393   // Sum contributions from SDigits 
394   // 
395   // Parameters:
396   //    sdigitsBranch Branch of SDigit data 
397   //
398   AliFMDDebug(3, ("Runnin our version of SumContributions"));
399
400   // Get a list of hits from the FMD manager 
401   TClonesArray *fmdSDigits = fFMD->SDigits();
402   
403   // Get number of entries in the tree 
404   Int_t nevents  = Int_t(sdigitsBranch->GetEntries());
405   
406   Int_t read = 0;
407   // Loop over the events in the 
408   for (Int_t event = 0; event < nevents; event++)  {
409     // Read in entry number `event' 
410     read += sdigitsBranch->GetEntry(event);
411     
412     // Get the number of sdigits 
413     Int_t nsdigits = fmdSDigits->GetEntries ();
414     AliFMDDebug(3, ("Got %5d SDigits", nsdigits));
415     for (Int_t sdigit = 0; sdigit < nsdigits; sdigit++) {
416       // Get the sdigit number `sdigit'
417       AliFMDSDigit* fmdSDigit = 
418         static_cast<AliFMDSDigit*>(fmdSDigits->UncheckedAt(sdigit));
419
420       AliFMDDebug(5, ("Adding contribution of %d tracks", 
421                       fmdSDigit->GetNTrack()));
422       AliFMDDebug(15, ("Contrib from FMD%d%c[%2d,%3d] (%s) from track %d", 
423                        fmdSDigit->Detector(),
424                        fmdSDigit->Ring(),
425                        fmdSDigit->Sector(),
426                        fmdSDigit->Strip(),
427                        fmdSDigit->GetName(), 
428                        fmdSDigit->GetTrack(0)));
429       
430       // Extract parameters 
431       AddContribution(fmdSDigit->Detector(),
432                       fmdSDigit->Ring(),
433                       fmdSDigit->Sector(),
434                       fmdSDigit->Strip(),
435                       fmdSDigit->Edep(), 
436                       kTRUE,
437                       fmdSDigit->GetNTrack(),
438                       fmdSDigit->GetTracks());
439     }  // sdigit loop
440   } // event loop
441
442
443   AliFMDDebug(3, ("Size of cache: %d bytes, read %d bytes", 
444                   int(sizeof(fEdep)), read));
445 }
446
447 //____________________________________________________________________
448 //
449 // EOF
450 // 
451
452
453
454