Adding Id to PWG3 classes for better tracking of the coverity defect fixes (Ivana)
[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 /* $Id$ */
18
19 ///////////////////////////////////////////////////////////////////////////
20 //                Add general description                                 //
21 //                                                                       //
22 //                                                                       //
23 /*
24 Detailed description
25
26
27 */
28 //                                                                       //
29 ///////////////////////////////////////////////////////////////////////////
30
31 #include <TMath.h>
32 #include <TVectorT.h>
33 #include <TH1.h>
34 #include <TROOT.h>
35 #include <TCanvas.h>
36
37 #include <AliCFContainer.h>
38 #include <AliCFGridSparse.h>
39
40 #include "AliDielectronSignalBase.h"
41 #include "AliDielectronSignalFunc.h"
42
43 #include "AliDielectronSpectrum.h"
44
45 ClassImp(AliDielectronSpectrum)
46
47 AliDielectronSpectrum::AliDielectronSpectrum() :
48   TNamed(),
49   fCFSignal(0x0),
50   fCFCorrection(0x0),
51   fCFSpectrum(0x0),
52   fCFCorrMatrix(0x0),
53   fStepSignal(kTRUE),
54   fStepSignificance(kFALSE),
55   fStepSOB(kFALSE),
56   fStepMass(kFALSE),
57   fStepMassWidth(kFALSE),
58   fSignalStep(-1),
59   fCorrNom(-1),
60   fCorrDenom(-1),
61   fSignalMethods(0),
62   fVariables("Pt"),
63   fOwnerSpectrum(kTRUE),
64   fVisualDebug(kFALSE),
65   fNvars(0),
66   fVars(0x0),
67   fNbins(0x0),
68   fCurrentBins(0x0),
69   fCurrentPositions(0x0)
70 {
71   //
72   // Default Constructor
73   //
74   fSignalMethods.SetOwner(kFALSE);
75 }
76
77 //______________________________________________
78 AliDielectronSpectrum::AliDielectronSpectrum(const char* name, const char* title) :
79   TNamed(name, title),
80   fCFSignal(0x0),
81   fCFCorrection(0x0),
82   fCFSpectrum(0x0),
83   fCFCorrMatrix(0x0),
84   fStepSignal(kTRUE),
85   fStepSignificance(kFALSE),
86   fStepSOB(kFALSE),
87   fStepMass(kFALSE),
88   fStepMassWidth(kFALSE),
89   fSignalStep(-1),
90   fCorrNom(-1),
91   fCorrDenom(-1),
92   fSignalMethods(0),
93   fVariables("Pt"),
94   fOwnerSpectrum(kTRUE),
95   fVisualDebug(kFALSE),
96   fNvars(0),
97   fVars(0x0),
98   fNbins(0x0),
99   fCurrentBins(0x0),
100   fCurrentPositions(0x0)
101 {
102   //
103   // Named Constructor
104   //
105   fSignalMethods.SetOwner(kFALSE);
106 }
107
108 //______________________________________________
109 AliDielectronSpectrum::~AliDielectronSpectrum()
110 {
111   //
112   // Default Destructor
113   //
114
115   if (fNbins) delete [] fNbins;
116   if (fVars) delete [] fVars;
117   if (fCFCorrMatrix) delete fCFCorrMatrix;
118   if ( fOwnerSpectrum && fCFSpectrum ) delete fCFSpectrum;
119 }
120
121 //______________________________________________
122 void AliDielectronSpectrum::Process()
123 {
124   //
125   // Extract signal and perform correction in the specified bins
126   //
127
128   //sanity checks
129   if (!fCFSignal){
130     AliError("Cannot perform signal extraction, no signal container set");
131     return;
132   }
133   
134   if (fSignalMethods.GetEntries()==0){
135     AliWarning("No Signal extraction method specified, using a default one");
136     AliDielectronSignalFunc *func=new AliDielectronSignalFunc("gaus+exp","Gauss for Signal and Exponential background");
137     func->SetDefaults(1);
138     fSignalMethods.Add(func);
139     fSignalMethods.SetOwner();
140   }
141
142   //setup configured variables
143   if (!SetupVariables()) return;
144   
145   //create container for the spectrum
146   CreateCFSpectrum();
147
148   if (!fCFSpectrum){
149     AliError("Could not create the Spectrum container");
150     return;
151   }
152
153   //get efficiency map if correction container is available
154   if (fCFCorrection&&!fCFCorrMatrix){
155     CreateCorrectionMatrix();
156   }
157
158   //loop over all configured bins and extract the signal
159   fCurrentBins=new Int_t[fNvars];
160   fCurrentPositions=new Double_t[fNvars];
161   ExtractSignalInBins();
162   delete [] fCurrentBins;
163   delete [] fCurrentPositions;
164   if (fSignalMethods.IsOwner()) {
165     fSignalMethods.Delete();
166     fSignalMethods.SetOwner(kFALSE);
167   }
168   
169 }
170
171 //______________________________________________
172 Bool_t AliDielectronSpectrum::SetupVariables()
173 {
174   //
175   // Setup the variables arrays
176   //
177   
178   TObjArray *arr=fVariables.Tokenize(":");
179   fNvars=arr->GetEntries();
180   fVars=new Int_t[fNvars];
181   fNbins=new Int_t[fNvars];
182   
183   for (Int_t iVar=0; iVar<fNvars; ++iVar){
184     fVars[iVar]=fCFSignal->GetVar(arr->UncheckedAt(iVar)->GetName());
185     if (fVars[iVar]==-1){
186       AliError(Form("Variable '%s' not found in Signal container!",arr->UncheckedAt(iVar)->GetName()));
187       delete [] fVars;
188       fVars=0x0;
189       delete [] fNbins;
190       fNbins=0x0;
191       delete arr;
192       return kFALSE;
193     }
194     
195     fNbins[iVar]=fCFSignal->GetNBins(fVars[iVar]);
196   }
197   delete arr;
198   return kTRUE;
199 }
200
201 //______________________________________________
202 void AliDielectronSpectrum::CreateCFSpectrum()
203 {
204   //
205   // Create CF container for the spectrum
206   //
207
208   Int_t nAddStep=0;
209   if (fStepSignal)       nAddStep+=2;
210   if (fStepSignificance) ++nAddStep;
211   if (fStepSOB)          ++nAddStep;
212   if (fStepMass)         ++nAddStep;
213   if (fStepMassWidth)    ++nAddStep;
214   
215   Int_t nStep=nAddStep*(fSignalMethods.GetEntries());
216   if (fSignalMethods.GetEntries()>1) nStep+=nAddStep;
217
218   fCFSpectrum = new AliCFContainer(GetName(), GetTitle(), nStep, fNvars, fNbins);
219
220   // initialize the variables and their bin limits
221   for (Int_t iVar=0; iVar<fNvars; iVar++) {
222     fCFSpectrum->SetBinLimits(iVar, fCFSignal->GetBinLimits(fVars[iVar]));
223     fCFSpectrum->SetVarTitle(iVar, fCFSignal->GetVarTitle(fVars[iVar]));
224   }
225
226   // setup step titles
227   Int_t steps=0;
228   for (Int_t iMethod=0; iMethod<fSignalMethods.GetEntries(); ++iMethod){
229     TString name(fSignalMethods.UncheckedAt(iMethod)->GetName());
230     if (fStepSignal){
231       fCFSpectrum->SetStepTitle(steps++,(name+" (Signal)").Data());
232       fCFSpectrum->SetStepTitle(steps++,(name+" (Corrected Signal)").Data());
233     }
234     if (fStepSignificance){
235       fCFSpectrum->SetStepTitle(steps++,(name+" (Significance)").Data());
236     }
237     if (fStepSOB){
238       fCFSpectrum->SetStepTitle(steps++,(name+" (S/B)").Data());
239     }
240     if (fStepMass){
241       fCFSpectrum->SetStepTitle(steps++,(name+" (Mass)").Data());
242     }
243     if (fStepMassWidth){
244       fCFSpectrum->SetStepTitle(steps++,(name+" (Mass width)").Data());
245     }
246   }
247
248   if (fSignalMethods.GetEntries()>1){
249     fCFSpectrum->SetStepTitle(steps++,"Mean of methods");
250     fCFSpectrum->SetStepTitle(steps++,"Mean of methods (Corrected)");
251   }
252
253   if (nStep!=steps){
254     AliError("Something went wrong in the step creation");
255     delete fCFSpectrum;
256     fCFSpectrum=0x0;
257   }
258 }
259
260 //______________________________________________
261 void AliDielectronSpectrum::CreateCorrectionMatrix()
262 {
263   //
264   // Get the correction matrix for the corresponding variables
265   //
266
267   if (!fCFCorrection) return;
268   
269   TObjArray *arr=fVariables.Tokenize(":");
270   Int_t nvars=arr->GetEntries();
271   Int_t *vars=new Int_t[nvars];
272   
273   for (Int_t iVar=0; iVar<fNvars; ++iVar){
274     vars[iVar]=fCFCorrection->GetVar(arr->UncheckedAt(iVar)->GetName());
275     if (vars[iVar]==-1){
276       AliError(Form("Variable '%s' not found in Correction container!",arr->UncheckedAt(iVar)->GetName()));
277       delete [] vars;
278       delete arr;
279       return;
280     }
281   }
282   delete arr;
283   
284   fCFCorrMatrix  =fCFCorrection->GetGrid(fCorrNom)->Project(nvars,vars,0,0);
285   AliCFGridSparse *hnDeNom=fCFCorrection->GetGrid(fCorrDenom)->Project(nvars,vars,0,0);
286   fCFCorrMatrix->Divide(hnDeNom);
287   delete hnDeNom;
288 }
289
290 //______________________________________________
291 void AliDielectronSpectrum::ExtractSignalInBins(Int_t variable)
292 {
293   //
294   // recursively loop over bins and extract signal
295   //
296
297   Int_t varPairType=fCFSignal->GetVar("PairType");
298   Int_t varMass=fCFSignal->GetVar("M");
299   
300   for (Int_t ibin=0; ibin<fNbins[variable]; ++ibin){
301     Int_t bin=ibin+1;
302     fCFSignal->GetGrid(fSignalStep)->GetGrid()->GetAxis(fVars[variable])->SetRange(bin,bin);
303     fCurrentBins[variable]=bin;
304     fCurrentPositions[variable]=fCFSignal->GetGrid(fSignalStep)->GetBinCenter(fVars[variable],bin);
305     
306     if (variable != fNvars-1) ExtractSignalInBins(variable+1);
307     
308     TObjArray arrPairTypes(10);
309     arrPairTypes.SetOwner();
310   
311     for (Int_t itype=0; itype<3; ++itype){
312 //     Int_t itype=1;
313       fCFSignal->SetRangeUser(varPairType,itype,itype,fSignalStep);
314       TH1 *h=fCFSignal->Project(varMass,fSignalStep);
315       h->SetDirectory(0);
316       arrPairTypes.AddAt(h,itype);
317     }
318     AliInfo(Form("Processing bin: %d (%.2f)",ibin, fCurrentPositions[variable]));
319     //loop over all signal extraction methods and retrieve signals
320     for (Int_t iMethod=0; iMethod<fSignalMethods.GetEntries(); ++iMethod){
321       AliDielectronSignalBase *signalMethod=(AliDielectronSignalBase*)fSignalMethods.At(iMethod);
322       signalMethod->Process(&arrPairTypes);
323       if (fVisualDebug){
324         TCanvas *c=(TCanvas*)gROOT->GetListOfCanvases()->FindObject("SpectrumVisualDebug");
325         if (!c) c=new TCanvas("SpectrumVisualDebug","SpectrumVisualDebug");
326         c->Clear();
327         TString binsProc;
328         for (Int_t ivar=0; ivar<fNvars; ++ivar) binsProc+=Form("_%02d",fCurrentBins[ivar]);
329         signalMethod->Draw("stat");
330         binsProc.Append(".png");
331         binsProc.Prepend("SpectrumVisualDebug");
332         c->Update();
333         c->Print(binsProc.Data());
334       }
335       const TVectorD &val=signalMethod->GetValues();
336       const TVectorD &err=signalMethod->GetErrors();
337       
338       //Fill sparse containers
339       Int_t step=0;
340       if (fStepSignal){
341         //signal
342         fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(0));
343         fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(0));
344         ++step;
345
346         //corrected signal
347         if (fCFCorrMatrix){
348           Float_t corrFactor = fCFCorrMatrix->GetElement(fCurrentPositions);
349           Float_t corrError  = fCFCorrMatrix->GetElementError(fCurrentPositions);
350
351           Float_t value=val(0)*corrFactor;
352           fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,value);
353           Float_t error=TMath::Sqrt( (err(0)/val(0))*(err(0)/val(0)) +
354                                      (corrError/corrFactor)*(corrError/corrFactor)
355                                    )*value;
356           fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,error);
357 //           printf("corrFactor: %f+-%f\n",value,error);
358         }
359         ++step;
360       }
361
362       if (fStepSignificance) {
363         fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(2));
364         fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(2));
365         ++step;
366       }
367       
368       if (fStepSOB) {
369         fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(3));
370         fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(3));
371         ++step;
372       }
373
374       if (fStepMass) {
375         fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(4));
376         fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(4));
377         ++step;
378       }
379       
380       if (fStepMassWidth) {
381         fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(5));
382         fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(5));
383         ++step;
384       }
385     }// end method loop
386     
387     arrPairTypes.Delete();
388     
389   }// end bin loop
390 }
391