RALICE/AliInvmass.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 // $Id$
  17
  18 ////////////////////////////////////////////////////////////////////////////////
  19 // Class AliInvmass
  20 // Construction of invariant mass and combinatorial background.
  21 //
  22 // Example :
  23 // ---------
  24 //
  25 // TObjArray* photons=new TObjArray(); // Array with photon tracks for pi0 rec.
  26 //
  27 // // Code to create some photon tracks from pi0 decays
  28 // Int_t ntracks=200;
  29 // for (Int_t i=0; i<ntracks; i++)
  30 // {
  31 //  photons->Add(new Alitrack);
  32 //  ...
  33 //  ...
  34 //  ...
  35 // }
  36 //
  37 // // Perform the invariant mass and comb. bkg. reconstruction
  38 //
  39 // TObjArray* allm=q.Invmass(photons,photons); // All reconstructed invariant masses
  40 //
  41 // TH1F* hall=new TH1F("hall","hall",200,0,2); // Histo with M_inv of all combinations
  42 //
  43 // Int_t nall=0;
  44 // if (allm) nall=allm->GetEntries();
  45 //
  46 // AliTrack* t;
  47 // Float_t minv;
  48 // for (Int_t j=0; j<nall; j++)
  49 // {
  50 //  t=(AliTrack*)allm->At(j);
  51 //  if (t)
  52 //  {
  53 //   minv=t->GetMass();
  54 //   hall->Fill(minv);
  55 //  }
  56 // }
  57 //
  58 // TObjArray* bkgm=q.CombBkg(photons,photons); // Reconstructed comb. background
  59 //
  60 // TH1F* hbkg=new TH1F("hbkg","hbkg",200,0,2); // Histo with M_inv. of comb. background
  61 //
  62 // Int_t nbkg=0;
  63 // if (bkgm) nbkg=bkgm->GetEntries();
  64 //
  65 // for (Int_t j=0; j<nbkg; j++)
  66 // {
  67 //  t=(AliTrack*)bkgm->At(j);
  68 //  if (t)
  69 //  {
  70 //   minv=t->GetMass();
  71 //   hbkg->Fill(minv);
  72 //  }
  73 // }
  74 //
  75 // TH1F* hsig=new TH1F("sig","sig",200,0,2);   // Histo with the bkg. subtracted signal
  76 // hsig->Sumw2();
  77 // hsig->Add(hall,hbkg,1,-1);
  78 //
  79 //
  80 // Note : By default the storage of the reconstructed information is performed
  81 //        in separate TObjArrays for the signal and comb. background resp.
  82 //        In order to limit the memory usage, AliInvmass::SetStorageMode(1) may be
  83 //        used to activate only a single TObjArray to store the reconstructed information.
  84 //        Consequently, the following statements
  85 //
  86 //         TObjArray* allm=q.Invmass(photons,photons);
  87 //         TObjArray* bkgm=q.CombBkg(photons,photons);
  88 //
  89 //        will result in the fact that after he invokation of CombBkg
  90 //        the information of "allm" is lost due to the fact that the storage is
  91 //        is re-used for "bkgm" in case the "single storage" option has been selected.
  92 //        Usage of the, in that case invalid, pointer "allm" may cause your
  93 //        program to crash.
  94 //
  95 //        * Thus : In case of single storage usage, all invokations of the returned
  96 //                 array pointer have to be completed before invoking any memberfunction
  97 //                 of the same AliInvmass object again.
  98 //
  99 //
 100 //
 101 //--- Author: Nick van Eijndhoven 12-apr-1999 UU-SAP Utrecht
 102 //- Modified: NvE $Date$ UU-SAP Utrecht
 103 ////////////////////////////////////////////////////////////////////////////////
 104
 105 #include "AliInvmass.h"
 106 #include "Riostream.h"
 107
 108 ClassImp(AliInvmass) // Class implementation to enable ROOT I/O
 109
 110 AliInvmass::AliInvmass()
 111 {
 112 // Creation of an AliInvmass object and initialisation of parameters
 113  fPi=acos(-1.);
 114  fMode=2;
 115  fBkg=0;
 116  fNewtheta=1;
 117  fNewphi=1;
 118  fMinv=0;
 119  fMbkg=0;
 120 }
 121 ////////////////////////////////////////////////////////////////////////////////
 122 AliInvmass::~AliInvmass()
 123 {
 124 // Destructor to delete dynamically allocated memory
 125  if (fMinv)
 126  {
 127   delete fMinv;
 128   fMinv=0;
 129  }
 130
 131  if (fMbkg)
 132  {
 133   delete fMbkg;
 134   fMbkg=0;
 135  }
 136 }
 137 ////////////////////////////////////////////////////////////////////////////////
 138 void AliInvmass::SetStorageMode(Int_t m)
 139 {
 140 // Set storage mode for the result arrays for inv. mass and comb. background
 141  fMode=2;
 142  if (m==1) fMode=1;
 143 }
 144 ////////////////////////////////////////////////////////////////////////////////
 145 void AliInvmass::SetThetaSwitch(Int_t i)
 146 {
 147 // Enable/Disable (1/0) switching of theta angle in comb. bkg. reconstruction.
 148 // Default : Switching of theta is enabled.
 149  fNewtheta=1;
 150  if (i==0) fNewtheta=0;
 151 }
 152 ////////////////////////////////////////////////////////////////////////////////
 153 void AliInvmass::SetPhiSwitch(Int_t i)
 154 {
 155 // Enable/Disable (1/0) switching of phi angle in comb. bkg. reconstruction.
 156 // Default : Switching of phi is enabled.
 157  fNewphi=1;
 158  if (i==0) fNewphi=0;
 159 }
 160 ////////////////////////////////////////////////////////////////////////////////
 161 Int_t AliInvmass::GetStorageMode()
 162 {
 163 // Provide mode of storage for the result arrays for inv. mass and comb. background
 164  return fMode;
 165 }
 166 ////////////////////////////////////////////////////////////////////////////////
 167 Int_t AliInvmass::GetThetaSwitch()
 168 {
 169 // Provide the theta switching flag
 170  return fNewtheta;
 171 }
 172 ////////////////////////////////////////////////////////////////////////////////
 173 Int_t AliInvmass::GetPhiSwitch()
 174 {
 175 // Provide the phi switching flag
 176  return fNewphi;
 177 }
 178 ////////////////////////////////////////////////////////////////////////////////
 179 TObjArray* AliInvmass::Invmass(TObjArray* a1,TObjArray* a2)
 180 {
 181 // Perform two-particle invariant mass reconstruction
 182  fBkg=0;
 183  Combine(a1,a2);
 184  return fMinv;
 185 }
 186 ////////////////////////////////////////////////////////////////////////////////
 187 TObjArray* AliInvmass::CombBkg(TObjArray* a1,TObjArray* a2)
 188 {
 189 // Perform two-particle combinatorial background reconstruction
 190  fBkg=1;
 191  Combine(a1,a2);
 192  if (fMode != 1)
 193  {
 194   return fMbkg;
 195  }
 196  else
 197  {
 198   return fMinv;
 199  }
 200 }
 201 ////////////////////////////////////////////////////////////////////////////////
 202 void AliInvmass::Combine(TObjArray* a1,TObjArray* a2)
 203 {
 204 // Perform two-particle invariant mass reconstruction
 205
 206  if ((!fBkg || fMode==1) && fMinv)
 207  {
 208   delete fMinv;
 209   fMinv=0;
 210  }
 211
 212  if (fBkg && (fMode !=1) && fMbkg)
 213  {
 214   delete fMbkg;
 215   fMbkg=0;
 216  }
 217
 218  Int_t isame; // Indicates whether both lists are identical
 219  isame=0;
 220  if (a1==a2) isame=1;
 221
 222  // Index i must loop over the shortest of a1 and a2
 223  TObjArray* listi=a1;
 224  TObjArray* listj=a2;
 225  Int_t ni=a1->GetEntries();
 226  Int_t nj=a2->GetEntries();
 227  if (nj < ni)
 228  {
 229   ni=a2->GetEntries();
 230   nj=a1->GetEntries();
 231   listi=a2;
 232   listj=a1;
 233  }
 234
 235  AliTrack* p1=0;
 236  AliTrack* p2=0;
 237  AliTrack* px=0;
 238  Ali4Vector ptot;
 239  AliTrack* t=0;
 240  Double_t v2[4],vx[4];
 241  Float_t q1,q2;
 242
 243  Int_t jmin; // Start index for list j
 244  Int_t jx;   // Index for randomly picked particle for comb. bkg. reconstruction
 245
 246  for (Int_t i=0; i<ni; i++) // Select first a particle from list i
 247  {
 248   p1=(AliTrack*)listi->At(i);
 249   p2=0;
 250
 251   if (!p1) continue;
 252
 253   jmin=0;
 254   if (isame) jmin=i+1;
 255   for (Int_t j=jmin; j<nj; j++) // Select also a particle from list j
 256   {
 257    p2=(AliTrack*)listj->At(j);
 258    if (p1==p2) p2=0; // Don't combine particle with itself
 259
 260    if (!p2) continue;
 261
 262    p2->GetVector(v2,"sph");
 263
 264    // Take theta and phi from randomly chosen other list j particle for bkg. reconstr.
 265    if (fBkg)
 266    {
 267     px=0;
 268     if ((!isame && nj>1) || (isame && nj>2))
 269     {
 270      jx=int(fRndm.Uniform(0,float(nj)));
 271      px=(AliTrack*)listj->At(jx);
 272
 273      while (!px || px==p2 || px==p1)
 274      {
 275       jx++;
 276       if (jx >= nj) jx=0;
 277       px=(AliTrack*)listj->At(jx);
 278      }
 279
 280      px->GetVector(vx,"sph");
 281      if (fNewtheta) v2[2]=vx[2]; // Replace the theta angle in the v2 vector
 282      if (fNewphi) v2[3]=vx[3]; // Replace the phi angle in the v2 vector
 283     }
 284    }
 285
 286    if ((!fBkg && p2) || (fBkg && px))
 287    {
 288     // Store the data of this two-particle combination
 289     ptot.SetVector(v2,"sph");
 290     ptot=(Ali4Vector)(ptot+(*p1));
 291     q1=p1->GetCharge();
 292     q2=p2->GetCharge();
 293     t=new AliTrack;
 294     t->Set4Momentum(ptot);
 295     t->SetCharge(q1+q2);
 296     if (!fBkg || fMode==1)
 297     {
 298      if (!fMinv)
 299      {
 300       fMinv=new TObjArray();
 301       fMinv->SetOwner();
 302      }
 303      fMinv->Add(t);
 304     }
 305     else
 306     {
 307      if (!fMbkg)
 308      {
 309       fMbkg=new TObjArray();
 310       fMbkg->SetOwner();
 311      }
 312      fMbkg->Add(t);
 313     }
 314    }
 315
 316   } // End of second particle loop
 317
 318  } // End of first particle loop
 319 }
 320 ////////////////////////////////////////////////////////////////////////////////