PWGDQ/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+" (Signal)").Data());
 230       fCFSpectrum->SetStepTitle(steps++,(name+" (Corrected Signal)").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