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