PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.cxx

   1 //
   2 // Class to do the sharing correction for MC data.
   3 //
   4 // Input:
   5 //    - AliESDFMD object  - from reconstruction
   6 //    - Kinematics
   7 //    - Track-References
   8 //
   9 // Output:
  10 //    - AliESDFMD object  - copy of input, but with signals merged
  11 //
  12 // Corrections used:
  13 //    - None
  14 //
  15 // Histograms:
  16 //    - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of
  17 //      signals before and after the filter.
  18 //    - For each ring (see above), an array of distributions of number of
  19 //      hit strips for each vertex bin (if enabled - see Init method)
  20 //
  21 #include "AliFMDMCSharingFilter.h"
  22 #include <AliESDFMD.h>
  23 #include <AliMCEvent.h>
  24 #include <AliTrackReference.h>
  25 #include <AliStack.h>
  26 #include <TAxis.h>
  27 #include <TList.h>
  28 #include <TH1.h>
  29 #include <TMath.h>
  30 #include "AliFMDStripIndex.h"
  31 #include <AliLog.h>
  32 #include <TROOT.h>
  33 #include <iostream>
  34 #include <iomanip>
  35
  36 ClassImp(AliFMDMCSharingFilter)
  37 #if 0
  38 ; // This is for Emacs
  39 #endif
  40
  41
  42 //____________________________________________________________________
  43 AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title)
  44   : AliFMDSharingFilter(title),
  45     fFMD1i(0),
  46     fFMD2i(0),
  47     fFMD2o(0),
  48     fFMD3i(0),
  49     fFMD3o(0),
  50     fSumEta(0)
  51 {
  52   //
  53   // Constructor
  54   //
  55   // Parameters:
  56   //    title Title of object  - not significant
  57   //
  58   fFMD1i = new TH2D("FMD1i_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
  59   fFMD2i = new TH2D("FMD2i_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
  60   fFMD2o = new TH2D("FMD2o_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
  61   fFMD3i = new TH2D("FMD3i_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
  62   fFMD3o = new TH2D("FMD3o_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
  63   fFMD1i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
  64   fFMD1i->SetXTitle("Hits (MC)");
  65   fFMD2i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
  66   fFMD2i->SetXTitle("Hits (MC)");
  67   fFMD2o->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
  68   fFMD2o->SetXTitle("Hits (MC)");
  69   fFMD3i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
  70   fFMD3i->SetXTitle("Hits (MC)");
  71   fFMD3o->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
  72   fFMD3o->SetXTitle("Hits (MC)");
  73   fFMD1i->SetDirectory(0);
  74   fFMD2i->SetDirectory(0);
  75   fFMD2o->SetDirectory(0);
  76   fFMD3i->SetDirectory(0);
  77   fFMD3o->SetDirectory(0);
  78   fSumEta = new TH1D("mcSumEta", "MC INEL Truth", 200, -4, 6);
  79   fSumEta->SetXTitle("#eta");
  80   fSumEta->SetYTitle("dN_{ch}/d#eta");
  81   fSumEta->SetDirectory(0);
  82   fSumEta->Sumw2();
  83   fSumEta->SetMarkerColor(kOrange+2);
  84   fSumEta->SetMarkerStyle(22);
  85   fSumEta->SetFillColor(0);
  86   fSumEta->SetFillStyle(0);
  87
  88 }
  89
  90 //____________________________________________________________________
  91 AliFMDMCSharingFilter::AliFMDMCSharingFilter(const AliFMDMCSharingFilter& o)
  92   : AliFMDSharingFilter(o),
  93     fFMD1i(o.fFMD1i),
  94     fFMD2i(o.fFMD2i),
  95     fFMD2o(o.fFMD2o),
  96     fFMD3i(o.fFMD3i),
  97     fFMD3o(o.fFMD3o),
  98     fSumEta(o.fSumEta)
  99 {
 100   //
 101   // Copy constructor
 102   //
 103   // Parameters:
 104   //    o Object to copy from
 105   //
 106 }
 107
 108 //____________________________________________________________________
 109 AliFMDMCSharingFilter::~AliFMDMCSharingFilter()
 110 {
 111   //
 112   // Destructor
 113   //
 114   if (fFMD1i)  delete fFMD1i;
 115   if (fFMD2i)  delete fFMD2i;
 116   if (fFMD2o)  delete fFMD2o;
 117   if (fFMD3i)  delete fFMD3i;
 118   if (fFMD3o)  delete fFMD3o;
 119   if (fSumEta) delete fSumEta;
 120 }
 121
 122 //____________________________________________________________________
 123 AliFMDMCSharingFilter&
 124 AliFMDMCSharingFilter::operator=(const AliFMDMCSharingFilter& o)
 125 {
 126   //
 127   // Assignment operator
 128   //
 129   // Parameters:
 130   //    o Object to assign from
 131   //
 132   // Return:
 133   //    Reference to this
 134   //
 135   AliFMDSharingFilter::operator=(o);
 136   return *this;
 137 }
 138
 139 //____________________________________________________________________
 140 void
 141 AliFMDMCSharingFilter::StoreParticle(UShort_t   d,
 142                                      Char_t     r,
 143                                      UShort_t   s,
 144                                      UShort_t   t,
 145                                      AliESDFMD& output) const
 146 {
 147   //
 148   // Store a particle hit in FMD<i>dr</i>[<i>s,t</i>] in @a output
 149   //
 150   // Parameters:
 151   //    d       Detector
 152   //    r       Ring
 153   //    s       Sector
 154   //    t       Strip
 155   //    output  Output ESD object
 156   //
 157   Double_t old = output.Multiplicity(d,r,s,t);
 158   if (old == AliESDFMD::kInvalidMult) old = 0;
 159   output.SetMultiplicity(d,r,s,t,old+1);
 160 }
 161
 162 //____________________________________________________________________
 163 Bool_t
 164 AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input,
 165                                 const AliMCEvent& event,
 166                                 Double_t          vz,
 167                                 AliESDFMD&        output,
 168                                 TH2D*             primary)
 169 {
 170   //
 171   // Filter the input kinematics and track references, using
 172   // some of the ESD information
 173   //
 174   // Parameters:
 175   //    input   Input ESD event
 176   //    event   Input MC event
 177   //    vz      Vertex position
 178   //    output  Output ESD-like object
 179   //    primary Per-event histogram of primaries
 180   //
 181   // Return:
 182   //    True on succes, false otherwise
 183   //
 184   output.Clear();
 185
 186   // Copy eta values to output
 187   for (UShort_t ed = 1; ed <= 3; ed++) {
 188     UShort_t nq = (ed == 1 ? 1 : 2);
 189     for (UShort_t eq = 0; eq < nq; eq++) {
 190       Char_t   er = (eq == 0 ? 'I' : 'O');
 191       UShort_t ns = (eq == 0 ?  20 :  40);
 192       UShort_t nt = (eq == 0 ? 512 : 256);
 193       for (UShort_t es = 0; es < ns; es++)
 194         for (UShort_t et = 0; et < nt; et++)
 195           output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
 196     }
 197   }
 198   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
 199   Int_t nTracks   = stack->GetNtrack();//event.GetNumberOfTracks();
 200   Int_t nPrim     = stack->GetNprimary();//event.GetNumberOfPrimary();
 201   for (Int_t iTr = 0; iTr < nTracks; iTr++) {
 202     AliMCParticle* particle =
 203       static_cast<AliMCParticle*>(event.GetTrack(iTr));
 204
 205     // Check the returned particle
 206     if (!particle) continue;
 207
 208     // Check if this charged and a primary
 209     Bool_t isCharged = particle->Charge() != 0;
 210     if (!isCharged) continue;
 211     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
 212
 213     // Fill 'dn/deta' histogram
 214     if (isPrimary && iTr < nPrim) {
 215       fSumEta->Fill(particle->Eta());
 216       primary->Fill(particle->Eta(), particle->Phi());
 217     }
 218
 219     Int_t    nTrRef  = particle->GetNumberOfTrackReferences();
 220     Int_t    longest = -1;
 221     Double_t angle   = 0;
 222     UShort_t oD = 0, oS = 1024, oT = 1024;
 223     Char_t   oR = '\0';
 224     for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
 225       AliTrackReference* ref = particle->GetTrackReference(iTrRef);
 226
 227       // Check existence
 228       if (!ref) continue;
 229
 230       // Check that we hit an FMD element
 231       if (ref->DetectorId() != AliTrackReference::kFMD)
 232         continue;
 233
 234       // Get the detector coordinates
 235       UShort_t d, s, t;
 236       Char_t r;
 237       AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
 238       if (oD > 0 && oR != '\0' && (d != oD || r != oR)) {
 239         longest = -1;
 240         StoreParticle(oD, oR, oS, oT, output);
 241       }
 242
 243       oD = d;
 244       oR = r;
 245       oS = s;
 246       oT = t;
 247
 248       // The longest passage is determined through the angle
 249       Double_t x    = ref->X();
 250       Double_t y    = ref->Y();
 251       Double_t z    = ref->Z()-vz;
 252       Double_t rr   = TMath::Sqrt(x*x+y*y);
 253       Double_t theta= TMath::ATan2(rr,z);
 254       Double_t ang  = TMath::Abs(TMath::Pi()-theta);
 255       if (ang > angle) {
 256         longest = iTrRef;
 257         angle   = ang;
 258       }
 259     } // Loop over track references
 260     if (longest < 0) continue;
 261
 262     // Get the reference corresponding to the longest path through the detector
 263     AliTrackReference* ref = particle->GetTrackReference(longest);
 264
 265     // Get the detector coordinates
 266     UShort_t d, s, t;
 267     Char_t r;
 268     AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
 269
 270     StoreParticle(d,r,s,t,output);
 271   } // Loop over tracks
 272   return kTRUE;
 273 }
 274
 275 //____________________________________________________________________
 276 void
 277 AliFMDMCSharingFilter::CompareResults(const AliESDFMD&  esd,
 278                                       const AliESDFMD&  mc)
 279 {
 280   //
 281   // Compare the result of merging to the monte-carlo truth.  This
 282   // fills the correlation histograms
 283   //
 284   // Parameters:
 285   //    esd  ESD after sharing correction
 286   //    mc   MC ESD
 287   //
 288
 289   // Copy eta values to output
 290   for (UShort_t d = 1; d <= 3; d++) {
 291     UShort_t nq = (d == 1 ? 1 : 2);
 292     for (UShort_t q = 0; q < nq; q++) {
 293       Char_t   r  = (q == 0 ? 'I' : 'O');
 294       UShort_t ns = (q == 0 ?  20 :  40);
 295       UShort_t nt = (q == 0 ? 512 : 256);
 296       TH2*     co = 0;
 297       switch (d) {
 298       case 1: co = fFMD1i; break;
 299       case 2: co = (q == 0 ? fFMD2i : fFMD2o); break;
 300       case 3: co = (q == 0 ? fFMD3i : fFMD3o); break;
 301       }
 302
 303       for (UShort_t s = 0; s < ns; s++) {
 304         for (UShort_t t = 0; t < nt; t++) {
 305           Float_t mEsd = esd.Multiplicity(d, r, s, t);
 306           Float_t mMc  = mc.Multiplicity(d, r, s, t);
 307
 308           co->Fill(mMc, mEsd);
 309         }
 310       }
 311     }
 312   }
 313 }
 314
 315 //____________________________________________________________________
 316 void
 317 AliFMDMCSharingFilter::DefineOutput(TList* dir)
 318 {
 319   //
 320   // Define the output histograms.  These are put in a sub list of the
 321   // passed list.   The histograms are merged before the parent task calls
 322   // AliAnalysisTaskSE::Terminate
 323   //
 324   // Parameters:
 325   //    dir Directory to add to
 326   //
 327   AliFMDSharingFilter::DefineOutput(dir);
 328   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 329   TList* cd = new TList;
 330   cd->SetName("esd_mc_comparion");
 331   d->Add(cd);
 332   cd->Add(fFMD1i);
 333   cd->Add(fFMD2i);
 334   cd->Add(fFMD2o);
 335   cd->Add(fFMD3i);
 336   cd->Add(fFMD3o);
 337   dir->Add(fSumEta);
 338 }
 339
 340 //____________________________________________________________________
 341 void
 342 AliFMDMCSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
 343 {
 344   //
 345   // Scale the histograms to the total number of events
 346   //
 347   // Parameters:
 348   //    dir     Where the output is
 349   //    nEvents Number of events
 350   //
 351   AliFMDSharingFilter::ScaleHistograms(dir, nEvents);
 352   TH1D* sumEta = static_cast<TH1D*>(dir->FindObject("mcSumEta"));
 353   if (!sumEta) {
 354     AliWarning(Form("No mcSumEta histogram found in output list"));
 355     return;
 356   }
 357   sumEta->Scale(1. / nEvents, "width");
 358 }
 359
 360 //____________________________________________________________________
 361 void
 362 AliFMDMCSharingFilter::Print(Option_t* option) const
 363 {
 364   //
 365   // Print information
 366   //
 367   // Parameters:
 368   //    option Not used
 369   //
 370   char ind[gROOT->GetDirLevel()+1];
 371   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 372   ind[gROOT->GetDirLevel()] = '\0';
 373   AliFMDSharingFilter::Print(option);
 374 }
 375
 376 //____________________________________________________________________
 377 //
 378 // EOF
 379 //