PWG2/FORWARD/analysis2/AliFMDSharingFilter.cxx

   1 //
   2 // Class to do the sharing correction.  That is, a filter that merges
   3 // adjacent strip signals presumably originating from a single particle
   4 // that impinges on the detector in such a way that it deposite energy
   5 // into two or more strips.
   6 //
   7 // Input:
   8 //    - AliESDFMD object  - from reconstruction
   9 //
  10 // Output:
  11 //    - AliESDFMD object  - copy of input, but with signals merged
  12 //
  13 // Corrections used:
  14 //    - AliFMDCorrELossFit
  15 //
  16 // Histograms:
  17 //    - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of
  18 //      signals before and after the filter.
  19 //    - For each ring (see above), an array of distributions of number of
  20 //      hit strips for each vertex bin (if enabled - see Init method)
  21 //
  22 //
  23 //
  24 #include "AliFMDSharingFilter.h"
  25 #include <AliESDFMD.h>
  26 #include <TAxis.h>
  27 #include <TList.h>
  28 #include <TH1.h>
  29 #include <TMath.h>
  30 #include "AliForwardCorrectionManager.h"
  31 #include <AliLog.h>
  32 #include <TROOT.h>
  33 #include <iostream>
  34 #include <iomanip>
  35
  36 ClassImp(AliFMDSharingFilter)
  37 #if 0
  38 ; // This is for Emacs
  39 #endif
  40
  41
  42 //____________________________________________________________________
  43 AliFMDSharingFilter::AliFMDSharingFilter()
  44   : TNamed(),
  45     fRingHistos(),
  46     fLowCut(0.),
  47     fCorrectAngles(kFALSE),
  48     fNXi(1),
  49     fDebug(0)
  50 {
  51   //
  52   // Default Constructor - do not use
  53   //
  54 }
  55
  56 //____________________________________________________________________
  57 AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
  58   : TNamed("fmdSharingFilter", title),
  59     fRingHistos(),
  60     fLowCut(0.),
  61     fCorrectAngles(kFALSE),
  62     fNXi(1),
  63     fDebug(0)
  64 {
  65   //
  66   // Constructor
  67   //
  68   // Parameters:
  69   //    title Title of object  - not significant
  70   //
  71   fRingHistos.Add(new RingHistos(1, 'I'));
  72   fRingHistos.Add(new RingHistos(2, 'I'));
  73   fRingHistos.Add(new RingHistos(2, 'O'));
  74   fRingHistos.Add(new RingHistos(3, 'I'));
  75   fRingHistos.Add(new RingHistos(3, 'O'));
  76 }
  77
  78 //____________________________________________________________________
  79 AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
  80   : TNamed(o),
  81     fRingHistos(),
  82     fLowCut(o.fLowCut),
  83     fCorrectAngles(o.fCorrectAngles),
  84     fNXi(o.fNXi),
  85     fDebug(o.fDebug)
  86 {
  87   //
  88   // Copy constructor
  89   //
  90   // Parameters:
  91   //    o Object to copy from
  92   //
  93   TIter    next(&o.fRingHistos);
  94   TObject* obj = 0;
  95   while ((obj = next())) fRingHistos.Add(obj);
  96 }
  97
  98 //____________________________________________________________________
  99 AliFMDSharingFilter::~AliFMDSharingFilter()
 100 {
 101   //
 102   // Destructor
 103   //
 104   fRingHistos.Delete();
 105 }
 106
 107 //____________________________________________________________________
 108 AliFMDSharingFilter&
 109 AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
 110 {
 111   //
 112   // Assignment operator
 113   //
 114   // Parameters:
 115   //    o Object to assign from
 116   //
 117   // Return:
 118   //    Reference to this
 119   //
 120   TNamed::operator=(o);
 121
 122   fLowCut        = o.fLowCut;
 123   fCorrectAngles = o.fCorrectAngles;
 124   fNXi           = o.fNXi;
 125   fDebug         = o.fDebug;
 126
 127   fRingHistos.Delete();
 128   TIter    next(&o.fRingHistos);
 129   TObject* obj = 0;
 130   while ((obj = next())) fRingHistos.Add(obj);
 131
 132   return *this;
 133 }
 134
 135 //____________________________________________________________________
 136 AliFMDSharingFilter::RingHistos*
 137 AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
 138 {
 139   //
 140   // Get the ring histogram container
 141   //
 142   // Parameters:
 143   //    d Detector
 144   //    r Ring
 145   //
 146   // Return:
 147   //    Ring histogram container
 148   //
 149   Int_t idx = -1;
 150   switch (d) {
 151   case 1: idx = 0; break;
 152   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 153   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 154   }
 155   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
 156
 157   return static_cast<RingHistos*>(fRingHistos.At(idx));
 158 }
 159
 160 //____________________________________________________________________
 161 Bool_t
 162 AliFMDSharingFilter::Filter(const AliESDFMD& input,
 163                             Bool_t           lowFlux,
 164                             AliESDFMD&       output)
 165 {
 166   //
 167   // Filter the input AliESDFMD object
 168   //
 169   // Parameters:
 170   //    input     Input
 171   //    lowFlux   If this is a low-flux event
 172   //    output    Output AliESDFMD object
 173   //
 174   // Return:
 175   //    True on success, false otherwise
 176   //
 177   output.Clear();
 178   TIter    next(&fRingHistos);
 179   RingHistos* o      = 0;
 180   Double_t    lowCut = GetLowCut();
 181   while ((o = static_cast<RingHistos*>(next())))
 182     o->Clear();
 183
 184   // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 185
 186   for(UShort_t d = 1; d <= 3; d++) {
 187     Int_t nRings = (d == 1 ? 1 : 2);
 188     for (UShort_t q = 0; q < nRings; q++) {
 189       Char_t      r    = (q == 0 ? 'I' : 'O');
 190       UShort_t    nsec = (q == 0 ?  20 :  40);
 191       UShort_t    nstr = (q == 0 ? 512 : 256);
 192
 193       RingHistos* histos = GetRingHistos(d, r);
 194
 195       for(UShort_t s = 0; s < nsec;  s++) {
 196         Bool_t usedThis = kFALSE;
 197         Bool_t usedPrev = kFALSE;
 198
 199         for(UShort_t t = 0; t < nstr; t++) {
 200           output.SetMultiplicity(d,r,s,t,0.);
 201           Float_t mult = SignalInStrip(input,d,r,s,t);
 202
 203           // Keep dead-channel information.
 204           if(mult == AliESDFMD::kInvalidMult)
 205             output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
 206
 207           // If no signal or dead strip, go on.
 208           if (mult == AliESDFMD::kInvalidMult || mult == 0) continue;
 209
 210           // Get the pseudo-rapidity
 211           Double_t eta = input.Eta(d,r,s,t);
 212
 213           // Fill the diagnostics histogram
 214           histos->fBefore->Fill(mult);
 215
 216           // Get next and previous signal - if any
 217           Double_t prevE = 0;
 218           Double_t nextE = 0;
 219           if (t != 0) {
 220             prevE = SignalInStrip(input,d,r,s,t-1);
 221             if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
 222           }
 223           if (t != nstr - 1) {
 224             nextE = SignalInStrip(input,d,r,s,t+1);
 225             if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
 226           }
 227
 228           Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
 229                                                       lowFlux,d,r,s,t,
 230                                                       usedPrev,usedThis);
 231           if (mergedEnergy > lowCut) histos->Incr();
 232           histos->fAfter->Fill(mergedEnergy);
 233
 234           output.SetMultiplicity(d,r,s,t,mergedEnergy);
 235           output.SetEta(d,r,s,t,eta);
 236         } // for strip
 237       } // for sector
 238     } // for ring
 239   } // for detector
 240
 241   next.Reset();
 242   while ((o = static_cast<RingHistos*>(next())))
 243     o->Finish();
 244
 245   return kTRUE;
 246 }
 247
 248 //_____________________________________________________________________
 249 Double_t
 250 AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
 251                                    UShort_t         d,
 252                                    Char_t           r,
 253                                    UShort_t         s,
 254                                    UShort_t         t) const
 255 {
 256   //
 257   // Get the signal in a strip
 258   //
 259   // Parameters:
 260   //    fmd   ESD object
 261   //    d     Detector
 262   //    r     Ring
 263   //    s     Sector
 264   //    t     Strip
 265   //
 266   // Return:
 267   //    The energy signal
 268   //
 269   Double_t mult = input.Multiplicity(d,r,s,t);
 270   if (mult == AliESDFMD::kInvalidMult || mult == 0) return mult;
 271
 272   if (fCorrectAngles && !input.IsAngleCorrected())
 273     mult = AngleCorrect(mult, input.Eta(d,r,s,t));
 274   else if (!fCorrectAngles && input.IsAngleCorrected())
 275     mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
 276   return mult;
 277 }
 278 //_____________________________________________________________________
 279 Double_t
 280 AliFMDSharingFilter::GetLowCut() const
 281 {
 282   //
 283   // Get the low cut.  Normally, the low cut is taken to be the lower
 284   // value of the fit range used when generating the energy loss fits.
 285   // However, if fLowCut is set (using SetLowCit) to a value greater
 286   // than 0, then that value is used.
 287   //
 288   if (fLowCut > 0) return fLowCut;
 289   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 290   AliFMDCorrELossFit* fits = fcm.GetELossFit();
 291   return fits->GetLowCut();
 292 }
 293
 294 //_____________________________________________________________________
 295 Double_t
 296 AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
 297 {
 298   //
 299   // Get the high cut.  The high cut is defined as the
 300   // most-probably-value peak found from the energy distributions, minus
 301   // 2 times the width of the corresponding Landau.
 302   //
 303   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 304
 305
 306   // Get the high cut.  The high cut is defined as the
 307   // most-probably-value peak found from the energy distributions, minus
 308   // 2 times the width of the corresponding Landau.
 309   AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
 310   Double_t delta = 1024;
 311   Double_t xi    = 1024;
 312   if (!fit) {
 313     AliError(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
 314   }
 315   else {
 316     delta = fit->fDelta;
 317     xi    = fit->fXi;
 318   }
 319
 320   if (delta > 100) {
 321     AliWarning(Form("FMD%d%c, eta=%f, Delta=%f, xi=%f", d, r, eta, delta, xi));
 322   }
 323   Double_t highCut = (delta - fNXi * xi);
 324   return highCut;
 325 }
 326
 327 //_____________________________________________________________________
 328 Double_t
 329 AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
 330                                          Double_t eta,
 331                                          Double_t prevE,
 332                                          Double_t nextE,
 333                                          Bool_t   lowFlux,
 334                                          UShort_t d,
 335                                          Char_t   r,
 336                                          UShort_t /*s*/,
 337                                          UShort_t /*t*/,
 338                                          Bool_t&  usedPrev,
 339                                          Bool_t&  usedThis)
 340 {
 341   //
 342   // The actual algorithm
 343   //
 344   // Parameters:
 345   //    mult      The unfiltered signal in the strip
 346   //    eta       Psuedo rapidity
 347   //    prevE     Previous strip signal (or 0)
 348   //    nextE     Next strip signal (or 0)
 349   //    lowFlux   Whether this is a low flux event
 350   //    d         Detector
 351   //    r         Ring
 352   //    s         Sector
 353   //    t         Strip
 354   //    usedPrev  Whether the previous strip was used in sharing or not
 355   //    usedThis  Wether this strip was used in sharing or not.
 356   //
 357   // Return:
 358   //    The filtered signal in the strip
 359   //
 360
 361   // IF the multiplicity is very large, or below the cut, reject it, and
 362   // flag it as candidate for sharing
 363   Double_t    lowCut = GetLowCut();
 364   if(mult > 12 || mult < lowCut) {
 365     usedThis      = kFALSE;
 366     usedPrev      = kFALSE;
 367     return 0;
 368   }
 369
 370   // If this was shared with the previous signal, return 0 and mark it as
 371   // not a candiate for sharing
 372   if(usedThis) {
 373     usedThis      = kFALSE;
 374     usedPrev      = kTRUE;
 375     return 0.;
 376   }
 377
 378   //analyse and perform sharing on one strip
 379
 380   // Get the high cut
 381   Double_t highCut = GetHighCut(d, r, eta);
 382
 383   // If this signal is smaller than the next, and the next signal is smaller
 384   // than then the high cut, and it's a low-flux event, then mark this and
 385   // the next signal as candidates for sharing
 386   //
 387   // This is the test if the signal is the smaller of two possibly
 388   // shared signals.
 389   if (mult  < nextE   &&
 390       nextE > highCut &&
 391       lowFlux) {
 392     usedThis      = kFALSE;
 393     usedPrev      = kFALSE;
 394     return 0;
 395   }
 396
 397   Double_t totalE  = mult;
 398
 399   // If the previous signal was larger than the low cut, and
 400   // the previous was smaller than high cut, and the previous signal
 401   // isn't marked as used, then add it's energy to this signal
 402   if (prevE > lowCut &&
 403       prevE < highCut &&
 404       !usedPrev)
 405     totalE += prevE;
 406
 407   // If the next signal is larger than the low cut, and
 408   // the next signal is smaller than the low cut, then add the next signal
 409   // to this, and marked the next signal as used
 410   if(nextE > lowCut && nextE < highCut ) {
 411     totalE      += nextE;
 412     usedThis =  kTRUE;
 413   }
 414
 415   // If we're processing on non-angle corrected data, we should do the
 416   // angle correction here
 417   if (!fCorrectAngles)
 418     totalE = AngleCorrect(totalE, eta);
 419
 420   // Fill post processing histogram
 421   // if(totalE > fLowCut)
 422   //   GetRingHistos(d,r)->fAfter->Fill(totalE);
 423
 424   // Return value
 425   Double_t mergedEnergy = 0;
 426
 427   if(totalE > 0) {
 428     // If we have a signal, then this is used (sharedPrev=true) and
 429     // the signal is set to the result
 430     mergedEnergy = totalE;
 431     usedPrev     = kTRUE;
 432   }
 433   else {
 434     // If we do not have a signal (too low), then this is not shared,
 435     // and the next strip is not shared either
 436     usedThis  = kFALSE;
 437     usedPrev  = kFALSE;
 438   }
 439
 440   return mergedEnergy;
 441 }
 442 //____________________________________________________________________
 443 Double_t
 444 AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
 445 {
 446   //
 447   // Angle correct the signal
 448   //
 449   // Parameters:
 450   //    mult Angle Un-corrected Signal
 451   //    eta  Pseudo-rapidity
 452   //
 453   // Return:
 454   //    Angle corrected signal
 455   //
 456   Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
 457   if (eta < 0) theta -= TMath::Pi();
 458   return mult * TMath::Cos(theta);
 459 }
 460 //____________________________________________________________________
 461 Double_t
 462 AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
 463 {
 464   //
 465   // Angle de-correct the signal
 466   //
 467   // Parameters:
 468   //    mult Angle corrected Signal
 469   //    eta  Pseudo-rapidity
 470   //
 471   // Return:
 472   //    Angle un-corrected signal
 473   //
 474   Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
 475   if (eta < 0) theta -= TMath::Pi();
 476   return mult / TMath::Cos(theta);
 477 }
 478
 479 //____________________________________________________________________
 480 void
 481 AliFMDSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
 482 {
 483   //
 484   // Scale the histograms to the total number of events
 485   //
 486   // Parameters:
 487   //    dir     Where the output is
 488   //    nEvents Number of events
 489   //
 490   if (nEvents <= 0) return;
 491   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 492   if (!d) return;
 493
 494   TIter    next(&fRingHistos);
 495   RingHistos* o = 0;
 496   while ((o = static_cast<RingHistos*>(next())))
 497     o->ScaleHistograms(d, nEvents);
 498 }
 499
 500 //____________________________________________________________________
 501 void
 502 AliFMDSharingFilter::DefineOutput(TList* dir)
 503 {
 504   //
 505   // Define the output histograms.  These are put in a sub list of the
 506   // passed list.   The histograms are merged before the parent task calls
 507   // AliAnalysisTaskSE::Terminate
 508   //
 509   // Parameters:
 510   //    dir Directory to add to
 511   //
 512   TList* d = new TList;
 513   d->SetName(GetName());
 514   dir->Add(d);
 515   TIter    next(&fRingHistos);
 516   RingHistos* o = 0;
 517   while ((o = static_cast<RingHistos*>(next()))) {
 518     o->Output(d);
 519   }
 520 }
 521 //____________________________________________________________________
 522 void
 523 AliFMDSharingFilter::Print(Option_t* /*option*/) const
 524 {
 525   //
 526   // Print information
 527   //
 528   // Parameters:
 529   //    option Not used
 530   //
 531   char ind[gROOT->GetDirLevel()+1];
 532   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 533   ind[gROOT->GetDirLevel()] = '\0';
 534   std::cout << ind << "AliFMDSharingFilter: " << GetName() << '\n'
 535             << ind << " Low cut:                " << fLowCut << '\n'
 536             << ind << " N xi factor:            " << fNXi    << '\n'
 537             << ind << " Use corrected angles:   "
 538             << (fCorrectAngles ? "yes" : "no") << std::endl;
 539 }
 540
 541 //====================================================================
 542 AliFMDSharingFilter::RingHistos::RingHistos()
 543   : AliForwardUtil::RingHistos(),
 544     fBefore(0),
 545     fAfter(0),
 546     fHits(0),
 547     fNHits(0)
 548 {
 549   //
 550   // Default CTOR
 551   //
 552
 553 }
 554
 555 //____________________________________________________________________
 556 AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
 557   : AliForwardUtil::RingHistos(d,r),
 558     fBefore(0),
 559     fAfter(0),
 560     fHits(0),
 561     fNHits(0)
 562 {
 563   //
 564   // Constructor
 565   //
 566   // Parameters:
 567   //    d detector
 568   //    r ring
 569   //
 570   fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()),
 571                      Form("Energy loss in %s (reconstruction)", fName.Data()),
 572                      1000, 0, 25);
 573   fAfter  = new TH1D(Form("%s_Ana_Eloss", fName.Data()),
 574                      Form("Energy loss in %s (sharing corrected)",
 575                           fName.Data()), 1000, 0, 25);
 576   fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
 577   fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
 578   fBefore->SetFillColor(kRed+1);
 579   fBefore->SetFillStyle(3001);
 580   // fBefore->Sumw2();
 581   fBefore->SetDirectory(0);
 582   fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
 583   fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
 584   fAfter->SetFillColor(kBlue+1);
 585   fAfter->SetFillStyle(3001);
 586   // fAfter->Sumw2();
 587   fAfter->SetDirectory(0);
 588
 589   fHits = new TH1D(Form("%s_Hits", fName.Data()),
 590                    Form("Number of hits in %s", fName.Data()),
 591                    200, 0, 200000);
 592   fHits->SetDirectory(0);
 593   fHits->SetXTitle("# of hits");
 594   fHits->SetFillColor(kGreen+1);
 595 }
 596 //____________________________________________________________________
 597 AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
 598   : AliForwardUtil::RingHistos(o),
 599     fBefore(o.fBefore),
 600     fAfter(o.fAfter),
 601     fHits(o.fHits),
 602     fNHits(o.fNHits)
 603 {
 604   //
 605   // Copy constructor
 606   //
 607   // Parameters:
 608   //    o Object to copy from
 609   //
 610 }
 611
 612 //____________________________________________________________________
 613 AliFMDSharingFilter::RingHistos&
 614 AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
 615 {
 616   //
 617   // Assignment operator
 618   //
 619   // Parameters:
 620   //    o Object to assign from
 621   //
 622   // Return:
 623   //    Reference to this
 624   //
 625   AliForwardUtil::RingHistos::operator=(o);
 626   fDet = o.fDet;
 627   fRing = o.fRing;
 628
 629   if (fBefore) delete  fBefore;
 630   if (fAfter)  delete  fAfter;
 631   if (fHits)   delete fHits;
 632
 633   fBefore = static_cast<TH1D*>(o.fBefore->Clone());
 634   fAfter  = static_cast<TH1D*>(o.fAfter->Clone());
 635   fHits  = static_cast<TH1D*>(o.fHits->Clone());
 636
 637   return *this;
 638 }
 639 //____________________________________________________________________
 640 AliFMDSharingFilter::RingHistos::~RingHistos()
 641 {
 642   //
 643   // Destructor
 644   //
 645   if (fBefore) delete fBefore;
 646   if (fAfter)  delete fAfter;
 647   if (fHits)   delete fHits;
 648 }
 649 //____________________________________________________________________
 650 void
 651 AliFMDSharingFilter::RingHistos::Finish()
 652 {
 653   //
 654   // Finish off
 655   //
 656   //
 657   fHits->Fill(fNHits);
 658 }
 659
 660 //____________________________________________________________________
 661 void
 662 AliFMDSharingFilter::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
 663 {
 664   //
 665   // Scale the histograms to the total number of events
 666   //
 667   // Parameters:
 668   //    nEvents Number of events
 669   //    dir     Where the output is
 670   //
 671   TList* l = GetOutputList(dir);
 672   if (!l) return;
 673
 674   TH1D* before = static_cast<TH1D*>(l->FindObject(Form("%s_ESD_ELoss",
 675                                                        fName.Data())));
 676   TH1D* after  = static_cast<TH1D*>(l->FindObject(Form("%s_Ana_ELoss",
 677                                                        fName.Data())));
 678   if (before) before->Scale(1./nEvents);
 679   if (after)  after->Scale(1./nEvents);
 680 }
 681
 682 //____________________________________________________________________
 683 void
 684 AliFMDSharingFilter::RingHistos::Output(TList* dir)
 685 {
 686   //
 687   // Make output
 688   //
 689   // Parameters:
 690   //    dir where to store
 691   //
 692   TList* d = DefineOutputList(dir);
 693   d->Add(fBefore);
 694   d->Add(fAfter);
 695   d->Add(fHits);
 696   dir->Add(d);
 697 }
 698
 699 //____________________________________________________________________
 700 //
 701 // EOF
 702 //