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 <TCanvas.h>
  20 #include <TStyle.h>
  21 #include <TROOT.h>
  22 #include <TFile.h>
  23 #include <iostream>
  24
  25 /** @class Poisson
  26     @brief Make a poisson reconstruction
  27     @code
  28     Root> .L Compile.C
  29     Root> Compile("Poisson.C")
  30     Root> Poisson c
  31     Root> c.Run();
  32     @endcode
  33     @ingroup FMD_script
  34  */
  35 class Poisson : public AliFMDInput
  36 {
  37 protected:
  38   TH2D*  fEmpty; // Histogram
  39   TH2D*  fTotal; // Histogram
  40   TH2D*  fMult;  // Histogram
  41   TFile* fFile;  // File
  42   Int_t  fEv;    // Event number
  43   Double_t fThreshold;
  44 public:
  45   /** Constructor
  46       @param threshold Threshold
  47       @param nEta      # of @f$ \eta@f$ bins
  48       @param minEta    minimum @f$ \eta@f$
  49       @param maxEta    maximum @f$ \eta@f$
  50       @param nPhi      # of @f$ \eta@f$ bins
  51       @param minPhi    minimum @f$ \varphi@f$
  52       @param maxPhi    maximum @f$ \varphi@f$ */
  53   Poisson(Double_t threshold=.3,
  54           Int_t nEta=120, Float_t minEta=-6, Float_t maxEta=6,
  55           Int_t nPhi=4,   Float_t minPhi=0,  Float_t maxPhi=2*TMath::Pi())
  56     : fFile(0), fEv(0), fThreshold(threshold)
  57   {
  58     AddLoad(kESD);
  59
  60     fEmpty = new TH2D("empty", "# of empty strips", nEta, minEta, maxEta,
  61                       nPhi, minPhi, maxPhi);
  62     fTotal = new TH2D("total", "Total # of strips", nEta, minEta, maxEta,
  63                       nPhi, minPhi, maxPhi);
  64     fMult = new TH2D("mult", "Multiplicity", nEta, minEta, maxEta,
  65                       nPhi, minPhi, maxPhi);
  66     fEmpty->SetXTitle("#eta");
  67     fEmpty->SetYTitle("#phi");
  68     fEmpty->SetZTitle("N");
  69     fTotal->SetXTitle("#eta");
  70     fTotal->SetYTitle("#phi");
  71     fTotal->SetZTitle("N");
  72     fMult->SetXTitle("#eta");
  73     fMult->SetYTitle("#phi");
  74     fMult->SetZTitle("<M_{ch}>");
  75
  76   }
  77   /** Initialize the analyser. Opens the output file.
  78       @return @c true on success. */
  79   virtual Bool_t Init()
  80   {
  81     if (!AliFMDInput::Init()) return kFALSE;
  82     fFile = TFile::Open("poisson.root", "RECREATE");
  83     if (!fFile) return kFALSE;
  84     return kTRUE;
  85   }
  86   /** Begining of event
  87       @param event Event number
  88       @return @c false on error */
  89   virtual Bool_t Begin(Int_t event)
  90   {
  91     if (!AliFMDInput::Begin(event)) return kFALSE;
  92     fEv = event;
  93     fEmpty->Clear();
  94     fTotal->Clear();
  95     fMult->Clear();
  96     return kTRUE;
  97   }
  98   /** Process ESD data.  For each strip, check if the
  99       psuedo-multiplicity is less than the threshold.  If it is, then
 100       count the strip as empty.
 101       @param esd ESD data
 102       @return @c true on success. */
 103   virtual Bool_t ProcessESD(AliESDFMD* esd)
 104   {
 105     for (UShort_t det = 1; det <= 3; det++) {
 106       for (UShort_t rng = 0; rng < 2; rng++) {
 107         Char_t ring = (rng == 0 ? 'I' : 'O');
 108         // Not covered channels
 109         for (UShort_t sec = 0; sec < 40; sec++) {
 110           for (UShort_t str = 0; str < 512; str++) {
 111             Float_t mult = esd->Multiplicity(det, ring, sec, str);
 112             Float_t eta  = esd->Eta(det, ring, sec, str);
 113             // Dead channels, or not covered.
 114             if (mult >= AliESDFMD::kInvalidMult) continue;
 115             if (eta  >= AliESDFMD::kInvalidEta) continue;
 116             Float_t phi;
 117             switch (ring) {
 118             case 'I':  phi = (sec + .5) * 2 * TMath::Pi() / 20; break;
 119             case 'O':  phi = (sec + .5) * 2 * TMath::Pi() / 40; break;
 120             }
 121             fTotal->Fill(eta, phi);
 122             if (mult < fThreshold) fEmpty->Fill(eta, phi);
 123           } // Loop over strips
 124         } // Loop over sectors
 125       } // Loop over rings
 126     } // Loop over detectors
 127     return kTRUE;
 128   }
 129   /** For each bin, reconstruct the charge particle multiplicity as
 130       @f[
 131       m = - N_{total} \log\left(\frac{N_{empty}}{N_{total}}\right)
 132       @f]
 133       where @f$ N_{total}@f$ is the total number of strips in the bin,
 134       and @f$ N_{empty}@f$ is the number of strips in the bin that did
 135       not fire.
 136       @return @c true  */
 137   virtual Bool_t End()
 138   {
 139     for (Int_t etaBin = 1; etaBin <= fEmpty->GetNbinsX(); etaBin++) {
 140       for (Int_t phiBin = 1; phiBin <= fEmpty->GetNbinsY(); phiBin++) {
 141         Double_t empty  = fEmpty->GetBinContent(etaBin, phiBin);
 142         Double_t total  = fTotal->GetBinContent(etaBin, phiBin);
 143         Double_t lambda = (empty > 0 ? - TMath::Log(empty / total) : 1);
 144         Double_t mult   = lambda * total;
 145         fMult->SetBinContent(etaBin, phiBin, mult);
 146       }
 147     }
 148     fFile->cd();
 149     fMult->Write(Form("mult%03d", fEv));
 150     if (!gROOT->IsBatch()) {
 151       gStyle->SetPalette(1);
 152       TCanvas* c = new TCanvas("poisson", "Poisson multiplicity");
 153       c->SetFillColor(0);
 154       fMult->Draw("colz");
 155     }
 156     return AliFMDInput::End();
 157   }
 158
 159   /** At end of run.  Write and close output file @c poisson.root
 160       @return @c true on success */
 161   virtual Bool_t Finish()
 162   {
 163     fFile->Write();
 164     fFile->Close();
 165     fFile = 0;
 166     return AliFMDInput::Finish();
 167   }
 168   ClassDef(Poisson,0);
 169 };
 170
 171 //____________________________________________________________________
 172 //
 173 // EOF
 174 //