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