PWG2/RESONANCES/AliRsnValue.cxx

   1
   2 //
   3 // Class AliRsnValue
   4 //
   5 // Definition of a single value which can be computed
   6 // from any of the defined input objects implemented
   7 // in the resonance package.
   8 //
   9
  10 #include "AliRsnEvent.h"
  11 #include "AliRsnDaughter.h"
  12 #include "AliRsnMother.h"
  13 #include "AliRsnPairDef.h"
  14
  15 #include "AliRsnValue.h"
  16
  17 ClassImp(AliRsnValue)
  18
  19 //_____________________________________________________________________________
  20 AliRsnValue::AliRsnValue() :
  21   TNamed(),
  22   fValue(0.0),
  23   fType(kValueTypes),
  24   fArray(0),
  25   fESDCuts()
  26 {
  27 //
  28 // Main constructor (version 1)
  29 // This can also be created without any argument.
  30 //
  31 }
  32
  33 //_____________________________________________________________________________
  34 AliRsnValue::AliRsnValue
  35 (const char *name, EValueType type, Int_t nbins, Double_t min, Double_t max) :
  36   TNamed(name, ""),
  37   fValue(0.0),
  38   fType(type),
  39   fArray(0),
  40   fESDCuts()
  41 {
  42 //
  43 // Main constructor (version 1)
  44 // This can also be created without any argument.
  45 //
  46
  47   SetBins(nbins, min, max);
  48 }
  49
  50 //_____________________________________________________________________________
  51 AliRsnValue::AliRsnValue
  52 (const char *name, EValueType type, Double_t min, Double_t max, Double_t step) :
  53   TNamed(name, ""),
  54   fValue(0.0),
  55   fType(type),
  56   fArray(0),
  57   fESDCuts()
  58 {
  59 //
  60 // Main constructor (version 2)
  61 //
  62
  63   SetBins(min, max, step);
  64 }
  65
  66 //_____________________________________________________________________________
  67 AliRsnValue::AliRsnValue
  68 (const char *name, EValueType type, Int_t nbins, Double_t *array) :
  69   TNamed(name, ""),
  70   fValue(0.0),
  71   fType(type),
  72   fArray(0),
  73   fESDCuts()
  74 {
  75 //
  76 // Main constructor (version 2)
  77 //
  78
  79   SetBins(nbins, array);
  80 }
  81
  82 //_____________________________________________________________________________
  83 AliRsnValue::AliRsnValue(const AliRsnValue& copy) :
  84   TNamed(copy),
  85   fValue(copy.fValue),
  86   fType(copy.fType),
  87   fArray(copy.fArray),
  88   fESDCuts(copy.fESDCuts)
  89 {
  90 //
  91 // Copy constructor
  92 //
  93 }
  94
  95 //_____________________________________________________________________________
  96 AliRsnValue& AliRsnValue::operator=(const AliRsnValue& copy)
  97 {
  98 //
  99 // Assignment operator
 100 //
 101
 102   SetName(copy.GetName());
 103
 104   fType = copy.fType;
 105   fValue = copy.fValue;
 106   fArray = copy.fArray;
 107   fESDCuts = copy.fESDCuts;
 108
 109   return (*this);
 110 }
 111
 112 //_____________________________________________________________________________
 113 void AliRsnValue::SetBins(Int_t nbins, Double_t min, Double_t max)
 114 {
 115 //
 116 // Set binning for the axis in equally spaced bins
 117 // where the number of bins, minimum and maximum are given.
 118 //
 119
 120   fArray.Set(nbins + 1);
 121
 122   Double_t mymax = TMath::Max(min, max);
 123   Double_t mymin = TMath::Min(min, max);
 124
 125   Int_t    k = 0;
 126   Double_t binSize = (mymax - mymin) / ((Double_t)nbins);
 127
 128   fArray[0] = mymin;
 129   for (k = 1; k <= nbins; k++) fArray[k] = fArray[k-1] + binSize;
 130   for (k = 0; k < fArray.GetSize() - 1; k++) AliDebug(AliLog::kDebug + 3, Form("Bin #%d: %f - %f", k, fArray[k], fArray[k+1]));
 131 }
 132
 133 //_____________________________________________________________________________
 134 void AliRsnValue::SetBins(Double_t min, Double_t max, Double_t step)
 135 {
 136 //
 137 // Set binning for the axis in equally spaced bins
 138 // where the bin size, minimum and maximum are given.
 139 //
 140
 141   Double_t dblNbins = TMath::Abs(max - min) / step;
 142   Int_t    intNbins = ((Int_t)dblNbins) + 1;
 143
 144   SetBins(intNbins, min, max);
 145 }
 146
 147 //_____________________________________________________________________________
 148 void AliRsnValue::SetBins(Int_t nbins, Double_t *array)
 149 {
 150 //
 151 // Set binning for the axis in unequally spaced bins
 152 // using the same way it is done in TAxis
 153 //
 154
 155   fArray.Adopt(nbins, array);
 156   for (Int_t k = 0; k < fArray.GetSize() - 1; k++) AliDebug(AliLog::kDebug + 3, Form("Bin #%d: %f - %f", k, fArray[k], fArray[k+1]));
 157 }
 158
 159 //_____________________________________________________________________________
 160 Bool_t AliRsnValue::Eval(AliRsnMother * const mother, AliRsnPairDef * const pairDef, AliRsnEvent * const event)
 161 {
 162 //
 163 // Evaluation of the required value.
 164 // Checks that the passed object is of the right type
 165 // and if this check is successful, returns the required value.
 166 // The output of the function tells if it was successful,
 167 // and the values must be taken with GetValue().
 168 //
 169
 170   // avoid segfaults
 171   if (!mother) return kFALSE;
 172   if (!pairDef) return kFALSE;
 173
 174   Double_t mass = pairDef->GetMotherMass();
 175
 176   switch (fType)
 177   {
 178     case kTrack1P:
 179       fValue = mother->GetDaughter(0)->P().Mag();
 180       break;
 181     case kTrack2P:
 182       fValue = mother->GetDaughter(1)->P().Mag();
 183       break;
 184     case kTrack1Pt:
 185       fValue = mother->GetDaughter(0)->P().Perp();
 186       break;
 187     case kTrack2Pt:
 188       fValue = mother->GetDaughter(1)->P().Perp();
 189       break;
 190     case kTrack1Px:
 191       fValue = mother->GetDaughter(0)->P().X();
 192       break;
 193     case kTrack1Py:
 194       fValue = mother->GetDaughter(0)->P().Y();
 195       break;
 196     case kTrack1Pz:
 197       fValue = mother->GetDaughter(0)->P().Z();
 198       break;
 199     case kTrack2Px:
 200       fValue = mother->GetDaughter(1)->P().X();
 201       break;
 202     case kTrack2Py:
 203       fValue = mother->GetDaughter(1)->P().Y();
 204       break;
 205     case kTrack2Pz:
 206       fValue = mother->GetDaughter(1)->P().Z();
 207       break;
 208     case kPairInvMass:
 209       fValue = mother->Sum().M();
 210       break;
 211     case kPairInvMassMC:
 212       fValue = mother->SumMC().M();
 213       break;
 214     case kPairInvMassRes:
 215       fValue = (mother->SumMC().M() - mother->Sum().M()) / mother->SumMC().M();
 216       break;
 217     case kPairPt:
 218       fValue = mother->Sum().Perp();
 219       break;
 220     case kPairEta:
 221       fValue = mother->Sum().Eta();
 222       break;
 223     case kPairMt:
 224       if (TMath::Abs(mass) < 1E-5) AliWarning(Form("Suspicious mass value specified: %f", mass));
 225       fValue = (TMath::Sqrt(mother->Sum().Perp2() + mass*mass) - mass);
 226       break;
 227     case kPairY:
 228       if (TMath::Abs(mass) < 1E-5) AliWarning(Form("Suspicious mass value specified: %f", mass));
 229       mother->SetDefaultMass(mass);
 230       fValue = mother->Ref().Rapidity();
 231       break;
 232     case kPairPhi:
 233       fValue = mother->Sum().Phi();
 234       break;
 235     case kPairPhiMC:
 236       fValue = mother->SumMC().Phi();
 237       break;
 238     case kPairPtRatio:
 239       fValue  = TMath::Abs(mother->GetDaughter(0)->P().Perp() - mother->GetDaughter(1)->P().Perp());
 240       fValue /= TMath::Abs(mother->GetDaughter(0)->P().Perp() + mother->GetDaughter(1)->P().Perp());
 241       break;
 242     case kPairDipAngle:
 243       fValue = mother->GetDaughter(0)->P().Angle(mother->GetDaughter(1)->P().Vect());
 244       fValue = TMath::Abs(TMath::ACos(fValue));
 245       break;
 246     case kPairCosThetaStar:
 247       fValue = mother->CosThetaStar();
 248       break;
 249     case kPairCosThetaStar1:
 250       //fValue = TMath::Cos(mother->ThetaStar(kTRUE, kFALSE));
 251       break;
 252     case kPairCosThetaStar2:
 253       //fValue = TMath::Cos(mother->ThetaStar(kFALSE, kFALSE));
 254       break;
 255     case kPairCosThetaStarMC1:
 256       //fValue = TMath::Cos(mother->ThetaStar(kTRUE, kTRUE));
 257       break;
 258     case kPairCosThetaStarMC2:
 259       //fValue = TMath::Cos(mother->ThetaStar(kFALSE, kTRUE));
 260       break;
 261     case kAngleToLeading:
 262       {
 263           int ID1 = (mother->GetDaughter(0))->GetID();
 264           int ID2 = (mother->GetDaughter(1))->GetID();
 265           int leadingID = event->SelectLeadingParticle(0);
 266           if(leadingID == ID1 || leadingID == ID2) return kFALSE;
 267           AliRsnDaughter  leadingPart = event->GetDaughter(leadingID);
 268           AliVParticle *ref = leadingPart.GetRef();
 269
 270           fValue = ref->Phi() - mother->Sum().Phi();
 271           //return angle w.r.t. leading particle in the range -pi/2, 3/2pi
 272           while(fValue >= TMath::Pi()) fValue -= 2*TMath::Pi();
 273           while(fValue < -0.5*TMath::Pi()) fValue += 2*TMath::Pi();
 274           //Printf("%g", fValue);
 275       }
 276       break;
 277     case kEventMult:
 278       if (!event)
 279       {
 280         fValue = 0.0;
 281         return kFALSE;
 282       }
 283       else fValue = (Double_t)event->GetMultiplicity();
 284       break;
 285     case kEventMultESDcuts:
 286       if (!event)
 287       {
 288         fValue = 0.0;
 289         return kFALSE;
 290       }
 291       else
 292       {
 293         AliESDEvent *esd = event->GetRefESD();
 294         if (!esd)
 295         {
 296           AliError("Cannot use method based on ESD cuts when input is not ESD.");
 297           fValue = 0.0;
 298           return kFALSE;
 299         }
 300         fValue = (Double_t)fESDCuts.CountAcceptedTracks(esd);
 301       }
 302     case kLeadingPt:
 303       if (!event)
 304       {
 305         fValue = 0.0;
 306         return kFALSE;
 307       }
 308       else
 309       {
 310           int leadingID = event->SelectLeadingParticle(0);
 311           if(leadingID >= 0) {
 312                   AliRsnDaughter leadingPart = event->GetDaughter(leadingID);
 313                   AliVParticle *ref = leadingPart.GetRef();
 314                   fValue = ref->Pt();
 315           }
 316           else fValue = 0;
 317       }
 318       break;
 319     case kQInv:
 320       {
 321         TLorentzVector diff = mother->GetDaughter(0)->P() - mother->GetDaughter(1)->P();
 322         fValue = diff.M();
 323       }
 324       break;
 325     default:
 326       AliWarning("Invalid value type");
 327       return kFALSE;
 328   }
 329
 330   return kTRUE;
 331 }
 332
 333 //_____________________________________________________________________________
 334 Bool_t AliRsnValue::Eval(AliRsnDaughter * const daughter, AliRsnEvent * const event)
 335 {
 336 //
 337 // Evaluation of the required value.
 338 // Checks that the passed object is of the right type
 339 // and if this check is successful, returns the required value.
 340 // The output of the function tells if it was successful,
 341 // and the values must be taken with GetValue().
 342 //
 343
 344   // avoid segfaults
 345   if (!daughter) return kFALSE;
 346
 347   switch (fType)
 348   {
 349     case kEventMult:
 350       if (!event)
 351       {
 352         fValue = 0.0;
 353         return kFALSE;
 354       }
 355       else fValue = (Double_t)event->GetMultiplicity();
 356       break;
 357     case kLeadingPt:
 358       if (!event)
 359       {
 360         fValue = 0.0;
 361         return kFALSE;
 362       }
 363       else
 364       {
 365           int leadingID = event->SelectLeadingParticle(0);
 366           if(leadingID >= 0) {
 367                   AliRsnDaughter leadingPart = event->GetDaughter(leadingID);
 368                   AliVParticle *ref = leadingPart.GetRef();
 369                   fValue = ref->Pt();
 370           }
 371           else fValue = 0;
 372       }
 373       break;
 374     default:
 375       AliWarning("Invalid value type");
 376       return kFALSE;
 377   }
 378
 379   return kTRUE;
 380 }
 381
 382 //_____________________________________________________________________________
 383 void AliRsnValue::Print(Option_t *) const
 384 {
 385 //
 386 // Print all bins
 387 //
 388
 389   Int_t   i, n = fArray.GetSize();
 390   TString msg("Array values: ");
 391
 392   for (i = 0; i < n; i++) msg += Form("%f, ", fArray[i]);
 393
 394   AliInfo(Form("Axis name: %s", GetName()));
 395   AliInfo(msg.Data());
 396 }