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