]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FMD/AliFMDDigitizer.cxx
Adding calibration object for the sharing efficiency
[u/mrichter/AliRoot.git] / FMD / AliFMDDigitizer.cxx
index 43a229c941f3435e7cde7a9c7abed7265a7e6528..c17a265073221a93ba75a8e3f49887dbe8a28a26 100644 (file)
@@ -1,5 +1,5 @@
 /**************************************************************************
- * Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights reserved. *
+ * Copyright(c) 2004, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
  * Author: The ALICE Off-line Project.                                    *
  * Contributors are mentioned in the code where appropriate.              *
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-
- //////////////////////////////////////////////////////////////////////////////
-//                                                                           //
-//  Forward Multiplicity Detector based on Silicon plates                    //
-//  This class contains the procedures simulation ADC  signal for            //
-//  the Forward Multiplicity detector  : hits -> digits                      //
-//  ADC signal consists                                                      //
-//   - number of detector;                                                   //
-//   - number of ring;                                                       //
-//   - number of sector;                                                     //
-//   - ADC signal in this channel                                            //
-//                                                                           //
- //////////////////////////////////////////////////////////////////////////////
-
-#include <TTree.h> 
-#include <TVector.h>
-#include <TObjArray.h>
+/* $Id$ */
+/** @file    AliFMDDigitizer.cxx
+    @author  Christian Holm Christensen <cholm@nbi.dk>
+    @date    Mon Mar 27 12:38:26 2006
+    @brief   FMD Digitizers implementation
+    @ingroup FMD_sim
+*/
+//////////////////////////////////////////////////////////////////////////////
+//
+//  This class contains the procedures simulation ADC  signal for the
+//  Forward Multiplicity detector  : SDigits->Digits
+// 
+//  Digits consists of
+//   - Detector #
+//   - Ring ID                                             
+//   - Sector #     
+//   - Strip #
+//   - ADC count in this channel                                  
+//
+//  Digits consists of
+//   - Detector #
+//   - Ring ID                                             
+//   - Sector #     
+//   - Strip #
+//   - Total energy deposited in the strip
+//   - ADC count in this channel                                  
+//
+// As the Digits and SDigits have so much in common, the classes
+// AliFMDDigitizer and AliFMDSDigitizer are implemented via a base
+// class AliFMDBaseDigitizer.
+//
+//                 +---------------------+
+//                 | AliFMDBaseDigitizer |
+//                 +---------------------+
+//                           ^
+//                           |
+//                +----------+---------+
+//                |                    |
+//      +-----------------+     +------------------+
+//      | AliFMDDigitizer |    | AliFMDSDigitizer |
+//      +-----------------+    +------------------+
+//                |
+//     +-------------------+
+//     | AliFMDSSDigitizer |
+//     +-------------------+
+//
+// These classes has several paramters: 
+//
+//     fPedestal
+//     fPedestalWidth
+//         (Only AliFMDDigitizer)
+//         Mean and width of the pedestal.  The pedestal is simulated
+//         by a Guassian, but derived classes my override MakePedestal
+//         to simulate it differently (or pick it up from a database).
+//
+//     fVA1MipRange
+//         The dymamic MIP range of the VA1_ALICE pre-amplifier chip 
+//
+//     fAltroChannelSize
+//         The largest number plus one that can be stored in one
+//         channel in one time step in the ALTRO ADC chip. 
+//
+//     fSampleRate
+//         How many times the ALTRO ADC chip samples the VA1_ALICE
+//         pre-amplifier signal.   The VA1_ALICE chip is read-out at
+//         10MHz, while it's possible to drive the ALTRO chip at
+//         25MHz.  That means, that the ALTRO chip can have time to
+//         sample each VA1_ALICE signal up to 2 times.  Although it's
+//         not certain this feature will be used in the production,
+//         we'd like have the option, and so it should be reflected in
+//         the code.
+//
+//
+// The shaping function of the VA1_ALICE is generally given by 
+//
+//      f(x) = A(1 - exp(-Bx))
+//
+// where A is the total charge collected in the pre-amp., and B is a
+// paramter that depends on the shaping time of the VA1_ALICE circut.
+// 
+// When simulating the shaping function of the VA1_ALICe
+// pre-amp. chip, we have to take into account, that the shaping
+// function depends on the previous value of read from the pre-amp. 
+//
+// That results in the following algorithm:
+//
+//    last = 0;
+//    FOR charge IN pre-amp. charge train DO 
+//      IF last < charge THEN 
+//        f(t) = (charge - last) * (1 - exp(-B * t)) + last
+//      ELSE
+//        f(t) = (last - charge) * exp(-B * t) + charge)
+//      ENDIF
+//      FOR i IN # samples DO 
+//        adc_i = f(i / (# samples))
+//      DONE
+//      last = charge
+//   DONE
+//
+// Here, 
+//
+//   pre-amp. charge train 
+//       is a series of 128 charges read from the VA1_ALICE chip
+//
+//   # samples
+//       is the number of times the ALTRO ADC samples each of the 128
+//       charges from the pre-amp. 
+//
+// Where Q is the total charge collected by the VA1_ALICE
+// pre-amplifier.   Q is then given by 
+//
+//           E S 
+//      Q =  - -
+//           e R
+//
+// where E is the total energy deposited in a silicon strip, R is the
+// dynamic range of the VA1_ALICE pre-amp (fVA1MipRange), e is the
+// energy deposited by a single MIP, and S ALTRO channel size in each
+// time step (fAltroChannelSize).  
+//
+// The energy deposited per MIP is given by 
+//
+//      e = M * rho * w 
+//
+// where M is the universal number 1.664, rho is the density of
+// silicon, and w is the depth of the silicon sensor. 
+//
+// The final ADC count is given by 
+//
+//      C' = C + P
+//
+// where P is the (randomized) pedestal (see MakePedestal)
+//
+// This class uses the class template AliFMDMap<Type> to make an
+// internal cache of the energy deposted of the hits.  The class
+// template is instantasized as 
+//
+//  typedef AliFMDMap<std::pair<Float_t, UShort_t> > AliFMDEdepMap;
+//
+// The first member of the values is the summed energy deposition in a
+// given strip, while the second member of the values is the number of
+// hits in a given strip.  Using the second member, it's possible to
+// do some checks on just how many times a strip got hit, and what
+// kind of error we get in our reconstructed hits.  Note, that this
+// information is currently not written to the digits tree.  I think a
+// QA (Quality Assurance) digit tree is better suited for that task.
+// However, the information is there to be used in the future. 
+//
+//
+// Latest changes by Christian Holm Christensen
+//
+//////////////////////////////////////////////////////////////////////////////
+
+//      /1
+//      |           A(-1 + B + exp(-B))
+//      | f(x) dx = ------------------- = 1
+//      |                    B
+//      / 0
+//
+// and B is the a parameter defined by the shaping time (fShapingTime).  
+//
+// Solving the above equation, for A gives
+//
+//                 B
+//      A = ----------------
+//          -1 + B + exp(-B)
+//
+// So, if we define the function g: [0,1] -> [0:1] by 
+//
+//               / v
+//               |              Bu + exp(-Bu) - Bv - exp(-Bv) 
+//      g(u,v) = | f(x) dx = -A -----------------------------
+//               |                            B
+//               / u
+//
+// we can evaluate the ALTRO sample of the VA1_ALICE pre-amp between
+// any two times (u, v), by 
+//       
+//
+//                                B        Bu + exp(-Bu) - Bv - exp(-Bv)
+//      C = Q g(u,v) = - Q ---------------- -----------------------------
+//                        -1 + B + exp(-B)              B                  
+//
+//               Bu + exp(-Bu) - Bv - exp(-Bv) 
+//        = -  Q -----------------------------
+//                    -1 + B + exp(-B)
+//
+
+#include <TTree.h>             // ROOT_TTree
 #include <TFile.h>
-#include <TDirectory.h>
-#include <TRandom.h>
-
-
-#include "AliFMDDigitizer.h"
-#include "AliFMD.h"
-#include "AliFMDhit.h"
-#include "AliFMDdigit.h"
-#include "AliRunDigitizer.h"
-
-#include "AliRun.h"
-#include "AliPDG.h"
-#include "AliLoader.h"
-#include "AliRunLoader.h"
-
-#include <stdlib.h>
-#include <Riostream.h>
-#include <Riostream.h>
-
+#include "AliFMDDebug.h"       // Better debug macros
+#include "AliFMDDigitizer.h"   // ALIFMDSSDIGITIZER_H
+#include "AliFMD.h"            // ALIFMD_H
+#include "AliFMDSDigit.h"      // ALIFMDDIGIT_H
+#include "AliFMDDigit.h"       // ALIFMDDIGIT_H
+#include "AliFMDParameters.h"   // ALIFMDPARAMETERS_H
+#include <AliRunDigitizer.h>   // ALIRUNDIGITIZER_H
+#include <AliRun.h>            // ALIRUN_H
+#include <AliLoader.h>         // ALILOADER_H
+#include <AliRunLoader.h>      // ALIRUNLOADER_H
+    
+//====================================================================
 ClassImp(AliFMDDigitizer)
+#if 0
+;
+#endif
 
-//___________________________________________
-  AliFMDDigitizer::AliFMDDigitizer()  :AliDigitizer()
-{
-// Default ctor - don't use it
-  ;
-}
-
-//___________________________________________
-AliFMDDigitizer::AliFMDDigitizer(AliRunDigitizer* manager) 
-    :AliDigitizer(manager) 
-{
-  // ctor which should be used
-  //  fDebug =0;
-  // if (GetDebug()>2)
-  //  cerr<<"AliFMDDigitizer::AliFMDDigitizer"
-  //     <<"(AliRunDigitizer* manager) was processed"<<endl;
-}
-
-//------------------------------------------------------------------------
-AliFMDDigitizer::~AliFMDDigitizer()
-{
-// Destructor
-}
-
- //------------------------------------------------------------------------
-Bool_t AliFMDDigitizer::Init()
+//____________________________________________________________________
+Bool_t
+AliFMDDigitizer::Init()
 {
-// Initialization
-// cout<<"AliFMDDigitizer::Init"<<endl;
- return kTRUE;
+  // Initialisation
+  if (!AliFMDBaseDigitizer::Init()) return kFALSE;
+  
+#if 0
+  // Get the AliRun object 
+  AliRun* run = fRunLoader->GetAliRun();
+  if (!run) { 
+    AliWarning("Loading gAlice");
+    fRunLoader->LoadgAlice();
+    if (!run) { 
+      AliError("Can not get Run from Run Loader");
+      return kFALSE;
+    }
+  }
+  
+  // Get the AliFMD object 
+  fFMD = static_cast<AliFMD*>(run->GetDetector("FMD"));
+  if (!fFMD) {
+    AliError("Can not get FMD from gAlice");
+    return kFALSE;
+  }  
+#endif
+  return kTRUE;
 }
 
-//---------------------------------------------------------------------
 
-void AliFMDDigitizer::Exec(Option_t * /*option*/)
+//____________________________________________________________________
+void
+AliFMDDigitizer::Exec(Option_t*)
 {
+  if (!fManager) { 
+    AliError("No digitisation manager defined");
+    return;
+  }
+
+  // Clear array of deposited energies 
+  fEdep.Reset();
+
+  AliRunLoader* runLoader = 0;
+  if (!gAlice) { 
+    TString folderName(fManager->GetInputFolderName(0));
+    runLoader = AliRunLoader::GetRunLoader(folderName.Data());
+    if (!runLoader) { 
+      AliError(Form("Failed at getting run loader from %s",
+                   folderName.Data()));
+      return;
+    }
+    if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
+    runLoader->GetAliRun();
+  }
+  if (!gAlice) { 
+    AliError("Can not get Run from Run Loader");
+    return;
+  }
+  
+  // Get the AliFMD object 
+  fFMD = static_cast<AliFMD*>(gAlice->GetDetector("FMD"));
+  if (!fFMD) {
+    AliError("Can not get FMD from gAlice");
+    return;
+  }  
+
+
+  // Loop over input files
+  Int_t nFiles= fManager->GetNinputs();
+  AliFMDDebug(1, (" Digitizing event number %d, got %d inputs",
+                 fManager->GetOutputEventNr(), nFiles));
+  for (Int_t inputFile = 0; inputFile < nFiles; inputFile++) {
+    AliFMDDebug(5, ("Now reading input # %d", inputFile));
+    // Get the current loader 
+    AliRunLoader* currentLoader = 
+      AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
+    if (!currentLoader) { 
+      Error("Exec", Form("no run loader for input file # %d", inputFile));
+      continue;
+    }
 
-  /*
-   Conver hits to digits:
-   - number of detector;
-   - number of ring;
-   - number of sector;
-   - ADC signal in this channel
-  */
+    // Cache contriutions 
+    AliFMDDebug(5, ("Now summing the contributions from input # %d",inputFile));
 
-  AliRunLoader *inRL, *outRL;//in and out Run Loaders
-  AliLoader *ingime, *outgime;// in and out ITSLoaders
+    // Get the FMD loader 
+    AliLoader* inFMD = currentLoader->GetLoader("FMDLoader");
+    // And load the summable digits
+    inFMD->LoadSDigits("READ");
+  
+    // Get the tree of summable digits
+    TTree* sdigitsTree = inFMD->TreeS();
+    if (!sdigitsTree)  {
+      AliError("No sdigit tree from manager");
+      continue;
+    }
+    if (AliLog::GetDebugLevel("FMD","") >= 10) {
+      TFile* file = sdigitsTree->GetCurrentFile();
+      if (!file) {
+       AliWarning("Input tree has no file!");
+      }
+      else { 
+       AliFMDDebug(10, ("Input tree file %s content:", file->GetName()));
+       file->ls();
+      }
+      // AliFMDDebug(5, ("Input tree %s file structure:", 
+      //                 sdigitsTree->GetName()));
+      // sdigitsTree->Print();
+    }
 
-  outRL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
-  outgime = outRL->GetLoader("FMDLoader");
+    // Get the FMD branch 
+    TBranch* sdigitsBranch = sdigitsTree->GetBranch("FMD");
+    if (!sdigitsBranch) {
+      AliError("Failed to get sdigit branch");
+      return;
+    }
 
-#ifdef DEBUG
-  cout<<"AliFMDDigitizer::>SDigits2Digits start...\n";
-#endif
+    // Set the branch addresses 
+    fFMD->SetTreeAddress();
+
+    // Sum contributions from the sdigits
+    AliFMDDebug(3, ("Will now sum contributions from SDigits"));
+    SumContributions(sdigitsBranch);
+
+    // Unload the sdigits
+    inFMD->UnloadSDigits();
+  }  
+
+  TString       outFolder(fManager->GetOutputFolderName());
+  AliRunLoader* out    = AliRunLoader::GetRunLoader(outFolder.Data());
+  AliLoader*    outFMD = out->GetLoader("FMDLoader");
+  if (!outFMD) { 
+    AliError("Cannot get the FMDLoader output folder");
+    return;
+  }
+  TTree* outTree = MakeOutputTree(outFMD);
+  if (!outTree) { 
+    AliError("Failed to get output tree");
+    return;
+  }
+  // Set the branch address 
+  fFMD->SetTreeAddress();
+  
+  // And digitize the cached data 
+  DigitizeHits();
+  
+  // Fill the tree
+  Int_t write = outTree->Fill();
+  AliFMDDebug(5, ("Wrote %d bytes to digit tree", write));
+  
+  // Store the digits
+  StoreDigits(outFMD);
+}
 
+//____________________________________________________________________
+void
+AliFMDDigitizer::SumContributions(TBranch* sdigitsBranch) 
+{
+  // Sum energy deposited contributions from each sdigits in a cache
+  // (fEdep).  
+  AliFMDDebug(3, ("Runnin our version of SumContributions"));
 
-  Int_t volume, sector, ring, charge;
-  Float_t e;
-  Float_t de[10][50][520];
-  Int_t hit;
-  Int_t digit[5];
-  Int_t ivol, iSector, iRing;
-  for (Int_t i=0; i<10; i++)
-    for(Int_t j=0; j<50; j++)
-      for(Int_t ij=0; ij<520; ij++)
-     de[i][j][ij]=0;
-  Int_t numberOfRings[5]=
-  {512,256,512,256,512};
-  Int_t numberOfSector[5]=
-  {20,40,20,40,20}; 
+  // Get a list of hits from the FMD manager 
+  TClonesArray *fmdSDigits = fFMD->SDigits();
   
-  AliFMDhit *fmdHit=0;
-  TTree *tH=0;
-  TBranch *brHits=0;
-  TBranch *brD=0;
-
-  inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
-  if (inRL == 0x0)
-    {
-      Error("Exec","Can not find Run Loader for input stream 0");
-      return;
-    }
-  //  Info("Exec","inRL->GetAliRun() %#x",inRL->GetAliRun());
-
-  if (!inRL->GetAliRun()) inRL->LoadgAlice();
-
-  AliFMD * fFMD = (AliFMD *) inRL->GetAliRun()->GetDetector("FMD");
-  //  Info("Exec","inRL->GetAliRun(): %#x, FMD: %#x, InRL %#x.",inRL->GetAliRun(),fFMD,inRL);
-  if (fFMD == 0x0)
-   {
-     Error("Exec","Can not get FMD from gAlice");
-     return;
-   }
-// Loop over files to digitize
-
-  Int_t nFiles=GetManager()->GetNinputs();
-  for (Int_t inputFile=0; inputFile<nFiles;inputFile++) 
-   {
-    cout<<" event "<<fManager->GetOutputEventNr()<<endl;
-    if (fFMD)
-     {
-
-      inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
-      ingime = inRL->GetLoader("FMDLoader");
-      ingime->LoadHits("READ");//probably it is necessary to load them before
-      
-      
-      tH = ingime->TreeH();
-      if (tH == 0x0)
-       {
-         ingime->LoadHits("read");
-         tH = ingime->TreeH();
-       }
-      brHits = tH->GetBranch("FMD");
-      if (brHits) {
-  //      brHits->SetAddress(&fHits);
-          fFMD->SetHitsAddressBranch(brHits);
-      }else{
-        Fatal("Exec","EXEC Branch FMD hit not found");
-      }
-      TClonesArray *fFMDhits = fFMD->Hits ();
-      
-      Int_t ntracks    = (Int_t) tH->GetEntries();
-      cout<<"Number of tracks TreeH"<<ntracks<<endl;
-      for (Int_t track = 0; track < ntracks; track++)
-       {
-         brHits->GetEntry(track);
-         Int_t nhits = fFMDhits->GetEntries ();
-         // if(nhits>0) cout<<"nhits "<<nhits<<endl;
-         for (hit = 0; hit < nhits; hit++)
-           {
-             fmdHit = (AliFMDhit *) fFMDhits->UncheckedAt(hit);
-
-             volume = fmdHit->Volume ();
-             sector = fmdHit->NumberOfSector ();
-             ring = fmdHit->NumberOfRing ();
-             e = fmdHit->Edep ();
-             de[volume][sector][ring] += e;
-            //   if (fManager->GetOutputEventNr()>1)
-                    //      cout<<" "<<volume<<" "<<sector<<" "<<ring<<endl;
-           }          //hit loop
-       }               //track loop
-    }               
-//if FMD
-
-  // Put noise and make ADC signal
-   Float_t mipI = 1.664 * 0.04 * 2.33 / 22400;     // = 6.923e-6;
-   for ( ivol=1; ivol<=5; ivol++){
-     for ( iSector=1; iSector<=numberOfSector[ivol-1]; iSector++){
-       for ( iRing=1; iRing<=numberOfRings[ivol-1]; iRing++){
-         digit[0]=ivol;
-         digit[1]=iSector;
-         digit[2]=iRing;
-         charge = Int_t (de[ivol][iSector][iRing] / mipI);
-         Int_t pedestal=Int_t(gRandom->Gaus(500,250));
-  //       digit[3]=PutNoise(charge);
-         digit[3]=charge + pedestal;
-         if(digit[3]<= 500) digit[3]=500; 
-    //dynamic range from MIP(0.155MeV) to 30MIP(4.65MeV)
-    //1024 ADC channels 
-         Float_t channelWidth=(22400*50)/1024;
-         digit[4]=Int_t(digit[3]/channelWidth);
-         if (digit[4]>1024) digit[4]=1024; 
-         fFMD->AddDigit(digit);
-       } //ivol
-     } //iSector
-   } //iRing
-
-   TTree* treeD = outgime->TreeD();
-   cout<<" treeD "<<treeD;
-   if (treeD == 0x0) {
-     outgime->MakeTree("D");
-     treeD = outgime->TreeD();
-     cout<<" After MakeTree "<<treeD<<endl;
-   }
-   cout<<" Before reset "<<treeD<<endl;
-   //   treeD->Clear();
-   treeD->Reset();
-   fFMD->MakeBranchInTreeD(treeD);
-   brD = treeD->GetBranch("FMD");
-   cout<<" Make branch "<<brD<<endl;
-
-   treeD->Fill();  //this operator does not work for events >1
-   //PH   treeD->Print();
-   outgime->WriteDigits("OVERWRITE");
+  // Get number of entries in the tree 
+  Int_t nevents  = Int_t(sdigitsBranch->GetEntries());
   
-   gAlice->ResetDigits();
-   }
+  Int_t read = 0;
+  // Loop over the events in the 
+  for (Int_t event = 0; event < nevents; event++)  {
+    // Read in entry number `event' 
+    read += sdigitsBranch->GetEntry(event);
+    
+    // Get the number of sdigits 
+    Int_t nsdigits = fmdSDigits->GetEntries ();
+    AliFMDDebug(3, ("Got %5d SDigits", nsdigits));
+    for (Int_t sdigit = 0; sdigit < nsdigits; sdigit++) {
+      // Get the sdigit number `sdigit'
+      AliFMDSDigit* fmdSDigit = 
+       static_cast<AliFMDSDigit*>(fmdSDigits->UncheckedAt(sdigit));
+
+      AliFMDDebug(5, ("Adding contribution of %d tracks", 
+                     fmdSDigit->GetNTrack()));
+      AliFMDDebug(15, ("Contrib from FMD%d%c[%2d,%3d] (%s) from track %d", 
+                      fmdSDigit->Detector(),
+                      fmdSDigit->Ring(),
+                      fmdSDigit->Sector(),
+                      fmdSDigit->Strip(),
+                      fmdSDigit->GetName(), 
+                      fmdSDigit->GetTrack(0)));
+      
+      // Extract parameters 
+      AddContribution(fmdSDigit->Detector(),
+                     fmdSDigit->Ring(),
+                     fmdSDigit->Sector(),
+                     fmdSDigit->Strip(),
+                     fmdSDigit->Edep(), 
+                     kTRUE,
+                     fmdSDigit->GetNTrack(),
+                     fmdSDigit->GetTracks());
+    }  // sdigit loop
+  } // event loop
+
+
+  AliFMDDebug(3, ("Size of cache: %d bytes, read %d bytes", 
+                  sizeof(fEdep), read));
 }
+
+//____________________________________________________________________
+//
+// EOF
+//