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 "AliMCEvent.h"
  19 #include "AliStack.h"
  20 #include "AliGenEventHeader.h"
  21 #include "AliAODEvent.h"
  22 #include "AliRsnCutPID.h"
  23 #include "AliESDtrackCuts.h"
  24
  25 #include "AliRsnEvent.h"
  26
  27 ClassImp(AliRsnEvent)
  28
  29 //_____________________________________________________________________________
  30 AliRsnEvent::AliRsnEvent(AliVEvent *ref, AliVEvent *refMC) :
  31   fRef(ref),
  32   fRefMC(refMC),
  33   fLeading(-1)
  34 {
  35 //
  36 // Default constructor.
  37 //
  38 }
  39
  40 //_____________________________________________________________________________
  41 AliRsnEvent::AliRsnEvent(const AliRsnEvent &event) :
  42   TObject(event),
  43   fRef(event.fRef),
  44   fRefMC(event.fRefMC),
  45   fLeading(event.fLeading)
  46 {
  47 //
  48 // Copy constructor.
  49 //
  50 }
  51
  52 //_____________________________________________________________________________
  53 AliRsnEvent& AliRsnEvent::operator= (const AliRsnEvent & event)
  54 {
  55 //
  56 // Works in the same way as the copy constructor.
  57 //
  58
  59   (TObject)(*this) = (TObject)event;
  60   fRef             = event.fRef;
  61   fRefMC           = event.fRefMC;
  62   fLeading         = event.fLeading;
  63
  64   return (*this);
  65 }
  66
  67 //_____________________________________________________________________________
  68 Bool_t AliRsnEvent::SetDaughter(AliRsnDaughter &out, Int_t i, AliRsnDaughter::ERefType type)
  69 {
  70 //
  71 // Using the second and third arguments, retrieves the i-th object in the
  72 // appropriate sample (tracks or V0s) and sets the first reference object
  73 // in order to point to that.
  74 // If a MonteCarlo information is provided, sets the useful informations from there,
  75 // and in case of a V0, sets the 'label' data member only when the two daughters
  76 // of the V0 point to the same mother.
  77 // Returns kFALSE whenever the operation fails (out of range, NULL references).
  78 //
  79
  80   Int_t result = 0;
  81
  82   if (IsESD() && type == AliRsnDaughter::kTrack) result = SetDaughterESDtrack(out, i);
  83   if (IsAOD() && type == AliRsnDaughter::kTrack) result = SetDaughterAODtrack(out, i);
  84   if (IsESD() && type == AliRsnDaughter::kV0   ) result = SetDaughterESDv0   (out, i);
  85   if (IsAOD() && type == AliRsnDaughter::kV0   ) result = SetDaughterAODv0   (out, i);
  86
  87   return (result == 0);
  88 }
  89
  90 //_____________________________________________________________________________
  91 Bool_t AliRsnEvent::SetDaughterMC(AliRsnDaughter &out, Int_t label)
  92 {
  93 //
  94 // Using the second argument, retrieves the i-th object in the
  95 // MC sample (if present) and assigns the track using only that,
  96 // so that it is considered both as main reference and MC reference.
  97 // (used for MC-only analysis).
  98 //
  99
 100   if (!fRefMC)
 101   {
 102     out.SetBad();
 103     return kFALSE;
 104   }
 105
 106   // try to convert into both types
 107   Int_t        imum;
 108   AliMCEvent  *esd = GetRefMCESD();
 109   AliAODEvent *aod = GetRefMCAOD();
 110
 111   // ESD
 112   if (esd)
 113   {
 114     // if the MC track exists, assign it
 115     AliMCParticle *track = (AliMCParticle*)fRef->GetTrack(label);
 116     if (!track)
 117     {
 118       out.SetBad();
 119       return kFALSE;
 120     }
 121     out.SetRef(track);
 122     out.SetRefMC(track);
 123     out.SetLabel(label);
 124     out.SetGood();
 125
 126     // search for its mother in stack
 127     imum = track->GetMother();
 128     if (imum >= 0 && imum < esd->Stack()->GetNtrack())
 129     {
 130       TParticle *mum = esd->Stack()->Particle(imum);
 131       if (mum) out.SetMotherPDG(TMath::Abs(mum->GetPdgCode()));
 132     }
 133   }
 134
 135   // AOD
 136   if (aod)
 137   {
 138     // checks that array of MC particles exists
 139     TClonesArray *mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 140     if(!mcArray)
 141     {
 142       out.SetBad();
 143       return kFALSE;
 144     }
 145
 146     // in this case one has to loop over the sample to find the good one
 147     TObjArrayIter next(mcArray);
 148     AliAODMCParticle *part = 0x0;
 149     while ( (part = (AliAODMCParticle*)next()) )
 150     {
 151       if (TMath::Abs(part->GetLabel()) == label)
 152       {
 153         // if the MC track exists, assign it
 154         out.SetRef(part);
 155         out.SetRefMC(part);
 156         out.SetLabel(label);
 157         out.SetGood();
 158
 159         // search for the mother
 160         imum = part->GetMother();
 161         if (imum >= 0 && imum < mcArray->GetEntriesFast())
 162         {
 163           AliAODMCParticle *mum = (AliAODMCParticle*)mcArray->At(imum);
 164           if (mum) out.SetMotherPDG(TMath::Abs(mum->GetPdgCode()));
 165         }
 166         break;
 167       }
 168     }
 169     return kTRUE;
 170   }
 171
 172   return kFALSE;
 173 }
 174
 175 //_____________________________________________________________________________
 176 AliRsnDaughter AliRsnEvent::GetDaughter(Int_t i, AliRsnDaughter::ERefType type)
 177 {
 178 //
 179 // Returns a daughter set using same criteria as SetDaughter
 180 //
 181
 182   AliRsnDaughter d;
 183   SetDaughter(d, i, type);
 184   return d;
 185 }
 186
 187 //_____________________________________________________________________________
 188 AliRsnDaughter AliRsnEvent::GetDaughterMC(Int_t i)
 189 {
 190 //
 191 // Returns a daughter set using same criteria as SetDaughterMC
 192 //
 193
 194   AliRsnDaughter d;
 195   SetDaughterMC(d, i);
 196   return d;
 197 }
 198
 199 //_____________________________________________________________________________
 200 Int_t AliRsnEvent::GetMultiplicity(AliESDtrackCuts *cuts)
 201 {
 202 //
 203 // Returns event multiplicity as the number of tracks.
 204 // If the argument is not NULL, returns instead the
 205 // number of tracks passing the cuts hereby defined.
 206 //
 207
 208   if (!fRef) return 0;
 209
 210   AliESDEvent *esd = GetRefESD();
 211   if (cuts && esd) return cuts->CountAcceptedTracks(esd);
 212   else return fRef->GetNumberOfTracks();
 213 }
 214
 215 //_____________________________________________________________________________
 216 Double_t AliRsnEvent::GetVz()
 217 {
 218 //
 219 // Return Z coord of primary vertex
 220 //
 221
 222   return fRef->GetPrimaryVertex()->GetZ();
 223 }
 224
 225 //_____________________________________________________________________________
 226 Int_t AliRsnEvent::SelectLeadingParticle
 227 (Double_t ptMin, AliRsnCutPID *cutPID)
 228 {
 229 //
 230 // Searches the collection of all particles with given PID type and charge,
 231 // and returns the one with largest momentum, provided that it is greater than 1st argument.
 232 // If one specifies AliRsnPID::kUnknown as type or AliRsnDaughter::kNoPID as method,
 233 // the check is done over all particles irrespectively of their PID.
 234 // If the sign argument is '+' or '-', the check is done over the particles of that charge,
 235 // otherwise it is done irrespectively of the charge.
 236 //
 237
 238   Int_t i, nTracks = fRef->GetNumberOfTracks();
 239   fLeading = -1;
 240   AliRsnDaughter leading;
 241   leading.SetBad();
 242
 243   for (i = 0; i < nTracks; i++)
 244   {
 245     AliRsnDaughter track = GetDaughter(i);
 246     if (cutPID) if (!cutPID->IsSelected(&track)) continue;
 247     const AliVParticle *ref = track.GetRef();
 248     if (ref->Pt() < ptMin) continue;
 249     //double pt = track.P().Perp();
 250     //Printf("track %d %g", i, pt);
 251     if (!leading.IsOK() || ref->Pt() > ptMin)
 252     {
 253       fLeading = i;
 254       //leading = track;
 255       ptMin = ref->Pt();
 256     }
 257   }
 258   return fLeading;
 259 }
 260
 261 //_________________________________________________________________________________________________
 262 Double_t AliRsnEvent::GetAverageMomentum(Int_t &count, AliRsnCutPID *cutPID)
 263 {
 264 //
 265 // Loops on the list of tracks and computes average total momentum.
 266 //
 267
 268   Int_t i, nTracks = fRef->GetNumberOfTracks();
 269   Double_t pmean = 0.0;
 270
 271   for (i = 0, count = 0; i < nTracks; i++) {
 272     AliRsnDaughter track = GetDaughter(i);
 273     if (cutPID) if (!cutPID->IsSelected(&track)) continue;
 274     pmean += track.P().Mag();
 275     count++;
 276   }
 277
 278   if (count > 0) pmean /= (Double_t)count;
 279   else pmean = 0.0;
 280
 281   return pmean;
 282 }
 283
 284 //_____________________________________________________________________________
 285 Bool_t AliRsnEvent::GetAngleDistr
 286 (Double_t &angleMean, Double_t &angleRMS, AliRsnDaughter leading)
 287 {
 288 //
 289 // Takes the leading particle and computes the mean and RMS
 290 // of the distribution of directions of all other tracks
 291 // with respect to the direction of leading particle.
 292 //
 293
 294   if (!leading.IsOK()) return kFALSE;
 295
 296   Int_t i, count, nTracks = fRef->GetNumberOfTracks();
 297   Double_t angle, angle2Mean = 0.0;
 298
 299   angleMean = angle2Mean = 0.0;
 300
 301   for (i = 0, count = 0; i < nTracks; i++) {
 302     AliRsnDaughter trk = GetDaughter(i);
 303     if (trk.GetID() == leading.GetID()) continue;
 304
 305     angle = leading.P().Angle(trk.P().Vect());
 306
 307     angleMean += angle;
 308     angle2Mean += angle * angle;
 309     count++;
 310   }
 311
 312   if (!count) return kFALSE;
 313
 314   angleMean /= (Double_t)count;
 315   angle2Mean /= (Double_t)count;
 316   angleRMS = TMath::Sqrt(angle2Mean - angleMean * angleMean);
 317
 318   return kTRUE;
 319 }
 320
 321 //_____________________________________________________________________________
 322 Int_t AliRsnEvent::SetDaughterESDtrack(AliRsnDaughter &out, Int_t i)
 323 {
 324 //
 325 // Setup the first argument to the track identified by the index.
 326 // When available, adds the MC information and references.
 327 // ---
 328 // Version #1: ESD tracks
 329 //
 330
 331   // check 1: index in good range
 332   if (i < 0 || i >= fRef->GetNumberOfTracks())
 333   {
 334     out.SetBad();
 335     return 1;
 336   }
 337
 338   // check 2: not NULL object
 339   AliVTrack *track = (AliVTrack*)fRef->GetTrack(i);
 340   if (!track)
 341   {
 342     out.SetBad();
 343     return 2;
 344   }
 345
 346   // assign references of reconstructed track
 347   Int_t label = TMath::Abs(track->GetLabel());
 348   out.SetRef(track);
 349   out.SetLabel(label);
 350   out.SetGood();
 351
 352   // search for MC track, if available
 353   AliMCEvent *mc = GetRefMCESD();
 354   if (!mc) return 0;
 355
 356   // loop on particles using the track label as reference
 357   // and then assign also the mother type, if found
 358   AliStack *stack = mc->Stack();
 359   if (label >= 0 && label < stack->GetNtrack())
 360   {
 361     TParticle *part = stack->Particle(label);
 362     if (part)
 363     {
 364       Int_t imum = part->GetFirstMother();
 365       if (imum >= 0 && imum <= stack->GetNtrack())
 366       {
 367         TParticle *mum = stack->Particle(imum);
 368         if (mum) out.SetMotherPDG(TMath::Abs(mum->GetPdgCode()));
 369       }
 370     }
 371     else
 372     {
 373       return 3;
 374     }
 375   }
 376
 377   return 0;
 378 }
 379
 380 //_____________________________________________________________________________
 381 Int_t AliRsnEvent::SetDaughterAODtrack(AliRsnDaughter &out, Int_t i)
 382 {
 383 //
 384 // Setup the first argument to the track identified by the index.
 385 // When available, adds the MC information and references.
 386 // ---
 387 // Version #2: AOD tracks
 388 //
 389
 390   // check 1: index in good range
 391   if (i < 0 || i >= fRef->GetNumberOfTracks())
 392   {
 393     out.SetBad();
 394     return 1;
 395   }
 396
 397   // check 2: not NULL object
 398   AliVTrack *track = (AliVTrack*)fRef->GetTrack(i);
 399   if (!track)
 400   {
 401     out.SetBad();
 402     return 2;
 403   }
 404
 405   // assign references of reconstructed track
 406   Int_t label = TMath::Abs(track->GetLabel());
 407   out.SetRef(track);
 408   out.SetLabel(label);
 409   out.SetGood();
 410
 411   // search for MC track, if available
 412   AliAODEvent *mc = GetRefMCAOD();
 413   if (!mc) return 0;
 414
 415   // loop on particles using the track label as reference
 416   // and then assign also the mother type, if found
 417   TClonesArray *mcArray = (TClonesArray*)mc->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 418   if(!mcArray) return 0;
 419   TObjArrayIter next(mcArray);
 420   AliAODMCParticle *part = 0x0;
 421   while ( (part = (AliAODMCParticle*)next()) )
 422   {
 423     if (TMath::Abs(part->GetLabel()) == label)
 424     {
 425       out.SetRefMC(part);
 426       Int_t imum = part->GetMother();
 427       if (imum >= 0 && imum <= mcArray->GetEntriesFast())
 428       {
 429         AliAODMCParticle *mum = (AliAODMCParticle*)mcArray->At(imum);
 430         if (mum) out.SetMotherPDG(TMath::Abs(mum->GetPdgCode()));
 431       }
 432       break;
 433     }
 434   }
 435
 436   return 0;
 437 }
 438
 439 //_____________________________________________________________________________
 440 Int_t AliRsnEvent::SetDaughterESDv0(AliRsnDaughter &out, Int_t i)
 441 {
 442 //
 443 // Setup the first argument to the track identified by the index.
 444 // When available, adds the MC information and references.
 445 // ---
 446 // Version #3: ESD v0
 447 //
 448
 449   // check 1: index in good range
 450   if (i > fRef->GetNumberOfV0s())
 451   {
 452     out.SetBad();
 453     return 1;
 454   }
 455
 456   // check 2: not NULL object
 457   AliESDEvent *ev = GetRefESD();
 458   AliESDv0    *v0 = ev->GetV0(i);
 459   if (!v0)
 460   {
 461     out.SetBad();
 462     return 2;
 463   }
 464
 465   // assign references of reconstructed track
 466   out.SetRef(v0);
 467   out.SetGood();
 468
 469   // this time, assigning label is not trivial,
 470   // it is done only if MC is present and both
 471   // daughters come from a true particle
 472   AliMCEvent  *mc = GetRefMCESD();
 473   AliESDtrack *tp = ev->GetTrack(v0->GetPindex());
 474   AliESDtrack *tn = ev->GetTrack(v0->GetNindex());
 475   if (mc && tp && tn)
 476   {
 477     Int_t        lp = TMath::Abs(tp->GetLabel());
 478     Int_t        ln = TMath::Abs(tn->GetLabel());
 479     TParticle   *pp = mc->Stack()->Particle(lp);
 480     TParticle   *pn = mc->Stack()->Particle(ln);
 481     // if their first mothers are the same, the V0 is true
 482     // otherwise no label can be assigned
 483     if (pp->GetFirstMother() == pn->GetFirstMother()) out.SetLabel(pp->GetFirstMother());
 484   }
 485
 486   return 0;
 487 }
 488
 489 //_____________________________________________________________________________
 490 Int_t AliRsnEvent::SetDaughterAODv0(AliRsnDaughter &out, Int_t i)
 491 {
 492 //
 493 // Setup the first argument to the track identified by the index.
 494 // When available, adds the MC information and references.
 495 // ---
 496 // Version #4: AOD v0
 497 //
 498
 499   // check 1: index in good range
 500   if (i > fRef->GetNumberOfV0s())
 501   {
 502     out.SetBad();
 503     return 1;
 504   }
 505
 506   // check 2: not NULL object
 507   AliAODEvent *ev = GetRefAOD();
 508   AliAODv0    *v0 = ev->GetV0(i);
 509   if (!v0)
 510   {
 511     out.SetBad();
 512     return 2;
 513   }
 514
 515   // assign references of reconstructed track
 516   out.SetRef(v0);
 517   out.SetGood();
 518   out.SetLabel(-1);
 519
 520   // this time, assigning label is not trivial,
 521   // it is done only if MC is present and both
 522   // daughters come from a true particle
 523   TClonesArray *mcArray = (TClonesArray*)ev->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 524   AliAODTrack  *tp = ev->GetTrack(v0->GetPosID());
 525   AliAODTrack  *tn = ev->GetTrack(v0->GetNegID());
 526   if (mcArray && tp && tn)
 527   {
 528     Int_t        lp = TMath::Abs(tp->GetLabel());
 529     Int_t        ln = TMath::Abs(tn->GetLabel());
 530     // loop on array to find MC daughters
 531     AliAODMCParticle *pp = 0x0, *pn = 0x0;
 532     TObjArrayIter next(mcArray);
 533     AliAODMCParticle *part = 0x0;
 534     while ( (part = (AliAODMCParticle*)next()) )
 535     {
 536       if (TMath::Abs(part->GetLabel()) == lp) pp = (AliAODMCParticle*)mcArray->IndexOf(part);
 537       if (TMath::Abs(part->GetLabel()) == ln) pn = (AliAODMCParticle*)mcArray->IndexOf(part);
 538     }
 539     // assign a MC reference and a label only to true V0s
 540     if (pp->GetMother() == pn->GetMother())
 541     {
 542       AliAODMCParticle *mcv0 = (AliAODMCParticle*)mcArray->At(pp->GetMother());
 543       out.SetRefMC(mcv0);
 544       out.SetLabel(pp->GetMother());
 545       Int_t imum = mcv0->GetMother();
 546       if (imum >= 0 && imum <= mcArray->GetEntriesFast())
 547       {
 548         AliAODMCParticle *mum = (AliAODMCParticle*)mcArray->At(imum);
 549         if (mum) out.SetMotherPDG(TMath::Abs(mum->GetPdgCode()));
 550       }
 551     }
 552   }
 553
 554   return 0;
 555 }