PWG2/FORWARD/analysis2/AliFMDSharingFilter.cxx

   1 //
   2 // Class to do the sharing correction of FMD ESD data
   3 //
   4 #include "AliFMDSharingFilter.h"
   5 #include <AliESDFMD.h>
   6 #include <TAxis.h>
   7 #include <TList.h>
   8 #include <TH1.h>
   9 #include <TMath.h>
  10 #include "AliFMDAnaParameters.h"
  11 #include <AliLog.h>
  12
  13 ClassImp(AliFMDSharingFilter)
  14 #if 0
  15 ; // This is for Emacs
  16 #endif
  17
  18
  19 //____________________________________________________________________
  20 AliFMDSharingFilter::AliFMDSharingFilter()
  21   : TNamed(),
  22     fRingHistos(),
  23     fLowCut(0.3),
  24     fCorrectAngles(kFALSE),
  25     fEtaCorr(0)
  26 {}
  27
  28 //____________________________________________________________________
  29 AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
  30   : TNamed("fmdSharingFilter", title),
  31     fRingHistos(),
  32     fLowCut(0.3),
  33     fCorrectAngles(kFALSE),
  34     fEtaCorr(0)
  35 {
  36   fRingHistos.Add(new RingHistos(1, 'I'));
  37   fRingHistos.Add(new RingHistos(2, 'I'));
  38   fRingHistos.Add(new RingHistos(2, 'O'));
  39   fRingHistos.Add(new RingHistos(3, 'I'));
  40   fRingHistos.Add(new RingHistos(3, 'O'));
  41   fEtaCorr = new TH2F("etaCorr", "Correction of #eta",
  42                       200, -4, 6,  200, -4, 6);
  43   fEtaCorr->SetXTitle("#eta from ESD");
  44   fEtaCorr->SetYTitle("#eta from AnaParameters");
  45 }
  46
  47 //____________________________________________________________________
  48 AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
  49   : TNamed(o),
  50     fRingHistos(),
  51     fLowCut(o.fLowCut),
  52     fCorrectAngles(o.fCorrectAngles),
  53     fEtaCorr(o.fEtaCorr)
  54 {
  55   TIter    next(&o.fRingHistos);
  56   TObject* obj = 0;
  57   while ((obj = next())) fRingHistos.Add(obj);
  58 }
  59
  60 //____________________________________________________________________
  61 AliFMDSharingFilter::~AliFMDSharingFilter()
  62 {
  63   fRingHistos.Delete();
  64 }
  65
  66 //____________________________________________________________________
  67 AliFMDSharingFilter&
  68 AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
  69 {
  70   SetName(o.GetName());
  71   SetTitle(o.GetTitle());
  72
  73   fLowCut        = o.fLowCut;
  74   fCorrectAngles = o.fCorrectAngles;
  75
  76   fRingHistos.Delete();
  77   TIter    next(&o.fRingHistos);
  78   TObject* obj = 0;
  79   while ((obj = next())) fRingHistos.Add(obj);
  80
  81   return *this;
  82 }
  83
  84 //____________________________________________________________________
  85 AliFMDSharingFilter::RingHistos*
  86 AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
  87 {
  88   Int_t idx = -1;
  89   switch (d) {
  90   case 1: idx = 0; break;
  91   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
  92   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
  93   }
  94   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
  95
  96   return static_cast<RingHistos*>(fRingHistos.At(idx));
  97 }
  98
  99 //____________________________________________________________________
 100 Bool_t
 101 AliFMDSharingFilter::Filter(const AliESDFMD& input,
 102                             Bool_t           lowFlux,
 103                             AliESDFMD&       output,
 104                             Double_t         vz)
 105 {
 106   output.Clear();
 107   TIter    next(&fRingHistos);
 108   RingHistos* o = 0;
 109   while ((o = static_cast<RingHistos*>(next())))
 110     o->Clear();
 111
 112   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 113
 114   for(UShort_t d = 1; d <= 3; d++) {
 115     Int_t nRings = (d == 1 ? 1 : 2);
 116     for (UShort_t q = 0; q < nRings; q++) {
 117       Char_t      r    = (q == 0 ? 'I' : 'O');
 118       UShort_t    nsec = (q == 0 ?  20 :  40);
 119       UShort_t    nstr = (q == 0 ? 512 : 256);
 120
 121       RingHistos* histos = GetRingHistos(d, r);
 122
 123       for(UShort_t s = 0; s < nsec;  s++) {
 124         Bool_t usedThis = kFALSE;
 125         Bool_t usedPrev = kFALSE;
 126
 127         for(UShort_t t = 0; t < nstr; t++) {
 128           output.SetMultiplicity(d,r,s,t,0.);
 129           Float_t mult = SignalInStrip(input,d,r,s,t);
 130
 131           // Keep dead-channel information.
 132           if(mult == AliESDFMD::kInvalidMult)
 133             output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
 134
 135           // If no signal or dead strip, go on.
 136           if (mult == AliESDFMD::kInvalidMult || mult == 0) continue;
 137
 138           // Get the pseudo-rapidity
 139           Double_t eta1 = input.Eta(d,r,s,t);
 140           Double_t eta2 = pars->GetEtaFromStrip(d,r,s,t,vz);
 141           fEtaCorr->Fill(eta1, eta2);
 142           Double_t eta = eta2;
 143
 144           // Fill the diagnostics histogram
 145           histos->fBefore->Fill(mult);
 146
 147           // Get next and previous signal - if any
 148           Double_t prevE = 0;
 149           Double_t nextE = 0;
 150           if (t != 0) {
 151             prevE = SignalInStrip(input,d,r,s,t-1);
 152             if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
 153           }
 154           if (t != nstr - 1) {
 155             nextE = SignalInStrip(input,d,r,s,t+1);
 156             if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
 157           }
 158
 159           Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
 160                                                       lowFlux,d,r,s,t,
 161                                                       usedPrev,usedThis);
 162           if (mergedEnergy > fLowCut) histos->Incr();
 163           histos->fAfter->Fill(mergedEnergy);
 164
 165           output.SetMultiplicity(d,r,s,t,mergedEnergy);
 166           output.SetEta(d,r,s,t,eta);
 167         } // for strip
 168       } // for sector
 169     } // for ring
 170   } // for detector
 171
 172   next.Reset();
 173   while ((o = static_cast<RingHistos*>(next())))
 174     o->Finish();
 175
 176   return kTRUE;
 177 }
 178
 179 //_____________________________________________________________________
 180 Double_t
 181 AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
 182                                    UShort_t         d,
 183                                    Char_t           r,
 184                                    UShort_t         s,
 185                                    UShort_t         t) const
 186 {
 187   Double_t mult = input.Multiplicity(d,r,s,t);
 188   if (mult == AliESDFMD::kInvalidMult || mult == 0) return mult;
 189
 190   if (fCorrectAngles && !input.IsAngleCorrected())
 191     mult = AngleCorrect(mult, input.Eta(d,r,s,t));
 192   else if (!fCorrectAngles && input.IsAngleCorrected())
 193     mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
 194   return mult;
 195 }
 196
 197 //_____________________________________________________________________
 198 Double_t
 199 AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
 200 {
 201   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 202
 203   // Get the high cut.  The high cut is defined as the
 204   // most-probably-value peak found from the energy distributions, minus
 205   // 2 times the width of the corresponding Landau.
 206   Double_t mpv     = pars->GetMPV(d,r,eta);
 207   Double_t w       = pars->GetSigma(d,r,eta);
 208   if (mpv > 100)
 209     AliWarning(Form("FMD%d%c, eta=%f, MPV=%f, W=%f", d, r, eta, mpv, w));
 210   Double_t highCut = (mpv - 2 * w);
 211   return highCut;
 212 }
 213
 214 //_____________________________________________________________________
 215 Double_t
 216 AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
 217                                          Double_t eta,
 218                                          Double_t prevE,
 219                                          Double_t nextE,
 220                                          Bool_t   lowFlux,
 221                                          UShort_t d,
 222                                          Char_t   r,
 223                                          UShort_t /*s*/,
 224                                          UShort_t /*t*/,
 225                                          Bool_t&  usedPrev,
 226                                          Bool_t&  usedThis)
 227 {
 228   // IF the multiplicity is very large, or below the cut, reject it, and
 229   // flag it as candidate for sharing
 230   if(mult > 12 || mult < fLowCut) {
 231     usedThis      = kFALSE;
 232     usedPrev      = kFALSE;
 233     return 0;
 234   }
 235
 236   // If this was shared with the previous signal, return 0 and mark it as
 237   // not a candiate for sharing
 238   if(usedThis) {
 239     usedThis      = kFALSE;
 240     usedPrev      = kTRUE;
 241     return 0.;
 242   }
 243
 244   //analyse and perform sharing on one strip
 245
 246   // Calculate the total energy loss
 247   Double_t highCut = GetHighCut(d, r, eta);
 248
 249   // If this signal is smaller than the next, and the next signal is smaller
 250   // than then the high cut, and it's a low-flux event, then mark this and
 251   // the next signal as candidates for sharing
 252   //
 253   // This is the test if the signal is the smaller of two possibly
 254   // shared signals.
 255   if (mult  < nextE   &&
 256       nextE > highCut &&
 257       lowFlux) {
 258     usedThis      = kFALSE;
 259     usedPrev      = kFALSE;
 260     return 0;
 261   }
 262
 263   Double_t totalE  = mult;
 264
 265   // If the previous signal was larger than the low cut, and
 266   // the previous was smaller than high cut, and the previous signal
 267   // isn't marked as used, then add it's energy to this signal
 268   if (prevE > fLowCut &&
 269       prevE < highCut &&
 270       !usedPrev)
 271     totalE += prevE;
 272
 273   // If the next signal is larger than the low cut, and
 274   // the next signal is smaller than the low cut, then add the next signal
 275   // to this, and marked the next signal as used
 276   if(nextE > fLowCut && nextE < highCut ) {
 277     totalE      += nextE;
 278     usedThis =  kTRUE;
 279   }
 280
 281   // If we're processing on non-angle corrected data, we should do the
 282   // angle correction here
 283   if (!fCorrectAngles)
 284     totalE = AngleCorrect(totalE, eta);
 285
 286   // Fill post processing histogram
 287   // if(totalE > fLowCut)
 288   //   GetRingHistos(d,r)->fAfter->Fill(totalE);
 289
 290   // Return value
 291   Double_t mergedEnergy = 0;
 292
 293   if(totalE > 0) {
 294     // If we have a signal, then this is used (sharedPrev=true) and
 295     // the signal is set to the result
 296     mergedEnergy = totalE;
 297     usedPrev     = kTRUE;
 298   }
 299   else {
 300     // If we do not have a signal (too low), then this is not shared,
 301     // and the next strip is not shared either
 302     usedThis  = kFALSE;
 303     usedPrev  = kFALSE;
 304   }
 305
 306   return mergedEnergy;
 307 }
 308 //____________________________________________________________________
 309 Double_t
 310 AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
 311 {
 312   Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
 313   if (eta < 0) theta -= TMath::Pi();
 314   return mult * TMath::Cos(theta);
 315 }
 316 //____________________________________________________________________
 317 Double_t
 318 AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
 319 {
 320   Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
 321   if (eta < 0) theta -= TMath::Pi();
 322   return mult / TMath::Cos(theta);
 323 }
 324
 325 //____________________________________________________________________
 326 void
 327 AliFMDSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
 328 {
 329   if (nEvents <= 0) return;
 330   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 331   if (!d) return;
 332
 333   TIter    next(&fRingHistos);
 334   RingHistos* o = 0;
 335   while ((o = static_cast<RingHistos*>(next())))
 336     o->ScaleHistograms(d, nEvents);
 337 }
 338
 339 //____________________________________________________________________
 340 void
 341 AliFMDSharingFilter::DefineOutput(TList* dir)
 342 {
 343   TList* d = new TList;
 344   d->SetName(GetName());
 345   dir->Add(d);
 346   d->Add(fEtaCorr);
 347   TIter    next(&fRingHistos);
 348   RingHistos* o = 0;
 349   while ((o = static_cast<RingHistos*>(next()))) {
 350     o->Output(d);
 351   }
 352 }
 353
 354 //====================================================================
 355 AliFMDSharingFilter::RingHistos::RingHistos()
 356   : AliForwardUtil::RingHistos(),
 357     fBefore(0),
 358     fAfter(0),
 359     fHits(0),
 360     fNHits(0)
 361 {}
 362
 363 //____________________________________________________________________
 364 AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
 365   : AliForwardUtil::RingHistos(d,r),
 366     fBefore(0),
 367     fAfter(0),
 368     fHits(0),
 369     fNHits(0)
 370 {
 371   fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()),
 372                      Form("Energy loss in %s (reconstruction)", fName.Data()),
 373                      1000, 0, 25);
 374   fAfter  = new TH1D(Form("%s_Ana_Eloss", fName.Data()),
 375                      Form("Energy loss in %s (sharing corrected)",
 376                           fName.Data()), 1000, 0, 25);
 377   fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
 378   fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
 379   fBefore->SetFillColor(kRed+1);
 380   fBefore->SetFillStyle(3001);
 381   // fBefore->Sumw2();
 382   fBefore->SetDirectory(0);
 383   fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
 384   fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
 385   fAfter->SetFillColor(kBlue+1);
 386   fAfter->SetFillStyle(3001);
 387   // fAfter->Sumw2();
 388   fAfter->SetDirectory(0);
 389
 390   fHits = new TH1D(Form("%s_Hits", fName.Data()),
 391                    Form("Number of hits in %s", fName.Data()),
 392                    200, 0, 200000);
 393   fHits->SetDirectory(0);
 394   fHits->SetXTitle("# of hits");
 395   fHits->SetFillColor(kGreen+1);
 396 }
 397 //____________________________________________________________________
 398 AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
 399   : AliForwardUtil::RingHistos(o),
 400     fBefore(o.fBefore),
 401     fAfter(o.fAfter),
 402     fHits(o.fHits),
 403     fNHits(o.fNHits)
 404 {}
 405
 406 //____________________________________________________________________
 407 AliFMDSharingFilter::RingHistos&
 408 AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
 409 {
 410   AliForwardUtil::RingHistos::operator=(o);
 411   fDet = o.fDet;
 412   fRing = o.fRing;
 413
 414   if (fBefore) delete  fBefore;
 415   if (fAfter)  delete  fAfter;
 416   if (fHits)   delete fHits;
 417
 418   fBefore = static_cast<TH1D*>(o.fBefore->Clone());
 419   fAfter  = static_cast<TH1D*>(o.fAfter->Clone());
 420   fHits  = static_cast<TH1D*>(o.fHits->Clone());
 421
 422   return *this;
 423 }
 424 //____________________________________________________________________
 425 AliFMDSharingFilter::RingHistos::~RingHistos()
 426 {
 427   if (fBefore) delete fBefore;
 428   if (fAfter)  delete fAfter;
 429   if (fHits)   delete fHits;
 430 }
 431 //____________________________________________________________________
 432 void
 433 AliFMDSharingFilter::RingHistos::Finish()
 434 {
 435   fHits->Fill(fNHits);
 436 }
 437
 438 //____________________________________________________________________
 439 void
 440 AliFMDSharingFilter::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
 441 {
 442   TList* l = GetOutputList(dir);
 443   if (!l) return;
 444
 445   TH1D* before = static_cast<TH1D*>(l->FindObject(Form("%s_ESD_ELoss",
 446                                                        fName.Data())));
 447   TH1D* after  = static_cast<TH1D*>(l->FindObject(Form("%s_Ana_ELoss",
 448                                                        fName.Data())));
 449   if (before) before->Scale(1./nEvents);
 450   if (after)  after->Scale(1./nEvents);
 451 }
 452
 453 //____________________________________________________________________
 454 void
 455 AliFMDSharingFilter::RingHistos::Output(TList* dir)
 456 {
 457   TList* d = DefineOutputList(dir);
 458   d->Add(fBefore);
 459   d->Add(fAfter);
 460   d->Add(fHits);
 461   dir->Add(d);
 462 }
 463
 464 //____________________________________________________________________
 465 //
 466 // EOF
 467 //