FMD/scripts/Poisson.C

   1 //____________________________________________________________________
   2 //
   3 // $Id$
   4 //
   5 // Script that contains a class to draw hits, using the
   6 // AliFMDInputHits class in the util library.
   7 //
   8 // It draws the energy loss versus the p/(mq^2).  It can be overlayed
   9 // with the Bethe-Bloc curve to show how the simulation behaves
  10 // relative to the expected.
  11 //
  12 // Use the script `Compile.C' to compile this class using ACLic.
  13 //
  14 #include <AliESDFMD.h>
  15 #include <AliFMDInput.h>
  16 #include <TH2D.h>
  17 #include <TStyle.h>
  18 #include <TArrayF.h>
  19 #include <iostream>
  20
  21 /** @class Poisson
  22     @brief Make a poisson reconstruction
  23     @code
  24     Root> .L Compile.C
  25     Root> Compile("Poisson.C")
  26     Root> Poisson c
  27     Root> c.Run();
  28     @endcode
  29     @ingroup FMD_script
  30  */
  31 class Poisson : public AliFMDInput
  32 {
  33 private:
  34   TH2D*  fEmpty; // Histogram
  35   TH2D*  fTotal; // Histogram
  36   TH2D*  fMult;  // Histogram
  37   TFile* fFile;  // File
  38   Int_t  fEv;    // Event number
  39 public:
  40   Poisson(Double_t threshold=.3,
  41           Int_t nEta=120, Float_t minEta=-6, Float_t maxEta=6,
  42           Int_t nPhi=4,   Float_t minPhi=0,  Float_t maxPhi=2*TMath::Pi())
  43     : fFile(0), fEv(0)
  44   {
  45     AddLoad(kESD);
  46
  47     fEmpty = new TH2D("empty", "# of empty strips", nEta, minEta, maxEta,
  48                       nPhi, minPhi, maxPhi);
  49     fTotal = new TH2D("total", "Total # of strips", nEta, minEta, maxEta,
  50                       nPhi, minPhi, maxPhi);
  51     fMult = new TH2D("mult", "Multiplicity", nEta, minEta, maxEta,
  52                       nPhi, minPhi, maxPhi);
  53     fEmpty->SetXTitle("#eta");
  54     fEmpty->SetYTitle("#phi");
  55     fEmpty->SetZTitle("N");
  56     fTotal->SetXTitle("#eta");
  57     fTotal->SetYTitle("#phi");
  58     fTotal->SetZTitle("N");
  59     fMult->SetXTitle("#eta");
  60     fMult->SetYTitle("#phi");
  61     fMult->SetZTitle("<M_{ch}>");
  62
  63   }
  64   Bool_t Init()
  65   {
  66     if (!AliFMDInput::Begin(event)) return kFALSE;
  67     fFile = TFile::Open("poisson.root");
  68     if (!fFile) return kFALSE;
  69     return kTRUE;
  70   }
  71   /** Begining of event
  72       @param event Event number
  73       @return @c false on error */
  74   Bool_t Begin(Int_t event)
  75   {
  76     if (!AliFMDInput::Begin(event)) return kFALSE;
  77     fEv = event;
  78     fEmpty->Clear();
  79     fTotal->Clear();
  80     fMult->Clear();
  81     return kTRUE;
  82   }
  83   Bool_t ProcessESD(AliESDFMDHit* esd)
  84   {
  85     for (UShort_t det = 1; det <= esd->MaxDetector(); det++) {
  86       for (UShort_t rng = 0; rng < esd->MaxRing(); rng++) {
  87         Char_t ring = (rng == 0 ? 'I' : 'O');
  88         // Not covered channels
  89         for (UShort_t sec = 0; sec < esd->MaxSector(); sec++) {
  90           for (UShort_t str = 0; str < esd->MaxStrip(); str++) {
  91             Float_t mult = esd->Multiplicity(det, ring, sec, str);
  92             Float_t eta  = esd->Eta(det, ring, sec, str);
  93             // Dead channels, or not covered.
  94             if (mult >= AliESDFMD::kInvalidMult) continue;
  95             if (esd  >= AliESDFMD::kInvalidEta) continue;
  96             Float_t phi;
  97             switch (ring) {
  98             case 'I':  phi = (sec + .5) * 2 * TMath::Pi() / 20; break;
  99             case 'O':  phi = (sec + .5) * 2 * TMath::Pi() / 40; break;
 100             }
 101             fTotal->Fill(eta, phi);
 102             if (mult < threshold) fEmpty->Fill(eta, phi);
 103           } // Loop over strips
 104         } // Loop over sectors
 105       } // Loop over rings
 106     } // Loop over detectors
 107   }
 108   /** For each bin, reconstruct the charge particle multiplicity as
 109       @f[
 110       m = - N_{total} \log\left(\frac{N_{empty}}{N_{total}}\right)
 111       @f]
 112       where @f$ N_{total}@f$ is the total number of strips in the bin,
 113       and @f$ N_{empty}@f$ is the number of strips in the bin that did
 114       not fire.
 115       @return @c true  */
 116   Bool_t End()
 117   {
 118     for (Int_t etaBin = 1; etaBin <= fEmpty->GetNbinsX(); etaBin++) {
 119       for (Int_t phiBin = 1; phiBin <= fEmpty->GetNbinsY(); phiBin++) {
 120         Double_t empty  = fEmpty->GetBinContent(etaBin, phiBin);
 121         Double_t total  = fTotal->GetBinContent(etaBin, phiBin);
 122         Double_t lambda = (empty > 0 ? - TMath::Log(empty / nTotal) : 1);
 123         Double_t mult   = lambda * nTotal;
 124         fMult->SetBinContent(etaBin, phiBin, mult);
 125       }
 126     }
 127     fFile->cd();
 128     fMult->Write(Form("mult%03d", fEv));
 129     return AliFMDInput::End();
 130   }
 131
 132   Bool_t Finish()
 133   {
 134     fFile->Write();
 135     fFile->Close();
 136     fFile = 0;
 137     return AliFMDInput::Finish();
 138   }
 139   ClassDef(Poisson,0);
 140 };
 141
 142 //____________________________________________________________________
 143 //
 144 // EOF
 145 //