PWG2/RESONANCES/AliRsnEvent.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //
  17 // *** Class AliRsnEvent ***
  18 //
  19 // A container for a collection of AliRsnDaughter objects from an event.
  20 // Contains also the primary vertex, useful for some cuts.
  21 // In order to retrieve easily the tracks which have been identified
  22 // as a specific type and charge, there is an array of indexes which
  23 // allows to avoid to loop on all tracks and have only the neede ones.
  24 //
  25 // authors: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
  26 //          M. Vala (email: martin.vala@cern.ch)
  27 //
  28
  29 #include <Riostream.h>
  30 #include <TH1.h>
  31
  32 #include "AliLog.h"
  33
  34 #include "AliRsnDaughter.h"
  35 #include "AliRsnEvent.h"
  36 #include "AliRsnMCInfo.h"
  37
  38 ClassImp(AliRsnEvent)
  39
  40 //_____________________________________________________________________________
  41 AliRsnEvent::AliRsnEvent() :
  42     TNamed("rsnEvent", ""),
  43     fPVx(0.0),
  44     fPVy(0.0),
  45     fPVz(0.0),
  46     fPhiMean(0.0),
  47     fMult(0),
  48     fTracks(0x0),
  49     fNoPID(0x0),
  50     fPerfectPID(0x0),
  51     fRealisticPID(0x0),
  52     fSelPIDType(AliRsnPID::kUnknown),
  53     fSelCharge('0'),
  54     fSelPIDMethod(AliRsnDaughter::kRealistic),
  55     fSelCuts(0x0)
  56 {
  57 //
  58 // Default constructor
  59 // (implemented but not recommended for direct use)
  60 //
  61 }
  62
  63 //_____________________________________________________________________________
  64 AliRsnEvent::AliRsnEvent(const AliRsnEvent &event) :
  65     TNamed(event),
  66     fPVx(event.fPVx),
  67     fPVy(event.fPVy),
  68     fPVz(event.fPVz),
  69     fPhiMean(event.fPhiMean),
  70     fMult(event.fMult),
  71     fTracks(0x0),
  72     fNoPID(0x0),
  73     fPerfectPID(0x0),
  74     fRealisticPID(0x0),
  75     fSelPIDType(AliRsnPID::kUnknown),
  76     fSelCharge('0'),
  77     fSelPIDMethod(AliRsnDaughter::kRealistic),
  78     fSelCuts(0x0)
  79 {
  80 //
  81 // Copy constructor.
  82 // Copies all the tracks from the argument's collection
  83 // to this' one, and then recreates the PID index arrays,
  84 // trusting on the PID informations in the copied tracks.
  85 //
  86
  87   // during track copy, counts how many faults happen
  88   Int_t errors = Fill(event.fTracks);
  89   if (errors) AliWarning(Form("%d errors occurred in copy", errors));
  90
  91   // fill PID index arrays
  92   // FillPIDArrays();
  93
  94   if (event.fNoPID) fNoPID = new AliRsnPIDIndex(* (event.fNoPID));
  95   if (event.fPerfectPID) fPerfectPID = new AliRsnPIDIndex(* (event.fPerfectPID));
  96   if (event.fRealisticPID) fRealisticPID = new AliRsnPIDIndex(* (event.fRealisticPID));
  97 }
  98
  99 //_____________________________________________________________________________
 100 AliRsnEvent& AliRsnEvent::operator= (const AliRsnEvent &event)
 101 {
 102 //
 103 // Works in the same way as the copy constructor.
 104 //
 105   // copy name and title
 106   SetName(event.GetName());
 107   SetTitle(event.GetTitle());
 108
 109   // copy primary vertex and initialize track counter to 0
 110   fPVx = event.fPVx;
 111   fPVy = event.fPVy;
 112   fPVz = event.fPVz;
 113
 114   // other data
 115   fPhiMean = event.fPhiMean;
 116   fMult = event.fMult;
 117
 118   // add tracks from array of argument
 119   Int_t errors = Fill(event.fTracks);
 120   if (errors) AliWarning(Form("%d errors occurred in copy", errors));
 121
 122   // fill PID arrays
 123   // FillPIDArrays();
 124   if (event.fNoPID)
 125   {
 126     if (!fNoPID) fNoPID = new AliRsnPIDIndex(* (event.fNoPID));
 127     else (*fNoPID) = * (event.fNoPID);
 128   }
 129   if (event.fPerfectPID)
 130   {
 131     if (!fPerfectPID) fPerfectPID = new AliRsnPIDIndex(* (event.fPerfectPID));
 132     else (*fPerfectPID) = * (event.fPerfectPID);
 133   }
 134   if (event.fRealisticPID)
 135   {
 136     if (!fRealisticPID) fRealisticPID = new AliRsnPIDIndex(* (event.fRealisticPID));
 137     else (*fRealisticPID) = * (event.fRealisticPID);
 138   }
 139
 140   // return this object
 141   return (*this);
 142 }
 143
 144 //_____________________________________________________________________________
 145 AliRsnEvent::~AliRsnEvent()
 146 {
 147 //
 148 // Destructor.
 149 // Deletes the TClonesArray, after clearing its content.
 150 // Other memory-allocating arrays are cleared by their
 151 // destructor, which is automatically called from here.
 152 //
 153
 154   Clear();
 155   if (fTracks) delete fTracks;
 156 }
 157
 158 //_____________________________________________________________________________
 159 void AliRsnEvent::Init()
 160 {
 161 //
 162 // Initialize TClonesArray data-member.
 163 //
 164
 165   fTracks = new TClonesArray("AliRsnDaughter", 1);
 166   //fTracks->BypassStreamer (kFALSE);
 167 }
 168
 169 //_____________________________________________________________________________
 170 void AliRsnEvent::Clear(Option_t* /*option*/)
 171 {
 172 //
 173 // Empties the collections (does not delete the objects).
 174 // The track collection is emptied only at the end.
 175 // Since some objects could be uninitialized, some
 176 // "if" statement are used.
 177 //
 178
 179   if (fTracks) fTracks->Delete();
 180   delete fNoPID;
 181   fNoPID = 0x0;
 182   delete fPerfectPID;
 183   fPerfectPID = 0x0;
 184   delete fRealisticPID;
 185   fRealisticPID = 0x0;
 186 }
 187
 188 //_____________________________________________________________________________
 189 AliRsnDaughter* AliRsnEvent::AddTrack(AliRsnDaughter track)
 190 {
 191 //
 192 // Stores a new track into the array and returns
 193 // a reference pointer to it (which is NULL in case of errors).
 194 //
 195
 196   Int_t nextIndex = fTracks->GetEntriesFast();
 197   TClonesArray &tracks = (*fTracks);
 198   AliRsnDaughter *copy = new(tracks[nextIndex]) AliRsnDaughter(track);
 199   return copy;
 200 }
 201
 202 //_____________________________________________________________________________
 203 AliRsnDaughter* AliRsnEvent::GetTrack(Int_t index)
 204 {
 205 //
 206 // Returns one track in the collection
 207 // given the absolute index in the global TClonesArray
 208 //
 209   return (AliRsnDaughter*) fTracks->UncheckedAt(index);
 210 }
 211
 212 //_____________________________________________________________________________
 213 TArrayI* AliRsnEvent::GetCharged(Char_t sign)
 214 {
 215 //
 216 // Returns an array with the indexes of all tracks with a given charge
 217 // (arg can be '+' or '-'), irrespective of its PID.
 218 // When the argument is wrong, a NULL pointer is returned.
 219 //
 220   if (fNoPID) return fNoPID->GetTracksArray(sign, AliRsnPID::kUnknown);
 221   return 0x0;
 222 }
 223
 224 //_____________________________________________________________________________
 225 TArrayI * AliRsnEvent::GetTracksArray
 226 (AliRsnDaughter::EPIDMethod pidtype, Char_t sign, AliRsnPID::EType type)
 227 {
 228 //
 229 // Returns an array of indexes of all tracks in this event
 230 // which match the charge sign and PID type in the arguments,
 231 // according to one of the allowed PID methods (perfect or realistic).
 232 // It retrieves this array from the AliRsnPIDIndex data members.
 233 // If the arguments are wrong a NULL pointer is returned.
 234 //
 235
 236   switch (pidtype)
 237   {
 238     case AliRsnDaughter::kRealistic:
 239       if (fRealisticPID)
 240       {
 241         return fRealisticPID->GetTracksArray(sign, type);
 242       }
 243       break;
 244     case AliRsnDaughter::kPerfect:
 245       if (fPerfectPID)
 246       {
 247         return fPerfectPID->GetTracksArray(sign, type);
 248       }
 249       break;
 250     case AliRsnDaughter::kNoPID:
 251       if (fNoPID)
 252       {
 253         return fNoPID->GetTracksArray(sign, AliRsnPID::kUnknown);
 254       }
 255       break;
 256     default:
 257       AliError("Handled PID methods here are only fNoPID,kPerfect and kRealistic. Nothing done.");
 258       return 0x0;
 259   }
 260   return 0x0;
 261 }
 262
 263 //_____________________________________________________________________________
 264 void AliRsnEvent::FillPIDArrays(Int_t arraySizeInit)
 265 {
 266 //
 267 // Initializes and fills the AliRsnPIDIndex objects containing
 268 // arrays of indexes for each possible charge and PID type.
 269 // This method is the unique way to do this, for safety reasons.
 270 //
 271
 272   if (fNoPID) delete fNoPID;
 273   if (fPerfectPID) delete fPerfectPID;
 274   if (fRealisticPID) delete fRealisticPID;
 275   fNoPID = new AliRsnPIDIndex(arraySizeInit);
 276   fPerfectPID = new AliRsnPIDIndex(arraySizeInit);
 277   fRealisticPID = new AliRsnPIDIndex(arraySizeInit);
 278
 279   // set the default type to Realistic
 280   AliRsnDaughter::SetPIDMethod(AliRsnDaughter::kRealistic);
 281
 282   // loop on tracks and create references
 283   Double_t prob;
 284   Int_t i, icharge, type;
 285   Short_t charge;
 286   AliRsnMCInfo *mcinfo = 0;
 287   AliRsnDaughter *track = 0;
 288   TObjArrayIter iter(fTracks);
 289   while ((track = (AliRsnDaughter*) iter.Next()))
 290   {
 291     charge = track->Charge();
 292     type = (Int_t) track->PIDType(prob);
 293     i = fTracks->IndexOf(track);
 294     mcinfo = track->GetMCInfo();
 295     if (charge > 0) icharge = 0;
 296     else if (charge < 0) icharge = 1;
 297     else
 298     {
 299       AliError("Found particle with ZERO charge!!!");
 300       continue;
 301     }
 302     // add to charged array
 303     fNoPID->AddIndex(i, icharge, (Int_t) AliRsnPID::kUnknown);
 304     // add to realistic PID array
 305     fRealisticPID->AddIndex(i, icharge, (Int_t) type);
 306     // add to perfect PID array (needs MCInfo present)
 307     if (mcinfo)
 308     {
 309       fPerfectPID->AddIndex(i, icharge, (Int_t) AliRsnPID::InternalType(mcinfo->PDG()));
 310     }
 311   }
 312
 313   // adjusts the size of arrays
 314   if (fNoPID) fNoPID->SetCorrectIndexSize();
 315   if (fPerfectPID) fPerfectPID->SetCorrectIndexSize();
 316   if (fRealisticPID) fRealisticPID->SetCorrectIndexSize();
 317 }
 318
 319 //_____________________________________________________________________________
 320 void AliRsnEvent::Print(Option_t *option) const
 321 {
 322 //
 323 // Lists the details of the event, and the ones of each
 324 // contained track.
 325 // The options are passed to AliRsnDaughter::Print().
 326 // Look at that method to understand option values.
 327 //
 328
 329   cout << "...Multiplicity     : " << fTracks->GetEntries() << endl;
 330   cout << "...Primary vertex   : " << fPVx << ' ' << fPVy << ' ' << fPVz << endl;
 331
 332   TObjArrayIter iter(fTracks);
 333   AliRsnDaughter *d = 0;
 334   while ((d = (AliRsnDaughter*) iter.Next()))
 335   {
 336     cout << "....Track #" << fTracks->IndexOf(d) << endl;
 337     d->Print(option);
 338   }
 339 }
 340
 341 //_____________________________________________________________________________
 342 void AliRsnEvent::MakeComputations()
 343 {
 344 //
 345 // Computes all required overall variables:
 346 // - multiplicity
 347 // - mean phi of tracks
 348 //
 349
 350   if (!fTracks)
 351   {
 352     fMult = 0;
 353     fPhiMean = 1000.0;
 354   }
 355   else
 356   {
 357     fMult = fTracks->GetEntries();
 358     if (fMult < 1) {
 359       fPhiMean = 1000.0;
 360     }
 361     else
 362     {
 363       fPhiMean = 0.0;
 364       AliRsnDaughter *d = 0;
 365       TObjArrayIter next(fTracks);
 366       while ( (d = (AliRsnDaughter*)next()) ) fPhiMean += d->Phi();
 367       fPhiMean /= (Double_t)fMult;
 368     }
 369   }
 370 }
 371
 372 //_____________________________________________________________________________
 373 void AliRsnEvent::CorrectByPrimaryVertex()
 374 {
 375 //
 376 // Corrects the X,Y,Z position of DCA vertex of all tracks
 377 // by the amount of stored primary vertex
 378 //
 379
 380     TObjArrayIter iter(fTracks);
 381     AliRsnDaughter *d = 0;
 382     while ((d = (AliRsnDaughter*) iter.Next())) {
 383         d->ShiftZero(fPVx, fPVy, fPVz);
 384     }
 385 }
 386
 387 /*
 388 //_____________________________________________________________________________
 389 Int_t AliRsnEvent::GetMultiplicity() const
 390 {
 391 //
 392 // Get number of all tracks
 393 //
 394
 395   if (!fTracks) return 0;
 396   return fTracks->GetEntries();
 397 }
 398 */
 399
 400 //_____________________________________________________________________________
 401 Int_t AliRsnEvent::GetNCharged(Char_t sign)
 402 {
 403 //
 404 // Get number of charged tracks
 405 //
 406
 407   Int_t icharge;
 408   icharge = ChargeIndex(sign);
 409   if (icharge < 0) return 0;
 410   TArrayI *charged = GetCharged(sign);
 411   if (!charged) return 0;
 412   return charged->GetSize();
 413 }
 414
 415 //_____________________________________________________________________________
 416 Int_t AliRsnEvent::Fill(TObjArray *array)
 417 {
 418 //
 419 // Fills the data-member TClonesArray of tracks with
 420 // the ones stored in the array passed as argument.
 421 // If this data-member is already present, it is cleared.
 422 // Returns the number of tracks which raised problems
 423 // while attempting to add them. Zero is the best.
 424 //
 425
 426   // clear the array if it is already instantiated,
 427   // create if otherwise
 428   if (fTracks) fTracks->Delete();
 429   else Init();
 430
 431   // copy argument entries into data-member
 432   Int_t errors = 0;
 433   AliRsnDaughter *track = 0;
 434   TObjArrayIter iter(array);
 435   while ((track = (AliRsnDaughter*) iter.Next()))
 436   {
 437     AliRsnDaughter *ref = AddTrack(*track);
 438     if (!ref)
 439     {
 440       AliWarning(Form("Problem occurred when copying track #%d from passed array", array->IndexOf(track)));
 441       errors++;
 442     }
 443   }
 444
 445   return errors;
 446 }
 447
 448 //_____________________________________________________________________________
 449 Int_t AliRsnEvent::ChargeIndex(Char_t sign) const
 450 //
 451 // Returns the array index corresponding to charge
 452 // 0 for positive, 1 for negative
 453 //
 454 {
 455   if (sign == '+') return 0;
 456   else if (sign == '-') return 1;
 457   else
 458   {
 459     AliError(Form("Character '%c' not recognized as charge sign", sign));
 460     return -1;
 461   }
 462 }
 463 //_____________________________________________________________________________
 464 inline void AliRsnEvent::SetSelection
 465 (AliRsnPID::EType pid, Char_t charge, AliRsnDaughter::EPIDMethod meth, AliRsnCutSet *cuts)
 466 {
 467 //
 468 // Set all selection parameters at once
 469 //
 470
 471   SetSelectionPIDType(pid);
 472   SetSelectionCharge(charge);
 473   SetSelectionPIDMethod(meth);
 474   SetSelectionTrackCuts(cuts);
 475 }
 476
 477 //_____________________________________________________________________________
 478 AliRsnDaughter* AliRsnEvent::GetLeadingParticle(Double_t ptMin)
 479 {
 480 //
 481 // Searches the collection of all particles with given PID type and charge,
 482 // and returns the one with largest momentum, provided that it is greater than 1st argument.
 483 // If one specifies AliRsnPID::kUnknown as type or AliRsnDaughter::kNoPID as method,
 484 // the check is done over all particles irrespectively of their PID.
 485 // If the sign argument is '+' or '-', the check is done over the particles of that charge,
 486 // otherwise it is done irrespectively of the charge.
 487 //
 488
 489   Int_t i;
 490   TArrayI *array = 0x0;
 491   AliRsnDaughter *track1 = 0x0, *track2 = 0x0, *leading = 0x0;
 492
 493   if (fSelCharge == '+' || fSelCharge == '-') {
 494     // if the charge '+' or '-' does simply the search
 495     array = GetTracksArray(fSelPIDMethod, fSelCharge, fSelPIDType);
 496     for (i = 0; i < array->GetSize(); i++) {
 497       track1 = (AliRsnDaughter *) fTracks->At(array->At(i));
 498       if (!track1) continue;
 499       if (track1->Pt() < ptMin) continue;
 500       if (!CutPass(track1)) continue;
 501       ptMin = track1->Pt();
 502       leading = track1;
 503     }
 504   }
 505   else {
 506     track1 = GetLeadingParticle(ptMin);
 507     track2 = GetLeadingParticle(ptMin);
 508     if (track1 && track2) {
 509       if (track1->Pt() > track2->Pt()) leading = track1;
 510       else leading = track2;
 511     }
 512     else if (track1) leading = track1;
 513     else if (track2) leading = track2;
 514     else leading = 0x0;
 515   }
 516
 517   return leading;
 518 }
 519
 520 //_____________________________________________________________________________
 521 Double_t AliRsnEvent::GetAverageMomentum(Int_t &count)
 522 {
 523 //
 524 // Loops on the list of tracks and computes average total momentum.
 525 //
 526
 527   Int_t i;
 528   Double_t pmean = 0.0;
 529   TArrayI *array = 0x0;
 530   AliRsnDaughter *d = 0x0;
 531
 532   if (fSelCharge == '+' || fSelCharge == '-') {
 533     // if the charge '+' or '-' does simply the search
 534     count = 0;
 535     array = GetTracksArray(fSelPIDMethod, fSelCharge, fSelPIDType);
 536     for (i = 0; i < array->GetSize(); i++) {
 537       d = (AliRsnDaughter *) fTracks->At(array->At(i));
 538       if (!d) continue;
 539       if (!CutPass(d)) continue;
 540       pmean += d->P();
 541       count++;
 542     }
 543     if (count > 0) pmean /= (Double_t)count;
 544     else pmean = 0.0;
 545   }
 546   else {
 547     Int_t countP, countM;
 548     Double_t pmeanP = GetAverageMomentum(countP);
 549     Double_t pmeanM = GetAverageMomentum(countM);
 550     if (countP && countM) {
 551       pmean = (pmeanP * (Double_t)countP + pmeanM * (Double_t)countM) / (countP + countM);
 552       count = countP + countM;
 553     }
 554     else if (countP) {
 555       pmean = pmeanP;
 556       count = countP;
 557     }
 558     else if (countM) {
 559       pmean = pmeanM;
 560       count = countM;
 561     }
 562     else {
 563       count = 0;
 564       pmean = 0.0;
 565     }
 566   }
 567
 568   return pmean;
 569 }
 570
 571 //_____________________________________________________________________________
 572 Bool_t AliRsnEvent::GetAngleDistrWRLeading
 573 (Double_t &angleMean, Double_t &angleRMS, Double_t ptMin)
 574 {
 575 //
 576 // Takes the leading particle and computes the mean and RMS
 577 // of the distribution of directions of all other tracks
 578 // with respect to the direction of leading particle.
 579 //
 580
 581   AliRsnDaughter *leading = GetLeadingParticle(ptMin);
 582   if (!leading) return kFALSE;
 583
 584   Int_t count = 0;
 585   Double_t angle, angle2Mean;
 586   AliRsnDaughter *trk = 0x0;
 587   TObjArrayIter next(fTracks);
 588
 589   angleMean = angle2Mean = 0.0;
 590
 591   while ( (trk = (AliRsnDaughter*)next()) )
 592   {
 593     if (trk == leading) continue;
 594
 595     angle = leading->AngleTo(trk);
 596
 597     angleMean += angle;
 598     angle2Mean += angle * angle;
 599     count++;
 600   }
 601
 602   if (!count) return kFALSE;
 603
 604   angleMean /= (Double_t)count;
 605   angle2Mean /= (Double_t)count;
 606   angleRMS = TMath::Sqrt(angle2Mean - angleMean*angleMean);
 607
 608   return kTRUE;
 609 }