PWG2/RESONANCES/AliRsnAnalysis.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 AliRsnAnalysis
  18 //            Reconstruction and analysis of a binary resonance
  19 // ........................................
  20 // ........................................
  21 // ........................................
  22 // ........................................
  23 // ........................................
  24 //
  25 // author: A. Pulvirenti             (email: alberto.pulvirenti@ct.infn.it)
  26 //-------------------------------------------------------------------------
  27
  28 #include <Riostream.h>
  29
  30 #include <TH1.h>
  31 #include <TTree.h>
  32 #include <TDatabasePDG.h>
  33
  34 #include "AliRsnDaughter.h"
  35 #include "AliRsnDaughterCut.h"
  36 #include "AliRsnDaughterCutPair.h"
  37 #include "AliRsnEvent.h"
  38 #include "AliRsnAnalysis.h"
  39
  40 ClassImp(AliRsnAnalysis)
  41
  42 //--------------------------------------------------------------------------------------------------------
  43 AliRsnAnalysis::AliRsnAnalysis() :
  44   TObject(),
  45   fRejectFakes(kFALSE),
  46   fNBins(0),
  47   fHistoMin(0.0),
  48   fHistoMax(0.0),
  49   fTrueMotherPDG(0),
  50   fMixPairDefs(0x0),
  51   fMixHistograms(0x0),
  52   fPairDefs(0x0),
  53   fHistograms(0x0),
  54   fPairCuts(0x0),
  55   fEventsTree(0x0)
  56 {
  57 //
  58 // Constructor
  59 // Initializes all pointers and collections to NULL.
  60 //
  61         Int_t i;
  62         for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i] = 0;
  63 }
  64 //--------------------------------------------------------------------------------------------------------
  65 AliRsnAnalysis::AliRsnAnalysis(const AliRsnAnalysis &copy) :
  66   TObject(copy),
  67   fRejectFakes(copy.fRejectFakes),
  68   fNBins(copy.fNBins),
  69   fHistoMin(copy.fHistoMin),
  70   fHistoMax(copy.fHistoMax),
  71   fTrueMotherPDG(copy.fTrueMotherPDG),
  72   fMixPairDefs(0x0),
  73   fMixHistograms(0x0),
  74   fPairDefs(0x0),
  75   fHistograms(0x0),
  76   fPairCuts(0x0),
  77   fEventsTree(0x0)
  78 {
  79 //
  80 // Copy constructor
  81 // Initializes all pointers and collections to NULL anyway,
  82 // but copies some settings from the argument.
  83 //
  84         Int_t i;
  85         for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i] = 0;
  86 }
  87 //--------------------------------------------------------------------------------------------------------
  88 void AliRsnAnalysis::AddCutPair(AliRsnDaughterCutPair *cut)
  89 {
  90 //
  91 // Add a cut on pairs.
  92 // This cut is global for all pairs.
  93 //
  94         if (!fPairCuts) fPairCuts = new TObjArray(0);
  95         fPairCuts->AddLast(cut);
  96 }
  97 //--------------------------------------------------------------------------------------------------------
  98 void AliRsnAnalysis::AddCutSingle(AliPID::EParticleType type, AliRsnDaughterCut *cut)
  99 {
 100 //
 101 // Add a cut on single particles.
 102 // This cut must be specified for each particle type.
 103 //
 104         if (type >= AliPID::kElectron && type <= AliPID::kProton) {
 105                 if (!fCuts[type]) fCuts[type] = new TObjArray(0);
 106                 fCuts[type]->AddLast(cut);
 107         }
 108 }
 109 //--------------------------------------------------------------------------------------------------------
 110 void AliRsnAnalysis::AddMixPairDef
 111 (AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2)
 112 {
 113 //
 114 // Adds a new pair definition to create a new histogram,
 115 // for the event mixing step.
 116 // If the pair defs array is NULL, it is initialized here.
 117 //
 118         if (!fMixPairDefs) fMixPairDefs = new TObjArray(0);
 119         fMixPairDefs->AddLast( new AliPairDef(p1, sign1, p2, sign2, 0, kFALSE) );
 120 }
 121 //--------------------------------------------------------------------------------------------------------
 122 void AliRsnAnalysis::AddPairDef
 123 (AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2, Bool_t onlyTrue)
 124 {
 125 //
 126 // Adds a new pair definition to create a new histogram,
 127 // for the signal evaluation (same event) step.
 128 // If the pair defs array is NULL, it is initialized here.
 129 //
 130         if (!fPairDefs) fPairDefs = new TObjArray(0);
 131         fPairDefs->AddLast( new AliPairDef(p1, sign1, p2, sign2, fTrueMotherPDG, onlyTrue) );
 132 }
 133 //--------------------------------------------------------------------------------------------------------
 134 void AliRsnAnalysis::Clear(Option_t* /* option */)
 135 {
 136 //
 137 // Clear heap
 138 //
 139         fHistograms->Clear("C");
 140         fPairDefs->Clear("C");
 141
 142         fMixHistograms->Clear("C");
 143         fMixPairDefs->Clear("C");
 144
 145         Int_t i;
 146         for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i]->Clear("C");
 147         fPairCuts->Clear("C");
 148 }
 149 //--------------------------------------------------------------------------------------------------------
 150 Stat_t AliRsnAnalysis::EventMix(Int_t nmix, Int_t multDiffMax, Double_t vzDiffMax, Bool_t compareTotMult)
 151 {
 152 //
 153 // Performs event mixing.
 154 // It takes the array of fMixPairDefs and stores results in fMixHistograms.
 155 // First parameter defines how many events must be mixed with each one.
 156 // Events to be mixed are chosen using the other parameters:
 157 //
 158 // - multDiffMax defines the maximum allowed difference in particle multiplicity,
 159 //   a) when last argument is kFALSE (default) the multiplicity comparison is made with respect
 160 //      to the particles of 'type 2' in the pair
 161 //   b) when last argument is kTRUE, the multiplicity comparison is made with respect of total
 162 //      particle multiplicity in the events
 163 //
 164 // - vzDiffMax defines maximum allowed difference in Z coordinate of prim. vertex.
 165 //
 166 // If one wants to exchange the particle types, one has to add both combinations of particles
 167 // as different pair defs.
 168 //
 169 // EXAMPLE:
 170 // analysis->AddMixPairDef(AliRsnEvent::kPion, '+', AliRsnEvent::kKaon, '-');
 171 // analysis->AddMixPairDef(AliRsnEvent::kKaon, '-', AliRsnEvent::kPion, '+');
 172 //
 173         // allocate the histograms array
 174         Int_t i, npairdefs = (Int_t)fMixPairDefs->GetEntries();
 175         fMixHistograms = new TObjArray(npairdefs);
 176         TObjArray &histos = *fMixHistograms;
 177         AliPairDef *pd = 0;
 178         for (i = 0; i < npairdefs; i++) {
 179                 pd = (AliPairDef*)fMixPairDefs->At(i);
 180                 histos[i] = new TH1D(Form("hmix_%s", pd->GetName()), pd->GetTitle(), fNBins, fHistoMin, fHistoMax);
 181         }
 182
 183         // Link events branch
 184         AliRsnEvent *event2 = new AliRsnEvent;
 185         fEventsTree->SetBranchAddress("events", &event2);
 186
 187         // define variables to store data about particles
 188         Stat_t nPairs = 0;
 189         Int_t iev, ievmix, nEvents = (Int_t)fEventsTree->GetEntries();
 190
 191         // loop on events
 192         Int_t nmixed, mult1, mult2, diffMult;
 193         Double_t vz1, vz2, diffVz;
 194         for (iev = 0; iev < nEvents; iev++) {
 195
 196                 // get event
 197                 event2->Clear("DELETE");
 198                 fEventsTree->GetEntry(iev);
 199
 200                 // copy this event into 'event 1'
 201                 AliRsnEvent *event1 = new AliRsnEvent(*event2);
 202
 203                 // get Z position of primary vertex
 204                 vz1 = event1->GetPrimaryVertexZ();
 205
 206                 // if total multiplicities must be used
 207                 // it is computed here
 208                 if (compareTotMult) {
 209                         mult1 = event1->GetMultiplicity();
 210                 }
 211                 else {
 212                         mult1 = 0;
 213                 }
 214
 215                 // message
 216                 if (iev % 10 == 0) cout << "\rMixing event " << iev << flush;
 217
 218                 // loop on pair definitions
 219                 for (i = 0; i < npairdefs; i++) {
 220
 221                         // get pair definition
 222                         pd = (AliPairDef*)fMixPairDefs->At(i);
 223
 224                         // get histogram reference
 225                         TH1D *histogram = (TH1D*)histos[i];
 226
 227                         // get multiplicity of particle type '2' in event '1'
 228                         if (!mult1) {
 229                                 mult1 = (Int_t)event1->GetTracks(pd->GetSign2(), pd->GetParticle2())->GetEntries();
 230                         }
 231
 232                         // loop on events for mixing
 233                         nmixed = 0;
 234                         ievmix = iev;
 235                         while (nmixed < nmix) {
 236
 237                                 // get other event (it starts from iev + 1, and
 238                                 // loops again to 0 if reachs end of tree)
 239                                 ievmix++;
 240                                 if (ievmix >= nEvents) ievmix -= nEvents;
 241
 242                                 // if event-2-mix is equal to 'iev', then
 243                                 // a complete loop has been done, and no more
 244                                 // mixing can be done, even if they are too few
 245                                 // then the loop breaks anyway
 246                                 if (ievmix == iev) break;
 247
 248                                 // get other event
 249                                 event2->Clear("DELETE");
 250                                 fEventsTree->GetEntry(ievmix);
 251
 252                                 // get event stats related to event 2
 253                                 vz2 = event2->GetPrimaryVertexZ();
 254                                 if (compareTotMult) {
 255                                         mult2 = event2->GetMultiplicity();
 256                                 }
 257                                 else {
 258                                         mult2 = event2->GetTracks(pd->GetSign2(), pd->GetParticle2())->GetEntries();
 259                                 }
 260
 261                                 // fill histogram if checks are passed
 262                                 diffMult = TMath::Abs(mult1 - mult2);
 263                                 diffVz = TMath::Abs(vz1 - vz2);
 264                                 if (diffVz <= vzDiffMax && diffMult <= multDiffMax) {
 265                                         //cout << ievmix << " " << flush;
 266                                         nPairs += Compute(pd, histogram, event1, event2);
 267                                         nmixed++;
 268                                 }
 269                         }
 270                 }
 271
 272                 event1->Clear("DELETE");
 273                 delete event1;
 274                 event1 = 0;
 275
 276         } // end loop on events
 277         cout << endl;
 278
 279         return nPairs;
 280 }
 281 //--------------------------------------------------------------------------------------------------------
 282 Stat_t AliRsnAnalysis::Process()
 283 {
 284 //
 285 // Reads the list 'fPairDefs', and builds an inv-mass histogram for each definition.
 286 // For each event, particle of 'type 1' are combined with particles of 'type 2' as
 287 // defined in each pair definition specified in the list, taking the list of both types
 288 // from the same event.
 289 // This method is supposed to evaluate the 'signal' of the resonance, containing the peak.
 290 // It can also be used when one wants to evaluate the 'background' with the 'wrong sign'
 291 // particle pairs.
 292 //
 293         // allocate the histograms array in order to contain
 294         // as many objects as the number of pair definitionss
 295         Int_t i, npairdefs = (Int_t)fPairDefs->GetEntries();
 296         fHistograms = new TObjArray(npairdefs);
 297         TObjArray &histos = *fHistograms;
 298
 299         // loop over pair defs in order to retrieve the particle species chosen
 300         // which are used to set an unique name for the output histogram
 301         // There will be a direct correspondence between the inder of pairdef in its array
 302         // and the corresponding histogram in the 'fHistograms' list.
 303         AliPairDef *pd = 0;
 304         for (i = 0; i < npairdefs; i++) {
 305                 pd = (AliPairDef*)fPairDefs->At(i);
 306                 histos[i] = new TH1D(Form("h_%s", pd->GetName()), pd->GetTitle(), fNBins, fHistoMin, fHistoMax);
 307         }
 308
 309         // Link events branch
 310         AliRsnEvent *event = new AliRsnEvent;
 311         fEventsTree->SetBranchAddress("events", &event);
 312
 313         // define variables to store data about particles
 314         Stat_t nPairs = 0;
 315         Int_t iev, nEvents = (Int_t)fEventsTree->GetEntries();
 316
 317         // loop on events
 318         for (iev = 0; iev < nEvents; iev++) {
 319
 320                 // get event
 321                 event->Clear("DELETE");
 322                 fEventsTree->GetEntry(iev);
 323                 if (iev % 1000 == 0) cout << "\rProcessing event " << iev << flush;
 324
 325                 // loop over the collection of pair defs
 326                 for (i = 0; i < npairdefs; i++) {
 327                         pd = (AliPairDef*)fPairDefs->At(i);
 328                         TH1D *histogram = (TH1D*)histos[i];
 329                         // invoke AliPairDef method to fill histogram
 330                         nPairs += Compute(pd, histogram, event, event);
 331                 }
 332
 333         } // end loop on events
 334         cout << endl;
 335
 336         return nPairs;
 337 }
 338 //--------------------------------------------------------------------------------------------------------
 339 void AliRsnAnalysis::WriteHistograms() const
 340 {
 341 //
 342 // Writes histograms in current directory
 343 //
 344         TH1D *histogram;
 345         TObjArrayIter iter(fHistograms);
 346         while ( (histogram = (TH1D*)iter.Next()) ) {
 347                 histogram->Write();
 348         }
 349
 350         if (fMixHistograms) {
 351                 TObjArrayIter iterMix(fMixHistograms);
 352                 while ( (histogram = (TH1D*)iterMix.Next()) ) {
 353                         histogram->Write();
 354                 }
 355         }
 356 }
 357 //--------------------------------------------------------------------------------------------------------
 358 Stat_t AliRsnAnalysis::Compute
 359 (AliPairDef *pd, TH1D* &histogram, AliRsnEvent *event1, AliRsnEvent *event2)
 360 {
 361 //
 362 // Adds to the specified histogram the invariant mass spectrum calculated taking
 363 // particles of type 1 from event 1, and particles of type 2 from event 2.
 364 // Events can be equal (signal) or different (background with event mixing).
 365 //
 366         // define two 'cursor' objects
 367         AliRsnDaughter *track1 = 0, *track2 = 0;
 368
 369         // define iterators for the two collections
 370         if (!event1 || !event1->GetTracks(pd->GetSign1(), pd->GetParticle1())) {
 371                 Error("Compute", "Null pointer for particle 1 collection");
 372                 return 0.0;
 373         }
 374         if (!event2 || !event2->GetTracks(pd->GetSign2(), pd->GetParticle2())) {
 375                 Error("Compute", "Null pointer for particle 2 collection");
 376                 return 0.0;
 377         }
 378         TObjArrayIter iter1(event1->GetTracks(pd->GetSign1(), pd->GetParticle1()));
 379         TObjArrayIter iter2(event2->GetTracks(pd->GetSign2(), pd->GetParticle2()));
 380
 381         // define temporary variables for better code readability
 382         Stat_t nPairs = 0;
 383
 384         // loop on particle of type 1 (in event 1)
 385         while ( (track1 = (AliRsnDaughter*)iter1.Next()) ) {
 386
 387                 if (fRejectFakes && (track1->GetLabel() < 0)) continue;
 388                 if (!SingleCutCheck(pd->GetParticle1(), track1)) continue;
 389
 390                 iter2.Reset();
 391
 392                 // loop on particles of type 2 (in event 2)
 393                 while ( (track2 = (AliRsnDaughter*)iter2.Next()) ) {
 394
 395                         if (fRejectFakes && (track2->GetLabel() < 0)) continue;
 396                         if (event1 == event2 && track1->GetIndex() == track2->GetIndex()) continue;
 397                         if (!SingleCutCheck(pd->GetParticle2(), track2)) continue;
 398                         if (!PairCutCheck(track1, track2)) continue;
 399
 400                         /*
 401                         // check
 402                         Char_t sign1 = (track1->GetSign() > 0) ? '+' : '-';
 403                         Char_t sign2 = (track2->GetSign() > 0) ? '+' : '-';
 404                         Info("Compute", "Particle 1: PDG = %d, sign = %c --- Particle 2: PDG = %d, sign = %c", track1->GetTruePDG(), sign1, track2->GetTruePDG(), sign2);
 405                         */
 406
 407                         // total 4-momentum
 408                         track1->SetMass(pd->GetMass1());
 409                         track2->SetMass(pd->GetMass2());
 410                         AliRsnDaughter sum = AliRsnDaughter::Sum(*track1, *track2);
 411
 412                         // if the choice to get ONLY true pairs is selected, a check is made on the mothers
 413                         if (pd->GetOnlyTrue()) {
 414                                 // the 'sum' AliRsnDaughter object is initialized with
 415                                 // the PDG code of the common mother (if it is the case)
 416                                 if (TMath::Abs(sum.GetMotherPDG()) == TMath::Abs(fTrueMotherPDG)) {
 417                                         histogram->Fill(sum.GetMass());
 418                                         nPairs++;
 419                                 }
 420                         }
 421                         else {
 422                                 histogram->Fill(sum.GetMass());
 423                                 nPairs++;
 424                         }
 425
 426                 } // end loop 2
 427
 428         } // end loop 1
 429
 430         return nPairs;
 431 }
 432 //--------------------------------------------------------------------------------------------------------
 433 Bool_t AliRsnAnalysis::SingleCutCheck(Int_t itype, AliRsnDaughter *track) const
 434 {
 435 //
 436 // Checks a track against single particle cuts (if defined)
 437 //
 438         if (!fCuts[itype]) return kTRUE;
 439
 440         TObjArrayIter iter(fCuts[itype]);
 441         AliRsnDaughterCut *cut = 0;
 442         while ( (cut = (AliRsnDaughterCut*)iter.Next()) ) {
 443                 if (!cut->Pass(track)) return kFALSE;
 444         }
 445
 446         return kTRUE;
 447 }
 448 //--------------------------------------------------------------------------------------------------------
 449 Bool_t AliRsnAnalysis::PairCutCheck(AliRsnDaughter *track1, AliRsnDaughter *track2) const
 450 {
 451 //
 452 // Checks a pair against pair cuts (if defined)
 453 //
 454         if (!fPairCuts) return kTRUE;
 455
 456         TObjArrayIter iter(fPairCuts);
 457         AliRsnDaughterCutPair *cut = 0;
 458         while ( (cut = (AliRsnDaughterCutPair*)iter.Next()) ) {
 459                 if (!cut->Pass(track1, track2)) return kFALSE;
 460         }
 461
 462         return kTRUE;
 463 }
 464 //--------------------------------------------------------------------------------------------------------
 465 AliRsnAnalysis::AliPairDef::AliPairDef
 466 (AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2, Int_t pdgMother, Bool_t onlyTrue) :
 467   TNamed(),
 468   fOnlyTrue(onlyTrue),
 469   fTrueMotherPDG(pdgMother),
 470   fMass1(0.0),
 471   fSign1(sign1),
 472   fParticle1(p1),
 473   fMass2(0.0),
 474   fSign2(sign2),
 475   fParticle2(p2)
 476 {
 477 //
 478 // Constructor for nested class
 479 //
 480         if (!fOnlyTrue) fTrueMotherPDG = 0;
 481
 482         // instance a PDG database to recovery true masses of particles
 483         TDatabasePDG *db = TDatabasePDG::Instance();
 484         Int_t pdg1 = AliPID::ParticleCode((Int_t)p1);
 485         Int_t pdg2 = AliPID::ParticleCode((Int_t)p2);
 486         fMass1 = db->GetParticle(pdg1)->Mass();
 487         fMass2 = db->GetParticle(pdg2)->Mass();
 488
 489         // set name according to the chosen particles
 490         fName.Append(Form("%s(%c)_%s(%c)", ParticleName(p1), sign1, ParticleName(p2), sign2));
 491         fTitle.Append(Form("Inv. mass for %s(%c) and %s(%c)", ParticleName(p1), sign1, ParticleName(p2), sign2));
 492         if (onlyTrue) {
 493                 fName.Append("_true");
 494                 fTitle.Append(" (true pairs)");
 495         }
 496 }
 497 //--------------------------------------------------------------------------------------------------------
 498 Text_t* AliRsnAnalysis::AliPairDef::ParticleName(AliPID::EParticleType part) const
 499 {
 500 //
 501 // [PRIVATE]
 502 // Returns the name of the particle in text format
 503 //
 504         if (part == AliPID::kElectron) return ("E");
 505         else if (part == AliPID::kMuon) return ("Mu");
 506         else if (part == AliPID::kPion) return ("Pi");
 507         else if (part == AliPID::kKaon) return ("K");
 508         else if (part == AliPID::kProton) return ("P");
 509         else {
 510                 Warning("ParticleName", "Unrecognized value of EParticle argument");
 511                 return ("???");
 512         }
 513 }
 514 //--------------------------------------------------------------------------------------------------------