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   // retrieve vertex and set impact parameters
 111   Double_t dx = out.Xv(), dy = out.Yv(), dz = out.Zv();
 112   const AliVVertex *v = fRef->GetPrimaryVertex();
 113   if (v) {
 114     dx -= v->GetX();
 115     dy -= v->GetY();
 116     dz -= v->GetZ();
 117   } else if (fRefMC) {
 118     // if reference is an MC event, no primary vertex is supplied
 119     // but it is possible to retrieve it from header
 120     TArrayF fvertex(3);
 121     fRefMC->GenEventHeader()->PrimaryVertex(fvertex);
 122     dx -= fvertex[0];
 123     dy -= fvertex[1];
 124     dz -= fvertex[2];
 125   }
 126   out.SetDr(TMath::Sqrt(dx*dx + dy*dy));
 127   out.SetDz(dz);
 128
 129   // compute PID probabilities by combining
 130   // the PID weights in the source with priors
 131   // and eventually using the PIDDefESD
 132   // (the AliRsnDaughter objec knows how to manage the latter)
 133   out.CombineWithPriors(fPrior, &fPIDDefESD);
 134
 135   // dynamic reference to true nature of referenced event
 136   // to get kink index
 137   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(fRef);
 138   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fRef);
 139
 140   if (esd) {
 141     AliESDtrack *esdTrack = esd->GetTrack(i);
 142     out.FindKinkIndex(esdTrack);
 143   } else if (aod) {
 144     out.FindKinkIndex(aod);
 145   }
 146 }
 147
 148 //_____________________________________________________________________________
 149 AliRsnDaughter AliRsnEvent::GetDaughter(Int_t i)
 150 {
 151 //
 152 // Return an AliRsnDaughter taken from this event,
 153 // with all additional data members well set.
 154 //
 155
 156   AliRsnDaughter out;
 157   SetDaughter(out, i);
 158
 159   return out;
 160 }
 161
 162 //_____________________________________________________________________________
 163 Int_t AliRsnEvent::GetMultiplicity()
 164 {
 165 //
 166 // Returns event multiplicity
 167 //
 168   AliDebug(AliLog::kDebug+2,"<-");
 169   if (!fRef) return 0;
 170   AliDebug(AliLog::kDebug+2,"->");
 171   return fRef->GetNumberOfTracks();
 172 }
 173
 174 //_____________________________________________________________________________
 175 Double_t AliRsnEvent::GetVz()
 176 {
 177 //
 178 // Return Z coord of primary vertex
 179 //
 180   AliDebug(AliLog::kDebug+2,"<-");
 181   return fRef->GetPrimaryVertex()->GetZ();
 182   AliDebug(AliLog::kDebug+2,"->");
 183 }
 184
 185 //_____________________________________________________________________________
 186 AliRsnDaughter AliRsnEvent::GetLeadingParticle
 187 (Double_t ptMin, AliPID::EParticleType type)
 188 {
 189 //
 190 // Searches the collection of all particles with given PID type and charge,
 191 // and returns the one with largest momentum, provided that it is greater than 1st argument.
 192 // If one specifies AliRsnPID::kUnknown as type or AliRsnDaughter::kNoPID as method,
 193 // the check is done over all particles irrespectively of their PID.
 194 // If the sign argument is '+' or '-', the check is done over the particles of that charge,
 195 // otherwise it is done irrespectively of the charge.
 196 //
 197
 198   Int_t i, nTracks = fRef->GetNumberOfTracks();
 199   AliRsnDaughter output;
 200
 201   for (i = 0; i < nTracks; i++) {
 202     AliRsnDaughter track = GetDaughter(i);
 203     if (!AcceptTrackPID(&track, type)) continue;
 204     if (track.Pt() < ptMin) continue;
 205     if (!output.IsOK() || track.Pt() > output.Pt()) {
 206       output = track;
 207       output.SetGood();
 208     }
 209   }
 210
 211   return output;
 212 }
 213
 214 //_________________________________________________________________________________________________
 215 Double_t AliRsnEvent::GetAverageMomentum(Int_t &count, AliPID::EParticleType type)
 216 {
 217 //
 218 // Loops on the list of tracks and computes average total momentum.
 219 //
 220
 221   Int_t i, nTracks = fRef->GetNumberOfTracks();
 222   Double_t pmean = 0.0;
 223
 224   for (i = 0, count = 0; i < nTracks; i++) {
 225     AliRsnDaughter track = GetDaughter(i);
 226     if (!AcceptTrackPID(&track, type)) continue;
 227     pmean += track.P();
 228     count++;
 229   }
 230
 231   if (count > 0) pmean /= (Double_t)count;
 232   else pmean = 0.0;
 233
 234   return pmean;
 235 }
 236
 237 //_____________________________________________________________________________
 238 Bool_t AliRsnEvent::GetAngleDistr
 239 (Double_t &angleMean, Double_t &angleRMS, AliRsnDaughter leading)
 240 {
 241 //
 242 // Takes the leading particle and computes the mean and RMS
 243 // of the distribution of directions of all other tracks
 244 // with respect to the direction of leading particle.
 245 //
 246
 247   if (!leading.IsOK()) return kFALSE;
 248
 249   Int_t i, count, nTracks = fRef->GetNumberOfTracks();
 250   Double_t angle, angle2Mean = 0.0;
 251
 252   angleMean = angle2Mean = 0.0;
 253
 254   for (i = 0, count = 0; i < nTracks; i++) {
 255     AliRsnDaughter trk = GetDaughter(i);
 256     if (trk.GetID() == leading.GetID()) continue;
 257
 258     angle = leading.AngleTo(trk);
 259
 260     angleMean += angle;
 261     angle2Mean += angle * angle;
 262     count++;
 263   }
 264
 265   if (!count) return kFALSE;
 266
 267   angleMean /= (Double_t)count;
 268   angle2Mean /= (Double_t)count;
 269   angleRMS = TMath::Sqrt(angle2Mean - angleMean*angleMean);
 270
 271   return kTRUE;
 272 }
 273
 274 //_____________________________________________________________________________
 275 void AliRsnEvent::SetPriorProbability(Double_t *const out)
 276 {
 277 //
 278 // Set all prior probabilities at once, using an assayr of values.
 279 //
 280
 281   Int_t i;
 282
 283   for (i = 0; i < AliPID::kSPECIES; i++)
 284   {
 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   {
 302     AliInfo(Form("Prior probability for %10s = %3.5f", AliPID::ParticleName((AliPID::EParticleType)i), fPrior[i]));
 303   }
 304 }
 305
 306 //_____________________________________________________________________________
 307 void AliRsnEvent::GetPriorProbability(Double_t *out) const
 308 {
 309 //
 310 // Stores in the passed argument all the values.
 311 //
 312
 313   Int_t i;
 314
 315   for (i = 0; i < AliPID::kSPECIES; i++)
 316   {
 317     out[i] = fPrior[i];
 318   }
 319
 320 }
 321
 322 //_____________________________________________________________________________
 323 Bool_t AliRsnEvent::AcceptTrackPID(AliRsnDaughter * const d, AliPID::EParticleType type)
 324 {
 325 //
 326 // [PRIVATE]
 327 // Checks if the track PID (according to method in use) corresponds
 328 // to the required identification species.
 329 // If the second argument is "kUnknown", answer of this method is always YES.
 330 //
 331
 332   if (type == AliPID::kUnknown) return kTRUE;
 333
 334   return (d->AssignedPID() == type);
 335 }