[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronHF.cxx

/*************************************************************************
* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

///////////////////////////////////////////////////////////////////////////
//                Dielectron HF                                  //
//                                                                       //
//                                                                       //
/*
Detailed description


*/
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#include <TVectorD.h>
#include <TH1.h>
#include <TAxis.h>
#include <AliVParticle.h>

#include <AliLog.h>

#include "AliDielectron.h"
#include "AliDielectronHelper.h"
#include "AliDielectronMC.h"
#include "AliDielectronPair.h"
#include "AliDielectronSignalMC.h"

#include "AliDielectronHF.h"

ClassImp(AliDielectronHF)

AliDielectronHF::AliDielectronHF() :
  TNamed(),
  fArrPairType(0x0),
  fPairType(kOSonly),
  fSignalsMC(0x0),
  fAxes(kMaxCuts),
  fVarBinLimits(0x0),
  fVar(0),
  fHasMC(kFALSE),
  fStepGenerated(kFALSE)
{
  //
  // Default Constructor
  //
  for (Int_t i=0; i<kMaxCuts; ++i){
    fVarCuts[i]=0;
    fVarCutType[i]=0;
    fBinType[i]=kStdBin;
  }
  fAxes.SetOwner(kTRUE);
}

//______________________________________________
AliDielectronHF::AliDielectronHF(const char* name, const char* title) :
  TNamed(name, title),
  fArrPairType(0x0),
  fPairType(kOSonly),
  fSignalsMC(0x0),
  fAxes(kMaxCuts),
  fVarBinLimits(0x0),
  fVar(0),
  fHasMC(kFALSE),
  fStepGenerated(kFALSE)
{
  //
  // Named Constructor
  //
  for (Int_t i=0; i<kMaxCuts; ++i){
    fVarCuts[i]=0;
    fVarCutType[i]=0;
    fBinType[i]=kStdBin;
  }
  fAxes.SetOwner(kTRUE);
}

//______________________________________________
AliDielectronHF::~AliDielectronHF()
{
  //
  // Default Destructor
  //
  fAxes.Delete();
}

//________________________________________________________________
void AliDielectronHF::SetVariable(AliDielectronVarManager::ValueTypes type,
					  Int_t nbins, Double_t min, Double_t max, Bool_t log)
{
  //
  // Set main variable for the histos
  //

  fVarBinLimits=0x0;
  if (!log) fVarBinLimits=AliDielectronHelper::MakeLinBinning(nbins,min,max);
  else fVarBinLimits=AliDielectronHelper::MakeLogBinning(nbins,min,max);
  if (!fVarBinLimits) return ;
 
  fVar=(UShort_t)type;

}

//________________________________________________________________
void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
                                             Int_t nbins, Double_t min, Double_t max, Bool_t log, Bool_t leg, EBinType btype)
{
  //
  // Add a variable to the mixing handler
  //

  // limit number of variables to kMaxCuts
  if (fAxes.GetEntriesFast()>=kMaxCuts) return;
  
  TVectorD *binLimits=0x0;
  if (!log) binLimits=AliDielectronHelper::MakeLinBinning(nbins,min,max);
  else binLimits=AliDielectronHelper::MakeLogBinning(nbins,min,max);
  if (!binLimits) return;

  Int_t size=fAxes.GetEntriesFast();
  fVarCuts[size]=(UShort_t)type;
  fVarCutType[size]=leg;
  fAxes.Add(binLimits->Clone());
  fBinType[size]=btype;
}

//________________________________________________________________
void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
                                             const char* binLimitStr, Bool_t leg, EBinType btype)
{
  //
  // Add a variable to the mixing handler with arbitrary binning
  //

  // limit number of variables to kMaxCuts
  if (fAxes.GetEntriesFast()>=kMaxCuts) return;
  
  TVectorD *binLimits=AliDielectronHelper::MakeArbitraryBinning(binLimitStr);
  if (!binLimits) return;
  
  Int_t size=fAxes.GetEntriesFast();
  fVarCuts[size]=(UShort_t)type;
  fVarCutType[size]=leg;
  fAxes.Add(binLimits);
  fBinType[size]=btype;
}

//________________________________________________________________
void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
                                             TVectorD * binLimits, Bool_t leg, EBinType btype)
{
  //
  // Add a variable to the mixing handler with a vector
  // the TVectorD is assumed to be surplus after the creation and will be deleted!!!
  //

  // limit number of variables to kMaxCuts
  if (fAxes.GetEntriesFast()>=kMaxCuts) return;
  
  if (!binLimits) return;
  
  Int_t size=fAxes.GetEntriesFast();
  fVarCuts[size]=(UShort_t)type;
  fVarCutType[size]=leg;
  fAxes.Add(binLimits);
  fBinType[size]=btype;
}

//______________________________________________
void AliDielectronHF::Fill(Int_t label1, Int_t label2, Int_t nSignal) 
{
  //
  // fill the pure MC part of the container starting from a pair of 2 particles (part1 and part2 are legs)
  //
  // fill only if we have asked for these steps
  if(!fStepGenerated) return;

  AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
  AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);

  AliDielectronMC* dieMC = AliDielectronMC::Instance();
  
  Int_t mLabel1 = dieMC->GetMothersLabel(label1);    // should work for both ESD and AOD
  Int_t mLabel2 = dieMC->GetMothersLabel(label2);

  // check the same mother option
  AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
  if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
  if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
    
  // fill the leg variables
  Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
  Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
  AliDielectronVarManager::Fill(part1,valuesLeg1);
  AliDielectronVarManager::Fill(part2,valuesLeg2);
    
  // fill the pair and event variables
  Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
  AliDielectronVarManager::Fill(dieMC->GetMCEvent(), valuesPair);
  AliDielectronVarManager::FillVarMCParticle2(part1,part2,valuesPair);

  if(part1->Charge()*part2->Charge()<0) {
    //valuesPair[AliDielectronVarManager::kPairType]=1;
//     for(Int_t i=0; i<fAxes.GetEntriesFast(); i++) {
//       printf("[D]: %s %f %f %f \n",
// 	     AliDielectronVarManager::GetValueName(fVarCuts[i]),
// 	     valuesLeg1[fVarCuts[i]], valuesLeg2[fVarCuts[i]], valuesPair[fVarCuts[i]]); 
//     }
    Fill(nSignal, valuesPair,  valuesLeg1, valuesLeg2);
  }
  // on OS at the moment
  // else if(part1->Charge()>0)
  //   valuesPair[AliDielectronVarManager::kPairType]=0;
  // else
  //   valuesPair[AliDielectronVarManager::kPairType]=2; // if one of the two particles is neutral, the pair will go here

  return;
}
//______________________________________________
void AliDielectronHF::Fill(Int_t pairIndex, const AliDielectronPair *particle)
{
  //
  // fill histograms for event, pair and daughter cuts and pair types
  //
  
  // only OS pairs in case of MC
  if(fHasMC && pairIndex!=AliDielectron::kEv1PM) return;

  // only selected pair types in case of data
  if(!IsPairTypeSelected(pairIndex)) return;

  // get event and pair variables
  Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
  AliDielectronVarManager::Fill(particle,valuesPair);

  // get leg variables
  Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues]={0};
  AliDielectronVarManager::Fill(particle->GetFirstDaughter(),valuesLeg1);
  Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues]={0};
  AliDielectronVarManager::Fill(particle->GetSecondDaughter(),valuesLeg2);

  // fill
  if(!fHasMC) { Fill(pairIndex, valuesPair,  valuesLeg1, valuesLeg2); }
  if(fHasMC && fSignalsMC) {
    for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
      if(AliDielectronMC::Instance()->IsMCTruth(particle, (AliDielectronSignalMC*)fSignalsMC->At(i))) 
	Fill(i, valuesPair,  valuesLeg1, valuesLeg2);
    }
  }
  
  return;
  
}

//______________________________________________
void AliDielectronHF::Fill(Int_t Index, Double_t * const valuesPair, Double_t * const valuesLeg1, Double_t * const valuesLeg2)
{
  //
  // main fill function using index and values as input
  //

  TObjArray *histArr = static_cast<TObjArray*>(fArrPairType.At(Index));
  if(!histArr) return;

  Int_t size  = GetNumberOfBins();
  // loop over all histograms
  for(Int_t ihist=0; ihist<size; ihist++) {

    Int_t sizeAdd   = 1; 
    Bool_t selected = kTRUE;
    
    // loop over all cut variables
    Int_t nvars = fAxes.GetEntriesFast();
    for(Int_t ivar=0; ivar<nvars; ivar++) {
      
      // get bin limits
      TVectorD *bins = static_cast<TVectorD*>(fAxes.At(ivar));
      Int_t nbins    = bins->GetNrows()-1;

      // bin limits for current ivar bin
      Int_t ibin   = (ihist/sizeAdd)%nbins; 
      Double_t lowEdge = (*bins)[ibin];
      Double_t upEdge  = (*bins)[ibin+1];
      switch(fBinType[ivar]) {
      case kStdBin:     upEdge=(*bins)[ibin+1];     break;
      case kBinToMax:   upEdge=(*bins)[nbins];      break;
      case kBinFromMin: lowEdge=(*bins)[0];         break;
      case kSymBin:     upEdge=(*bins)[nbins-ibin]; 
	if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins]; // to avoid low>up
	break;
      }

      // leg variable
      if(fVarCutType[ivar]) {
	if( (valuesLeg1[fVarCuts[ivar]] < lowEdge || valuesLeg1[fVarCuts[ivar]] >= upEdge) ||
	    (valuesLeg2[fVarCuts[ivar]] < lowEdge || valuesLeg2[fVarCuts[ivar]] >= upEdge) ) {
	  selected=kFALSE;
	  break;
	}
      } 
      else { // pair and event variables
	if( (valuesPair[fVarCuts[ivar]] < lowEdge || valuesPair[fVarCuts[ivar]] >= upEdge) ) {
	  selected=kFALSE;
	  break;
	}
      }
      
      sizeAdd*=nbins;
    } //end of var cut loop
    
    // do not fill the histogram
    if(!selected) continue;
    
    // fill the histogram
    TH1F *tmp=static_cast<TH1F*>(histArr->At(ihist));
    TString title = tmp->GetTitle();
    AliDebug(10,title.Data());
    
    AliDebug(10,Form("Fill var %d %s value %f in %s \n",fVar,AliDielectronVarManager::GetValueName(fVar),valuesPair[fVar],tmp->GetName()));
    tmp->Fill(valuesPair[fVar]);
  
  } //end of hist loop
  
}

//______________________________________________
void AliDielectronHF::Init()
{
  //
  // initialise event buffers
  //

  // has MC signals
  fHasMC=AliDielectronMC::Instance()->HasMC();
  Int_t steps = 0;
  if(fHasMC) steps=fSignalsMC->GetEntries();
  if(fStepGenerated) steps*=2; 

  // init pair type array
  fArrPairType.SetName(Form("%s_HF",GetName()));
  if(fHasMC) fArrPairType.Expand(steps);
  else fArrPairType.Expand(AliDielectron::kEv1PMRot+1);

  TH1F *hist  = 0x0;
  Int_t size  = GetNumberOfBins();
  AliDebug(10,Form("Creating a histo array with size %d \n",size));

  Int_t sizeAdd  = 1; 

  // fill object array with the histograms
  TObjArray *histArr = new TObjArray();
  histArr->Expand(size);

  for(Int_t ihist=0; ihist<size; ihist++) {
    hist=new TH1F(Form("histarr%04d",ihist),"",fVarBinLimits->GetNrows()-1,fVarBinLimits->GetMatrixArray());
    hist->SetXTitle(Form("%s",AliDielectronVarManager::GetValueName(fVar)));
    hist->SetYTitle("#pairs");
    histArr->AddAt(hist,ihist);
  }

  // loop over all cut variables
  Int_t nvars = fAxes.GetEntriesFast();
  for(Int_t ivar=0; ivar<nvars; ivar++) {
    
    // get bin limits
    TVectorD *bins = static_cast<TVectorD*>(fAxes.At(ivar));
    Int_t nbins    = bins->GetNrows()-1;


    // loop over all histograms an set unique titles
    for(Int_t ihist=0; ihist<size; ihist++) {

      // get the lower limit for current ivar bin
      Int_t ibin   = (ihist/sizeAdd)%nbins; 
      Double_t lowEdge = (*bins)[ibin];
      Double_t upEdge  = (*bins)[ibin+1];
      switch(fBinType[ivar]) {
      case kStdBin:     upEdge=(*bins)[ibin+1];     break;
      case kBinToMax:   upEdge=(*bins)[nbins];      break;
      case kBinFromMin: lowEdge=(*bins)[0];         break;
      case kSymBin:     upEdge=(*bins)[nbins-ibin];
	if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins]; // to avoid low>up
	break;
      }

      TH1F *tmp=static_cast<TH1F*>(histArr->At(ihist));
      TString title = tmp->GetTitle();
      if(ivar!=0)           title+=":";
      if(fVarCutType[ivar]) title+="Leg";
      title+=AliDielectronVarManager::GetValueName(fVarCuts[ivar]);
      title+=Form("#%.2f#%.2f",lowEdge,upEdge);
      tmp->SetNameTitle(title.Data(),title.Data());
      AliDebug(10,title.Data());
    }
    sizeAdd*=nbins;
  }

  // copy array to the selected pair types/ MC sources
  if(fHasMC) {
    for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
	TString title = Form("(Signal: %s)",fSignalsMC->At(i)->GetTitle());
	fArrPairType[i]=(TObjArray*)histArr->Clone(title.Data());
	if(fStepGenerated)  {
	  title+=" MC truth";
	  fArrPairType[i+fSignalsMC->GetEntries()]=(TObjArray*)histArr->Clone(title.Data());
	}
      }
   }
   else {
    for(Int_t i=0; i<AliDielectron::kEv1PMRot+1; i++) {
      if(IsPairTypeSelected(i)) fArrPairType[i]=(TObjArray*)histArr->Clone(Form("%s",AliDielectron::PairClassName(i)));
      else fArrPairType[i]=0x0;
    }
  }
  
  // clean up
  if(histArr) {
    delete histArr;
    histArr=0;
  }

}

//______________________________________________
Int_t AliDielectronHF::GetNumberOfBins() const
{
  //
  // return the number of bins this mixing handler has
  //
  Int_t size=1;
  for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i)
    size*=((static_cast<TVectorD*>(fAxes.At(i)))->GetNrows()-1);
  return size;
}

//______________________________________________
Bool_t AliDielectronHF::IsPairTypeSelected(Int_t itype)
{
  //
  // check whether a pair type was selected
  //
  
  Bool_t selected = kFALSE;

  // fill all
  if(fPairType==kAll) selected = kTRUE;

  switch(itype) {
  case AliDielectron::kEv1PP: 
  case AliDielectron::kEv1MM:
    if(fPairType==kOSandLS)   selected = kTRUE;
    break;
  case AliDielectron::kEv1PM: selected = kTRUE;
    break;
  case AliDielectron::kEv1PEv2P:
  case AliDielectron::kEv1MEv2P:
  case AliDielectron::kEv1PEv2M:
  case AliDielectron::kEv1MEv2M:
    if(fPairType==kOSandMIX)  selected = kTRUE;
    break;
  case AliDielectron::kEv2PP:
  case AliDielectron::kEv2PM: 
  case AliDielectron::kEv2MM:
    selected = kFALSE;
    break;
  case AliDielectron::kEv1PMRot:
    if(fPairType==kOSandROT)  selected = kTRUE;
    break;
  }
  
  return selected;

}
Commit	Line	Data
5e2cf960	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 <TAxis.h>
443a091c	31	#include <AliVParticle.h>
5e2cf960	32
	33	#include <AliLog.h>
	34
	35	#include "AliDielectron.h"
	36	#include "AliDielectronHelper.h"
443a091c	37	#include "AliDielectronMC.h"
	38	#include "AliDielectronPair.h"
	39	#include "AliDielectronSignalMC.h"
5e2cf960	40
	41	#include "AliDielectronHF.h"
	42
	43	ClassImp(AliDielectronHF)
	44
	45	AliDielectronHF::AliDielectronHF() :
	46	TNamed(),
	47	fArrPairType(0x0),
	48	fPairType(kOSonly),
443a091c	49	fSignalsMC(0x0),
5e2cf960	50	fAxes(kMaxCuts),
5e2cf960	51	fVarBinLimits(0x0),
443a091c	52	fVar(0),
4d7704c5	53	fHasMC(kFALSE),
4d7704c5	54	fStepGenerated(kFALSE)
5e2cf960	55	{
	56	//
	57	// Default Constructor
	58	//
	59	for (Int_t i=0; i<kMaxCuts; ++i){
	60	fVarCuts[i]=0;
	61	fVarCutType[i]=0;
	62	fBinType[i]=kStdBin;
	63	}
	64	fAxes.SetOwner(kTRUE);
	65	}
	66
	67	//______________________________________________
	68	AliDielectronHF::AliDielectronHF(const char* name, const char* title) :
	69	TNamed(name, title),
	70	fArrPairType(0x0),
	71	fPairType(kOSonly),
443a091c	72	fSignalsMC(0x0),
5e2cf960	73	fAxes(kMaxCuts),
5e2cf960	74	fVarBinLimits(0x0),
443a091c	75	fVar(0),
4d7704c5	76	fHasMC(kFALSE),
4d7704c5	77	fStepGenerated(kFALSE)
5e2cf960	78	{
	79	//
	80	// Named Constructor
	81	//
	82	for (Int_t i=0; i<kMaxCuts; ++i){
	83	fVarCuts[i]=0;
	84	fVarCutType[i]=0;
	85	fBinType[i]=kStdBin;
	86	}
	87	fAxes.SetOwner(kTRUE);
	88	}
	89
	90	//______________________________________________
	91	AliDielectronHF::~AliDielectronHF()
	92	{
	93	//
	94	// Default Destructor
	95	//
	96	fAxes.Delete();
	97	}
	98
	99	//________________________________________________________________
	100	void AliDielectronHF::SetVariable(AliDielectronVarManager::ValueTypes type,
	101	Int_t nbins, Double_t min, Double_t max, Bool_t log)
	102	{
	103	//
	104	// Set main variable for the histos
	105	//
	106
	107	fVarBinLimits=0x0;
	108	if (!log) fVarBinLimits=AliDielectronHelper::MakeLinBinning(nbins,min,max);
	109	else fVarBinLimits=AliDielectronHelper::MakeLogBinning(nbins,min,max);
	110	if (!fVarBinLimits) return ;
	111
	112	fVar=(UShort_t)type;
	113
	114	}
	115
	116	//________________________________________________________________
	117	void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
	118	Int_t nbins, Double_t min, Double_t max, Bool_t log, Bool_t leg, EBinType btype)
	119	{
	120	//
	121	// Add a variable to the mixing handler
	122	//
	123
	124	// limit number of variables to kMaxCuts
	125	if (fAxes.GetEntriesFast()>=kMaxCuts) return;
	126
	127	TVectorD *binLimits=0x0;
	128	if (!log) binLimits=AliDielectronHelper::MakeLinBinning(nbins,min,max);
	129	else binLimits=AliDielectronHelper::MakeLogBinning(nbins,min,max);
	130	if (!binLimits) return;
	131
	132	Int_t size=fAxes.GetEntriesFast();
	133	fVarCuts[size]=(UShort_t)type;
	134	fVarCutType[size]=leg;
	135	fAxes.Add(binLimits->Clone());
	136	fBinType[size]=btype;
	137	}
	138
	139	//________________________________________________________________
	140	void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
	141	const char* binLimitStr, Bool_t leg, EBinType btype)
142	{
143	//
144	// Add a variable to the mixing handler with arbitrary binning
145	//
146
147	// limit number of variables to kMaxCuts
148	if (fAxes.GetEntriesFast()>=kMaxCuts) return;
149
150	TVectorD *binLimits=AliDielectronHelper::MakeArbitraryBinning(binLimitStr);
151	if (!binLimits) return;
152
153	Int_t size=fAxes.GetEntriesFast();
154	fVarCuts[size]=(UShort_t)type;
155	fVarCutType[size]=leg;
156	fAxes.Add(binLimits);
157	fBinType[size]=btype;
158	}
159
160	//________________________________________________________________
161	void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
162	TVectorD * binLimits, Bool_t leg, EBinType btype)
163	{
164	//
165	// Add a variable to the mixing handler with a vector
166	// the TVectorD is assumed to be surplus after the creation and will be deleted!!!
167	//
168
169	// limit number of variables to kMaxCuts
170	if (fAxes.GetEntriesFast()>=kMaxCuts) return;
171
172	if (!binLimits) return;
173
174	Int_t size=fAxes.GetEntriesFast();
175	fVarCuts[size]=(UShort_t)type;
176	fVarCutType[size]=leg;
177	fAxes.Add(binLimits);
178	fBinType[size]=btype;
179	}
180
181	//______________________________________________
443a091c	182	void AliDielectronHF::Fill(Int_t label1, Int_t label2, Int_t nSignal)
5e2cf960	183	{
5e2cf960	184	//
443a091c	185	// fill the pure MC part of the container starting from a pair of 2 particles (part1 and part2 are legs)
5e2cf960	186	//
4d7704c5	187	// fill only if we have asked for these steps
4d7704c5	188	if(!fStepGenerated) return;
5e2cf960	189
443a091c	190	AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
443a091c	191	AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
5e2cf960	192
443a091c	193	AliDielectronMC* dieMC = AliDielectronMC::Instance();
	194
	195	Int_t mLabel1 = dieMC->GetMothersLabel(label1); // should work for both ESD and AOD
	196	Int_t mLabel2 = dieMC->GetMothersLabel(label2);
5e2cf960	197
443a091c	198	// check the same mother option
	199	AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
	200	if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
	201	if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
	202
	203	// fill the leg variables
	204	Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
	205	Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
	206	AliDielectronVarManager::Fill(part1,valuesLeg1);
	207	AliDielectronVarManager::Fill(part2,valuesLeg2);
	208
	209	// fill the pair and event variables
	210	Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
	211	AliDielectronVarManager::Fill(dieMC->GetMCEvent(), valuesPair);
	212	AliDielectronVarManager::FillVarMCParticle2(part1,part2,valuesPair);
	213
	214	if(part1->Charge()*part2->Charge()<0) {
88204efa	215	//valuesPair[AliDielectronVarManager::kPairType]=1;
	216	// for(Int_t i=0; i<fAxes.GetEntriesFast(); i++) {
	217	// printf("[D]: %s %f %f %f \n",
	218	// AliDielectronVarManager::GetValueName(fVarCuts[i]),
	219	// valuesLeg1[fVarCuts[i]], valuesLeg2[fVarCuts[i]], valuesPair[fVarCuts[i]]);
	220	// }
443a091c	221	Fill(nSignal, valuesPair, valuesLeg1, valuesLeg2);
	222	}
	223	// on OS at the moment
	224	// else if(part1->Charge()>0)
	225	// valuesPair[AliDielectronVarManager::kPairType]=0;
	226	// else
	227	// valuesPair[AliDielectronVarManager::kPairType]=2; // if one of the two particles is neutral, the pair will go here
	228
	229	return;
	230	}
	231	//______________________________________________
	232	void AliDielectronHF::Fill(Int_t pairIndex, const AliDielectronPair *particle)
	233	{
	234	//
	235	// fill histograms for event, pair and daughter cuts and pair types
	236	//
5e2cf960	237
443a091c	238	// only OS pairs in case of MC
	239	if(fHasMC && pairIndex!=AliDielectron::kEv1PM) return;
	240
	241	// only selected pair types in case of data
	242	if(!IsPairTypeSelected(pairIndex)) return;
5e2cf960	243
	244	// get event and pair variables
	245	Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
	246	AliDielectronVarManager::Fill(particle,valuesPair);
	247
	248	// get leg variables
	249	Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues]={0};
	250	AliDielectronVarManager::Fill(particle->GetFirstDaughter(),valuesLeg1);
	251	Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues]={0};
	252	AliDielectronVarManager::Fill(particle->GetSecondDaughter(),valuesLeg2);
	253
443a091c	254	// fill
	255	if(!fHasMC) { Fill(pairIndex, valuesPair, valuesLeg1, valuesLeg2); }
	256	if(fHasMC && fSignalsMC) {
	257	for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
	258	if(AliDielectronMC::Instance()->IsMCTruth(particle, (AliDielectronSignalMC*)fSignalsMC->At(i)))
	259	Fill(i, valuesPair, valuesLeg1, valuesLeg2);
	260	}
	261	}
	262
	263	return;
	264
	265	}
	266
	267	//______________________________________________
	268	void AliDielectronHF::Fill(Int_t Index, Double_t * const valuesPair, Double_t * const valuesLeg1, Double_t * const valuesLeg2)
	269	{
	270	//
	271	// main fill function using index and values as input
	272	//
	273
	274	TObjArray histArr = static_cast<TObjArray>(fArrPairType.At(Index));
	275	if(!histArr) return;
	276
	277	Int_t size = GetNumberOfBins();
5e2cf960	278	// loop over all histograms
	279	for(Int_t ihist=0; ihist<size; ihist++) {
	280
	281	Int_t sizeAdd = 1;
	282	Bool_t selected = kTRUE;
	283
	284	// loop over all cut variables
	285	Int_t nvars = fAxes.GetEntriesFast();
	286	for(Int_t ivar=0; ivar<nvars; ivar++) {
	287
	288	// get bin limits
	289	TVectorD bins = static_cast<TVectorD>(fAxes.At(ivar));
	290	Int_t nbins = bins->GetNrows()-1;
	291
	292	// bin limits for current ivar bin
	293	Int_t ibin = (ihist/sizeAdd)%nbins;
	294	Double_t lowEdge = (*bins)[ibin];
	295	Double_t upEdge = (*bins)[ibin+1];
	296	switch(fBinType[ivar]) {
	297	case kStdBin: upEdge=(*bins)[ibin+1]; break;
	298	case kBinToMax: upEdge=(*bins)[nbins]; break;
	299	case kBinFromMin: lowEdge=(*bins)[0]; break;
	300	case kSymBin: upEdge=(*bins)[nbins-ibin];
	301	if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins]; // to avoid low>up
	302	break;
	303	}
443a091c	304
5e2cf960	305	// leg variable
443a091c	306	if(fVarCutType[ivar]) {
5e2cf960	307	if( (valuesLeg1[fVarCuts[ivar]] < lowEdge \|\| valuesLeg1[fVarCuts[ivar]] >= upEdge) \|\|
	308	(valuesLeg2[fVarCuts[ivar]] < lowEdge \|\| valuesLeg2[fVarCuts[ivar]] >= upEdge) ) {
	309	selected=kFALSE;
	310	break;
	311	}
	312	}
	313	else { // pair and event variables
	314	if( (valuesPair[fVarCuts[ivar]] < lowEdge \|\| valuesPair[fVarCuts[ivar]] >= upEdge) ) {
	315	selected=kFALSE;
	316	break;
	317	}
	318	}
	319
	320	sizeAdd*=nbins;
	321	} //end of var cut loop
	322
	323	// do not fill the histogram
	324	if(!selected) continue;
	325
	326	// fill the histogram
	327	TH1F tmp=static_cast<TH1F>(histArr->At(ihist));
	328	TString title = tmp->GetTitle();
	329	AliDebug(10,title.Data());
	330
443a091c	331	AliDebug(10,Form("Fill var %d %s value %f in %s \n",fVar,AliDielectronVarManager::GetValueName(fVar),valuesPair[fVar],tmp->GetName()));
5e2cf960	332	tmp->Fill(valuesPair[fVar]);
	333
	334	} //end of hist loop
	335
	336	}
	337
	338	//______________________________________________
	339	void AliDielectronHF::Init()
	340	{
	341	//
	342	// initialise event buffers
	343	//
	344
443a091c	345	// has MC signals
443a091c	346	fHasMC=AliDielectronMC::Instance()->HasMC();
4d7704c5	347	Int_t steps = 0;
	348	if(fHasMC) steps=fSignalsMC->GetEntries();
	349	if(fStepGenerated) steps*=2;
443a091c	350
5e2cf960	351	// init pair type array
5e2cf960	352	fArrPairType.SetName(Form("%s_HF",GetName()));
4d7704c5	353	if(fHasMC) fArrPairType.Expand(steps);
443a091c	354	else fArrPairType.Expand(AliDielectron::kEv1PMRot+1);
5e2cf960	355
	356	TH1F *hist = 0x0;
	357	Int_t size = GetNumberOfBins();
	358	AliDebug(10,Form("Creating a histo array with size %d \n",size));
	359
	360	Int_t sizeAdd = 1;
	361
	362	// fill object array with the histograms
	363	TObjArray *histArr = new TObjArray();
	364	histArr->Expand(size);
	365
	366	for(Int_t ihist=0; ihist<size; ihist++) {
	367	hist=new TH1F(Form("histarr%04d",ihist),"",fVarBinLimits->GetNrows()-1,fVarBinLimits->GetMatrixArray());
	368	hist->SetXTitle(Form("%s",AliDielectronVarManager::GetValueName(fVar)));
	369	hist->SetYTitle("#pairs");
	370	histArr->AddAt(hist,ihist);
	371	}
	372
	373	// loop over all cut variables
	374	Int_t nvars = fAxes.GetEntriesFast();
	375	for(Int_t ivar=0; ivar<nvars; ivar++) {
	376
	377	// get bin limits
	378	TVectorD bins = static_cast<TVectorD>(fAxes.At(ivar));
	379	Int_t nbins = bins->GetNrows()-1;
	380
	381
	382	// loop over all histograms an set unique titles
	383	for(Int_t ihist=0; ihist<size; ihist++) {
	384
	385	// get the lower limit for current ivar bin
	386	Int_t ibin = (ihist/sizeAdd)%nbins;
	387	Double_t lowEdge = (*bins)[ibin];
	388	Double_t upEdge = (*bins)[ibin+1];
	389	switch(fBinType[ivar]) {
	390	case kStdBin: upEdge=(*bins)[ibin+1]; break;
	391	case kBinToMax: upEdge=(*bins)[nbins]; break;
	392	case kBinFromMin: lowEdge=(*bins)[0]; break;
	393	case kSymBin: upEdge=(*bins)[nbins-ibin];
	394	if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins]; // to avoid low>up
	395	break;
	396	}
	397
	398	TH1F tmp=static_cast<TH1F>(histArr->At(ihist));
	399	TString title = tmp->GetTitle();
	400	if(ivar!=0) title+=":";
	401	if(fVarCutType[ivar]) title+="Leg";
	402	title+=AliDielectronVarManager::GetValueName(fVarCuts[ivar]);
	403	title+=Form("#%.2f#%.2f",lowEdge,upEdge);
	404	tmp->SetNameTitle(title.Data(),title.Data());
	405	AliDebug(10,title.Data());
	406	}
	407	sizeAdd*=nbins;
	408	}
	409
443a091c	410	// copy array to the selected pair types/ MC sources
	411	if(fHasMC) {
	412	for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
4d7704c5	413	TString title = Form("(Signal: %s)",fSignalsMC->At(i)->GetTitle());
	414	fArrPairType[i]=(TObjArray*)histArr->Clone(title.Data());
	415	if(fStepGenerated) {
	416	title+=" MC truth";
	417	fArrPairType[i+fSignalsMC->GetEntries()]=(TObjArray*)histArr->Clone(title.Data());
	418	}
	419	}
	420	}
	421	else {
443a091c	422	for(Int_t i=0; i<AliDielectron::kEv1PMRot+1; i++) {
	423	if(IsPairTypeSelected(i)) fArrPairType[i]=(TObjArray*)histArr->Clone(Form("%s",AliDielectron::PairClassName(i)));
	424	else fArrPairType[i]=0x0;
	425	}
5e2cf960	426	}
	427
	428	// clean up
	429	if(histArr) {
	430	delete histArr;
	431	histArr=0;
	432	}
	433
	434	}
	435
	436	//______________________________________________
	437	Int_t AliDielectronHF::GetNumberOfBins() const
	438	{
	439	//
	440	// return the number of bins this mixing handler has
	441	//
	442	Int_t size=1;
	443	for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i)
	444	size=((static_cast<TVectorD>(fAxes.At(i)))->GetNrows()-1);
	445	return size;
	446	}
	447
	448	//______________________________________________
	449	Bool_t AliDielectronHF::IsPairTypeSelected(Int_t itype)
	450	{
	451	//
	452	// check whether a pair type was selected
	453	//
	454
	455	Bool_t selected = kFALSE;
	456
	457	// fill all
	458	if(fPairType==kAll) selected = kTRUE;
	459
	460	switch(itype) {
	461	case AliDielectron::kEv1PP:
	462	case AliDielectron::kEv1MM:
	463	if(fPairType==kOSandLS) selected = kTRUE;
	464	break;
	465	case AliDielectron::kEv1PM: selected = kTRUE;
	466	break;
	467	case AliDielectron::kEv1PEv2P:
	468	case AliDielectron::kEv1MEv2P:
	469	case AliDielectron::kEv1PEv2M:
	470	case AliDielectron::kEv1MEv2M:
	471	if(fPairType==kOSandMIX) selected = kTRUE;
	472	break;
	473	case AliDielectron::kEv2PP:
	474	case AliDielectron::kEv2PM:
	475	case AliDielectron::kEv2MM:
	476	selected = kFALSE;
	477	break;
	478	case AliDielectron::kEv1PMRot:
	479	if(fPairType==kOSandROT) selected = kTRUE;
	480	break;
	481	}
	482
	483	return selected;
	484
	485	}
	486
	487