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