FMD/AliFMDMultPoisson.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
  18 //____________________________________________________________________
  19 //
  20 // Reconstruct charged particle multiplicity in the FMD
  21 //
  22 // [See also the AliFMDReconstructor class]
  23 //
  24 // This class reconstructs the muliplicity in regions based on the
  25 // ratio of empty to full strips.
  26 //
  27 #include "AliFMD.h"                     // ALIFMD_H
  28 #include "AliFMDGeometry.h"             // ALIFMDGEOMETRY_H
  29 #include "AliFMDDetector.h"             // ALIFMDDETECTOR_H
  30 #include "AliFMDRing.h"                 // ALIFMDRING_H
  31 #include "AliFMDMultPoisson.h"          // ALIFMDMULTPOISSON_H
  32 #include "AliFMDMultRegion.h"           // ALIFMDMULTREGION_H
  33 #include "AliFMDDigit.h"                // ALIFMDDIGIT_H
  34 #include "AliLog.h"                     // ALILOG_H
  35 #include <TClonesArray.h>               // ROOT_TClonesArray
  36 #include <TTree.h>                      // ROOT_TTree
  37
  38 //____________________________________________________________________
  39 ClassImp(AliFMDMultPoisson)
  40 #if 0
  41   ; // This is here to keep Emacs for indenting the next line
  42 #endif
  43
  44 //____________________________________________________________________
  45 AliFMDMultPoisson::AliFMDMultPoisson()
  46   : AliFMDMultAlgorithm("Poisson", "Poisson"),
  47     fDeltaEta(0),
  48     fDeltaPhi(0),
  49     fThreshold(0)
  50 {
  51   SetDeltaEta();
  52   SetDeltaPhi();
  53   SetThreshold();
  54   fMult = new TClonesArray("AliFMDMultRegion", 1000);
  55 }
  56
  57 //____________________________________________________________________
  58 void
  59 AliFMDMultPoisson::PreEvent(TTree* tree, Float_t ipZ)
  60 {
  61   // Reset internal data
  62   AliFMDMultAlgorithm::PreEvent(tree, ipZ);
  63   fCurrentVertexZ = ipZ;
  64   fEmpty.Reset(kFALSE);
  65
  66   // Make a branch in the reconstruction tree.
  67   const Int_t kBufferSize = 16000;
  68   fTreeR->Branch("FMDPoisson", &fMult, kBufferSize);
  69
  70 }
  71
  72 //____________________________________________________________________
  73 void
  74 AliFMDMultPoisson::ProcessDigit(AliFMDDigit*  digit,
  75                                 Float_t       /* eta */,
  76                                 Float_t       /* phi */,
  77                                 UShort_t      count)
  78 {
  79   // Process one digit.
  80   //
  81   // Parameters:
  82   //
  83   //   digit            Digit to process
  84   //   ipZ              Z--coordinate of the primary interaction
  85   //                    vertex of this event
  86   //
  87   if (!digit) return;
  88   if (count < fThreshold) fEmpty(digit->Detector() - 1,
  89                                  digit->Ring(),
  90                                  digit->Sector(),
  91                                  digit->Strip()) = kTRUE;
  92 }
  93
  94 //____________________________________________________________________
  95 void
  96 AliFMDMultPoisson::PostEvent()
  97 {
  98   // Fill the branch
  99   // Based on the information in the cache, do the reconstruction.
 100
 101   // Loop over the detectors
 102   for (Int_t i = 1; i <= 3; i++) {
 103     AliFMDGeometry* fmd = AliFMDGeometry::Instance();
 104     AliFMDDetector* sub = fmd->GetDetector(i);
 105     if (!sub) continue;
 106
 107     // Loop over the rings in the detector
 108     for (Int_t j = 0; j < 2; j++) {
 109       AliFMDRing* r  = sub->GetRing((j == 0 ? 'I' : 'O'));
 110       Float_t     rZ = sub->GetRingZ((j == 0 ? 'I' : 'O'));
 111       if (!r) continue;
 112
 113       // Calculate low/high theta and eta
 114       // FIXME: Is this right?
 115       Float_t realZ    = fCurrentVertexZ + rZ;
 116       Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
 117       Float_t thetaIn  = TMath::ATan2(r->GetLowR(), realZ);
 118       Float_t etaOut   = - TMath::Log(TMath::Tan(thetaOut / 2));
 119       Float_t etaIn    = - TMath::Log(TMath::Tan(thetaIn / 2));
 120       if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
 121         Float_t tmp = etaIn;
 122         etaIn       = etaOut;
 123         etaOut      = tmp;
 124       }
 125
 126       //-------------------------------------------------------------
 127       //
 128       // Here starts poisson method
 129       //
 130       // Calculate eta step per strip, number of eta steps, number of
 131       // phi steps, and check the sign of the eta increment
 132       Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
 133       Int_t   nEta     = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta);
 134       Int_t   nPhi     = Int_t(360. / fDeltaPhi);
 135       Float_t sign     = TMath::Sign(Float_t(1.), etaIn);
 136
 137       AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
 138                         sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
 139
 140       // Loop over relevant phi values
 141       for (Int_t p = 0; p < nPhi; p++) {
 142         Float_t  minPhi    = p * fDeltaPhi;
 143         Float_t  maxPhi    = minPhi + fDeltaPhi;
 144         UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
 145         UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
 146
 147         AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
 148                           minPhi, maxPhi, minSector, maxSector));
 149         // Loop over relevant eta values
 150         for (Int_t e = nEta; e >= 0; --e) {
 151           Float_t  maxEta   = etaIn  - sign * e * fDeltaEta;
 152           Float_t  minEta   = maxEta - sign * fDeltaEta;
 153           if (sign > 0)  minEta = TMath::Max(minEta, etaOut);
 154           else           minEta = TMath::Min(minEta, etaOut);
 155           Float_t  theta1   = 2 * TMath::ATan(TMath::Exp(-minEta));
 156           Float_t  theta2   = 2 * TMath::ATan(TMath::Exp(-maxEta));
 157           Float_t  minR     = TMath::Abs(realZ * TMath::Tan(theta2));
 158           Float_t  maxR     = TMath::Abs(realZ * TMath::Tan(theta1));
 159           UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
 160           UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
 161
 162           AliDebug(10, Form("  Now in eta range %f, %f (strips %d, %d)\n"
 163                             "    [radii %f, %f, thetas %f, %f, sign %d]",
 164                             minEta, maxEta, minStrip, maxStrip,
 165                             minR, maxR, theta1, theta2, sign));
 166
 167           // Count number of empty strips
 168           Int_t   emptyStrips = 0;
 169           for (Int_t sector = minSector; sector < maxSector; sector++)
 170             for (Int_t strip = minStrip; strip < maxStrip; strip++)
 171               if (fEmpty(sub->GetId() - 1, r->GetId(), sector, strip))
 172                 emptyStrips++;
 173
 174           // The total number of strips
 175           Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
 176
 177           // Log ratio of empty to total number of strips
 178           AliDebug(10, Form("Lambda= %d / %d = %f",
 179                             emptyStrips, nTotal,
 180                             Float_t(emptyStrips) / nTotal));
 181
 182           Double_t lambda = (emptyStrips > 0 ?
 183                              - TMath::Log(Double_t(emptyStrips) / nTotal) :
 184                              1);
 185
 186           // The reconstructed number of particles is then given by
 187           Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
 188
 189           // Add a AliFMDMultRegion to the reconstruction tree.
 190           AliFMDMultRegion* m = new((*fMult)[fNMult])
 191             AliFMDMultRegion(sub->GetId(), r->GetId(),
 192                              minSector, maxSector, minStrip, maxStrip,
 193                              minEta, maxEta, minPhi, maxPhi,
 194                              reconstructed, AliFMDMultRegion::kPoission);
 195           (void)m;
 196           fNMult++;
 197         } // phi
 198       } // eta
 199     } // ring
 200   } // detector
 201 }
 202
 203
 204 //____________________________________________________________________
 205 //
 206 // EOF
 207 //