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