PWG2/RESONANCES/AliRsnEvent.cxx

   1 //
   2 // *** Class AliRsnEvent ***
   3 //
   4 // A container for a collection of AliRsnDaughter objects from an event.
   5 // Contains also the primary vertex, useful for some cuts.
   6 // In order to retrieve easily the tracks which have been identified
   7 // as a specific type and charge, there is an array of indexes which
   8 // allows to avoid to loop on all tracks and have only the neede ones.
   9 //
  10 // authors: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
  11 //          M. Vala (email: martin.vala@cern.ch)
  12 //
  13
  14 #include <TArrayF.h>
  15
  16 #include "AliLog.h"
  17 #include "AliVEvent.h"
  18 #include "AliESDEvent.h"
  19 #include "AliAODEvent.h"
  20 #include "AliMCEvent.h"
  21 #include "AliStack.h"
  22 #include "AliGenEventHeader.h"
  23
  24 #include "AliRsnEvent.h"
  25
  26 ClassImp(AliRsnEvent)
  27
  28 //_____________________________________________________________________________
  29 AliRsnEvent::AliRsnEvent(AliVEvent *ref, AliMCEvent *refMC) :
  30     fRef(ref),
  31     fRefMC(refMC),
  32     fPIDDefESD()
  33 {
  34 //
  35 // Default constructor.
  36 // Set the prior probabilities to some default values
  37 //
  38
  39   fPrior[0] = 0.02;
  40   fPrior[1] = 0.02;
  41   fPrior[2] = 0.83;
  42   fPrior[3] = 0.07;
  43   fPrior[4] = 0.06;
  44 }
  45
  46 //_____________________________________________________________________________
  47 AliRsnEvent::AliRsnEvent(const AliRsnEvent &event) :
  48     TObject(event),
  49     fRef(event.fRef),
  50     fRefMC(event.fRefMC),
  51     fPIDDefESD(event.fPIDDefESD)
  52 {
  53 //
  54 // Copy constructor.
  55 //
  56 }
  57
  58 //_____________________________________________________________________________
  59 AliRsnEvent& AliRsnEvent::operator= (const AliRsnEvent & event)
  60 {
  61 //
  62 // Works in the same way as the copy constructor.
  63 //
  64
  65   (TObject)(*this) = (TObject)event;
  66   fRef = event.fRef;
  67   fRefMC = event.fRefMC;
  68   fPIDDefESD = event.fPIDDefESD;
  69
  70   return (*this);
  71 }
  72
  73 //_____________________________________________________________________________
  74 AliRsnEvent::~AliRsnEvent()
  75 {
  76 //
  77 // Destructor.
  78 //
  79 }
  80
  81 //_____________________________________________________________________________
  82 void AliRsnEvent::SetDaughter(AliRsnDaughter &out, Int_t i)
  83 {
  84 //
  85 // Return a track stored here in format of AliRsnDaughter.
  86 // and finds in the reference event the informations to set
  87 // the proprietary data members of AliRsnDaughter
  88 //
  89
  90   // retrieve reference particle from reference event
  91   // if it is found, by defaul track can be used (good)
  92   AliVParticle *ref = (AliVParticle*)fRef->GetTrack(i);
  93   if (!ref) return;
  94   out.SetRef(ref);
  95   out.SetGood();
  96
  97   // if MC info is present, retrieve from it
  98   TParticle *refMC = 0;
  99   if (fRefMC) {
 100     Int_t label = TMath::Abs(ref->GetLabel());
 101     refMC = fRefMC->Stack()->Particle(label);
 102     out.SetParticle(refMC);
 103     out.FindMotherPDG(fRefMC->Stack());
 104   }
 105
 106   // if fRef is MC event return
 107   AliMCEvent *mc = dynamic_cast<AliMCEvent *>(fRef);
 108   if (mc) return;
 109
 110   // dynamic reference to true nature of referenced event
 111   // to get kink index
 112   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(fRef);
 113   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fRef);
 114
 115   if (esd) {
 116     // retrieve vertex and set impact parameters
 117     Double_t dx = out.Xv(), dy = out.Yv(), dz = out.Zv();
 118     const AliVVertex *v = fRef->GetPrimaryVertex();
 119     if (v) {
 120       dx -= v->GetX();
 121       dy -= v->GetY();
 122       dz -= v->GetZ();
 123     } else if (fRefMC) {
 124       // if reference is an MC event, no primary vertex is supplied
 125       // but it is possible to retrieve it from header
 126       TArrayF fvertex(3);
 127       fRefMC->GenEventHeader()->PrimaryVertex(fvertex);
 128       dx -= fvertex[0];
 129       dy -= fvertex[1];
 130       dz -= fvertex[2];
 131     }
 132     out.SetDr(TMath::Sqrt(dx*dx + dy*dy));
 133     out.SetDz(dz);
 134   }
 135   // compute PID probabilities by combining
 136   // the PID weights in the source with priors
 137   // and eventually using the PIDDefESD
 138   // (the AliRsnDaughter objec knows how to manage the latter)
 139   out.CombineWithPriors(fPrior, &fPIDDefESD);
 140
 141   if (esd) {
 142     AliESDtrack *esdTrack = esd->GetTrack(i);
 143     out.FindKinkIndex(esdTrack);
 144   } else if (aod) {
 145     out.FindKinkIndex(aod);
 146   }
 147 }
 148
 149 //_____________________________________________________________________________
 150 AliRsnDaughter AliRsnEvent::GetDaughter(Int_t i)
 151 {
 152 //
 153 // Return an AliRsnDaughter taken from this event,
 154 // with all additional data members well set.
 155 //
 156
 157   AliRsnDaughter out;
 158   SetDaughter(out, i);
 159
 160   return out;
 161 }
 162
 163 //_____________________________________________________________________________
 164 Int_t AliRsnEvent::GetMultiplicity()
 165 {
 166 //
 167 // Returns event multiplicity
 168 //
 169   AliDebug(AliLog::kDebug+2,"<-");
 170   if (!fRef) return 0;
 171   AliDebug(AliLog::kDebug+2,"->");
 172   return fRef->GetNumberOfTracks();
 173 }
 174
 175 //_____________________________________________________________________________
 176 Double_t AliRsnEvent::GetVz()
 177 {
 178 //
 179 // Return Z coord of primary vertex
 180 //
 181   AliDebug(AliLog::kDebug+2,"<-");
 182   return fRef->GetPrimaryVertex()->GetZ();
 183   AliDebug(AliLog::kDebug+2,"->");
 184 }
 185
 186 //_____________________________________________________________________________
 187 AliRsnDaughter AliRsnEvent::GetLeadingParticle
 188 (Double_t ptMin, AliPID::EParticleType type)
 189 {
 190 //
 191 // Searches the collection of all particles with given PID type and charge,
 192 // and returns the one with largest momentum, provided that it is greater than 1st argument.
 193 // If one specifies AliRsnPID::kUnknown as type or AliRsnDaughter::kNoPID as method,
 194 // the check is done over all particles irrespectively of their PID.
 195 // If the sign argument is '+' or '-', the check is done over the particles of that charge,
 196 // otherwise it is done irrespectively of the charge.
 197 //
 198
 199   Int_t i, nTracks = fRef->GetNumberOfTracks();
 200   AliRsnDaughter output;
 201
 202   for (i = 0; i < nTracks; i++) {
 203     AliRsnDaughter track = GetDaughter(i);
 204     if (!AcceptTrackPID(&track, type)) continue;
 205     if (track.Pt() < ptMin) continue;
 206     if (!output.IsOK() || track.Pt() > output.Pt()) {
 207       output = track;
 208       output.SetGood();
 209     }
 210   }
 211
 212   return output;
 213 }
 214
 215 //_________________________________________________________________________________________________
 216 Double_t AliRsnEvent::GetAverageMomentum(Int_t &count, AliPID::EParticleType type)
 217 {
 218 //
 219 // Loops on the list of tracks and computes average total momentum.
 220 //
 221
 222   Int_t i, nTracks = fRef->GetNumberOfTracks();
 223   Double_t pmean = 0.0;
 224
 225   for (i = 0, count = 0; i < nTracks; i++) {
 226     AliRsnDaughter track = GetDaughter(i);
 227     if (!AcceptTrackPID(&track, type)) continue;
 228     pmean += track.P();
 229     count++;
 230   }
 231
 232   if (count > 0) pmean /= (Double_t)count;
 233   else pmean = 0.0;
 234
 235   return pmean;
 236 }
 237
 238 //_____________________________________________________________________________
 239 Bool_t AliRsnEvent::GetAngleDistr
 240 (Double_t &angleMean, Double_t &angleRMS, AliRsnDaughter leading)
 241 {
 242 //
 243 // Takes the leading particle and computes the mean and RMS
 244 // of the distribution of directions of all other tracks
 245 // with respect to the direction of leading particle.
 246 //
 247
 248   if (!leading.IsOK()) return kFALSE;
 249
 250   Int_t i, count, nTracks = fRef->GetNumberOfTracks();
 251   Double_t angle, angle2Mean = 0.0;
 252
 253   angleMean = angle2Mean = 0.0;
 254
 255   for (i = 0, count = 0; i < nTracks; i++) {
 256     AliRsnDaughter trk = GetDaughter(i);
 257     if (trk.GetID() == leading.GetID()) continue;
 258
 259     angle = leading.AngleTo(trk);
 260
 261     angleMean += angle;
 262     angle2Mean += angle * angle;
 263     count++;
 264   }
 265
 266   if (!count) return kFALSE;
 267
 268   angleMean /= (Double_t)count;
 269   angle2Mean /= (Double_t)count;
 270   angleRMS = TMath::Sqrt(angle2Mean - angleMean * angleMean);
 271
 272   return kTRUE;
 273 }
 274
 275 //_____________________________________________________________________________
 276 void AliRsnEvent::SetPriorProbability(Double_t *const out)
 277 {
 278 //
 279 // Set all prior probabilities at once, using an assayr of values.
 280 //
 281
 282   Int_t i;
 283
 284   for (i = 0; i < AliPID::kSPECIES; i++) {
 285     fPrior[i] = out[i];
 286   }
 287 }
 288
 289 //_____________________________________________________________________________
 290 void AliRsnEvent::DumpPriors()
 291 {
 292 //
 293 // Print all prior probabilities.
 294 // Printout is done using AliInfo, so this will not appear when
 295 // the GlobalLogLevel is set to higher level errors.
 296 //
 297
 298   Int_t i;
 299
 300   for (i = 0; i < AliPID::kSPECIES; i++) {
 301     AliInfo(Form("Prior probability for %10s = %3.5f", AliPID::ParticleName((AliPID::EParticleType)i), fPrior[i]));
 302   }
 303 }
 304
 305 //_____________________________________________________________________________
 306 void AliRsnEvent::GetPriorProbability(Double_t *out) const
 307 {
 308 //
 309 // Stores in the passed argument all the values.
 310 //
 311
 312   Int_t i;
 313
 314   for (i = 0; i < AliPID::kSPECIES; i++) {
 315     out[i] = fPrior[i];
 316   }
 317
 318 }
 319
 320 //_____________________________________________________________________________
 321 Bool_t AliRsnEvent::AcceptTrackPID(AliRsnDaughter * const d, AliPID::EParticleType type)
 322 {
 323 //
 324 // [PRIVATE]
 325 // Checks if the track PID (according to method in use) corresponds
 326 // to the required identification species.
 327 // If the second argument is "kUnknown", answer of this method is always YES.
 328 //
 329
 330   if (type == AliPID::kUnknown) return kTRUE;
 331
 332   return (d->AssignedPID() == type);
 333 }