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