PWGDQ/dielectron/AliDielectronHF.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, 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
  16 ///////////////////////////////////////////////////////////////////////////
  17 //                Dielectron HF                                  //
  18 //                                                                       //
  19 //                                                                       //
  20 /*
  21 Detailed description
  22
  23
  24 */
  25 //                                                                       //
  26 ///////////////////////////////////////////////////////////////////////////
  27
  28 #include <TVectorD.h>
  29 #include <TH1.h>
  30 #include <TH1F.h>
  31 #include <TH2.h>
  32 #include <TH3.h>
  33 #include <TProfile.h>
  34 #include <TProfile2D.h>
  35 #include <TProfile3D.h>
  36 #include <THnSparse.h>
  37 #include <TAxis.h>
  38 #include <TString.h>
  39 #include <TObjString.h>
  40 #include <TObjArray.h>
  41
  42 #include <AliVParticle.h>
  43 #include <AliLog.h>
  44
  45 #include "AliDielectron.h"
  46 #include "AliDielectronHelper.h"
  47 #include "AliDielectronMC.h"
  48 #include "AliDielectronPair.h"
  49 #include "AliDielectronSignalMC.h"
  50
  51 #include "AliDielectronHistos.h"
  52 #include "AliDielectronHF.h"
  53
  54 ClassImp(AliDielectronHF)
  55
  56 AliDielectronHF::AliDielectronHF() :
  57   TNamed(),
  58   fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
  59   fArrPairType(),
  60   fPairType(kSeOnlyOS),
  61   fSignalsMC(0x0),
  62   fVarCutType(new TBits(kMaxCuts)),
  63   fAxes(kMaxCuts),
  64   fHasMC(kFALSE),
  65   fStepGenerated(kFALSE),
  66   fEventArray(kFALSE),
  67   fRefObj(1)
  68 {
  69   //
  70   // Default Constructor
  71   //
  72   for (Int_t i=0; i<kMaxCuts; ++i){
  73     fVarCuts[i]=0;
  74     //    fVarCutType[i]=0;
  75     fBinType[i]=kStdBin;
  76   }
  77   fAxes.SetOwner(kTRUE);
  78   fRefObj.SetOwner(kTRUE);
  79   fArrPairType.SetOwner(kTRUE);
  80 }
  81
  82 //______________________________________________
  83 AliDielectronHF::AliDielectronHF(const char* name, const char* title) :
  84   TNamed(name, title),
  85   fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
  86   fArrPairType(),
  87   fPairType(kSeOnlyOS),
  88   fSignalsMC(0x0),
  89   fVarCutType(new TBits(kMaxCuts)),
  90   fAxes(kMaxCuts),
  91   fHasMC(kFALSE),
  92   fStepGenerated(kFALSE),
  93   fEventArray(kFALSE),
  94   fRefObj(1)
  95 {
  96   //
  97   // Named Constructor
  98   //
  99   for (Int_t i=0; i<kMaxCuts; ++i){
 100     fVarCuts[i]=0;
 101     //    fVarCutType[i]=0;
 102     fBinType[i]=kStdBin;
 103   }
 104   fAxes.SetOwner(kTRUE);
 105   fRefObj.SetOwner(kTRUE);
 106   fArrPairType.SetOwner(kTRUE);
 107 }
 108
 109 //______________________________________________
 110 AliDielectronHF::~AliDielectronHF()
 111 {
 112   //
 113   // Default Destructor
 114   //
 115   if(fUsedVars)   delete fUsedVars;
 116   if(fVarCutType) delete fVarCutType;
 117   fAxes.Delete();
 118   fRefObj.Delete();
 119   fArrPairType.Delete();
 120 }
 121
 122 //_____________________________________________________________________________
 123 void AliDielectronHF::UserProfile(const char* histClass, UInt_t valTypeP,
 124                                       const TVectorD * const binsX,
 125                                       UInt_t valTypeX, TString option, UInt_t valTypeW)
 126 {
 127   //
 128   // Histogram creation 1D case with arbitraty binning X
 129   // the TVectorD is assumed to be surplus after the creation and will be deleted!!!
 130   //
 131
 132   TH1 *hist=0x0;
 133   if(valTypeP==AliDielectronHistos::kNoProfile)
 134     hist=new TH1F("","",binsX->GetNrows()-1,binsX->GetMatrixArray());
 135   else {
 136     TString opt=""; Double_t pmin=0., pmax=0.;
 137     if(!option.IsNull()) {
 138       TObjArray *arr=option.Tokenize(";");
 139       arr->SetOwner();
 140       opt=((TObjString*)arr->At(0))->GetString();
 141       if(arr->GetEntriesFast()>1) pmin=(((TObjString*)arr->At(1))->GetString()).Atof();
 142       if(arr->GetEntriesFast()>2) pmax=(((TObjString*)arr->At(2))->GetString()).Atof();
 143       delete arr;
 144     }
 145     hist=new TProfile("","",binsX->GetNrows()-1,binsX->GetMatrixArray(),pmin,pmax,opt.Data());
 146     //      printf(" name %s PROFILE options: pmin %.1f pmax %.1f err %s \n",name,((TProfile*)hist)->GetYmin(),((TProfile*)hist)->GetYmax(),((TProfile*)hist)->GetErrorOption() );
 147   }
 148
 149   // store variales in axes
 150   UInt_t valType[4] = {0};
 151   valType[0]=valTypeX;     valType[1]=valTypeP;
 152   AliDielectronHistos::StoreVariables(hist, valType);
 153   hist->SetUniqueID(valTypeW); // store weighting variable
 154
 155   for(Int_t i=0; i<4; i++)   fUsedVars->SetBitNumber(valType[i],kTRUE);
 156   fUsedVars->SetBitNumber(valTypeW,kTRUE);
 157
 158   // adapt the name and title of the histogram in case they are empty
 159   AliDielectronHistos::AdaptNameTitle(hist, histClass);
 160   hist->SetName(Form("HF_%s",hist->GetName()));
 161
 162   fRefObj.AddLast(hist);
 163   delete binsX;
 164 }
 165
 166 //_____________________________________________________________________________
 167 void AliDielectronHF::UserProfile(const char* histClass, UInt_t valTypeP,
 168                                       const TVectorD * const binsX, const TVectorD * const binsY,
 169                                       UInt_t valTypeX, UInt_t valTypeY, TString option, UInt_t valTypeW)
 170 {
 171   //
 172   // Histogram creation 2D case with arbitraty binning X and Y
 173   // the TVectorD is assumed to be surplus after the creation and will be deleted!!!
 174   //
 175
 176   TH1 *hist=0x0;
 177   if(valTypeP==AliDielectronHistos::kNoProfile) {
 178     hist=new TH2F("","",
 179                   binsX->GetNrows()-1,binsX->GetMatrixArray(),
 180                   binsY->GetNrows()-1,binsY->GetMatrixArray());
 181   }
 182   else  {
 183     TString opt=""; Double_t pmin=0., pmax=0.;
 184     if(!option.IsNull()) {
 185       TObjArray *arr=option.Tokenize(";");
 186       arr->SetOwner();
 187       opt=((TObjString*)arr->At(0))->GetString();
 188       if(arr->GetEntriesFast()>1) pmin=(((TObjString*)arr->At(1))->GetString()).Atof();
 189       if(arr->GetEntriesFast()>2) pmax=(((TObjString*)arr->At(2))->GetString()).Atof();
 190       delete arr;
 191     }
 192     hist=new TProfile2D("","",
 193                         binsX->GetNrows()-1,binsX->GetMatrixArray(),
 194                         binsY->GetNrows()-1,binsY->GetMatrixArray());
 195     ((TProfile2D*)hist)->BuildOptions(pmin,pmax,opt.Data());
 196   }
 197
 198   // store variales in axes
 199   UInt_t valType[4] = {0};
 200   valType[0]=valTypeX;     valType[1]=valTypeY; valType[2]=valTypeP;
 201   AliDielectronHistos::StoreVariables(hist, valType);
 202   hist->SetUniqueID(valTypeW); // store weighting variable
 203
 204   for(Int_t i=0; i<4; i++)   fUsedVars->SetBitNumber(valType[i],kTRUE);
 205   fUsedVars->SetBitNumber(valTypeW,kTRUE);
 206
 207   // adapt the name and title of the histogram in case they are empty
 208   AliDielectronHistos::AdaptNameTitle(hist, histClass);
 209   hist->SetName(Form("HF_%s",hist->GetName()));
 210
 211   fRefObj.AddLast(hist);
 212   delete binsX;
 213   delete binsY;
 214 }
 215
 216 //_____________________________________________________________________________
 217 void AliDielectronHF::UserProfile(const char* histClass, UInt_t valTypeP,
 218                                       const TVectorD * const binsX, const TVectorD * const binsY, const TVectorD * const binsZ,
 219                                       UInt_t valTypeX, UInt_t valTypeY, UInt_t valTypeZ, TString option, UInt_t valTypeW)
 220 {
 221   //
 222   // Histogram creation 3D case with arbitraty binning X, Y, Z
 223   // the TVectorD is assumed to be surplus after the creation and will be deleted!!!
 224   //
 225   TH1 *hist=0x0;
 226   if(valTypeP==AliDielectronHistos::kNoProfile) {
 227     hist=new TH3F("","",
 228                   binsX->GetNrows()-1,binsX->GetMatrixArray(),
 229                   binsY->GetNrows()-1,binsY->GetMatrixArray(),
 230                   binsZ->GetNrows()-1,binsZ->GetMatrixArray());
 231   }
 232   else {
 233     TString opt=""; Double_t pmin=0., pmax=0.;
 234     if(!option.IsNull()) {
 235       TObjArray *arr=option.Tokenize(";");
 236       arr->SetOwner();
 237       opt=((TObjString*)arr->At(0))->GetString();
 238       if(arr->GetEntriesFast()>1) pmin=(((TObjString*)arr->At(1))->GetString()).Atof();
 239       if(arr->GetEntriesFast()>2) pmax=(((TObjString*)arr->At(2))->GetString()).Atof();
 240       delete arr;
 241     }
 242     hist=new TProfile3D("","",
 243                         binsX->GetNrows()-1,binsX->GetMatrixArray(),
 244                         binsY->GetNrows()-1,binsY->GetMatrixArray(),
 245                         binsZ->GetNrows()-1,binsZ->GetMatrixArray());
 246     ((TProfile3D*)hist)->BuildOptions(pmin,pmax,opt.Data());
 247   }
 248
 249   // store variales in axes
 250   UInt_t valType[4] = {0};
 251   valType[0]=valTypeX;     valType[1]=valTypeY;     valType[2]=valTypeZ;     valType[3]=valTypeP;
 252   AliDielectronHistos::StoreVariables(hist, valType);
 253   hist->SetUniqueID(valTypeW); // store weighting variable
 254
 255   for(Int_t i=0; i<4; i++)   fUsedVars->SetBitNumber(valType[i],kTRUE);
 256   fUsedVars->SetBitNumber(valTypeW,kTRUE);
 257
 258   // adapt the name and title of the histogram in case they are empty
 259   AliDielectronHistos::AdaptNameTitle(hist, histClass);
 260   hist->SetName(Form("HF_%s",hist->GetName()));
 261
 262   fRefObj.AddLast(hist);
 263   delete binsX;
 264   delete binsY;
 265   delete binsZ;
 266 }
 267
 268 //_____________________________________________________________________________
 269 void AliDielectronHF::UserSparse(const char* histClass, Int_t ndim, TObjArray *limits, UInt_t *vars, UInt_t valTypeW)
 270 {
 271   //
 272   // THnSparse creation with non-linear binning
 273   //
 274
 275   THnSparseF *hist=0;
 276   Int_t bins[ndim];
 277   // get number of bins
 278   for(Int_t idim=0 ;idim<ndim; idim++) {
 279     TVectorD *vec = (TVectorD*) limits->At(idim);
 280     bins[idim]=vec->GetNrows()-1;
 281   }
 282
 283   hist=new THnSparseF("",histClass, ndim, bins, 0x0, 0x0);
 284
 285   // set binning
 286   for(Int_t idim=0 ;idim<ndim; idim++) {
 287     TVectorD *vec = (TVectorD*) limits->At(idim);
 288     hist->SetBinEdges(idim,vec->GetMatrixArray());
 289   }
 290
 291   // store variales in axes
 292   AliDielectronHistos::StoreVariables(hist, vars);
 293   hist->SetUniqueID(valTypeW); // store weighting variable
 294
 295   // store which variables are used
 296   for(Int_t i=0; i<20; i++)   fUsedVars->SetBitNumber(vars[i],kTRUE);
 297   fUsedVars->SetBitNumber(valTypeW,kTRUE);
 298
 299   // adapt the name and title of the histogram in case they are empty
 300   TString name;
 301   for(Int_t iv=0; iv < ndim; iv++) name+=Form("%s_",AliDielectronVarManager::GetValueName(vars[iv]));
 302   name.Resize(name.Length()-1);
 303   hist->SetName(Form("HF_%s",name.Data()));
 304
 305   fRefObj.AddLast(hist);
 306   delete limits;
 307
 308 }
 309
 310 //________________________________________________________________
 311 void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
 312                                      Int_t nbins, Double_t min, Double_t max, Bool_t log, Bool_t leg, EBinType btype)
 313 {
 314   //
 315   // Add a variable to the mixing handler
 316   //
 317
 318   // limit number of variables to kMaxCuts
 319   if (fAxes.GetEntriesFast()>=kMaxCuts) return;
 320
 321   TVectorD *binLimits=0x0;
 322   if (!log) binLimits=AliDielectronHelper::MakeLinBinning(nbins,min,max);
 323   else binLimits=AliDielectronHelper::MakeLogBinning(nbins,min,max);
 324   if (!binLimits) return;
 325
 326   Int_t size=fAxes.GetEntriesFast();
 327   fVarCuts[size]=(UShort_t)type;
 328   //  fVarCutType[size]=leg;
 329   fVarCutType->SetBitNumber(size,leg);
 330   fAxes.Add(binLimits->Clone());
 331   fBinType[size]=btype;
 332   fUsedVars->SetBitNumber(type,kTRUE);
 333 }
 334
 335 //________________________________________________________________
 336 void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
 337                                              const char* binLimitStr, Bool_t leg, EBinType btype)
 338 {
 339   //
 340   // Add a variable to the mixing handler with arbitrary binning
 341   //
 342
 343   // limit number of variables to kMaxCuts
 344   if (fAxes.GetEntriesFast()>=kMaxCuts) return;
 345
 346   TVectorD *binLimits=AliDielectronHelper::MakeArbitraryBinning(binLimitStr);
 347   if (!binLimits) return;
 348
 349   Int_t size=fAxes.GetEntriesFast();
 350   fVarCuts[size]=(UShort_t)type;
 351   //  fVarCutType[size]=leg;
 352   fVarCutType->SetBitNumber(size,leg);
 353   fAxes.Add(binLimits);
 354   fBinType[size]=btype;
 355   fUsedVars->SetBitNumber(type,kTRUE);
 356 }
 357
 358 //________________________________________________________________
 359 void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
 360                                              TVectorD * binLimits, Bool_t leg, EBinType btype)
 361 {
 362   //
 363   // Add a variable to the mixing handler with a vector
 364   // the TVectorD is assumed to be surplus after the creation and will be deleted!!!
 365   //
 366
 367   // limit number of variables to kMaxCuts
 368   if (fAxes.GetEntriesFast()>=kMaxCuts) return;
 369
 370   if (!binLimits) return;
 371
 372   Int_t size=fAxes.GetEntriesFast();
 373   fVarCuts[size]=(UShort_t)type;
 374   //  fVarCutType[size]=leg;
 375   fVarCutType->SetBitNumber(size,leg);
 376   fAxes.Add(binLimits);
 377   fBinType[size]=btype;
 378   fUsedVars->SetBitNumber(type,kTRUE);
 379 }
 380
 381 //______________________________________________
 382 void AliDielectronHF::Fill(Int_t label1, Int_t label2, Int_t nSignal)
 383 {
 384   //
 385   // fill the pure MC part of the container starting from a pair of 2 particles (part1 and part2 are legs)
 386   //
 387   // fill only if we have asked for these steps
 388   if(!fStepGenerated || fEventArray) return;
 389
 390   AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
 391   AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
 392   if(!part1 || !part2) return;
 393
 394   AliDielectronMC* dieMC = AliDielectronMC::Instance();
 395
 396   Int_t mLabel1 = dieMC->GetMothersLabel(label1);    // should work for both ESD and AOD
 397   Int_t mLabel2 = dieMC->GetMothersLabel(label2);
 398
 399   // check the same mother option
 400   AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
 401   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
 402   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
 403
 404   AliDielectronVarManager::SetFillMap(fUsedVars);
 405   // fill the leg variables
 406   Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
 407   Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
 408   AliDielectronVarManager::Fill(part1,valuesLeg1);
 409   AliDielectronVarManager::Fill(part2,valuesLeg2);
 410
 411   // fill the pair and event variables
 412   Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
 413   AliDielectronVarManager::Fill(dieMC->GetMCEvent(), valuesPair);
 414   AliDielectronVarManager::FillVarMCParticle2(part1,part2,valuesPair);
 415
 416   // if pair types are filled, fill mc sources at the end
 417   Int_t istep=0;
 418   if(fPairType!=kMConly) istep=AliDielectron::kEv1PMRot+1;
 419
 420   // only OS at the moment
 421   if(part1->Charge()*part2->Charge()<0) {
 422     Fill(istep+nSignal+fSignalsMC->GetEntries(), valuesPair,  valuesLeg1, valuesLeg2);
 423   }
 424
 425   return;
 426 }
 427 //______________________________________________
 428 void AliDielectronHF::Fill(Int_t pairIndex, const AliDielectronPair *particle)
 429 {
 430   //
 431   // fill histograms for event, pair and daughter cuts and pair types
 432   //
 433
 434   // only OS pairs in case of MC
 435   //////////////////////////////  if(fHasMC && pairIndex!=AliDielectron::kEv1PM) return;
 436
 437   // only selected pair types in case of data
 438   if(!IsPairTypeSelected(pairIndex) || fEventArray) return;
 439
 440   // get event and pair variables
 441   Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
 442   AliDielectronVarManager::SetFillMap(fUsedVars);
 443   AliDielectronVarManager::Fill(particle,valuesPair);
 444
 445   // get leg variables (TODO: do not fill for the moment since leg cuts are not opened)
 446   Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues]={0};
 447   if(fVarCutType->CountBits())  AliDielectronVarManager::Fill(particle->GetFirstDaughterP(),valuesLeg1);
 448   Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues]={0};
 449   if(fVarCutType->CountBits())  AliDielectronVarManager::Fill(particle->GetSecondDaughterP(),valuesLeg2);
 450
 451   // fill
 452
 453   // if pair types are filled, fill mc sources at the end
 454   Int_t istep = 0;
 455   if(fPairType!=kMConly) istep=AliDielectron::kEv1PMRot+1;
 456
 457   // mc source steps (only OS SE pairs)
 458   if(fHasMC && fSignalsMC && pairIndex==AliDielectron::kEv1PM) {
 459     for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 460       if(AliDielectronMC::Instance()->IsMCTruth(particle, (AliDielectronSignalMC*)fSignalsMC->At(i)))
 461         Fill(istep+i, valuesPair,  valuesLeg1, valuesLeg2);
 462     }
 463   }
 464
 465   // all pair types w/o use of mc information
 466   if(fPairType==kMConly) return;
 467
 468   // remove comments
 469   //// select correct step if we are looking at signals too
 470   ////  if(fHasMC && fSignalsMC) pairIndex += ( fSignalsMC->GetEntries() * (fStepGenerated ? 2 : 1) );
 471   Fill(pairIndex, valuesPair,  valuesLeg1, valuesLeg2);
 472
 473   return;
 474 }
 475
 476 //______________________________________________
 477 void AliDielectronHF::Fill(Int_t index, Double_t * const valuesPair, Double_t * const valuesLeg1, Double_t * const valuesLeg2)
 478 {
 479   //
 480   // main fill function using index and values as input
 481   //
 482
 483   TObjArray *histArr = static_cast<TObjArray*>(fArrPairType.At(index));
 484   if(!histArr) return;
 485
 486   Int_t size  = GetNumberOfBins();
 487   // loop over all histograms
 488   for(Int_t ihist=0; ihist<size; ihist++) {
 489
 490     Int_t sizeAdd   = 1;
 491     Bool_t selected = kTRUE;
 492
 493     // loop over all cut variables
 494     Int_t nvars = fAxes.GetEntriesFast();
 495     for(Int_t ivar=0; ivar<nvars; ivar++) {
 496
 497       // get bin limits
 498       TVectorD *bins = static_cast<TVectorD*>(fAxes.At(ivar));
 499       Int_t nbins    = bins->GetNrows()-1;
 500
 501       // bin limits for current ivar bin
 502       Int_t ibin   = (ihist/sizeAdd)%nbins;
 503       Double_t lowEdge = (*bins)[ibin];
 504       Double_t upEdge  = (*bins)[ibin+1];
 505       switch(fBinType[ivar]) {
 506       case kStdBin:     upEdge=(*bins)[ibin+1];     break;
 507       case kBinToMax:   upEdge=(*bins)[nbins];      break;
 508       case kBinFromMin: lowEdge=(*bins)[0];         break;
 509       case kSymBin:     upEdge=(*bins)[nbins-ibin];
 510         if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins]; // to avoid low>up
 511         break;
 512       }
 513
 514       // leg variable
 515       if(fVarCutType->TestBitNumber(ivar)) {
 516         if( (valuesLeg1[fVarCuts[ivar]] < lowEdge || valuesLeg1[fVarCuts[ivar]] >= upEdge) ||
 517             (valuesLeg2[fVarCuts[ivar]] < lowEdge || valuesLeg2[fVarCuts[ivar]] >= upEdge) ) {
 518           selected=kFALSE;
 519           break;
 520         }
 521       }
 522       else { // pair and event variables
 523         if( (valuesPair[fVarCuts[ivar]] < lowEdge || valuesPair[fVarCuts[ivar]] >= upEdge) ) {
 524           selected=kFALSE;
 525           break;
 526         }
 527       }
 528
 529       sizeAdd*=nbins;
 530     } //end of var cut loop
 531
 532     // do not fill the histogram
 533     if(!selected) continue;
 534
 535     // fill the object with Pair and event values
 536     TObjArray *tmp = (TObjArray*) histArr->At(ihist);
 537     TString title = tmp->GetName();
 538     AliDebug(10,title.Data());
 539     for(Int_t i=0; i<tmp->GetEntriesFast(); i++) {
 540       AliDielectronHistos::FillValues(tmp->At(i), valuesPair);
 541     }
 542     //    AliDebug(10,Form("Fill var %d %s value %f in %s \n",fVar,AliDielectronVarManager::GetValueName(fVar),valuesPair[fVar],tmp->GetName()));
 543   } //end of hist loop
 544
 545 }
 546
 547 //______________________________________________
 548 void AliDielectronHF::Init()
 549 {
 550   //
 551   // initialise event buffers
 552   //
 553
 554   // has MC signals
 555   fHasMC=AliDielectronMC::Instance()->HasMC();
 556   Int_t steps = 0;
 557   if(fHasMC) steps=fSignalsMC->GetEntries();
 558   if(fStepGenerated) steps*=2;
 559   if(fEventArray) steps=1;
 560
 561   // init pair type array
 562   fArrPairType.SetName(Form("%s_HF",GetName()));
 563   if( (fHasMC && fPairType==kMConly) || fEventArray) fArrPairType.Expand(steps);
 564   else fArrPairType.Expand(AliDielectron::kEv1PMRot+1+steps);
 565
 566   Int_t size  = GetNumberOfBins();
 567   AliDebug(10,Form("Creating a histo array with size %d \n",size));
 568
 569   Int_t sizeAdd  = 1;
 570
 571   // fill object array with the array of bin cells
 572   TObjArray *histArr = new TObjArray(0);
 573   if(!histArr) return;
 574   histArr->SetOwner(kTRUE);
 575   histArr->Expand(size);
 576
 577   //  printf("fRefObj %p \n",fRefObj);
 578   // array of histograms to each bin cell
 579   for(Int_t ihist=0; ihist<size; ihist++) {
 580     histArr->AddAt(fRefObj.Clone(""), ihist);
 581     //histArr->AddAt(fRefObj.Clone(Form("h%04d",ihist)), ihist);
 582   }
 583
 584   // loop over all cut variables and do the naming according to its bin cell
 585   Int_t nvars = fAxes.GetEntriesFast();
 586   for(Int_t ivar=0; ivar<nvars; ivar++) {
 587
 588     // get bin limits
 589     TVectorD *bins = static_cast<TVectorD*>(fAxes.At(ivar));
 590     Int_t nbins    = bins->GetNrows()-1;
 591
 592
 593     // loop over all bin cells an set unique titles
 594     for(Int_t ihist=0; ihist<size; ihist++) {
 595
 596       // get the lower limit for current ivar bin
 597       Int_t ibin   = (ihist/sizeAdd)%nbins;
 598       Double_t lowEdge = (*bins)[ibin];
 599       Double_t upEdge  = (*bins)[ibin+1];
 600       switch(fBinType[ivar]) {
 601       case kStdBin:     upEdge=(*bins)[ibin+1];     break;
 602       case kBinToMax:   upEdge=(*bins)[nbins];      break;
 603       case kBinFromMin: lowEdge=(*bins)[0];         break;
 604       case kSymBin:     upEdge=(*bins)[nbins-ibin];
 605         if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins]; // to avoid low>up
 606         break;
 607       }
 608
 609       TObjArray *tmp= (TObjArray*) histArr->At(ihist);
 610       TString title = tmp->GetName();
 611       if(!ivar)             title ="";
 612       if( ivar)             title+=":";
 613       if(fVarCutType->TestBitNumber(ivar)) title+="Leg";
 614       title+=AliDielectronVarManager::GetValueName(fVarCuts[ivar]);
 615       title+=Form("#%.2f#%.2f",lowEdge,upEdge);
 616       tmp->SetName(title.Data());
 617       AliDebug(10,title.Data());
 618     } // end: array of bin cell
 619     sizeAdd*=nbins;
 620   } //end: cut loop
 621
 622   // copy array to the selected event,  pair types/ MC sources
 623   Int_t istep=0;
 624
 625   ////////////////// only event array
 626   if(fEventArray) {
 627     // add a deep copy of the array
 628     fArrPairType[istep]=(TObjArray*)histArr->Clone("Event");
 629     ((TObjArray*)fArrPairType[istep])->SetOwner();
 630   }
 631   else {
 632     /////////////// pair types
 633     if(fPairType != kMConly) {
 634       for(istep=0; istep<AliDielectron::kEv1PMRot+1; istep++) {
 635
 636         // pair type should be filled
 637         if(IsPairTypeSelected(istep)) {
 638           // add a deep copy of the array
 639           fArrPairType[istep]=(TObjArray*)histArr->Clone(AliDielectron::PairClassName(istep));
 640           ((TObjArray*)fArrPairType[istep])->SetOwner();
 641         }
 642         else { //empty array
 643           fArrPairType[istep]=new TObjArray(0);
 644           ((TObjArray*)fArrPairType[istep])->SetOwner();
 645           ((TObjArray*)fArrPairType[istep])->SetName(AliDielectron::PairClassName(istep));
 646         }
 647       } //end: loop over pair types
 648     }
 649
 650     // mc sources
 651     if(fHasMC) {
 652       for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 653         TString title = Form("(Signal: %s)",fSignalsMC->At(i)->GetTitle());
 654         fArrPairType[istep+i]=(TObjArray*)histArr->Clone(title.Data());
 655         if(fStepGenerated)  {
 656           title+=" MC truth";
 657           fArrPairType[istep+i+fSignalsMC->GetEntries()]=(TObjArray*)histArr->Clone(title.Data());
 658         }
 659       } // end: loop over sources
 660     } //end: hasMC
 661   } //end: pair type array
 662
 663   // clean up
 664   if(histArr) {
 665     delete histArr;
 666     histArr=0;
 667   }
 668
 669 }
 670
 671 //______________________________________________
 672 Int_t AliDielectronHF::GetNumberOfBins() const
 673 {
 674   //
 675   // return the number of bins this histogram grid has
 676   //
 677   Int_t size=1;
 678   for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i)
 679     size*=((static_cast<TVectorD*>(fAxes.At(i)))->GetNrows()-1);
 680   return size;
 681 }
 682
 683 //______________________________________________
 684 Bool_t AliDielectronHF::IsPairTypeSelected(Int_t itype)
 685 {
 686   //
 687   // check whether a pair type was selected
 688   // TODO: cross check or replace by mixinghandlers processsing
 689
 690   Bool_t selected = kFALSE;
 691
 692   // fill all
 693   if(fPairType==kAll) return kTRUE;
 694
 695   switch(itype) {
 696   case AliDielectron::kEv1PP:
 697   case AliDielectron::kEv1MM:
 698     if(fPairType==kSeAll || fPairType==kSeMeAll || fPairType==kSeReAll )   selected = kTRUE;
 699     break;
 700   case AliDielectron::kEv1PM:
 701     if(fPairType!=kMeOnlyOS)  selected = kTRUE;
 702     break;
 703   case AliDielectron::kEv1PEv2P:
 704   case AliDielectron::kEv1MEv2M:
 705     if(fPairType==kMeAll || fPairType==kSeMeAll)   selected = kTRUE;
 706     break;
 707   case AliDielectron::kEv1PEv2M:
 708     if(fPairType==kMeAll || fPairType==kSeMeAll) selected = kTRUE;
 709     break;
 710   case AliDielectron::kEv1MEv2P:
 711     if(fPairType==kMeAll || fPairType==kSeMeAll || fPairType==kMeOnlyOS || fPairType==kSeMeOnlyOS)   selected = kTRUE;
 712     break;
 713   case AliDielectron::kEv2PP:
 714   case AliDielectron::kEv2PM:
 715   case AliDielectron::kEv2MM:
 716     selected = kFALSE;
 717     break;
 718   case AliDielectron::kEv1PMRot:
 719     if(fPairType==kSeReAll || fPairType==kSeReOnlyOS)   selected = kTRUE;
 720     break;
 721   }
 722
 723   return selected;
 724
 725 }
 726
 727