PWG/FLOW/Base/AliFlowMSPHistograms.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 #include <iostream>
  16
  17 #include "AliFlowMSPHistograms.h"
  18
  19 #include <TObjString.h>
  20 #include <TList.h>
  21
  22 AliFlowMSPHistograms::AliFlowMSPHistograms()
  23    : TNamed(),fVarNames(),fPAvg(),fPProd(),fXName()
  24 {
  25    // Default constructor for root IO
  26    //std::cerr << "Default constructor called" << std::endl;
  27 }
  28
  29 AliFlowMSPHistograms::AliFlowMSPHistograms(const AliFlowMSPHistograms &x)
  30    : TNamed(),fVarNames(),fPAvg(x.fPAvg),fPProd(x.fPProd),fXName(x.fXName)
  31 {
  32    // Copy constructor, clone all the TProfiles from the other object
  33    // The copy constructor made a shallow copy, i.e. all entries point to the originial profiles.
  34    // Therefore we need to Clone each entry.
  35    for( int i=0; i<fVarNames.GetEntries(); ++i ) {
  36       TObjString *org=(TObjString *)(fVarNames[i]);
  37       fVarNames[i]=new TObjString(*org);
  38    }
  39    fVarNames.SetOwner(kTRUE);
  40
  41    for( int i=0; i<fPAvg.GetEntries(); ++i ) {
  42       TProfile *org=(TProfile *)(fPAvg[i]);
  43       fPAvg[i]=org->Clone();
  44    }
  45    fPAvg.SetOwner(kTRUE);
  46
  47    for( int i=0; i<fPProd.GetEntries(); ++i ) {
  48       TProfile *org=(TProfile *)(fPProd[i]);
  49       fPProd[i]=org->Clone();
  50    }
  51    fPProd.SetOwner(kTRUE);
  52 }
  53
  54 AliFlowMSPHistograms::AliFlowMSPHistograms(const int dim, const char *name, const int nBins, const double xmin, const double xmax)
  55    : TNamed(name,"MSP histograms"),fVarNames(),fPAvg(),fPProd(),fXName()
  56 {
  57    Init(dim, name, nBins, xmin, xmax);
  58 }
  59
  60 AliFlowMSPHistograms::AliFlowMSPHistograms(const int dim, const char *name, const int nBins, const double *xBins)
  61    : TNamed(name,"MSP histograms"),fVarNames(),fPAvg(),fPProd(),fXName()
  62 {
  63    Init(dim, name, nBins, xBins);
  64 }
  65
  66 AliFlowMSPHistograms::~AliFlowMSPHistograms()
  67 {
  68    // Destructor. All histograms are owned by this object and will be deleted here.
  69    fPAvg.Clear();
  70    fPProd.Clear();
  71 }
  72
  73 void AliFlowMSPHistograms::Init(const int dim, const char *name, const int nBins, const double *xBins)
  74 {
  75    // Initialize all profiles after deleting the existing profiles
  76    // Variable binning is used for the xaxis
  77
  78    fVarNames.SetOwner(kTRUE);
  79    fPAvg.SetOwner(kTRUE);
  80    fPProd.SetOwner(kTRUE);
  81    fVarNames.Clear();
  82    fPAvg.Clear();
  83    fPProd.Clear();
  84
  85    fXName="";    // Name of x-axis
  86
  87    TString profileName;
  88
  89    for(int i=0; i<dim; ++i) {
  90       TString vname;
  91       vname+="<";
  92       vname+=i;
  93       vname+=">";
  94       fVarNames.AddAtAndExpand(new TObjString(vname),i);
  95    }
  96
  97    for(int i=0; i<dim; ++i) {
  98       profileName=name; profileName+=i;
  99       fPAvg[i] = new TProfile(profileName.Data(), "", nBins, xBins, "s");
 100       Avg(i)->SetXTitle(name);
 101       Avg(i)->SetDirectory(0);
 102       Avg(i)->Sumw2();
 103    }
 104
 105    for(int i=0; i<dim; ++i) for(int j=0; j<=i; ++j) {
 106       profileName=name; profileName+=i; profileName+=j;
 107       fPProd[Idx(i,j)]=new TProfile(profileName.Data(), "", nBins, xBins, "s");
 108       Prod(i,j)->SetXTitle(name);
 109       Prod(i,j)->SetDirectory(0);
 110       Prod(i,j)->Sumw2();
 111    }
 112
 113    UpdateTitles();
 114 }
 115
 116 void AliFlowMSPHistograms::Init(const int dim, const char *name, const int nBins, const double xmin, const double xmax)
 117 {
 118    // Initialize all profiles after deleting the existing profiles
 119    // Fixed size binning is used
 120
 121    fVarNames.SetOwner(kTRUE);
 122    fPAvg.SetOwner(kTRUE);
 123    fPProd.SetOwner(kTRUE);
 124    fVarNames.Clear();
 125    fPAvg.Clear();
 126    fPProd.Clear();
 127
 128    fXName="";    // Name of x-axis
 129
 130    TString profileName;
 131
 132    for(int i=0; i<dim; ++i) {
 133       TString vname;
 134       vname+="<";
 135       vname+=i;
 136       vname+=">";
 137       fVarNames.AddAtAndExpand(new TObjString(vname),i);
 138    }
 139
 140    for(int i=0; i<dim; ++i) {
 141       profileName=name; profileName+=i;
 142       fPAvg[i] = new TProfile(profileName.Data(), "", nBins, xmin, xmax, "s");
 143       Avg(i)->SetXTitle(name);
 144       Avg(i)->SetDirectory(0);
 145       Avg(i)->Sumw2();
 146    }
 147
 148    for(int i=0; i<dim; ++i) for(int j=0; j<=i; ++j) {
 149       profileName=name; profileName+=i; profileName+=j;
 150       fPProd[Idx(i,j)]=new TProfile(profileName.Data(), "", nBins, xmin, xmax, "s");
 151       Prod(i,j)->SetXTitle(name);
 152       Prod(i,j)->SetDirectory(0);
 153       Prod(i,j)->Sumw2();
 154    }
 155
 156    UpdateTitles();
 157 }
 158
 159 void  AliFlowMSPHistograms::SetVarName(const char *name, const int i)
 160 {
 161    if( i<0 || i>NVars() ) return;
 162
 163    TObjString *old=(TObjString *)(fVarNames.At(i));
 164    if( old ) {
 165       delete old;
 166       fVarNames.RemoveAt(i);
 167    }
 168    fVarNames.AddAtAndExpand(new TObjString(name),i);
 169    UpdateTitles(i);
 170 }
 171
 172 void  AliFlowMSPHistograms::SetVarName(const TString name, const int i)
 173 {
 174    if( i<0 || i>=NVars() ) return;
 175
 176    TObjString *old=(TObjString *)fVarNames.At(i);
 177    if( old ) {
 178       delete old;
 179       fVarNames.RemoveAt(i);
 180    }
 181    fVarNames.AddAtAndExpand(new TObjString(name),i);
 182    UpdateTitles(i);
 183 }
 184
 185 const char * AliFlowMSPHistograms::VarName(const int i) const
 186 {
 187    if( i<0 || i>=NVars() ) return 0;
 188    return ((TObjString *)fVarNames.At(i))->String().Data();
 189 }
 190
 191 void AliFlowMSPHistograms::Fill(const double x, const double *par, const double *wpar)
 192 {
 193    // Fill all profiles in the bin corresponding to x with the averages needed for the covariance matrix
 194    int dim=NVars();
 195    for(int i=0; i<dim; ++i) {
 196       Avg(i)->Fill(x,par[i],wpar[i]);
 197       for(int j=0; j<=i; ++j) {
 198          Prod(i,j)->Fill(x,par[i]*par[j],wpar[i]*wpar[j]);
 199       }
 200    }
 201 }
 202
 203 double AliFlowMSPHistograms::Average(const int bin, const int i)const
 204 {
 205    if( bin<=0 ) { // Average over all Xbins
 206       double stats[6];
 207       Avg(i)->GetStats(stats);
 208       double w=stats[0];            // sum of weights
 209       double y=stats[4];            // sum of w*y
 210       return (w>0?y/w:0);
 211    }
 212    return Avg(i)->GetBinContent(bin);
 213 }
 214
 215 double AliFlowMSPHistograms::Covariance(const int bin, const int i, const int j)const
 216 {
 217    double x,wx,y,wy,xy,wxy;
 218
 219    if( bin<=0 ) {  // Special case: sum over all bins
 220       double stats[6];
 221       Avg(i)->GetStats(stats);
 222       wx=stats[0];                        // sum of weights from TProfile
 223       x=(wx>0?stats[4]/wx:0);             // weighted mean from TProfile
 224
 225       Avg(j)->GetStats(stats);
 226       wy=stats[0];
 227       y=(wy>0?stats[4]/wy:0);
 228
 229       Prod(i,j)->GetStats(stats);
 230       wxy=stats[0];
 231       xy=(wxy>0?stats[4]/wxy:0);
 232    }else{
 233       x=Avg(i)->GetBinContent(bin);       // weighted average per bin from TProfile
 234       wx=Avg(i)->GetBinEntries(bin);      // sum of weights per bin from TProfile
 235       y=Avg(j)->GetBinContent(bin);
 236       wy=Avg(j)->GetBinEntries(bin);
 237       xy=Prod(i,j)->GetBinContent(bin);
 238       wxy=Prod(i,j)->GetBinEntries(bin);
 239    }
 240    return Covariance(x, y, xy, wx, wy, wxy);
 241 }
 242
 243 Long64_t AliFlowMSPHistograms::Merge(TCollection *c)
 244 {
 245    // Merge all MSPHistograms in the collection with this
 246    // The operator+= is used for this
 247    // No compatibility checks are done here, only the number of variables has to be compattible
 248    // The profiles are added, assuming the TProfile::Add function does more checks
 249    if( !c ) return 0;
 250    if( c->IsEmpty() ) return 0;
 251
 252    // Loop over given collection of AliFlowMSPHistograms
 253    AliFlowMSPHistograms *toMerge=0;
 254    Long64_t count=0;
 255    TIter next(c);
 256    while ( (toMerge=dynamic_cast<AliFlowMSPHistograms *>(next())) ) {
 257       *this+=*toMerge;
 258       ++count;
 259    }
 260    return count;
 261 }
 262
 263 AliFlowMSPHistograms &AliFlowMSPHistograms::operator+=(const AliFlowMSPHistograms &x)
 264 {
 265    // Merge this with another AliFlowMSPHistograms object.
 266    // The variable names are neither checked nor changed, only the number of variables is checked
 267    // The profiles are added using the TProfile::Add() function
 268
 269    int nvars=NVars();
 270    if( x.NVars() != nvars ) {
 271       std::cout << "Trying to merge incompattible MSPHistograms. NVars: " << x.NVars() << "!=" << nvars << std::endl;
 272       return *this;
 273    }
 274
 275    for(int i=0; i<nvars; ++i) {  // Loop over averages histograms
 276       Avg(i)->Add(x.Avg(i));
 277       for(int j=0; j<=i; ++j) {  // Loop over products histograms
 278          Prod(i,j)->Add(x.Prod(i,j));
 279       }
 280    }
 281    return *this;
 282 }
 283
 284 AliFlowMSPHistograms &AliFlowMSPHistograms::operator=(const AliFlowMSPHistograms &x)
 285 {
 286    // Assignment, clone all the TProfiles from the other object
 287    // We need to Clone each entry and delete existing ones.
 288    for( int i=0; i<fVarNames.GetEntries(); ++i ) {
 289       TObjString *org=(TObjString *)(x.fVarNames[i]);
 290       delete fVarNames[i];
 291       fVarNames[i]=new TObjString(*org);
 292    }
 293    fVarNames.SetOwner(kTRUE);
 294
 295    for( int i=0; i<fPAvg.GetEntries(); ++i ) {
 296       TProfile *org=(TProfile *)(x.fPAvg[i]);
 297       delete fPAvg[i];
 298       fPAvg[i]=org->Clone();
 299    }
 300    fPAvg.SetOwner(kTRUE);
 301
 302    for( int i=0; i<fPProd.GetEntries(); ++i ) {
 303       TProfile *org=(TProfile *)(x.fPProd[i]);
 304       delete fPProd[i];
 305       fPProd[i]=org->Clone();
 306    }
 307    fPProd.SetOwner(kTRUE);
 308    return *this;
 309 }
 310
 311 TObject *AliFlowMSPHistograms::Clone(const char *n) const
 312 {
 313    AliFlowMSPHistograms *x=new AliFlowMSPHistograms(*this);
 314    if(n)x->SetName(n);
 315    return x;
 316 }
 317
 318 // Private functions...................................................................................
 319
 320 double AliFlowMSPHistograms::Covariance(double x, double y, double xy, double wx, double wy, double wxy)const
 321 {
 322    // Calculate the estimated covariance between two variables
 323    // inputs are the measured averages of x, y and xy and the corresponding sums of weights
 324    // These are taken from TProfile::GetBinContent(ibin) and TProfile::GetBinEntries(ibin)
 325    // See also equations C12-C16 in Ante's thesis
 326
 327    double enumerator=(xy-x*y)*wxy;
 328    double denominator= wx*wy-wxy;
 329
 330    if(TMath::Abs(enumerator)<=1e-50*TMath::Abs(denominator)  ) { // TODO Get smallest Double_t from root ??
 331       return 0;
 332    }
 333    if( TMath::Abs(denominator)<=1e-50*TMath::Abs(enumerator) ) { // protect against infinity
 334       return 1e50;
 335    }
 336    double result=enumerator/denominator;
 337    return result;
 338 }
 339
 340 void AliFlowMSPHistograms::UpdateTitles(int k)
 341 {
 342    // Update the titles of all (k<0) or the profiles for variable k using the variable and parameter names
 343    // Parameter and x variable names should be set before this function is called
 344
 345    int dim=NVars();
 346
 347    for(int i=0; i<dim; ++i) {
 348       if( i==k || k<0 ) {
 349          Avg(i)->SetXTitle(XName());
 350          TString profileTitle="<"; profileTitle+=VarName(i); profileTitle+=">";
 351          Avg(i)->SetYTitle(profileTitle.Data());
 352          profileTitle+=" vs "; profileTitle+=XName();
 353          Avg(i)->SetTitle(profileTitle.Data());
 354       }
 355    }
 356
 357    for(int i=0; i<dim; ++i) for(int j=0; j<=i; ++j) {
 358       if( i==k || j==k || k<0 ) {
 359          Prod(i,j)->SetXTitle(XName());
 360          TString profileTitle="<"; profileTitle+=VarName(i); profileTitle+="*"; profileTitle+=VarName(j); profileTitle+=">";
 361          Prod(i,j)->SetYTitle(profileTitle.Data());
 362          profileTitle+=" vs "; profileTitle+=XName();
 363          Prod(i,j)->SetTitle(profileTitle.Data());
 364       }
 365    }
 366 }
 367