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