PWGLF/FORWARD/analysis2/AliFMDESDFixer.cxx

   1 #include "AliFMDESDFixer.h"
   2 #include "AliESDFMD.h"
   3 #include "AliFMDStripIndex.h"
   4 #include "AliForwardUtil.h"
   5 #include "AliForwardCorrectionManager.h"
   6 #include "AliFMDCorrNoiseGain.h"
   7 #include <TH1.h>
   8 #include <TList.h>
   9 #include <TObjArray.h>
  10 #include <TROOT.h>
  11 #include <TMath.h>
  12 #include <iostream>
  13 #include <iomanip>
  14
  15
  16 //____________________________________________________________________
  17 AliFMDESDFixer::AliFMDESDFixer()
  18   : TObject(),
  19     fRecoFactor(1),
  20     fMaxNoiseCorr(0.05),
  21     fRecalculateEta(true),
  22     fXtraDead(0),
  23     fHasXtraDead(false),
  24     fInvalidIsEmpty(false),
  25     fNoiseChange(0),
  26     fEtaChange(0),
  27     fDeadChange(0)
  28 {}
  29
  30 //____________________________________________________________________
  31 AliFMDESDFixer::AliFMDESDFixer(const char*)
  32   : TObject(),
  33     fRecoFactor(1),
  34     fMaxNoiseCorr(0.05),
  35     fRecalculateEta(false),
  36     fXtraDead(AliFMDStripIndex::Pack(3,'O',19,511)+1),
  37     fHasXtraDead(false),
  38     fInvalidIsEmpty(false),
  39     fNoiseChange(0),
  40     fEtaChange(0),
  41     fDeadChange(0)
  42 {}
  43
  44
  45 //____________________________________________________________________
  46 AliFMDESDFixer::AliFMDESDFixer(const AliFMDESDFixer& o)
  47   : TObject(),
  48     fRecoFactor(o.fRecoFactor),
  49     fMaxNoiseCorr(o.fMaxNoiseCorr),
  50     fRecalculateEta(o.fRecalculateEta),
  51     fXtraDead(o.fXtraDead),
  52     fHasXtraDead(o.fHasXtraDead),
  53     fInvalidIsEmpty(o.fInvalidIsEmpty),
  54     fNoiseChange(0),
  55     fEtaChange(0),
  56     fDeadChange(0)
  57 {}
  58
  59 //____________________________________________________________________
  60 AliFMDESDFixer&
  61 AliFMDESDFixer::operator=(const AliFMDESDFixer& o)
  62 {
  63   if (&o == this) return *this;
  64
  65   fRecoFactor = o.fRecoFactor;
  66   fMaxNoiseCorr = o.fMaxNoiseCorr;
  67   fRecalculateEta = o.fRecalculateEta;
  68   fXtraDead = o.fXtraDead;
  69   fHasXtraDead = o.fHasXtraDead;
  70   fInvalidIsEmpty = o.fInvalidIsEmpty;
  71   fNoiseChange = 0;
  72   fEtaChange  = 0;
  73   fDeadChange = 0;
  74
  75   return *this;
  76 };
  77
  78 //____________________________________________________________________
  79 void
  80 AliFMDESDFixer::CreateOutputObjects(TList* l)
  81 {
  82   TList* d = new TList;
  83   d->SetName(GetName());
  84   l->Add(d);
  85
  86   d->Add(AliForwardUtil::MakeParameter("recoFactor", fRecoFactor));
  87   d->Add(AliForwardUtil::MakeParameter("recalcEta",  fRecalculateEta));
  88   d->Add(AliForwardUtil::MakeParameter("invalidIsEmpty",  fInvalidIsEmpty));
  89
  90   fNoiseChange = new TH1D("noiseChange", "#delta#Delta due to noise",30,0,.3);
  91   fNoiseChange->SetDirectory(0);
  92   fNoiseChange->SetFillColor(kYellow+1);
  93   fNoiseChange->SetFillStyle(3001);
  94   d->Add(fNoiseChange);
  95
  96   fEtaChange = new TH1D("etaChange", "#delta#eta",40,-1,1);
  97   fEtaChange->SetDirectory(0);
  98   fEtaChange->SetFillColor(kCyan+1);
  99   fEtaChange->SetFillStyle(3001);
 100   d->Add(fEtaChange);
 101
 102   fDeadChange = new TH1D("deadChange", "#deltaN_{dead}",100,0,51200);
 103   fDeadChange->SetDirectory(0);
 104   fDeadChange->SetFillColor(kMagenta+1);
 105   fDeadChange->SetFillStyle(3001);
 106   d->Add(fDeadChange);
 107
 108   TObjArray* extraDead = new TObjArray;
 109   extraDead->SetOwner();
 110   extraDead->SetName("extraDead");
 111   fXtraDead.Compact();
 112   UInt_t firstBit = fXtraDead.FirstSetBit();
 113   UInt_t nBits    = fXtraDead.GetNbits();
 114   for (UInt_t i = firstBit; i < nBits; i++) {
 115     if (!fXtraDead.TestBitNumber(i)) continue;
 116     UShort_t dd, s, t;
 117     Char_t   r;
 118     AliFMDStripIndex::Unpack(i, dd, r, s, t);
 119     extraDead->Add(AliForwardUtil::MakeParameter(Form("FMD%d%c[%02d,%03d]",
 120                                                       dd, r, s, t),
 121                                                  UShort_t(i)));
 122   }
 123   d->Add(extraDead);
 124   fHasXtraDead = nBits > 0;
 125
 126 }
 127
 128 //____________________________________________________________________
 129 void
 130 AliFMDESDFixer::AddDead(UShort_t d, Char_t r, UShort_t s, UShort_t t)
 131 {
 132   if (d < 1 || d > 3) {
 133     Warning("AddDead", "Invalid detector FMD%d", d);
 134     return;
 135   }
 136   Bool_t inner = (r == 'I' || r == 'i');
 137   if (d == 1 && !inner) {
 138     Warning("AddDead", "Invalid ring FMD%d%c", d, r);
 139     return;
 140   }
 141   if ((inner && s >= 20) || (!inner && s >= 40)) {
 142     Warning("AddDead", "Invalid sector FMD%d%c[%02d]", d, r, s);
 143     return;
 144   }
 145   if ((inner && t >= 512) || (!inner && t >= 256)) {
 146     Warning("AddDead", "Invalid strip FMD%d%c[%02d,%03d]", d, r, s, t);
 147     return;
 148   }
 149
 150   Int_t id = AliFMDStripIndex::Pack(d, r, s, t);
 151   // Int_t i  = 0;
 152   fXtraDead.SetBitNumber(id, true);
 153 }
 154 //____________________________________________________________________
 155 void
 156 AliFMDESDFixer::AddDeadRegion(UShort_t d,  Char_t r,
 157                                    UShort_t s1, UShort_t s2,
 158                                    UShort_t t1, UShort_t t2)
 159 {
 160   // Add a dead region spanning from FMD<d><r>[<s1>,<t1>] to
 161   // FMD<d><r>[<s2>,<t2>] (both inclusive)
 162   for (Int_t s = s1; s <= s2; s++)
 163     for (Int_t t = t1; t <= t2; t++)
 164       AddDead(d, r, s, t);
 165 }
 166 //____________________________________________________________________
 167 void
 168 AliFMDESDFixer::AddDead(const Char_t* script)
 169 {
 170   if (!script || script[0] == '\0') return;
 171
 172   gROOT->Macro(Form("%s((AliFMDESDFixer*)%p);", script, this));
 173 }
 174
 175 //____________________________________________________________________
 176 Bool_t
 177 AliFMDESDFixer::IsDead(UShort_t d, Char_t r, UShort_t s, UShort_t t) const
 178 {
 179   Int_t id = AliFMDStripIndex::Pack(d, r, s, t);
 180   return fXtraDead.TestBitNumber(id);
 181 }
 182
 183 //____________________________________________________________________
 184 Int_t
 185 AliFMDESDFixer::FindTargetNoiseFactor(const AliESDFMD& esd, Bool_t check) const
 186 {
 187   if (!IsUseNoiseCorrection())
 188     // If the reconstruction factor was high (4 or more), do nothing
 189     return 0;
 190
 191   Int_t target = 0;
 192   if (AliESDFMD::Class_Version() < 4) {
 193     // IF we running with older STEER - we fix it here
 194     target = 4;
 195   } else {
 196 #if 1
 197     if (!esd.TestBit(1 << 14)) {
 198       // If the bit isn't set, do nothing
 199       return 0;
 200     }
 201 #else
 202     // Uncommented until Peter commits patch to STEER/ESD
 203     if (!esd.NeedNoiseFix()) {
 204       // If the bit isn't set, do nothing
 205       return 0;
 206     }
 207 #endif
 208     target = Int_t(esd.GetNoiseFactor());
 209   }
 210   // Get the target factor - even thought the method below returns a
 211   // floating point value, we know that the noise factor is always
 212   // integer, so we coerce it to be the same here.
 213   target -= fRecoFactor;
 214
 215   // If the target factor is the same or smaller than the assumed
 216   // factor, we have nothing to do here, and we return immediately
 217   if (target <= 0) return 0;
 218
 219   // Get the scaled noise from the correction mananger
 220   if (check && !AliForwardCorrectionManager::Instance().GetNoiseGain())
 221     return 0;
 222
 223   return target;
 224 }
 225
 226 //____________________________________________________________________
 227 #define ETA2COS(ETA)                                            \
 228   TMath::Cos(2*TMath::ATan(TMath::Exp(-TMath::Abs(ETA))))
 229
 230 //____________________________________________________________________
 231 void
 232 AliFMDESDFixer::Fix(AliESDFMD& esd, Double_t zvtx)
 233 {
 234
 235   const AliFMDCorrNoiseGain* ng  = 0;
 236   Int_t tgtFactor = FindTargetNoiseFactor(esd, false);
 237   if (tgtFactor > 0)
 238     ng  = AliForwardCorrectionManager::Instance().GetNoiseGain();
 239
 240   if (!ng && !fHasXtraDead && !fRecalculateEta && !fInvalidIsEmpty)
 241     // We have nothing to do!
 242     return;
 243
 244   UShort_t nDead = 0;
 245   for (UShort_t d = 1; d <= 3; d++) {
 246     UShort_t nQ = d == 1 ? 1 : 2;
 247
 248     for (UShort_t q = 0; q < nQ; q++) {
 249       Char_t   r  = (q == 0 ? 'I' : 'O');
 250       UShort_t nS = (q == 0 ?  20 :  40);
 251       UShort_t nT = (q == 0 ? 512 : 256);
 252
 253       for (UShort_t s = 0; s < nS; s++) {
 254         for (UShort_t t = 0; t < nT; t++) {
 255           Double_t mult     = esd.Multiplicity(d,r,s,t);
 256           Double_t eta      = esd.Eta(d,r,s,t);
 257           Double_t cosTheta = 0;
 258
 259           if (CheckDead(d,r,s,t,mult)) nDead++;
 260
 261           // Possibly re-calculate eta
 262           if (fRecalculateEta) RecalculateEta(d,r,s,t,zvtx,eta,mult,cosTheta);
 263
 264           // Possibly correct for poor treatment of ZS in reconstruction.
 265           if (ng && mult != AliESDFMD::kInvalidMult) {
 266             if (cosTheta <= 0) cosTheta = ETA2COS(eta);
 267             if (!NoiseCorrect(tgtFactor,ng->Get(d,r,s,t), cosTheta, mult))
 268               nDead++;
 269           }
 270
 271           // Write out final values to object
 272           if (mult >= AliESDFMD::kInvalidMult) mult = AliESDFMD::kInvalidMult;
 273           esd.SetMultiplicity(d,r,s,t,mult);
 274           esd.SetEta(d,r,s,t,eta);
 275         } // for t
 276       } // for s
 277     } // for q
 278   } // for d
 279   fDeadChange->Fill(nDead);
 280 }
 281
 282 //____________________________________________________________________
 283 Bool_t
 284 AliFMDESDFixer::CheckDead(UShort_t d, Char_t r, UShort_t s, UShort_t t,
 285                           Double_t& mult)
 286 {
 287   // Correct for zero's being flagged as invalid
 288   if (mult == AliESDFMD::kInvalidMult && fInvalidIsEmpty) mult = 0;
 289
 290   // Take into account what we're defined as dead
 291   if (IsDead(d,r,s,t)) {
 292     mult = AliESDFMD::kInvalidMult;
 293     return true;
 294   }
 295   return false;
 296 }
 297
 298 //____________________________________________________________________
 299 void
 300 AliFMDESDFixer::RecalculateEta(UShort_t d, Char_t r, UShort_t s, UShort_t t,
 301                                Double_t zvtx, Double_t& eta, Double_t& mult,
 302                                Double_t& cosTheta)
 303 {
 304   Double_t oldEta = eta;
 305   Double_t newEta = AliForwardUtil::GetEtaFromStrip(d,r,s,t, zvtx);
 306   eta             = newEta;
 307
 308   fEtaChange->Fill(newEta-oldEta);
 309
 310   if (mult == AliESDFMD::kInvalidMult) return;
 311
 312   Double_t newCos = ETA2COS(newEta);
 313   Double_t oldCos = ETA2COS(oldEta);
 314   Double_t corr   = newCos / oldCos;
 315   cosTheta        = newCos;
 316   mult            *= corr;
 317 }
 318 //____________________________________________________________________
 319 Bool_t
 320 AliFMDESDFixer::NoiseCorrect(Int_t target, Double_t corr, Double_t cosTheta,
 321                              Double_t& mult)
 322 {
 323   if (corr > fMaxNoiseCorr || corr <= 0) {
 324     mult = AliESDFMD::kInvalidMult;
 325     return false;
 326   }
 327   Double_t add = corr * target * cosTheta;
 328   fNoiseChange->Fill(add);
 329
 330   mult += add;
 331   return true;
 332 }
 333
 334 #define PF(N,V,...)                                     \
 335   AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
 336 #define PFB(N,FLAG)                             \
 337   do {                                                                  \
 338     AliForwardUtil::PrintName(N);                                       \
 339     std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
 340   } while(false)
 341 #define PFV(N,VALUE)                                    \
 342   do {                                                  \
 343     AliForwardUtil::PrintName(N);                       \
 344     std::cout << (VALUE) << std::endl; } while(false)
 345
 346 //____________________________________________________________________
 347 void
 348 AliFMDESDFixer::Print(Option_t*) const
 349 {
 350   AliForwardUtil::PrintTask(*this);
 351   gROOT->IncreaseDirLevel();
 352   PFB("Consider invalid null", fInvalidIsEmpty);
 353   PFB("Has extra dead", fHasXtraDead || fXtraDead.GetNbits() > 0);
 354   PFV("Reco noise factor", fRecoFactor);
 355   PFV("Max noise corr", fMaxNoiseCorr);
 356   PFB("Recalc. eta", fRecalculateEta);
 357   gROOT->DecreaseDirLevel();
 358 }
 359
 360 //____________________________________________________________________
 361 //
 362 // EOF
 363 //