Updates and additions: Classes for signal and spectrum extraction; saving of
[u/mrichter/AliRoot.git] / PWG3 / dielectron / AliDielectronSpectrum.cxx
1
2 /*************************************************************************
3 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 *                                                                        *
5 * Author: The ALICE Off-line Project.                                    *
6 * Contributors are mentioned in the code where appropriate.              *
7 *                                                                        *
8 * Permission to use, copy, modify and distribute this software and its   *
9 * documentation strictly for non-commercial purposes is hereby granted   *
10 * without fee, provided that the above copyright notice appears in all   *
11 * copies and that both the copyright notice and this permission notice   *
12 * appear in the supporting documentation. The authors make no claims     *
13 * about the suitability of this software for any purpose. It is          *
14 * provided "as is" without express or implied warranty.                  *
15 **************************************************************************/
16
17 ///////////////////////////////////////////////////////////////////////////
18 //                Add general description                                 //
19 //                                                                       //
20 //                                                                       //
21 /*
22 Detailed description
23
24
25 */
26 //                                                                       //
27 ///////////////////////////////////////////////////////////////////////////
28
29 #include <TMath.h>
30 #include <TVectorT.h>
31 #include <TH1.h>
32
33 #include <AliCFContainer.h>
34 #include <AliCFGridSparse.h>
35
36 #include "AliDielectronSignalBase.h"
37 #include "AliDielectronSignalFunc.h"
38
39 #include "AliDielectronSpectrum.h"
40
41 ClassImp(AliDielectronSpectrum)
42
43 AliDielectronSpectrum::AliDielectronSpectrum() :
44   TNamed(),
45   fCFSignal(0x0),
46   fCFCorrection(0x0),
47   fCFSpectrum(0x0),
48   fCFCorrMatrix(0x0),
49   fStepSignal(kTRUE),
50   fStepSignificance(kFALSE),
51   fStepSOB(kFALSE),
52   fSignalStep(-1),
53   fCorrNom(-1),
54   fCorrDenom(-1),
55   fSignalMethods(0),
56   fVariables("Pt"),
57   fOwnerSpectrum(kTRUE),
58   fNvars(0),
59   fVars(0x0),
60   fNbins(0x0),
61   fCurrentBins(0x0),
62   fCurrentPositions(0x0)
63 {
64   //
65   // Default Constructor
66   //
67   fSignalMethods.SetOwner(kFALSE);
68 }
69
70 //______________________________________________
71 AliDielectronSpectrum::AliDielectronSpectrum(const char* name, const char* title) :
72   TNamed(name, title),
73   fCFSignal(0x0),
74   fCFCorrection(0x0),
75   fCFSpectrum(0x0),
76   fCFCorrMatrix(0x0),
77   fStepSignal(kTRUE),
78   fStepSignificance(kFALSE),
79   fStepSOB(kFALSE),
80   fSignalStep(-1),
81   fCorrNom(-1),
82   fCorrDenom(-1),
83   fSignalMethods(0),
84   fVariables("Pt"),
85   fOwnerSpectrum(kTRUE),
86   fNvars(0),
87   fVars(0x0),
88   fNbins(0x0),
89   fCurrentBins(0x0),
90   fCurrentPositions(0x0)
91 {
92   //
93   // Named Constructor
94   //
95   fSignalMethods.SetOwner(kFALSE);
96 }
97
98 //______________________________________________
99 AliDielectronSpectrum::~AliDielectronSpectrum()
100 {
101   //
102   // Default Destructor
103   //
104
105   if (fNbins) delete [] fNbins;
106   if (fVars) delete [] fVars;
107   if (fCFCorrMatrix) delete fCFCorrMatrix;
108   if ( fOwnerSpectrum && fCFSpectrum ) delete fCFSpectrum;
109 }
110
111 //______________________________________________
112 void AliDielectronSpectrum::Process()
113 {
114   //
115   // Extract signal and perform correction in the specified bins
116   //
117
118   //sanity checks
119   if (!fCFSignal){
120     AliError("Cannot perform signal extraction, no signal container set");
121     return;
122   }
123   
124   if (fSignalMethods.GetEntries()==0){
125     AliWarning("No Signal extraction method specified, using a default one");
126     AliDielectronSignalFunc *func=new AliDielectronSignalFunc("gaus+exp","Gauss for Signal and Exponential background");
127     func->SetDefaults(1);
128     fSignalMethods.Add(func);
129     fSignalMethods.SetOwner();
130   }
131
132   //setup configured variables
133   if (!SetupVariables()) return;
134   
135   //create container for the spectrum
136   CreateCFSpectrum();
137
138   if (!fCFSpectrum){
139     AliError("Could not create the Spectrum container");
140     return;
141   }
142
143   //get efficiency map if correction container is available
144   if (fCFCorrection&&!fCFCorrMatrix){
145     CreateCorrectionMatrix();
146   }
147
148   //loop over all configured bins and extract the signal
149   fCurrentBins=new Int_t[fNvars];
150   fCurrentPositions=new Double_t[fNvars];
151   ExtractSignalInBins();
152   delete [] fCurrentBins;
153   delete [] fCurrentPositions;
154   if (fSignalMethods.IsOwner()) {
155     fSignalMethods.Delete();
156     fSignalMethods.SetOwner(kFALSE);
157   }
158   
159 }
160
161 //______________________________________________
162 Bool_t AliDielectronSpectrum::SetupVariables()
163 {
164   //
165   // Setup the variables arrays
166   //
167   
168   TObjArray *arr=fVariables.Tokenize(":");
169   fNvars=arr->GetEntries();
170   fVars=new Int_t[fNvars];
171   fNbins=new Int_t[fNvars];
172   
173   for (Int_t iVar=0; iVar<fNvars; ++iVar){
174     fVars[iVar]=fCFSignal->GetVar(arr->UncheckedAt(iVar)->GetName());
175     if (fVars[iVar]==-1){
176       AliError(Form("Variable '%s' not found in Signal container!",arr->UncheckedAt(iVar)->GetName()));
177       delete [] fVars;
178       fVars=0x0;
179       delete [] fNbins;
180       fNbins=0x0;
181       delete arr;
182       return kFALSE;
183     }
184     
185     fNbins[iVar]=fCFSignal->GetNBins(fVars[iVar]);
186   }
187   delete arr;
188   return kTRUE;
189 }
190
191 //______________________________________________
192 void AliDielectronSpectrum::CreateCFSpectrum()
193 {
194   //
195   // Create CF container for the spectrum
196   //
197
198   Int_t nAddStep=0;
199   if (fStepSignal)       nAddStep+=2;
200   if (fStepSignificance) ++nAddStep;
201   if (fStepSOB)          ++nAddStep;
202   
203   Int_t nStep=nAddStep*(fSignalMethods.GetEntries());
204   if (fSignalMethods.GetEntries()>1) nStep+=nAddStep;
205
206   fCFSpectrum = new AliCFContainer(GetName(), GetTitle(), nStep, fNvars, fNbins);
207
208   // initialize the variables and their bin limits
209   for (Int_t iVar=0; iVar<fNvars; iVar++) {
210     fCFSpectrum->SetBinLimits(iVar, fCFSignal->GetBinLimits(fVars[iVar]));
211     fCFSpectrum->SetVarTitle(iVar, fCFSignal->GetVarTitle(fVars[iVar]));
212   }
213
214   // setup step titles
215   Int_t steps=0;
216   for (Int_t iMethod=0; iMethod<fSignalMethods.GetEntries(); ++iMethod){
217     TString name(fSignalMethods.UncheckedAt(iMethod)->GetName());
218     if (fStepSignal){
219       fCFSpectrum->SetStepTitle(steps++,(name+" (Siganl)").Data());
220       fCFSpectrum->SetStepTitle(steps++,(name+" (Corrected Siganl)").Data());
221     }
222     if (fStepSignificance){
223       fCFSpectrum->SetStepTitle(steps++,(name+" (Significance)").Data());
224     }
225     if (fStepSOB){
226       fCFSpectrum->SetStepTitle(steps++,(name+" (S/B)").Data());
227     }
228   }
229
230   if (fSignalMethods.GetEntries()>1){
231     fCFSpectrum->SetStepTitle(steps++,"Mean of methods");
232     fCFSpectrum->SetStepTitle(steps++,"Mean of methods (Corrected)");
233   }
234
235   if (nStep!=steps){
236     AliError("Something went wrong in the step creation");
237     delete fCFSpectrum;
238     fCFSpectrum=0x0;
239   }
240 }
241
242 //______________________________________________
243 void AliDielectronSpectrum::CreateCorrectionMatrix()
244 {
245   //
246   // Get the correction matrix for the corresponding variables
247   //
248
249   if (!fCFCorrection) return;
250   
251   TObjArray *arr=fVariables.Tokenize(":");
252   Int_t nvars=arr->GetEntries();
253   Int_t *vars=new Int_t[nvars];
254   
255   for (Int_t iVar=0; iVar<fNvars; ++iVar){
256     vars[iVar]=fCFCorrection->GetVar(arr->UncheckedAt(iVar)->GetName());
257     if (vars[iVar]==-1){
258       AliError(Form("Variable '%s' not found in Correction container!",arr->UncheckedAt(iVar)->GetName()));
259       delete [] vars;
260       delete arr;
261       return;
262     }
263   }
264   delete arr;
265   
266   fCFCorrMatrix  =fCFCorrection->GetGrid(fCorrNom)->Project(nvars,vars,0,0);
267   AliCFGridSparse *hnDeNom=fCFCorrection->GetGrid(fCorrDenom)->Project(nvars,vars,0,0);
268   fCFCorrMatrix->Divide(hnDeNom);
269   delete hnDeNom;
270 }
271
272 //______________________________________________
273 void AliDielectronSpectrum::ExtractSignalInBins(Int_t variable)
274 {
275   //
276   // recursively loop over bins and extract signal
277   //
278
279   Int_t varPairType=fCFSignal->GetVar("PairType");
280   Int_t varMass=fCFSignal->GetVar("M");
281   
282   for (Int_t ibin=0; ibin<fNbins[variable]; ++ibin){
283     Int_t bin=ibin+1;
284     fCFSignal->GetGrid(fSignalStep)->GetGrid()->GetAxis(fVars[variable])->SetRange(bin,bin);
285     fCurrentBins[variable]=bin;
286     fCurrentPositions[variable]=fCFSignal->GetGrid(fSignalStep)->GetBinCenter(fVars[variable],bin);
287     
288     if (variable != fNvars-1) ExtractSignalInBins(variable+1);
289     
290     TObjArray arrPairTypes(10);
291     arrPairTypes.SetOwner();
292   
293     for (Int_t itype=0; itype<3; ++itype){
294 //     Int_t itype=1;
295       fCFSignal->SetRangeUser(varPairType,itype,itype,fSignalStep);
296       TH1 *h=fCFSignal->Project(varMass,fSignalStep);
297       h->SetDirectory(0);
298       arrPairTypes.AddAt(h,itype);
299     }
300     AliInfo(Form("Processing bin: %d (%.2f)",ibin, fCurrentPositions[variable]));
301     //loop over all signal extraction methods and retrieve signals
302     for (Int_t iMethod=0; iMethod<fSignalMethods.GetEntries(); ++iMethod){
303       AliDielectronSignalBase *signalMethod=(AliDielectronSignalBase*)fSignalMethods.At(iMethod);
304       signalMethod->Process(&arrPairTypes);
305       const TVectorD &val=signalMethod->GetValues();
306       const TVectorD &err=signalMethod->GetErrors();
307       
308       //Fill sparse containers
309       Int_t step=0;
310       if (fStepSignal){
311         //signal
312         fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(0));
313         fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(0));
314         ++step;
315
316         //corrected signal
317         if (fCFCorrMatrix){
318           Float_t corrFactor = fCFCorrMatrix->GetElement(fCurrentPositions);
319           Float_t corrError  = fCFCorrMatrix->GetElementError(fCurrentPositions);
320
321           Float_t value=val(0)*corrFactor;
322           fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,value);
323           Float_t error=TMath::Sqrt( (err(0)/val(0))*(err(0)/val(0)) +
324                                      (corrError/corrFactor)*(corrError/corrFactor)
325                                    )*value;
326           fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,error);
327 //           printf("corrFactor: %f+-%f\n",value,error);
328         }
329         ++step;
330       }
331
332       if (fStepSignificance) {
333         fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(2));
334         fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(2));
335         ++step;
336       }
337       
338       if (fStepSOB) {
339         fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(3));
340         fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(3));
341         ++step;
342       }
343     }// end method loop
344     
345     arrPairTypes.Delete();
346     
347   }// end bin loop
348 }
349