08-mar-2003 NvE Compiler option /GR introduced for MSVC++ in mklibs.bat to explicitly...
[u/mrichter/AliRoot.git] / 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 ////////////////////////////////////////////////////////////////////////////////