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