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 : public AliFMDInput
  22 {
  23 private:
  24   TH2D*  fEmpty; // Histogram
  25   TH2D*  fTotal; // Histogram
  26   TH2D*  fMult;  // Histogram
  27   TFile* fFile;  // File
  28   Int_t  fEv;    // Event number
  29 public:
  30   Poisson(Double_t threshold=.3,
  31           Int_t nEta=120, Float_t minEta=-6, Float_t maxEta=6,
  32           Int_t nPhi=4,   Float_t minPhi=0,  Float_t maxPhi=2*TMath::Pi())
  33     : fFile(0), fEv(0)
  34   {
  35     AddLoad(kESD);
  36
  37     fEmpty = new TH2D("empty", "# of empty strips", nEta, minEta, maxEta,
  38                       nPhi, minPhi, maxPhi);
  39     fTotal = new TH2D("total", "Total # of strips", nEta, minEta, maxEta,
  40                       nPhi, minPhi, maxPhi);
  41     fMult = new TH2D("mult", "Multiplicity", nEta, minEta, maxEta,
  42                       nPhi, minPhi, maxPhi);
  43     fEmpty->SetXTitle("#eta");
  44     fEmpty->SetYTitle("#phi");
  45     fEmpty->SetZTitle("N");
  46     fTotal->SetXTitle("#eta");
  47     fTotal->SetYTitle("#phi");
  48     fTotal->SetZTitle("N");
  49     fMult->SetXTitle("#eta");
  50     fMult->SetYTitle("#phi");
  51     fMult->SetZTitle("<M_{ch}>");
  52
  53   }
  54   Bool_t Init()
  55   {
  56     if (!AliFMDInput::Begin(event)) return kFALSE;
  57     fFile = TFile::Open("poisson.root");
  58     if (!fFile) return kFALSE;
  59     return kTRUE;
  60   }
  61   Bool_t Begin(Int_t event)
  62   {
  63     if (!AliFMDInput::Begin(event)) return kFALSE;
  64     fEv = event;
  65     fEmpty->Clear();
  66     fTotal->Clear();
  67     fMult->Clear();
  68     return kTRUE;
  69   }
  70   Bool_t ProcessESD(AliESDFMDHit* esd)
  71   {
  72     for (UShort_t det = 1; det <= esd->MaxDetector(); det++) {
  73       for (UShort_t rng = 0; rng < esd->MaxRing(); rng++) {
  74         Char_t ring = (rng == 0 ? 'I' : 'O');
  75         // Not covered channels
  76         for (UShort_t sec = 0; sec < esd->MaxSector(); sec++) {
  77           for (UShort_t str = 0; str < esd->MaxStrip(); str++) {
  78             Float_t mult = esd->Multiplicity(det, ring, sec, str);
  79             Float_t eta  = esd->Eta(det, ring, sec, str);
  80             // Dead channels, or not covered.
  81             if (mult >= AliESDFMD::kInvalidMult) continue;
  82             if (esd  >= AliESDFMD::kInvalidEta) continue;
  83             Float_t phi;
  84             switch (ring) {
  85             case 'I':  phi = (sec + .5) * 2 * TMath::Pi() / 20; break;
  86             case 'O':  phi = (sec + .5) * 2 * TMath::Pi() / 40; break;
  87             }
  88             fTotal->Fill(eta, phi);
  89             if (mult < threshold) fEmpty->Fill(eta, phi);
  90           } // Loop over strips
  91         } // Loop over sectors
  92       } // Loop over rings
  93     } // Loop over detectors
  94   }
  95   Bool_t End()
  96   {
  97     for (Int_t etaBin = 1; etaBin <= fEmpty->GetNbinsX(); etaBin++) {
  98       for (Int_t phiBin = 1; phiBin <= fEmpty->GetNbinsY(); phiBin++) {
  99         Double_t empty  = fEmpty->GetBinContent(etaBin, phiBin);
 100         Double_t total  = fTotal->GetBinContent(etaBin, phiBin);
 101         Double_t lambda = (empty > 0 ? - TMath::Log(empty / nTotal) : 1);
 102         Double_t mult   = lambda * nTotal;
 103         fMult->SetBinContent(etaBin, phiBin, mult);
 104       }
 105     }
 106     fFile->cd();
 107     fMult->Write(Form("mult%03d", fEv));
 108     return AliFMDInput::End();
 109   }
 110
 111   Bool_t Finish()
 112   {
 113     fFile->Write();
 114     fFile->Close();
 115     fFile = 0;
 116     return AliFMDInput::Finish();
 117   }
 118   ClassDef(Poisson,0);
 119 };
 120
 121 //____________________________________________________________________
 122 //
 123 // EOF
 124 //