PWG2/RESONANCES/AliRsnFunction.cxx

   1 //
   2 // Class AliRsnFunction
   3 //
   4 // This class defines a base classe to implement a function
   5 // which uses the internal RSN package event format (AliRsnEvent).
   6 // It contains some default flags which turn out to be useful:
   7 //  - a flag to select only the "true" pairs (tracks from same resonance)
   8 //  - a flag to know if the computation is done over two events (mixing)
   9 //
  10 // Any kind of analysis object should be implemented as inheriting from this
  11 // because the AliRsnAnalyzer which executes the analysis will accept a collection
  12 // of such objects, in order to have a unique format of processing method
  13 //
  14 // The user who implements a kind of computation type should inherit from
  15 // this class and override the virtual functions defined in it, which
  16 // initialize the final output histogram and define how to process data.
  17 //
  18 //
  19 // author: A. Pulvirenti             (email: alberto.pulvirenti@ct.infn.it)
  20 //
  21
  22 #include <Riostream.h>
  23
  24 #include <TH1.h>
  25 #include <TList.h>
  26 #include <TString.h>
  27
  28 #include "AliLog.h"
  29
  30 #include "AliRsnMCInfo.h"
  31 #include "AliRsnDaughter.h"
  32 #include "AliRsnEvent.h"
  33 #include "AliRsnPairDef.h"
  34 #include "AliRsnCut.h"
  35
  36 #include "AliRsnFunction.h"
  37
  38 ClassImp(AliRsnFunction)
  39
  40 //________________________________________________________________________________________
  41 AliRsnFunction::AliRsnFunction() :
  42     fFcnType(kFcnTypes),
  43     fRotAngle(0.0),
  44     fUseBins(kFALSE),
  45     fSkipOutsideInterval(kFALSE),
  46     fBins(0),
  47     fBinningCut(),
  48     fBinningCutType(AliRsnCut::kLastCutType),
  49     fHistoDef(0x0)
  50 {
  51   //
  52   // Constructor.
  53   // The histogram data member cannot be passed externally,
  54   // its initialization MUST be defined inside the Init() method,
  55   // which must be overridden in any derivate implementation.
  56   //
  57
  58   Int_t i;
  59   for (i = 0; i < 100; i++)
  60   {
  61     fHisto[i] = 0x0;
  62   }
  63 }
  64
  65 //________________________________________________________________________________________
  66 AliRsnFunction::AliRsnFunction
  67 (EFcnType type, AliRsnHistoDef *hd, Bool_t skipOut) :
  68     fFcnType(type),
  69     fRotAngle(0.0),
  70     fUseBins(kFALSE),
  71     fSkipOutsideInterval(skipOut),
  72     fBins(0),
  73     fBinningCut(),
  74     fBinningCutType(AliRsnCut::kLastCutType),
  75     fHistoDef(hd)
  76 {
  77   //
  78   // Constructor.
  79   // The histogram data member cannot be passed externally,
  80   // its initialization MUST be defined inside the Init() method,
  81   // which must be overridden in any derivate implementation.
  82   //
  83
  84   Int_t i;
  85   for (i = 0; i < 100; i++)
  86   {
  87     fHisto[i] = 0x0;
  88   }
  89 }
  90
  91 //________________________________________________________________________________________
  92 AliRsnFunction::AliRsnFunction(const AliRsnFunction &copy) :
  93     TObject(copy),
  94     fFcnType(copy.fFcnType),
  95     fRotAngle(copy.fRotAngle),
  96     fUseBins(copy.fUseBins),
  97     fSkipOutsideInterval(copy.fSkipOutsideInterval),
  98     fBins(0),
  99     fBinningCut(),
 100     fBinningCutType(AliRsnCut::kLastCutType),
 101     fHistoDef(copy.fHistoDef)
 102 {
 103   //
 104   // Copy constructor.
 105   // Calls the function to define binning.
 106   //
 107
 108   Int_t i, n = 100;
 109   for (i = 0; i < n; i++)
 110   {
 111     fHisto[i] = 0x0;
 112   }
 113
 114   if (fUseBins)
 115   {
 116     n = copy.fBins.GetSize();
 117     Double_t *array = new Double_t[n];
 118     for (i = 0; i < n; i++) array[i] = copy.fBins[i];
 119     SetBinningCut(copy.fBinningCutType, copy.fBins.GetSize(), array);
 120     delete [] array;
 121   }
 122 }
 123 //________________________________________________________________________________________
 124 const AliRsnFunction& AliRsnFunction::operator=(const AliRsnFunction& /*copy*/)
 125 {
 126   //
 127   // Assignment operator.
 128   // Behaves like copy constructor.
 129   // Also in this case, the histogram is not copied, and,
 130   // if it was present, it is destroyed and will need to be recreated.
 131   //
 132
 133   return (*this);
 134 }
 135 //________________________________________________________________________________________
 136 void AliRsnFunction::Clear(Option_t* /*option*/)
 137 {
 138   //
 139   // Clear arrays and histogram.
 140   // For the sake of security, all pointers are also set explicitly to NULL.
 141   //
 142
 143   Int_t i;
 144   for (i = 0; i < 100; i++)
 145   {
 146     delete fHisto[i];
 147     fHisto[i] = 0x0;
 148   }
 149 }
 150
 151 //________________________________________________________________________________________
 152 TList* AliRsnFunction::Init(const char *histoName, const char *histoTitle)
 153 {
 154   //
 155   // Initialization function.
 156   // By default, it initializes the owned histogram using the method
 157   // from AliRsnHistoDef class, giving the same name and title of this.
 158   // A user can override this behaviour, if necessary.
 159   // Before creating, the HistoDef is checked for proper initialization.
 160   //
 161
 162   Clear();
 163
 164   Int_t i, ibin, nbins = fHistoDef->GetNBins();
 165   Double_t min = fHistoDef->GetMin(), max = fHistoDef->GetMax();
 166
 167   // list is created and named after the general
 168   // settings used for the contained histograms
 169   TList *histos = new TList;
 170   histos->SetName(Form("%s", GetFcnName().Data()));
 171
 172   // a general histogram is always added,
 173   // which overrides the binning and collects everything
 174   fHisto[0] = new TH1D(histoName, histoTitle, nbins, min, max);
 175   histos->AddLast(fHisto[0]);
 176
 177   // if requested a binning w.r. to some cut variable, histograms are added
 178   // for that in this part of the method (one per each bin)
 179   Char_t hName[255];
 180   Char_t hTitle[255];
 181   if (fUseBins)
 182   {
 183     for (ibin = 0, i = 1; ibin < fBins.GetSize() - 1; ibin++, i++)
 184     {
 185       sprintf(hName, "%s[%.2f-%.2f]", histoName, fBins[ibin], fBins[ibin+1]);
 186       sprintf(hTitle, "%s [%.2f-%.2f]", histoTitle, fBins[ibin], fBins[ibin+1]);
 187       fHisto[i] = new TH1D(hName, hTitle, nbins, min, max);
 188       histos->AddLast(fHisto[i]);
 189     }
 190   }
 191
 192   // returns the full list at the end
 193   return histos;
 194 }
 195
 196 //________________________________________________________________________________________
 197 void AliRsnFunction::Init(const char *histoName, const char *histoTitle, TList *histos)
 198 {
 199   //
 200   // Initialization function.
 201   // By default, it initializes the owned histogram using the method
 202   // from AliRsnHistoDef class, giving the same name and title of this.
 203   // A user can override this behaviour, if necessary.
 204   // Before creating, the HistoDef is checked for proper initialization.
 205   //
 206
 207   Clear();
 208
 209   Int_t i, ibin, nbins = fHistoDef->GetNBins();
 210   Double_t min = fHistoDef->GetMin(), max = fHistoDef->GetMax();
 211
 212   // list is created and named after the general
 213   // settings used for the contained histograms
 214   if (!histos)
 215   {
 216     AliError("NULL target list!");
 217     return;
 218   }
 219
 220   // a general histogram is always added,
 221   // which overrides the binning and collects everything
 222   fHisto[0] = new TH1D(histoName, histoTitle, nbins, min, max);
 223   histos->AddLast(fHisto[0]);
 224
 225   // if requested a binning w.r. to some cut variable, histograms are added
 226   // for that in this part of the method (one per each bin)
 227   Char_t hName[255];
 228   Char_t hTitle[255];
 229   if (fUseBins)
 230   {
 231     for (ibin = 0, i = 1; ibin < fBins.GetSize() - 1; ibin++, i++)
 232     {
 233       sprintf(hName, "%s[%.2f-%.2f]", histoName, fBins[ibin], fBins[ibin+1]);
 234       sprintf(hTitle, "%s [%.2f-%.2f]", histoTitle, fBins[ibin], fBins[ibin+1]);
 235       fHisto[i] = new TH1D(hName, hTitle, nbins, min, max);
 236       histos->AddLast(fHisto[i]);
 237     }
 238   }
 239 }
 240
 241 //________________________________________________________________________________________
 242 void AliRsnFunction::SetBinningCut
 243 (AliRsnCut::EType type, Double_t min, Double_t max, Double_t step)
 244 {
 245   //
 246   // Set fixed bins
 247   //
 248
 249   fUseBins = kTRUE;
 250
 251   Int_t i, nBins = (Int_t)((max - min) / step) + 1;
 252   fBinningCutType = type;
 253   fBins.Set(nBins);
 254   for (i = 0; i < nBins; i++)
 255   {
 256     fBins[i] = min + (Double_t)i * step;
 257   }
 258 }
 259
 260 //________________________________________________________________________________________
 261 void AliRsnFunction::SetBinningCut
 262 (AliRsnCut::EType type, Int_t nbins, Double_t *bins)
 263 {
 264   //
 265   // Set variable bins
 266   //
 267
 268   fUseBins = kTRUE;
 269
 270   Int_t i;
 271   fBinningCutType = type;
 272   fBins.Set(nbins);
 273   for (i = 0; i < nbins; i++)
 274   {
 275     fBins[i] = bins[i];
 276   }
 277 }
 278
 279 //________________________________________________________________________________________
 280 TString AliRsnFunction::GetFcnName()
 281 {
 282   //
 283   // Return a string which names the function type
 284   //
 285
 286   TString text("Undef");
 287
 288   switch (fFcnType)
 289   {
 290     case kInvMass:
 291       text = "IM";
 292       break;
 293     case kInvMassMC:
 294       text = "IM_MC";
 295       break;
 296     case kInvMassRotated:
 297       text = Form("IMR%.2f", fRotAngle);
 298       break;
 299     case kResolution:
 300       text = "RES";
 301       break;
 302     case kPtSpectrum:
 303       text = "PT";
 304       break;
 305     case kEtaSpectrum:
 306       text = "ETA";
 307       break;
 308     case kAngleSpectrum:
 309       text = "ANGLE";
 310     default:
 311       AliError("Type not defined");
 312   }
 313
 314   return text;
 315 }
 316
 317 //________________________________________________________________________________________
 318 TString AliRsnFunction::GetFcnTitle()
 319 {
 320   //
 321   // Return a string which names the function type
 322   //
 323
 324   TString text("Undef");
 325
 326   switch (fFcnType)
 327   {
 328     case kInvMass:
 329       text = "Invariant mass";
 330       break;
 331     case kInvMassMC:
 332       text = "Invariant mass (MC)";
 333       break;
 334     case kResolution:
 335       text = "Resolution";
 336       break;
 337     case kPtSpectrum:
 338       text = "p_{#perp} distribution";
 339       break;
 340     case kEtaSpectrum:
 341       text = "#eta distribution";
 342       break;
 343     case kAngleSpectrum:
 344       text = "angle distribution";
 345       break;
 346     default:
 347       AliError("Type not defined");
 348   }
 349
 350   return text;
 351 }
 352
 353 //________________________________________________________________________________________
 354 Bool_t AliRsnFunction::Fill(AliRsnPairParticle *pair, AliRsnPairDef *ref, Double_t weight)
 355 {
 356   //
 357   // Fillse the histogram with data contained in a defined pair.
 358   // This method must be overidden by an appropriate definition in each inheriting class.
 359   //
 360
 361   Double_t value = FcnValue(pair, ref);
 362   if (fSkipOutsideInterval)
 363   {
 364     if (value < fHistoDef->GetMin()) return kFALSE;
 365     if (value > fHistoDef->GetMax()) return kFALSE;
 366   }
 367
 368   // fill global histogram
 369   if (weight == 0.0) fHisto[0]->Fill(value);
 370   else fHisto[0]->Fill(value, weight);
 371
 372   // if bins are allocated, find right one and fill it
 373   if (fUseBins)
 374   {
 375     Int_t i, ibin;
 376     for (ibin = 0, i = 1; ibin < fBins.GetSize() - 1; ibin++, i++)
 377     {
 378       if (!fHisto[i]) continue;
 379       fBinningCut.SetCutValues(fBinningCutType, fBins[ibin], fBins[ibin+1]);
 380       if (fBinningCut.IsSelected(AliRsnCut::kPair, pair))
 381       {
 382         if (weight == 0.0) fHisto[i]->Fill(value);
 383         else fHisto[i]->Fill(value, weight);
 384         break;
 385       }
 386     }
 387   }
 388
 389   return kTRUE;
 390 }
 391
 392 //________________________________________________________________________________________
 393 Double_t AliRsnFunction::FcnValue(AliRsnPairParticle *pair, AliRsnPairDef *ref)
 394 {
 395   //
 396   // This method must be overridden in all inheritin functions.
 397   // It computes the value which must be used to fill the histogram.
 398   //
 399
 400   switch (fFcnType)
 401   {
 402     case kInvMass:
 403       return pair->GetInvMass(ref->GetMass(0), ref->GetMass(1));
 404     case kInvMassMC:
 405       return pair->GetInvMassMC(ref->GetMass(0), ref->GetMass(1));
 406     case kInvMassRotated:
 407       //AliInfo(Form("*** ROTATION ANGLE = %f ***", fRotAngle));
 408       //AliInfo(Form("UNROTATED INV MASS = %f", pair->GetInvMass(ref->GetMass(0), ref->GetMass(1))));
 409       //pair->GetDaughter(1)->Print("P");
 410       pair->GetDaughter(1)->RotateP(fRotAngle * TMath::DegToRad());
 411       pair->ResetPair();
 412       //AliInfo(Form("  ROTATED INV MASS = %f", pair->GetInvMass(ref->GetMass(0), ref->GetMass(1))));
 413       //pair->GetDaughter(1)->Print("P");
 414       return pair->GetInvMass(ref->GetMass(0), ref->GetMass(1));
 415     case kResolution:
 416       return FcnResolution(pair, ref);
 417     case kPtSpectrum:
 418       return pair->GetPt();
 419     case kEtaSpectrum:
 420       return pair->GetEta();
 421     case kAngleSpectrum:
 422       return pair->GetAngle();
 423     default:
 424       AliError("Type not defined");
 425   }
 426
 427   return 0.0;
 428 }
 429
 430 //________________________________________________________________________________________
 431 inline Double_t AliRsnFunction::FcnResolution(AliRsnPairParticle *pair, AliRsnPairDef *ref)
 432 {
 433   //
 434   // Invariant mass resolution (compared between reconstructed and montecarlo)
 435   //
 436
 437   Double_t recInvMass = pair->GetInvMass(ref->GetMass(0), ref->GetMass(1));
 438   Double_t simInvMass = pair->GetInvMassMC(ref->GetMass(0), ref->GetMass(1));
 439
 440   return (simInvMass - recInvMass) / simInvMass;
 441 }