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