PWG2/FORWARD/analysis2/AliFMDDensityCalculator.cxx

   1 #include "AliFMDDensityCalculator.h"
   2 #include <AliESDFMD.h>
   3 #include <TAxis.h>
   4 #include <TList.h>
   5 #include <TMath.h>
   6 #include "AliFMDAnaParameters.h"
   7 #include "AliLog.h"
   8 #include <TH2D.h>
   9
  10 ClassImp(AliFMDDensityCalculator)
  11 #if 0
  12 ; // For Emacs
  13 #endif
  14
  15 //____________________________________________________________________
  16 AliFMDDensityCalculator::AliFMDDensityCalculator()
  17   : TNamed(),
  18     fRingHistos(),
  19     fMultCut(0.3)
  20 {}
  21
  22 //____________________________________________________________________
  23 AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
  24   : TNamed("fmdDensityCalculator", title),
  25     fRingHistos(),
  26     fMultCut(0.3)
  27 {
  28   fRingHistos.SetName(GetName());
  29   fRingHistos.Add(new RingHistos(1, 'I'));
  30   fRingHistos.Add(new RingHistos(2, 'I'));
  31   fRingHistos.Add(new RingHistos(2, 'O'));
  32   fRingHistos.Add(new RingHistos(3, 'I'));
  33   fRingHistos.Add(new RingHistos(3, 'O'));
  34 }
  35
  36 //____________________________________________________________________
  37 AliFMDDensityCalculator::AliFMDDensityCalculator(const
  38                                                  AliFMDDensityCalculator& o)
  39   : TNamed(o),
  40     fRingHistos(),
  41     fMultCut(o.fMultCut)
  42 {
  43   TIter    next(&o.fRingHistos);
  44   TObject* obj = 0;
  45   while ((obj = next())) fRingHistos.Add(obj);
  46 }
  47
  48 //____________________________________________________________________
  49 AliFMDDensityCalculator::~AliFMDDensityCalculator()
  50 {
  51   fRingHistos.Delete();
  52 }
  53
  54 //____________________________________________________________________
  55 AliFMDDensityCalculator&
  56 AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
  57 {
  58   SetName(o.GetName());
  59   SetTitle(o.GetTitle());
  60
  61   fMultCut = o.fMultCut;
  62
  63   fRingHistos.Delete();
  64   TIter    next(&o.fRingHistos);
  65   TObject* obj = 0;
  66   while ((obj = next())) fRingHistos.Add(obj);
  67
  68   return *this;
  69 }
  70
  71 //____________________________________________________________________
  72 AliFMDDensityCalculator::RingHistos*
  73 AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
  74 {
  75   Int_t idx = -1;
  76   switch (d) {
  77   case 1: idx = 0; break;
  78   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
  79   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
  80   }
  81   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
  82
  83   return static_cast<RingHistos*>(fRingHistos.At(idx));
  84 }
  85
  86 //____________________________________________________________________
  87 Bool_t
  88 AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
  89                                    AliForwardUtil::Histos& hists,
  90                                    Int_t                   vtxbin,
  91                                    Bool_t                  lowFlux)
  92 {
  93   for (UShort_t d=1; d<=3; d++) {
  94     UShort_t nr = (d == 1 ? 1 : 2);
  95     for (UShort_t q=0; q<nr; q++) {
  96       Char_t      r = (q == 0 ? 'I' : 'O');
  97       UShort_t    ns= (q == 0 ?  20 :  40);
  98       UShort_t    nt= (q == 0 ? 512 : 256);
  99       TH2D*       h = hists.Get(d,r);
 100       RingHistos* rh= GetRingHistos(d,r);
 101
 102       for (UShort_t s=0; s<ns; s++) {
 103         for (UShort_t t=0; t<nt; t++) {
 104           Float_t mult = fmd.Multiplicity(d,r,s,t);
 105
 106           if (mult == 0 || mult > 20) continue;
 107
 108           Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
 109           Float_t eta = fmd.Eta(d,r,s,t);
 110
 111           Float_t n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
 112           rh->fEvsN->Fill(mult,n);
 113
 114           Float_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
 115           if (c > 0) n /= c;
 116           rh->fEvsM->Fill(mult,n);
 117
 118           h->Fill(eta,phi,n);
 119           rh->fDensity->Fill(eta,phi,n);
 120         } // for t
 121       } // for s
 122     } // for q
 123   } // for d
 124
 125   return kTRUE;
 126 }
 127
 128 //_____________________________________________________________________
 129 Float_t
 130 AliFMDDensityCalculator::NParticles(Float_t  mult,
 131                                     UShort_t d,
 132                                     Char_t   r,
 133                                     UShort_t /*s*/,
 134                                     UShort_t /*t*/,
 135                                     Int_t    /*v*/,
 136                                     Float_t  eta,
 137                                     Bool_t   lowFlux) const
 138 {
 139   if (mult <= fMultCut) return 0;
 140   if (lowFlux) return 1;
 141
 142   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 143
 144   Float_t mpv  = pars->GetMPV(d,r,eta);
 145   Float_t w    = pars->GetSigma(d,r,eta);
 146   Float_t w2   = pars->Get2MIPWeight(d,r,eta);
 147   Float_t w3   = pars->Get3MIPWeight(d,r,eta);
 148   Float_t mpv2 = 2*mpv+2*w*TMath::Log(2);
 149   Float_t mpv3 = 3*mpv+3*w*TMath::Log(3);
 150
 151   Float_t sum  = (TMath::Landau(mult,mpv,w,kTRUE) +
 152                   w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) +
 153                   w3  * TMath::Landau(mult,mpv3,3*w,kTRUE));
 154   Float_t wsum = (TMath::Landau(mult,mpv,w,kTRUE) +
 155                   2*w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) +
 156                   3*w3 * TMath::Landau(mult,mpv3,3*w,kTRUE));
 157
 158   return (sum > 0) ? wsum / sum : 1;
 159 }
 160
 161 //_____________________________________________________________________
 162 Float_t
 163 AliFMDDensityCalculator::Correction(UShort_t d, Char_t r, UShort_t /*s*/,
 164                                     UShort_t t, Int_t v, Float_t eta,
 165                                     Bool_t lowFlux) const
 166 {
 167   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 168
 169   Float_t correction = AcceptanceCorrection(r,t);
 170   if (lowFlux) {
 171     TH1F* dblHitCor = pars->GetDoubleHitCorrection(d,r);
 172     if (dblHitCor) {
 173       Float_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
 174       if (dblC > 0)
 175         correction *= dblC;
 176     }
 177     else
 178       AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
 179   }
 180
 181   TH1F* deadCor = pars->GetFMDDeadCorrection(v);
 182   if (deadCor) {
 183     Float_t deadC = deadCor->GetBinContent(deadCor->FindBin(eta));
 184     if (deadC > 0)
 185       correction *= deadC;
 186   }
 187
 188   return correction;
 189 }
 190
 191
 192 //_____________________________________________________________________
 193 Float_t
 194 AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
 195 {
 196   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 197
 198   //AliFMDRing fmdring(ring);
 199   //fmdring.Init();
 200   Float_t   rad       = pars->GetMaxR(r)-pars->GetMinR(r);
 201   Float_t   slen      = pars->GetStripLength(r,t);
 202   Float_t   sblen     = pars->GetBaseStripLength(r,t);
 203   Float_t   nstrips   = (r == 'I' ? 512 : 256);
 204   Float_t   segment   = rad / nstrips;
 205   Float_t   radius    = pars->GetMinR(r) + segment*t;
 206
 207   Float_t   basearea1 = 0.5*sblen*TMath::Power(radius,2);
 208   Float_t   basearea2 = 0.5*sblen*TMath::Power((radius-segment),2);
 209   Float_t   basearea  = basearea1 - basearea2;
 210
 211   Float_t   area1     = 0.5*slen*TMath::Power(radius,2);
 212   Float_t   area2     = 0.5*slen*TMath::Power((radius-segment),2);
 213   Float_t   area      = area1 - area2;
 214
 215   return area/basearea;
 216 }
 217
 218 //____________________________________________________________________
 219 void
 220 AliFMDDensityCalculator::ScaleHistograms(Int_t nEvents)
 221 {
 222   if (nEvents <= 0) return;
 223
 224   TIter    next(&fRingHistos);
 225   RingHistos* o = 0;
 226   while ((o = static_cast<RingHistos*>(next()))) {
 227     o->fDensity->Scale(1. / nEvents);
 228   }
 229 }
 230
 231 //____________________________________________________________________
 232 void
 233 AliFMDDensityCalculator::Output(TList* dir)
 234 {
 235   TList* d = new TList;
 236   d->SetName(GetName());
 237   dir->Add(d);
 238   TIter    next(&fRingHistos);
 239   RingHistos* o = 0;
 240   while ((o = static_cast<RingHistos*>(next()))) {
 241     o->Output(d);
 242   }
 243 }
 244
 245 //====================================================================
 246 AliFMDDensityCalculator::RingHistos::RingHistos()
 247   : fDet(0),
 248     fRing('\0'),
 249     fEvsN(0),
 250     fEvsM(0),
 251     fDensity(0)
 252 {}
 253
 254 //____________________________________________________________________
 255 AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
 256   : fDet(d),
 257     fRing(r),
 258     fEvsN(0),
 259     fEvsM(0),
 260     fDensity(0)
 261 {
 262   fEvsN = new TH2D(Form("FMD%d%c_Eloss_N_nocorr", d, r),
 263                    Form("Energy loss vs uncorrected inclusive "
 264                         "N_{ch} in FMD%d%c", d, r),
 265                    100, -.5, 24.5, 100, -.5, 24.5);
 266   fEvsM = new TH2D(Form("FMD%d%c_Eloss_N_corr", d, r),
 267                    Form("Energy loss vs corrected inclusive N_{ch} in FMD%d%c",
 268                         d, r), 100, -.5, 24.5, 100, -.5, 24.5);
 269   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
 270   fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
 271   fEvsN->Sumw2();
 272   fEvsN->SetDirectory(0);
 273   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
 274   fEvsN->SetYTitle("Inclusive N_{ch} (corrected)");
 275   fEvsM->Sumw2();
 276   fEvsM->SetDirectory(0);
 277
 278   fDensity = new TH2D(Form("FMD%d%c_Incl_Density", d, r),
 279                       Form("in FMD%d%c", d, r),
 280                       200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
 281                       0, 2*TMath::Pi());
 282   fDensity->SetDirectory(0);
 283   fDensity->SetXTitle("#eta");
 284   fDensity->SetYTitle("#phi [radians]");
 285   fDensity->SetZTitle("Inclusive N_{ch} density");
 286 }
 287 //____________________________________________________________________
 288 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
 289   : TObject(o),
 290     fDet(o.fDet),
 291     fRing(o.fRing),
 292     fEvsN(o.fEvsN),
 293     fEvsM(o.fEvsM),
 294     fDensity(o.fDensity)
 295 {}
 296
 297 //____________________________________________________________________
 298 AliFMDDensityCalculator::RingHistos&
 299 AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
 300 {
 301   fDet  = o.fDet;
 302   fRing = o.fRing;
 303
 304   if (fEvsN)    delete  fEvsN;
 305   if (fEvsM)    delete  fEvsM;
 306   if (fDensity) delete fDensity;
 307
 308   fEvsN    = static_cast<TH2D*>(o.fEvsN->Clone());
 309   fEvsM    = static_cast<TH2D*>(o.fEvsM->Clone());
 310   fDensity = static_cast<TH2D*>(o.fDensity->Clone());
 311
 312   return *this;
 313 }
 314 //____________________________________________________________________
 315 AliFMDDensityCalculator::RingHistos::~RingHistos()
 316 {
 317   if (fEvsN)    delete fEvsN;
 318   if (fEvsM)    delete fEvsM;
 319   if (fDensity) delete fDensity;
 320 }
 321
 322 //____________________________________________________________________
 323 void
 324 AliFMDDensityCalculator::RingHistos::Output(TList* dir)
 325 {
 326   TList* d = new TList;
 327   d->SetName(Form("FMD%d%c", fDet, fRing));
 328   d->Add(fEvsN);
 329   d->Add(fEvsM);
 330   d->Add(fDensity);
 331   dir->Add(d);
 332 }
 333
 334 //____________________________________________________________________
 335 //
 336 // EOF
 337 //
 338
 339
 340