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