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