PWGDQ/dielectron/AliDielectronCF.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 //       Dielectron Correction framework manager                         //
  19 //                                                                       //
  20 /*
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30 */
  31 //                                                                       //
  32 ///////////////////////////////////////////////////////////////////////////
  33
  34 #include <TList.h>
  35 #include <TObjArray.h>
  36 #include <TVectorD.h>
  37 #include <TString.h>
  38 #include <TObjString.h>
  39
  40 #include <AliCFContainer.h>
  41 #include <AliAnalysisFilter.h>
  42 #include <AliAnalysisCuts.h>
  43 #include <AliVParticle.h>
  44 #include <AliLog.h>
  45
  46 #include "AliDielectronCF.h"
  47 #include "AliDielectronMC.h"
  48 #include "AliDielectronPair.h"
  49 #include "AliDielectronSignalMC.h"
  50
  51 ClassImp(AliDielectronCF)
  52
  53 AliDielectronCF::AliDielectronCF() :
  54   TNamed("DielectronCF","DielectronCF"),
  55   fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
  56   fNSteps(0),
  57   fNVars(0),
  58   fVarBinLimits(0x0),
  59   fNVarsLeg(0),
  60   fVarBinLimitsLeg(0x0),
  61   fNCuts(0),
  62   fValues(0x0),
  63   fIsMCTruth(0x0),
  64   fStepForMCtruth(kFALSE),
  65   fStepForNoCutsMCmotherPid(kFALSE),
  66   fStepForAfterAllCuts(kTRUE),
  67   fStepForPreFilter(kFALSE),
  68   fStepsForEachCut(kFALSE),
  69   fStepsForCutsIncreasing(kFALSE),
  70   fStepsForSignal(kTRUE),
  71   fStepsForBackground(kFALSE),
  72   fStepsForMCtruthOnly(kFALSE),
  73   fNStepMasks(0),
  74   fPdgMother(-1),
  75   fSignalsMC(0x0),
  76   fCfContainer(0x0),
  77   fHasMC(kFALSE),
  78   fNAddSteps(0)
  79 {
  80   //
  81   // Default constructor
  82   //
  83   for (Int_t i=0; i<AliDielectronVarManager::kNMaxValues; ++i){
  84     fVariables[i]=0;
  85     fVariablesLeg[i]=0;
  86   }
  87
  88   for (Int_t i=0; i<kNmaxAddSteps; ++i){
  89     fStepMasks[i]=0xFFFFFF;
  90   }
  91 }
  92
  93 //________________________________________________________________
  94 AliDielectronCF::AliDielectronCF(const char* name, const char* title) :
  95   TNamed(name, title),
  96   fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
  97   fNSteps(0),
  98   fNVars(0),
  99   fVarBinLimits(0x0),
 100   fNVarsLeg(0),
 101   fVarBinLimitsLeg(0x0),
 102   fNCuts(0),
 103   fValues(0x0),
 104   fIsMCTruth(0x0),
 105   fStepForMCtruth(kFALSE),
 106   fStepForNoCutsMCmotherPid(kFALSE),
 107   fStepForAfterAllCuts(kTRUE),
 108   fStepForPreFilter(kFALSE),
 109   fStepsForEachCut(kFALSE),
 110   fStepsForCutsIncreasing(kFALSE),
 111   fStepsForSignal(kTRUE),
 112   fStepsForBackground(kFALSE),
 113   fStepsForMCtruthOnly(kFALSE),
 114   fNStepMasks(0),
 115   fPdgMother(-1),
 116   fSignalsMC(0x0),
 117   fCfContainer(0x0),
 118   fHasMC(kFALSE),
 119   fNAddSteps(0)
 120 {
 121   //
 122   // Named constructor
 123   //
 124   for (Int_t i=0; i<AliDielectronVarManager::kNMaxValues; ++i){
 125     fVariables[i]=0;
 126     fVariablesLeg[i]=0;
 127   }
 128
 129   for (Int_t i=0; i<kNmaxAddSteps; ++i){
 130     fStepMasks[i]=0xFFFFFF;
 131   }
 132 }
 133
 134 //________________________________________________________________
 135 AliDielectronCF::~AliDielectronCF()
 136 {
 137   //
 138   // Destructor
 139   //
 140   if (fUsedVars) delete fUsedVars;
 141   if (fValues) delete [] fValues;
 142   if (fIsMCTruth) delete [] fIsMCTruth;
 143   if (fVarBinLimits) delete fVarBinLimits;
 144   if (fVarBinLimitsLeg) delete fVarBinLimitsLeg;
 145 }
 146
 147 //________________________________________________________________
 148 void AliDielectronCF::AddVariable(AliDielectronVarManager::ValueTypes type, Int_t nbins,
 149                                   Double_t min, Double_t max, Bool_t leg, Bool_t log)
 150 {
 151   //
 152   // Add a variable to the CF configuration
 153   // if leg is true it will add the variables of the leg
 154   // if log is true log binning will be created
 155   //
 156
 157   TVectorD *binLimits=0x0;
 158   if (!log) binLimits=MakeLinBinning(nbins,min,max);
 159   else binLimits=MakeLogBinning(nbins,min,max);
 160   AddVariable(type,binLimits,leg);
 161 }
 162
 163 //________________________________________________________________
 164 void AliDielectronCF::AddVariable(AliDielectronVarManager::ValueTypes type, const char* binLimitStr, Bool_t leg/*=kFALSE*/)
 165 {
 166   //
 167   // Add a variable to the CF configuration
 168   // specify arbitrary binning in a string.
 169   // Bin limits need to be separated by a ","
 170   //
 171   TString limits(binLimitStr);
 172   if (limits.IsNull()){
 173     AliError(Form("Bin Limit string is empty, cannot add the variable '%s'",AliDielectronVarManager::GetValueName(type)));
 174     return;
 175   }
 176
 177   TObjArray *arr=limits.Tokenize(",");
 178   Int_t nLimits=arr->GetEntries();
 179   if (nLimits<2){
 180     AliError(Form("Need at leas 2 bin limits, cannot add the variable '%s'",AliDielectronVarManager::GetValueName(type)));
 181     delete arr;
 182     return;
 183   }
 184
 185   TVectorD *binLimits=new TVectorD(nLimits);
 186   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 187     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 188   }
 189
 190   delete arr;
 191   AddVariable(type,binLimits,leg);
 192 }
 193
 194 //________________________________________________________________
 195 void AliDielectronCF::AddVariable(AliDielectronVarManager::ValueTypes type, TVectorD *binLimits, Bool_t leg/*=kFALSE*/)
 196 {
 197   //
 198   // Add variable with the binning given in the TVectorD
 199   //
 200   if (!leg){
 201     if (!fVarBinLimits){
 202       fVarBinLimits=new TObjArray;
 203       fVarBinLimits->SetOwner();
 204     }
 205     fVarBinLimits->Add(binLimits);
 206     fVariables[fNVars]  = (UInt_t)type;
 207     ++fNVars;
 208   } else {
 209     if (!fVarBinLimitsLeg){
 210       fVarBinLimitsLeg=new TObjArray;
 211       fVarBinLimitsLeg->SetOwner();
 212     }
 213     fVarBinLimitsLeg->Add(binLimits);
 214     fVariablesLeg[fNVarsLeg]  = (UInt_t)type;
 215     fUsedVars->SetBitNumber(type,kTRUE);
 216     ++fNVarsLeg;
 217   }
 218 }
 219
 220 //________________________________________________________________
 221 void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter)
 222 {
 223   //
 224   // Initialise container based on the cuts in the analysis filter
 225   //
 226
 227   fNCuts=filter.GetCuts()->GetEntries();
 228
 229   fHasMC=AliDielectronMC::Instance()->HasMC();
 230   fNAddSteps=1;
 231   if (fHasMC){
 232     if (fStepsForSignal && fSignalsMC) fNAddSteps+=fSignalsMC->GetEntries();
 233     if (fStepsForBackground) ++fNAddSteps;
 234     if (fStepsForMCtruthOnly) --fNAddSteps; // No Step for Pair information
 235   } else {
 236     //if
 237     fStepForMCtruth=kFALSE;
 238     fStepForNoCutsMCmotherPid=kFALSE;
 239     fStepsForSignal=kFALSE;
 240     fStepsForBackground=kFALSE;
 241   }
 242   // consitency checks to not duplicate steps
 243   if (fStepsForCutsIncreasing)     fStepForAfterAllCuts=kFALSE;
 244   if (fStepsForEachCut&&fNCuts==1) fStepForAfterAllCuts=kFALSE;
 245
 246   fNSteps=0;
 247   if (fStepForMCtruth && fSignalsMC)           fNSteps+=fSignalsMC->GetEntries();
 248   if (fStepForNoCutsMCmotherPid && fSignalsMC) fNSteps+=fSignalsMC->GetEntries();
 249   if (fStepForAfterAllCuts)                    fNSteps+=fNAddSteps;
 250
 251   if (fStepsForEachCut)         fNSteps+=(fNAddSteps*fNCuts); //one step for each cut + Signal (MC)
 252   if (fStepsForCutsIncreasing)  fNSteps+=(fNAddSteps*fNCuts); //one step for the increasing cuts + Signal (MC)
 253   // e.g. cut1, cut1&cut2, cut1&cut2&cut3, ...
 254
 255   fNSteps+=(fNAddSteps*fNStepMasks);                              // cuts for the additional cut masks
 256
 257   if (fStepForPreFilter) fNSteps+=fNAddSteps; //Add at the end for Prefilter (maxcutmask+1)
 258
 259   // create the container
 260
 261   Int_t *nbins=new Int_t[fNVars+2*fNVarsLeg];
 262   for (Int_t i=0;i<fNVars;++i) {
 263     Int_t nBins=(static_cast<TVectorD*>(fVarBinLimits->At(i)))->GetNrows()-1;
 264     nbins[i]=nBins;
 265   }
 266   for (Int_t i=0;i<fNVarsLeg;++i){
 267     Int_t nBins=(static_cast<TVectorD*>(fVarBinLimitsLeg->At(i)))->GetNrows()-1;
 268     nbins[i+fNVars]=nBins;
 269     nbins[i+fNVars+fNVarsLeg]=nBins;
 270   }
 271
 272   fCfContainer = new AliCFContainer(GetName(), GetTitle(), fNSteps, fNVars+2*fNVarsLeg, nbins);
 273   delete [] nbins;
 274
 275   // initialize the variables and their bin limits
 276   for (Int_t iVar=0; iVar<fNVars; iVar++) {
 277     UInt_t type=fVariables[iVar];
 278     Double_t *binLim = (static_cast<TVectorD*>(fVarBinLimits->At(iVar)))->GetMatrixArray();
 279
 280     fCfContainer->SetBinLimits(iVar, binLim);
 281     fCfContainer->SetVarTitle(iVar, AliDielectronVarManager::GetValueName(type));
 282   }
 283
 284   // initialize the variables and their bin limits for the Legs
 285   for (Int_t iVar=0; iVar<fNVarsLeg; iVar++) {
 286     UInt_t type=fVariablesLeg[iVar];
 287     Double_t *binLim=(static_cast<TVectorD*>(fVarBinLimitsLeg->At(iVar)))->GetMatrixArray();
 288
 289     //Leg1
 290     fCfContainer->SetBinLimits(iVar+fNVars, binLim);
 291     fCfContainer->SetVarTitle(iVar+fNVars, Form("Leg1_%s",AliDielectronVarManager::GetValueName(type)));
 292
 293     //Leg2
 294     fCfContainer->SetBinLimits(iVar+fNVars+fNVarsLeg, binLim);
 295     fCfContainer->SetVarTitle(iVar+fNVars+fNVarsLeg, Form("Leg2_%s",AliDielectronVarManager::GetValueName(type)));
 296   }
 297
 298   // array for storing values
 299   fValues = new Double_t[fNVars+2*fNVarsLeg];
 300
 301   // array for storing MC info
 302   if (fHasMC && fSignalsMC && fSignalsMC->GetEntries()>0) fIsMCTruth=new Bool_t[fSignalsMC->GetEntries()];
 303   //=================//
 304   // Set step titles //
 305   //=================//
 306   Int_t step=0;
 307
 308   //Pure MC truth
 309   if(fStepForMCtruth && fSignalsMC) {
 310     for(Int_t i=0; i<fSignalsMC->GetEntries(); i++)
 311       fCfContainer->SetStepTitle(step++, Form("MC truth (Signal: %s)", fSignalsMC->At(i)->GetTitle()));
 312   }
 313
 314   //before cuts (MC truth)
 315   if (fStepForNoCutsMCmotherPid && fSignalsMC){
 316     for(Int_t i=0; i<fSignalsMC->GetEntries(); i++)
 317       fCfContainer->SetStepTitle(step++,Form("No cuts (Signal: %s)",fSignalsMC->At(i)->GetTitle()));
 318   }
 319
 320   TString cutName;
 321   //Steps for each of the cuts
 322   if (fStepsForEachCut){
 323     for (Int_t iCut=0; iCut<fNCuts;++iCut) {
 324       cutName=filter.GetCuts()->At(iCut)->GetName(); //TODO: User GetTitle???
 325       if (!fStepsForMCtruthOnly) {
 326         fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
 327       }
 328       if (fHasMC){
 329         if (fStepsForSignal && fSignalsMC) {
 330           for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 331             fCfContainer->SetStepTitle(step++, Form("%s (Signal: %s)", cutName.Data(), fSignalsMC->At(i)->GetTitle())); //Step for  the cut with MC truth
 332           }
 333         }
 334         if (fStepsForBackground)
 335           fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
 336       }
 337     }
 338   }
 339
 340   //Steps for increasing cut match
 341   if (fStepsForCutsIncreasing){
 342     cutName=""; //TODO: User GetTitle???
 343     for (Int_t iCut=0; iCut<fNCuts;++iCut) {
 344       if (!cutName.IsNull()) cutName+="&";
 345       cutName+=filter.GetCuts()->At(iCut)->GetName();
 346       if (!fStepsForMCtruthOnly) {
 347         fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
 348       }
 349       if (fHasMC){
 350         if (fStepsForSignal && fSignalsMC)
 351           for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 352             fCfContainer->SetStepTitle(step++, Form("%s (Signal: %s)", cutName.Data(), fSignalsMC->At(i)->GetTitle())); //Step for the cut with MC truth
 353           }
 354           if (fStepsForBackground)
 355             fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
 356       }
 357     }
 358   }
 359
 360   //Steps of user defined cut combinations
 361   for (UInt_t iComb=0; iComb<fNStepMasks; ++iComb){
 362     cutName="";
 363     UInt_t mask=fStepMasks[iComb];
 364     for (Int_t iCut=0; iCut<fNCuts;++iCut) {
 365       if (mask&(1<<iCut)){
 366         if (cutName.IsNull()){
 367           cutName=filter.GetCuts()->At(iCut)->GetName();
 368         }else{
 369           cutName+="&";
 370           cutName+=filter.GetCuts()->At(iCut)->GetName();
 371         }
 372       }
 373     }
 374
 375     fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
 376
 377     if (fHasMC){
 378       if (fStepsForSignal && fSignalsMC)
 379         for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 380           fCfContainer->SetStepTitle(step++, Form("%s (Signal: %s)", cutName.Data(), fSignalsMC->At(i)->GetTitle())); //Step for the cut with MC truth
 381         }
 382       if (fStepsForBackground)
 383         fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
 384     }
 385   }
 386
 387   //After All cuts
 388   if (fStepForAfterAllCuts){
 389     cutName="No pair cuts";
 390     if (filter.GetCuts()->At(0)){
 391       cutName=filter.GetCuts()->At(0)->GetName(); //TODO: User GetTitle???
 392       for (Int_t iCut=1; iCut<fNCuts;++iCut) {
 393         cutName+="&";
 394         cutName+=filter.GetCuts()->At(iCut)->GetName();
 395       }
 396     }
 397     if (!fStepsForMCtruthOnly) {
 398       fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
 399     }
 400     if (fHasMC){
 401       if (fStepsForSignal && fSignalsMC)
 402         for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 403           fCfContainer->SetStepTitle(step++, Form("%s (Signal: %s)", cutName.Data(), fSignalsMC->At(i)->GetTitle())); //Step for the cut with MC truth
 404         }
 405         if (fStepsForBackground)
 406           fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
 407     }
 408   }
 409
 410   //Additional Step for result after PreFilter
 411   if (fStepForPreFilter){
 412     cutName="PreFilter";
 413     fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
 414     if (fHasMC){
 415       if (fStepsForSignal && fSignalsMC)
 416         for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 417           fCfContainer->SetStepTitle(step++, Form("%s (Signal %s)", cutName.Data(), fSignalsMC->At(i)->GetTitle())); //Step for the cut with MC truth
 418         }
 419         if (fStepsForBackground)
 420           fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
 421     }
 422   }
 423
 424
 425
 426   if (step!=fNSteps) {
 427     AliError(Form("Something went wrong in the naming of the steps!!! (%d != %d)",step,fNSteps));
 428   }
 429 }
 430
 431 //________________________________________________________________
 432 void AliDielectronCF::Fill(UInt_t mask, const AliDielectronPair *particle)
 433 {
 434   //
 435   // Fill the containers
 436   //
 437
 438   // Check the MC truths
 439   if(fIsMCTruth) {
 440     for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) fIsMCTruth[i]=kFALSE;
 441   }
 442
 443   //TODO: for the moment don't fill truth information for mixed event paris. No valid MC info is available
 444   //      in the mixing handler
 445   Bool_t isMixedPair=(particle->GetType()>2&&particle->GetType()<10);
 446
 447   Bool_t isBackground = kFALSE;
 448   if(fIsMCTruth && !isMixedPair) {
 449     for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
 450       fIsMCTruth[i] = AliDielectronMC::Instance()->IsMCTruth(particle, (AliDielectronSignalMC*)fSignalsMC->At(i));
 451       isBackground = (isBackground || fIsMCTruth[i]);
 452     }
 453     // background is considered that pair which does not fulfill any of the signals
 454     isBackground = !isBackground;
 455   }
 456
 457   Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
 458   AliDielectronVarManager::SetFillMap(fUsedVars);
 459   AliDielectronVarManager::Fill(particle,valuesPair);
 460
 461   for (Int_t iVar=0; iVar<fNVars; ++iVar){
 462     Int_t var=fVariables[iVar];
 463     fValues[iVar]=valuesPair[var];
 464   }
 465
 466   if (fNVarsLeg>0){
 467     Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues]={0};
 468     AliDielectronVarManager::Fill(particle->GetFirstDaughterP(),valuesLeg1);
 469     Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues]={0};
 470     AliDielectronVarManager::Fill(particle->GetSecondDaughterP(),valuesLeg2);
 471
 472     for (Int_t iVar=0; iVar<fNVarsLeg; ++iVar){
 473       Int_t var=fVariablesLeg[iVar];
 474       fValues[iVar+fNVars]=valuesLeg1[var];
 475       fValues[iVar+fNVars+fNVarsLeg]=valuesLeg2[var];
 476     }
 477   }
 478
 479   UInt_t selectedMask=(1<<fNCuts)-1;
 480
 481   //============//
 482   // Fill steps //
 483   //============//
 484   // Pure MC steps are handled in FillMC
 485   Int_t step=0;
 486   if (fStepForMCtruth && fIsMCTruth) step+=fSignalsMC->GetEntries();
 487
 488   //No cuts (MC truth)
 489   if (fStepForNoCutsMCmotherPid && fIsMCTruth){
 490     for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
 491       if(fIsMCTruth[i]) {
 492         fCfContainer->Fill(fValues,step);
 493       }
 494       ++step;
 495     }
 496   }
 497
 498   //Steps for each of the cuts
 499   if (fStepsForEachCut){
 500     for (Int_t iCut=0; iCut<fNCuts;++iCut) {
 501       UInt_t cutMask=1<<iCut;
 502       if ((mask&cutMask)==cutMask) {
 503         if(!fStepsForMCtruthOnly) {
 504           fCfContainer->Fill(fValues,step);
 505           ++step;
 506         }
 507         if (fHasMC){
 508           if ( fStepsForSignal && fIsMCTruth){
 509             for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
 510               if(fIsMCTruth[i]) {
 511                 fCfContainer->Fill(fValues,step);
 512               }
 513               ++step;
 514             }
 515           }
 516           if ( fStepsForBackground ){
 517             if (isBackground) fCfContainer->Fill(fValues,step);
 518             ++step;
 519           }
 520         }
 521       } else {
 522         step+=fNAddSteps;
 523       }
 524     }
 525   }
 526
 527
 528   //Steps for increasing cut match
 529   if (fStepsForCutsIncreasing&&fNCuts>2){
 530     for (Int_t iCut=0; iCut<fNCuts;++iCut) {
 531       UInt_t cutMask=(1<<(iCut+1))-1;
 532       if ((mask&cutMask)==cutMask) {
 533         if(!fStepsForMCtruthOnly) {
 534           fCfContainer->Fill(fValues,step);
 535           ++step;
 536         }
 537
 538         if (fHasMC){
 539           if ( fStepsForSignal && fIsMCTruth){
 540             for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
 541               if(fIsMCTruth[i]) {
 542                 fCfContainer->Fill(fValues,step);
 543               }
 544               ++step;
 545             }
 546           }
 547           if ( fStepsForBackground ){
 548             if (isBackground) fCfContainer->Fill(fValues,step);
 549             ++step;
 550           }
 551         }
 552       } else {
 553         step+=fNAddSteps;
 554       }
 555     }
 556   }
 557
 558   //Steps of user defined cut combinations
 559   for (UInt_t iComb=0; iComb<fNStepMasks; ++iComb){
 560     UInt_t userMask=fStepMasks[iComb];
 561     if ((mask&userMask)==userMask) {
 562       if(!fStepsForMCtruthOnly) {
 563         fCfContainer->Fill(fValues,step);
 564         ++step;
 565       }
 566       if (fHasMC){
 567         if ( fStepsForSignal && fIsMCTruth){
 568           for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
 569             if(fIsMCTruth[i]) {
 570               fCfContainer->Fill(fValues,step);
 571             }
 572             ++step;
 573           }
 574         }
 575         if ( fStepsForBackground ){
 576           if (isBackground) fCfContainer->Fill(fValues,step);
 577           ++step;
 578         }
 579       }
 580     } else {
 581       step+=fNAddSteps;
 582     }
 583   }
 584
 585   //All cuts
 586   if (fStepForAfterAllCuts){
 587     if (mask == selectedMask){
 588       if(!fStepsForMCtruthOnly) {
 589         fCfContainer->Fill(fValues,step);
 590         ++step;
 591       }
 592
 593       if (fHasMC){
 594         if ( fStepsForSignal && fIsMCTruth){
 595           for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
 596             if(fIsMCTruth[i]) {
 597               fCfContainer->Fill(fValues,step);
 598             }
 599             ++step;
 600           }
 601         }
 602         if ( fStepsForBackground ){
 603           if (isBackground) fCfContainer->Fill(fValues,step);
 604           ++step;
 605         }
 606       }
 607     } else {
 608       step+=fNAddSteps;
 609     }
 610   }
 611
 612   //prefilter
 613   if (fStepForPreFilter) {
 614     if (mask&(1<<fNCuts)) {
 615       if(!fStepsForMCtruthOnly) {
 616         fCfContainer->Fill(fValues,step);
 617         ++step;
 618       }
 619       if (fHasMC){
 620         if ( fStepsForSignal && fIsMCTruth){
 621           for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
 622             if(fIsMCTruth[i]) {
 623               fCfContainer->Fill(fValues,step);
 624             }
 625             ++step;
 626           }
 627         }
 628         if ( fStepsForBackground ){
 629           if (isBackground) fCfContainer->Fill(fValues,step);
 630           ++step;
 631         }
 632       }
 633     }
 634     else {
 635       step+=fNAddSteps;
 636     }
 637   }
 638
 639   if (step!=fNSteps) {
 640     AliError("Something went wrong in the step filling!!!");
 641   }
 642 }
 643
 644 //________________________________________________________________
 645 void AliDielectronCF::FillMC(const TObject *particle)
 646 {
 647   //
 648   // fill MC part of the Container
 649   //
 650   if (!fStepForMCtruth) return;
 651
 652   Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
 653   AliDielectronVarManager::SetFillMap(fUsedVars);
 654   AliDielectronVarManager::Fill(particle,valuesPair);
 655
 656   AliVParticle *d1=0x0;
 657   AliVParticle *d2=0x0;
 658   AliDielectronMC::Instance()->GetDaughters(particle,d1,d2);
 659
 660   //TODO: temporary solution, set manually the pair type to 1: unlikesign SE
 661   valuesPair[AliDielectronVarManager::kPairType]=1;
 662
 663   for (Int_t iVar=0; iVar<fNVars; ++iVar){
 664     Int_t var=fVariables[iVar];
 665     fValues[iVar]=valuesPair[var];
 666   }
 667
 668   if (fNVarsLeg>0){
 669     Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
 670     Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
 671     if (d1->Pt()>d2->Pt()){
 672       AliDielectronVarManager::Fill(d1,valuesLeg1);
 673       AliDielectronVarManager::Fill(d2,valuesLeg2);
 674     } else {
 675       AliDielectronVarManager::Fill(d2,valuesLeg1);
 676       AliDielectronVarManager::Fill(d1,valuesLeg2);
 677     }
 678
 679     for (Int_t iVar=0; iVar<fNVarsLeg; ++iVar){
 680       Int_t var=fVariablesLeg[iVar];
 681       fValues[iVar+fNVars]=valuesLeg1[var];
 682       fValues[iVar+fNVars+fNVarsLeg]=valuesLeg2[var];
 683     }
 684   }
 685
 686   fCfContainer->Fill(fValues,0);
 687 }
 688
 689
 690 //________________________________________________________________
 691 void AliDielectronCF::FillMC(Int_t label1, Int_t label2, Int_t nSignal) {
 692   //
 693   // fill the pure MC part of the container starting from a pair of 2 particles (part1 and part2 are legs)
 694   //
 695   if (!fStepForMCtruth) return;
 696
 697   AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
 698   AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
 699   if(!part1 || !part2) return;
 700
 701   AliDielectronMC* dieMC = AliDielectronMC::Instance();
 702
 703   Int_t mLabel1 = dieMC->GetMothersLabel(label1);    // should work for both ESD and AOD
 704   Int_t mLabel2 = dieMC->GetMothersLabel(label2);
 705   // check the same mother option
 706   AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
 707   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
 708   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
 709
 710   AliDielectronVarManager::SetFillMap(fUsedVars);
 711   // fill the leg variables
 712   if (fNVarsLeg>0){
 713     Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
 714     Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
 715     if (part1->Pt()>part2->Pt()){
 716       AliDielectronVarManager::Fill(part1,valuesLeg1);
 717       AliDielectronVarManager::Fill(part2,valuesLeg2);
 718     } else {
 719       AliDielectronVarManager::Fill(part2,valuesLeg1);
 720       AliDielectronVarManager::Fill(part1,valuesLeg2);
 721     }
 722
 723     for (Int_t iVar=0; iVar<fNVarsLeg; ++iVar){
 724       Int_t var=fVariablesLeg[iVar];
 725       fValues[iVar+fNVars]=valuesLeg1[var];
 726       fValues[iVar+fNVars+fNVarsLeg]=valuesLeg2[var];
 727     }
 728   }
 729
 730   Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
 731   AliDielectronVarManager::Fill(dieMC->GetMCEvent(), valuesPair);
 732   AliDielectronVarManager::FillVarMCParticle2(part1,part2,valuesPair);
 733
 734   if(part1->Charge()*part2->Charge()<0)
 735     valuesPair[AliDielectronVarManager::kPairType]=1;
 736   else if(part1->Charge()>0)
 737     valuesPair[AliDielectronVarManager::kPairType]=0;
 738   else
 739     valuesPair[AliDielectronVarManager::kPairType]=2; // if one of the two particles is neutral, the pair will go here
 740
 741   for(Int_t iVar=0; iVar<fNVars; ++iVar){
 742     Int_t var=fVariables[iVar];
 743     fValues[iVar]=valuesPair[var];
 744   }
 745
 746   fCfContainer->Fill(fValues,nSignal);
 747 }
 748
 749
 750 //_____________________________________________________________________________
 751 TVectorD* AliDielectronCF::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const
 752 {
 753   //
 754   // Make logarithmic binning
 755   // the user has to delete the array afterwards!!!
 756   //
 757
 758   //check limits
 759   if (xmin<1e-20 || xmax<1e-20){
 760     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
 761     return MakeLinBinning(nbinsX, xmin, xmax);
 762   }
 763   if (xmax<xmin){
 764     Double_t tmp=xmin;
 765     xmin=xmax;
 766     xmax=tmp;
 767   }
 768   TVectorD *binLim=new TVectorD(nbinsX+1);
 769   Double_t first=xmin;
 770   Double_t last=xmax;
 771   Double_t expMax=TMath::Log(last/first);
 772   for (Int_t i=0; i<nbinsX+1; ++i){
 773     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
 774   }
 775   return binLim;
 776 }
 777
 778 //_____________________________________________________________________________
 779 TVectorD* AliDielectronCF::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const
 780 {
 781   //
 782   // Make logarithmic binning
 783   // the user has to delete the array afterwards!!!
 784   //
 785   if (xmax<xmin){
 786     Double_t tmp=xmin;
 787     xmin=xmax;
 788     xmax=tmp;
 789   }
 790   TVectorD *binLim=new TVectorD(nbinsX+1);
 791   Double_t first=xmin;
 792   Double_t last=xmax;
 793   Double_t binWidth=(last-first)/nbinsX;
 794   for (Int_t i=0; i<nbinsX+1; ++i){
 795     (*binLim)[i]=first+binWidth*(Double_t)i;
 796   }
 797   return binLim;
 798 }
 799