1 /**************************************************************************
2 * Copyright(c) 2004, 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 //////////////////////////////////////////////////////////////////////////////
20 // This class contains the procedures simulation ADC signal for the
21 // Forward Multiplicity detector : Hits->Digits and Hits->SDigits
28 // - ADC count in this channel
35 // - Total energy deposited in the strip
36 // - ADC count in this channel
38 // As the Digits and SDigits have so much in common, the classes
39 // AliFMDDigitizer and AliFMDSDigitizer are implemented via a base
40 // class AliFMDBaseDigitizer.
42 // +---------------------+
43 // | AliFMDBaseDigitizer |
44 // +---------------------+
47 // +----------+---------+
49 // +-----------------+ +------------------+
50 // | AliFMDDigitizer | | AliFMDSDigitizer |
51 // +-----------------+ +------------------+
53 // These classes has several paramters:
57 // (Only AliFMDDigitizer)
58 // Mean and width of the pedestal. The pedestal is simulated
59 // by a Guassian, but derived classes my override MakePedestal
60 // to simulate it differently (or pick it up from a database).
63 // The dymamic MIP range of the VA1_ALICE pre-amplifier chip
66 // The largest number plus one that can be stored in one
67 // channel in one time step in the ALTRO ADC chip.
70 // How many times the ALTRO ADC chip samples the VA1_ALICE
71 // pre-amplifier signal. The VA1_ALICE chip is read-out at
72 // 10MHz, while it's possible to drive the ALTRO chip at
73 // 25MHz. That means, that the ALTRO chip can have time to
74 // sample each VA1_ALICE signal up to 2 times. Although it's
75 // not certain this feature will be used in the production,
76 // we'd like have the option, and so it should be reflected in
79 // The shaping function of the VA1_ALICE is given by
81 // f(x) = A(1 - exp(-Bx))
83 // Where A is a normalization constant, tuned so that the integral
86 // | A(-1 + B + exp(-B))
87 // | f(x) dx = ------------------- = 1
91 // and B is the a parameter defined by the shaping time (fShapingTime).
93 // Solving the above equation, for A gives
96 // A = ----------------
99 // So, if we define the function g: [0,1] -> [0:1] by
102 // | Bu + exp(-Bu) - Bv - exp(-Bv)
103 // g(u,v) = | f(x) dx = -A -----------------------------
107 // we can evaluate the ALTRO sample of the VA1_ALICE pre-amp between
108 // any two times (u, v), by
111 // B Bu + exp(-Bu) - Bv - exp(-Bv)
112 // C = Q g(u,v) = - Q ---------------- -----------------------------
113 // -1 + B + exp(-B) B
115 // Bu + exp(-Bu) - Bv - exp(-Bv)
116 // = - Q -----------------------------
119 // Where Q is the total charge collected by the VA1_ALICE
120 // pre-amplifier. Q is then given by
126 // where E is the total energy deposited in a silicon strip, R is the
127 // dynamic range of the VA1_ALICE pre-amp (fVA1MipRange), e is the
128 // energy deposited by a single MIP, and S ALTRO channel size in each
129 // time step (fAltroChannelSize).
131 // The energy deposited per MIP is given by
135 // where M is the universal number 1.664, rho is the density of
136 // silicon, and w is the depth of the silicon sensor.
138 // The final ADC count is given by
142 // where P is the (randomized) pedestal (see MakePedestal)
144 // This class uses the class template AliFMDMap<Type> to make an
145 // internal cache of the energy deposted of the hits. The class
146 // template is instantasized as
148 // typedef AliFMDMap<std::pair<Float_t, UShort_t> > AliFMDEdepMap;
150 // The first member of the values is the summed energy deposition in a
151 // given strip, while the second member of the values is the number of
152 // hits in a given strip. Using the second member, it's possible to
153 // do some checks on just how many times a strip got hit, and what
154 // kind of error we get in our reconstructed hits. Note, that this
155 // information is currently not written to the digits tree. I think a
156 // QA (Quality Assurance) digit tree is better suited for that task.
157 // However, the information is there to be used in the future.
160 // Latest changes by Christian Holm Christensen
162 //////////////////////////////////////////////////////////////////////////////
168 # include <TRandom.h>
173 #ifndef ALIFMDDIGITIZER_H
174 # include "AliFMDDigitizer.h"
180 # include "AliFMDHit.h"
182 #ifndef ALIFMDDIGIT_H
183 # include "AliFMDDigit.h"
185 #ifndef ALIFMDDIGIT_H
186 # include "AliFMDSDigit.h"
188 #ifndef ALIRUNDIGITIZER_H
189 # include "AliRunDigitizer.h"
195 # include "AliLoader.h"
197 #ifndef ALIRUNLOADER_H
198 # include "AliRunLoader.h"
201 //____________________________________________________________________
202 ClassImp(AliFMDEdepMap);
204 //====================================================================
205 ClassImp(AliFMDBaseDigitizer);
207 //____________________________________________________________________
208 AliFMDBaseDigitizer::AliFMDBaseDigitizer()
211 // Default ctor - don't use it
214 //____________________________________________________________________
215 AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager)
216 : AliDigitizer(manager, "AliFMDBaseDigitizer", "FMD Digitizer base class"),
218 fEdep(kMaxDetectors, kMaxRings, kMaxSectors, kMaxStrips)
221 AliDebug(1," processed");
223 SetAltroChannelSize();
228 //____________________________________________________________________
229 AliFMDBaseDigitizer::AliFMDBaseDigitizer(const Char_t* name,
231 : AliDigitizer(name, title),
233 fEdep(kMaxDetectors, kMaxRings, kMaxSectors, kMaxStrips)
236 AliDebug(1," processed");
238 SetAltroChannelSize();
243 //____________________________________________________________________
244 AliFMDBaseDigitizer::~AliFMDBaseDigitizer()
249 //____________________________________________________________________
251 AliFMDBaseDigitizer::Init()
258 //____________________________________________________________________
260 AliFMDBaseDigitizer::SumContributions(AliFMD* fmd)
262 // Sum energy deposited contributions from each hit in a cache
265 Fatal("SumContributions", "no run loader");
267 // Clear array of deposited energies
270 // Get the FMD loader
271 AliLoader* inFMD = fRunLoader->GetLoader("FMDLoader");
273 inFMD->LoadHits("READ");
275 // Get the tree of hits
276 TTree* hitsTree = inFMD->TreeH();
279 inFMD->LoadHits("READ");
280 hitsTree = inFMD->TreeH();
283 // Get the FMD branch
284 TBranch* hitsBranch = hitsTree->GetBranch("FMD");
285 if (hitsBranch) fmd->SetHitsAddressBranch(hitsBranch);
286 else AliFatal("Branch FMD hit not found");
288 // Get a list of hits from the FMD manager
289 TClonesArray *fmdHits = fmd->Hits();
291 // Get number of entries in the tree
292 Int_t ntracks = Int_t(hitsTree->GetEntries());
294 // Loop over the tracks in the
295 for (Int_t track = 0; track < ntracks; track++) {
296 // Read in entry number `track'
297 hitsBranch->GetEntry(track);
299 // Get the number of hits
300 Int_t nhits = fmdHits->GetEntries ();
301 for (Int_t hit = 0; hit < nhits; hit++) {
302 // Get the hit number `hit'
304 static_cast<AliFMDHit*>(fmdHits->UncheckedAt(hit));
306 // Extract parameters
307 UShort_t detector = fmdHit->Detector();
308 Char_t ring = fmdHit->Ring();
309 UShort_t sector = fmdHit->Sector();
310 UShort_t strip = fmdHit->Strip();
311 Float_t edep = fmdHit->Edep();
312 if (fEdep(detector, ring, sector, strip).first != 0)
313 AliDebug(1, Form("Double hit in %d%c(%d,%d)",
314 detector, ring, sector, strip));
316 fEdep(detector, ring, sector, strip).first += edep;
317 fEdep(detector, ring, sector, strip).second += 1;
318 // Add this to the energy deposited for this strip
323 //____________________________________________________________________
325 AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
327 // For the stored energy contributions in the cache (fEdep), convert
328 // the energy signal to ADC counts, and store the created digit in
329 // the digits array (AliFMD::fDigits)
332 for (UShort_t detector=1; detector <= 3; detector++) {
333 // Get pointer to subdetector
334 AliFMDSubDetector* det = 0;
336 case 1: det = fmd->GetFMD1(); break;
337 case 2: det = fmd->GetFMD2(); break;
338 case 3: det = fmd->GetFMD3(); break;
341 for (UShort_t ringi = 0; ringi <= 1; ringi++) {
342 // Get pointer to Ring
345 case 0: if (det->GetInner()) r = det->GetInner(); break;
346 case 1: if (det->GetOuter()) r = det->GetOuter(); break;
350 // Get number of sectors
351 UShort_t nSectors = UShort_t(360. / r->GetTheta());
352 // Loop over the number of sectors
353 for (UShort_t sector = 0; sector < nSectors; sector++) {
354 // Get number of strips
355 UShort_t nStrips = r->GetNStrips();
356 // Loop over the stips
357 for (UShort_t strip = 0; strip < nStrips; strip++) {
359 Float_t edep = fEdep(detector, r->GetId(), sector, strip).first;
360 ConvertToCount(edep, r->GetSiThickness(), fmd->GetSiDensity(),
362 AddDigit(fmd, detector, r->GetId(), sector, strip,
363 edep, UShort_t(counts[0]),
364 Short_t(counts[1]), Short_t(counts[2]));
366 // This checks if the digit created will give the `right'
367 // number of particles when reconstructed, using a naiive
368 // approach. It's here only as a quality check - nothing
370 CheckDigit(fEdep(detector, r->GetId(), sector, strip).first,
371 fEdep(detector, r->GetId(), sector, strip).second,
380 //____________________________________________________________________
382 AliFMDBaseDigitizer::ShapeIntegral(Float_t u, Float_t v) const
384 // Calculates the integral
387 // | Bu + exp(-Bu) - Bv - exp(-Bv)
388 // | f(x) dx = - A -----------------------------
392 // of the shaping function of the VA1_ALICE between times u and v
394 // f(x) = A(1 - exp(-Bx))
396 // where A is a normalization constant, tuned so that the integral
400 // | f(x) dx = 1 => A = ----------------
401 // | -1 + B + exp(-B)
404 // and B is the a parameter defined by the shaping time (fShapingTime).
406 // That is, the function return the value
408 // Bu + exp(-Bu) - Bv - exp(-Bv)
409 // - -----------------------------
412 // u,v should lie in the interval [0,1], and u < v
413 if (u == 0 && v == 1) return 1;
414 Float_t B = fShapingTime;
416 // Calculate the integral
417 Float_t res = - ((B * u + TMath::Exp(-B * u) - B * v - TMath::Exp(-B * v)) /
418 (-1 + B + TMath::Exp(-B)));
422 //____________________________________________________________________
424 AliFMDBaseDigitizer::ConvertToCount(Float_t edep,
427 TArrayI& counts) const
429 // Put noise and make ADC signal
430 // This is calculated as the product
432 // DeltaEmip * SiThickness * SiDensity / Number
436 // DeltaEmip is the energy loss of a MIP
437 // SiThickness is the thickness of the silicon
438 // SiDensity is the Silicon density
439 // Number is # of e^- per MIP
441 // Note: Need to check this is correct.
443 const Float_t mipI = 1.664 * siThickness * siDensity;
444 // const Float_t mipI = 1.664 * 0.04 * 2.33 / 22400; // = 6.923e-6;
447 UShort_t ped = MakePedestal();
449 Float_t convf = 1 / mipI * Float_t(fAltroChannelSize) / fVA1MipRange;
450 Int_t n = fSampleRate;
451 for (Ssiz_t i = 0; i < n; i++) {
452 Float_t w = ShapeIntegral(Float_t(i)/n, Float_t(i+1)/n);
453 counts[i] = UShort_t(TMath::Min(w * edep * convf + ped,
454 Float_t(fAltroChannelSize)));
459 //====================================================================
460 ClassImp(AliFMDDigitizer);
462 //____________________________________________________________________
463 AliFMDDigitizer::AliFMDDigitizer()
464 : AliFMDBaseDigitizer()
466 // Default ctor - don't use it
469 //____________________________________________________________________
470 AliFMDDigitizer::AliFMDDigitizer(AliRunDigitizer* manager)
471 : AliFMDBaseDigitizer(manager)
474 AliDebug(1," processed");
478 //____________________________________________________________________
480 AliFMDDigitizer::Exec(Option_t*)
482 // Get the output manager
483 TString outFolder(fManager->GetOutputFolderName());
485 AliRunLoader::GetRunLoader(outFolder.Data());
486 // Get the FMD output manager
487 AliLoader* outFMD = out->GetLoader("FMDLoader");
489 // Get the input loader
490 TString inFolder(fManager->GetInputFolderName(0));
492 AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
494 AliError("Can not find Run Loader for input stream 0");
497 // Get the AliRun object
498 if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
500 // Get the AliFMD object
501 AliFMD* fmd = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
503 AliError("Can not get FMD from gAlice");
507 Int_t nFiles= fManager->GetNinputs();
508 for (Int_t inputFile = 0; inputFile < nFiles; inputFile++) {
509 AliDebug(1,Form(" Digitizing event number %d",
510 fManager->GetOutputEventNr()));
511 // Get the current loader
513 AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
514 if (!fRunLoader) Fatal("Exec", "no run loader");
515 // Cache contriutions
516 SumContributions(fmd);
518 // Digitize the event
521 // Load digits from the tree
522 outFMD->LoadDigits("update");
523 // Get the tree of digits
524 TTree* digitTree = outFMD->TreeD();
526 outFMD->MakeTree("D");
527 digitTree = outFMD->TreeD();
530 // Make a branch in the tree
531 TClonesArray* digits = fmd->Digits();
532 fmd->MakeBranchInTree(digitTree, fmd->GetName(), &(digits), 4000, 0);
533 // TBranch* digitBranch = digitTree->GetBranch(fmd->GetName());
537 // Write the digits to disk
538 outFMD->WriteDigits("OVERWRITE");
539 outFMD->UnloadHits();
540 outFMD->UnloadDigits();
542 // Reset the digits in the AliFMD object
547 //____________________________________________________________________
549 AliFMDDigitizer::MakePedestal() const
551 return UShort_t(TMath::Max(gRandom->Gaus(fPedestal, fPedestalWidth), 0.));
554 //____________________________________________________________________
556 AliFMDDigitizer::AddDigit(AliFMD* fmd,
564 Short_t count3) const
566 fmd->AddDigit(detector, ring, sector, strip, count1, count2, count3);
569 //____________________________________________________________________
571 AliFMDDigitizer::CheckDigit(Float_t /* edep */,
573 const TArrayI& counts)
575 Int_t integral = counts[0];
576 if (counts[1] >= 0) integral += counts[1];
577 if (counts[2] >= 0) integral += counts[2];
578 integral -= Int_t(fPedestal + 2 * fPedestalWidth);
579 if (integral < 0) integral = 0;
581 Float_t convf = Float_t(fVA1MipRange) / fAltroChannelSize;
582 Float_t mips = integral * convf;
583 if (mips > Float_t(nhits) + .5 || mips < Float_t(nhits) - .5)
584 Warning("CheckDigit", "Digit -> %4.2f MIPS != %d +/- .5 hits",
588 //====================================================================
589 ClassImp(AliFMDSDigitizer);
591 //____________________________________________________________________
592 AliFMDSDigitizer::AliFMDSDigitizer()
594 // Default ctor - don't use it
597 //____________________________________________________________________
598 AliFMDSDigitizer::AliFMDSDigitizer(const Char_t* headerFile,
599 const Char_t* /* sdigfile */)
600 : AliFMDBaseDigitizer("FMDSDigitizer", "FMD SDigitizer")
603 AliDebug(1," processed");
605 fRunLoader = AliRunLoader::GetRunLoader(); // Open(headerFile);
607 Fatal("AliFMDSDigitizer", "cannot open session, header file '%s'",
609 AliLoader* loader = fRunLoader->GetLoader("FMDLoader");
611 Fatal("AliFMDSDigitizer", "cannot find FMD loader in specified event");
613 // Add task to tasks folder
614 loader->PostSDigitizer(this);
618 //____________________________________________________________________
619 AliFMDSDigitizer::~AliFMDSDigitizer()
621 AliLoader* loader = fRunLoader->GetLoader("FMDLoader");
622 loader->CleanSDigitizer();
625 //____________________________________________________________________
627 AliFMDSDigitizer::Exec(Option_t*)
629 // Get the output manager
631 Error("Exec", "Run loader is not set");
634 if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
635 if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
637 AliLoader* fmdLoader = fRunLoader->GetLoader("FMDLoader");
638 if (!fmdLoader) Fatal("Exec", "no FMD loader");
640 // Get the AliFMD object
642 static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
644 AliError("Can not get FMD from gAlice");
648 Int_t nEvents = Int_t(fRunLoader->TreeE()->GetEntries());
649 for (Int_t event = 0; event < nEvents; event++) {
650 AliDebug(1,Form(" Digitizing event number %d", event));
651 // Get the current loader
652 fRunLoader->GetEvent(event);
654 if (!fmdLoader->TreeS()) fmdLoader->MakeTree("S");
656 fmd->MakeBranch("S");
658 // Cache contriutions
659 SumContributions(fmd);
661 // Digitize the event
664 fmdLoader->TreeS()->Reset();
665 fmdLoader->TreeS()->Fill();
666 fmdLoader->WriteSDigits("OVERWRITE");
670 //____________________________________________________________________
672 AliFMDSDigitizer::AddDigit(AliFMD* fmd,
680 Short_t count3) const
682 fmd->AddSDigit(detector, ring, sector, strip, edep, count1, count2, count3);
687 //____________________________________________________________________