a54c867e316d2f97fc335fe0d2bd5303521ac8ec
[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 /*
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 ////////////////////////////////////////////////////////////////////////////////