[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronHFhelper.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 Histogram framework helper                           //
//                                                                       //
/*

A helper class to extract objects(histograms and/or profiles) from
a AliDielctronHF array of objects.


How to use it:

  AliDielectronHFhelper *hf = new AliDielectronHFhelper("path/to/the/output/file.root", "ConfigName");
  // print the structure
  hf->Print();

  //apply some cuts and print them
  hf->SetRangeUser("cut1name",cutmin1,cutmax1);
  hf->SetRangeUser(AliDielectronVarManager::kPt,ptmin,ptmax);
  hf->PrintCuts();

  // collect 1-,2- or 3-dim histograms or profiles with error option (default:"")
  TObjArray *arrHists = hf->CollectHistos(AliDielectronVarManager::kM);
  TObjArray *arrProfs = hf->CollectProfiles("",AliDielectronVarManager::kM,AliDielectronVarManager::kPt);

  // then you are left with an array of histograms for all pair types or MC signals

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

#include <TObjArray.h>
#include <TKey.h>
#include <TList.h>
#include <TClass.h>
#include <TObject.h>
#include <TFile.h>
#include <TString.h>
#include <TObjString.h>
#include <TVectorD.h>
#include <TMath.h>
#include <TH1.h>

#include <AliLog.h>

#include "AliDielectron.h"
#include "AliDielectronHFhelper.h"
#include "AliDielectronHF.h"

ClassImp(AliDielectronHFhelper)

//const char* AliDielectronHFhelper::fCutVars[AliDielectronHFhelper::kMaxCuts] = {"","","","","","","","","",""};

//________________________________________________________________
AliDielectronHFhelper::AliDielectronHFhelper(const char* filename, const char* container) :
  TNamed(),
  fMainArr(0x0),
  fCutVars(0x0),
  fCutLowLimits(0),
  fCutUpLimits(0)
{
  //
  // get HF container(s) from file 'filename'
  //
  SetHFArray(filename, container);

}

//________________________________________________________________
AliDielectronHFhelper::~AliDielectronHFhelper()
{
  //
  // dtor
  //
  if(fMainArr) delete fMainArr;
  if(fCutVars)     delete fCutVars;

}

//________________________________________________________________
void AliDielectronHFhelper::SetHFArray(const char* filename, const char* container)
{
  //
  // get HF container from file
  //

  TFile *f = TFile::Open(filename);

  TList *l=f->GetListOfKeys();
  TIter nextKey(l);
  TKey *k=0x0;
  while ( (k=static_cast<TKey*>(nextKey())) ){

    TObject *o=k->ReadObj();
    if (o->IsA()==TList::Class()){

      TList *tlist=(TList*)o;

      TIter next(tlist);
      TObject *obj=0x0;
      while ((obj = next())) {
	TString objname(obj->GetName());

	if( objname.Contains(Form("%s_HF",container)) && obj->IsA()==TObjArray::Class()) {
	  fMainArr = new TObjArray( *(dynamic_cast<TObjArray*>(obj)) );
	  fMainArr->SetOwner();
	  //	  fMainArr->Print();
	  return;
	}
      }
    }
  }

}
//________________________________________________________________
void AliDielectronHFhelper::SetRangeUser(const char *varname, Double_t min, Double_t max, Bool_t leg)
{
  //
  // Set range from variable name
  //

  Int_t size=fCutLowLimits.GetNrows();

  // check if cut is already set
  for(Int_t icut=0; icut<size; icut++) {
    TString cutName = fCutVars->At(icut)->GetName();
    if(!cutName.CompareTo(Form("%s%s",(leg?"Leg":""),varname))) {
      UnsetRangeUser(varname,leg);
      SetRangeUser(varname, min, max, leg);
      return;
    }
  }

  if(size>=kMaxCuts) return;

  // arrays
  if(!fCutVars) {
    fCutVars = new TObjArray();
    fCutVars->SetOwner();
  }
  fCutLowLimits.ResizeTo(size+1);
  fCutUpLimits.ResizeTo(size+1);

  // fill
  TObjString *str = new TObjString(Form("%s%s",(leg?"Leg":""),varname));
  fCutVars->Add(str);
  fCutLowLimits(size) = min;
  fCutUpLimits(size)  = max;
  AliWarning(Form(" %s [%.2f,%.2f]",fCutVars->At(size)->GetName(),fCutLowLimits(size),fCutUpLimits(size)));
}

//________________________________________________________________
void AliDielectronHFhelper::SetRangeUser(AliDielectronVarManager::ValueTypes type, Double_t min, Double_t max, Bool_t leg)
{
  //
  // Set range from AliDielectronVarManager
  //
  SetRangeUser(AliDielectronVarManager::GetValueName(type), min, max, leg);
}

//________________________________________________________________
void AliDielectronHFhelper::UnsetRangeUser(const char *varname, Bool_t leg)
{
  //
  // unset range from variable name
  //

  Int_t size=fCutLowLimits.GetNrows();
  TVectorD newlow;
  TVectorD newup;

  // find cut and build new vectors w/o it
  Int_t ientries = 0;
  for(Int_t icut=0; icut<size; icut++) {

    TString cutName = fCutVars->At(icut)->GetName();
    if(!cutName.CompareTo(Form("%s%s",(leg?"Leg":""),varname))) {
      fCutVars->AddAt(0x0,icut);
      continue;
    }

    // fill new vectors
    newlow.ResizeTo(ientries+1);
    newup.ResizeTo(ientries+1);
    newlow(ientries) = fCutLowLimits(icut);
    newup(ientries)  = fCutUpLimits(icut);

    ientries++;
  }

  // adapt new arrays/vectors
  fCutVars->Compress();

  fCutLowLimits.ResizeTo(ientries);
  fCutUpLimits.ResizeTo(ientries);
  for(Int_t icut=0; icut<ientries; icut++) {
    fCutLowLimits(icut) = newlow(icut);
    fCutUpLimits(icut)  = newup(icut);
  }
  // PrintCuts();

}

//________________________________________________________________
void AliDielectronHFhelper::UnsetRangeUser(AliDielectronVarManager::ValueTypes type, Bool_t leg)
{
  //
  // Unset range from AliDielectronVarManager
  //
  UnsetRangeUser(AliDielectronVarManager::GetValueName(type), leg);
}

//________________________________________________________________
TObjArray* AliDielectronHFhelper::CollectProfiles(TString option,
						  AliDielectronVarManager::ValueTypes varx,
						  AliDielectronVarManager::ValueTypes vary,
						  AliDielectronVarManager::ValueTypes varz,
						  AliDielectronVarManager::ValueTypes vart)
{
  //
  // collect 1-3 dimensional TProfiles for all kind of pair types or sources
  //

  // reconstruct the histogram/profile name resp. key
  Int_t dim    = 0;
  if(varx < AliDielectronVarManager::kNMaxValues) dim++;
  if(vary < AliDielectronVarManager::kNMaxValues) dim++;
  if(varz < AliDielectronVarManager::kNMaxValues) dim++;
  if(vart < AliDielectronVarManager::kNMaxValues) dim++;
  Bool_t bPairClass=0;
  if( varx < AliDielectronVarManager::kPairMax ||
      vary < AliDielectronVarManager::kPairMax ||
      varz < AliDielectronVarManager::kPairMax ||
      vart < AliDielectronVarManager::kPairMax  ) bPairClass=kTRUE;
  Bool_t bwght= (option.Contains("weight",TString::kIgnoreCase));   // weighted histogram
  Bool_t bprf = !(option.Contains("hist",TString::kIgnoreCase));     // tprofile
  Bool_t bStdOpt=kTRUE;
  if(option.Contains("S",TString::kIgnoreCase) || option.Contains("rms",TString::kIgnoreCase)) bStdOpt=kFALSE;

  TString key = "";
  if(bwght) dim--;
  if(bprf)  dim--;
  switch(dim) {
  case 3:
    key+=Form("%s_",AliDielectronVarManager::GetValueName(varx));
    key+=Form("%s_",AliDielectronVarManager::GetValueName(vary));
    key+=Form("%s",AliDielectronVarManager::GetValueName(varz));
    if(bprf) key+=Form("-%s%s",AliDielectronVarManager::GetValueName(vart),(bStdOpt ? "avg" : "rms"));
    break;
  case 2:
    key+=Form("%s_",AliDielectronVarManager::GetValueName(varx));
    key+=Form("%s",AliDielectronVarManager::GetValueName(vary));
    if(bprf) key+=Form("-%s%s",AliDielectronVarManager::GetValueName(varz),(bStdOpt ? "avg" : "rms"));
    break;
  case 1:
    key+=Form("%s",AliDielectronVarManager::GetValueName(varx));
    if(bprf) key+=Form("-%s%s",AliDielectronVarManager::GetValueName(vary),(bStdOpt ? "avg" : "rms"));
    break;
  }
  // to differentiate btw. leg and pair histos
  if(bPairClass) key.Prepend("p");
  // prepend HF
  key.Prepend("HF_");

  // add the weighted part to the key
  if(bwght) {
    if(bprf) dim++;
    switch(dim) {
    case 3:   key+=Form("-wght%s",AliDielectronVarManager::GetValueName(vart)); break;
    case 2:   key+=Form("-wght%s",AliDielectronVarManager::GetValueName(varz)); break;
    case 1:   key+=Form("-wght%s",AliDielectronVarManager::GetValueName(vary)); break;
    case 4:   AliError(Form(" NO weighted 3D profiles supported by the framework"));
    }
  }

  //printf("--------> KEY: %s \n",key.Data());
  // init array of histograms of interest
  TObjArray *collection = new TObjArray(fMainArr->GetEntriesFast());
  collection->SetOwner(kTRUE);

  TObjArray *cloneArr = new TObjArray( *(dynamic_cast<TObjArray*>(fMainArr)) );
  if(!cloneArr) return 0x0;
  cloneArr->SetOwner(kTRUE);

  // loop over all pair types or sources
  for(Int_t i=0; i<cloneArr->GetEntriesFast(); i++) {
    if(!cloneArr->At(i))                             continue;
    if(!((TObjArray*)cloneArr->At(i))->GetEntries()) continue;

    // Get signal/source array with its histograms of interest
    AliDebug(1,Form(" Looking into step %s selected",cloneArr->At(i)->GetName()));
    TObjArray *arr = FindObjects((TObjArray*)cloneArr->At(i));
    if(arr) {

      // find requested histogram
      collection->AddAt(arr->FindObject(key.Data()), i);
      if(!collection->At(i)) { AliError(Form("Object %s does not exist",key.Data())); return 0x0; }

      // modify the signal/source name if needed (MConly)
      TString stepName = cloneArr->At(i)->GetName();
      stepName.ReplaceAll("(","");
      stepName.ReplaceAll(")","");
      stepName.ReplaceAll(": ","");
      stepName.ReplaceAll(" ","_");
      stepName.ReplaceAll("Signal","");
      ((TH1*)collection->At(i))->SetName(Form("%s_%s",key.Data(),stepName.Data()));
    }
    else
      AliError(Form("Step %d not found",i));

  }

  //clean up
  //delete cloneArr;
  //cloneArr=0;

  return collection;
}

//________________________________________________________________
TObjArray* AliDielectronHFhelper::FindObjects(TObjArray *stepArr)
{
  // rename DoCuts, return values is a tobjarray
  // apply cuts and exclude objects from the array for merging (CUT selection)
  //

  // debug
  // TString title    = stepArr->At(0)->GetTitle();
  // TObjArray* vars  = title.Tokenize(":");
  // AliDebug(1,Form(" number of cuts/vars: %d/%d",fCutLowLimits.GetNrows(),vars->GetEntriesFast()));
  if(!stepArr) { AliError("step is empty"); return 0x0;}

  // check for missing cuts
  CheckCuts(stepArr);

  // loop over all cuts
  for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++) {

    Bool_t bFndBin      = kFALSE; // bin with exact limits found
    const char *cutvar  = fCutVars->At(icut)->GetName(); // cut variable
    Double_t min        = fCutLowLimits(icut);           // lower limit
    Double_t max        = fCutUpLimits(icut);            // upper limit
    AliDebug(5,Form(" Cut %d: %s [%.2f,%.2f]",icut,cutvar,min,max));

    // loop over the full grid of the signalArray (pair,source)
    for(Int_t i=0; i<stepArr->GetEntriesFast(); i++) {

      // continue if already empty
      if(!stepArr->At(i)) continue;

      // collect bins from the name
      TString title    = stepArr->At(i)->GetName();
      if(title.IsNull()) continue;
      AliDebug(5,Form(" %03d object name: %s",i,title.Data()));

      // loop over all variables
      TObjArray *vars  = title.Tokenize(":");
      for(Int_t ivar=0; ivar<vars->GetEntriesFast(); ivar++) {
	TString binvar = vars->At(ivar)->GetName();
	AliDebug(10,Form(" --> %d check bin var %s",ivar,binvar.Data()));

	// check cuts set by the user, and compare to the bin limits
	if(binvar.Contains(cutvar)) {
	  // bin limits
	  TObjArray *limits = binvar.Tokenize("#");
	  if(!limits) continue;
	  Double_t binmin = atof(limits->At(1)->GetName()); // lower bin limit
	  Double_t binmax = atof(limits->At(2)->GetName()); // upper bin limit
	  AliDebug(10,Form(" bin %s var %s [%.2f,%.2f]",binvar.Data(),limits->At(0)->GetName(),binmin,binmax));
	  delete limits;

	  // cut and remove objects from the array
	  if(binmin < min || binmax < min || binmin > max || binmax > max ) {
	    AliDebug(10,Form(" removed, out of range! lower bin,cut: %.2f,%.2f or upper bin,cut: %.2f>%.2f",binmin,min,binmax,max));
	    stepArr->AddAt(0x0,i);
	  }
	  if(bFndBin && !(binmin == min && binmax == max)) {
	    stepArr->AddAt(0x0,i);
	    AliDebug(10,Form(" removed, within range! lower bin,cut: %.2f,%.2f or upper bin,cut: %.2f,%.2f",binmin,min,binmax,max));
	  }

	  // did we found a bin with exact the cut ranges
	  // this can happen only once per variable
	  if(binmin==min && binmax==max) bFndBin=kTRUE;

	}

      }
      // clean up
      if(vars)   delete vars;

    } //end loop: pair step

  } //end: loop cuts

  // compress the array by removing all empty entries
  stepArr->Compress();
  AliDebug(1,Form(" Compression: %d arrays left",stepArr->GetEntriesFast()));

  // merge left objects
  TObjArray* hist = Merge(stepArr);
  if(hist) AliDebug(1,Form(" final array has %d objects",hist->GetEntries()));
  return hist;
}

//________________________________________________________________
TObjArray* AliDielectronHFhelper::Merge(TObjArray *arr)
{
  //
  // merge left objects into a single one (LAST step)
  //

  if(arr->GetEntriesFast()<1) { AliError(" No more objects left!"); return 0x0; }

  TObjArray *final=(TObjArray*) arr->At(0)->Clone();
  if(!final) return 0x0;
  final->SetOwner();

  TList listH;
  TString listHargs;
  listHargs.Form("((TCollection*)0x%lx)", (ULong_t)&listH);
  Int_t error = 0;

  //  final->Reset("CE");
  //  final->SetTitle(""); //TODO: change in future
  for(Int_t i=1; i<arr->GetEntriesFast(); i++) {
    listH.Add(arr->At(i));
    //   printf("%d: ent %.0f \n",i,((TH1*)((TObjArray*)arr->At(i))->At(0))->GetEntries());
    //    final->Add((TH1*)arr->At(i));
  }
  //  arr->Clear();

  final->Execute("Merge", listHargs.Data(), &error);
  // returns a tobjarray of histograms/profiles
  return final;
}

//________________________________________________________________
void AliDielectronHFhelper::CheckCuts(TObjArray *arr)
{
  //
  // Compare binning and cut variables. Add necessary cuts (full range, no exclusion)
  //

  // build array with bin variable, minimum and maximum bin values
  TString titleFIRST       = arr->First()->GetName();
  TString titleLAST        = arr->Last()->GetName();
  TObjArray* binvarsF  = titleFIRST.Tokenize(":#");
  TObjArray* binvarsL  = titleLAST.Tokenize(":#");
  Double_t binmin[kMaxCuts]= {0.0};
  Double_t binmax[kMaxCuts]= {0.0};
  if(!binvarsF) { delete binvarsL; return; }
  if(!binvarsL) { delete binvarsF; return; }
  for(Int_t ivar=0; ivar<binvarsF->GetEntriesFast(); ivar++) {

    TString elementF=binvarsF->At(ivar)->GetName();
    TString elementL=binvarsL->At(ivar)->GetName();
    AliDebug(1,Form(" binvar %d: %s,%s",ivar,elementF.Data(),elementL.Data()));

    switch(ivar%3) {
    case 0: continue; break;
    case 1: binmin[(int)ivar/3]=atof(elementF.Data()); break;
    case 2: binmax[(int)ivar/3]=atof(elementL.Data()); break;
    }

    binvarsF->AddAt(0x0,ivar);
  }
  binvarsF->Compress();

  // loop over all vars and cuts, check for missing stuff
  for(Int_t ivar=0; ivar<binvarsF->GetEntriesFast(); ivar++) {

    TString binvar=binvarsF->At(ivar)->GetName();
    Bool_t selected=kFALSE;

    AliDebug(1,Form(" check cuts %d %s [%.2f,%.2f]",ivar,binvar.Data(),binmin[ivar],binmax[ivar]));
    // loop over all cuts and check for missing stuff
    for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++) {
      if(binvar.Contains(fCutVars->At(icut)->GetName())) { selected=kTRUE; break; }
      //      else break;
    }

    // add missing cut with max limits
    if(!selected) {
      AliWarning(Form(" Bin variable %s not covered. Add cut!",binvar.Data()));
      Bool_t leg = binvar.BeginsWith("Leg");
      if(leg) binvar.Remove(0,3);
      SetRangeUser(binvar.Data(),binmin[ivar],binmax[ivar], leg);
    }

  }

  // clean up
  delete binvarsF;
  delete binvarsL;
}

//________________________________________________________________
void AliDielectronHFhelper::Print(const Option_t* /*option*/) const
{

  //
  // Print out object contents
  //
  AliInfo(Form(" Container:               %s",fMainArr->GetName()));

  // pairtypes, steps and sources
  Int_t stepLast=0;
  AliInfo(Form(" Number of filled steps:  %d",fMainArr->GetEntries()));
  for(Int_t istep=0; istep<fMainArr->GetEntriesFast(); istep++) {
    if(!fMainArr->At(istep))                             continue;
    if(!((TObjArray*)fMainArr->At(istep))->GetEntries()) continue;
    AliInfo(Form(" step %d: %s",istep,fMainArr->At(istep)->GetName()));
    stepLast=istep;
  }

  AliInfo(Form(" Number of objects:    %d",
	       ((TObjArray*) ((TObjArray*)fMainArr->At(stepLast)) ->First())->GetEntriesFast()));

  TString title       = ((TObjArray*)fMainArr->At(stepLast))->First()->GetName();
  TObjArray* binvars  = title.Tokenize(":");
  AliInfo(Form(" Number of variables:     %d",binvars->GetEntriesFast()));
  delete binvars;

  /* REACTIVATE
  // check for missing cuts
  CheckCuts(((TObjArray*)fMainArr->At(stepLast)));
  // loop over all cuts
  for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++) {
    const char *cutvar  = fCutVars->At(icut)->GetName(); // cut variable
    Double_t min        = fCutLowLimits(icut);           // lower limit
    Double_t max        = fCutUpLimits(icut);            // upper limit
    AliInfo(Form(" variable %d: %s [%.2f,%.2f]",icut,cutvar,min,max));
  }
  */
  TObjArray* binvars2  = title.Tokenize(":#");
  for(Int_t ivar=0; ivar<binvars2->GetEntriesFast(); ivar++) {
    if(ivar%3) continue;
    AliInfo(Form(" variable %.0f: %s",((Double_t)ivar)/3+1,binvars2->At(ivar)->GetName()));
  }
  delete binvars2;

}

//________________________________________________________________
void AliDielectronHFhelper::PrintCuts()
{

  //
  // Print cuts
  //

  // loop over all cuts
  AliInfo(" Selected cuts:");
  for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++)
    AliInfo(Form(" %d: %s [%.2f,%.2f]",icut,fCutVars->At(icut)->GetName(),fCutLowLimits(icut),fCutUpLimits(icut)));

}
Commit	Line	Data
	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 Histogram framework helper //
	18	// //
	19	/*
	20
	21	A helper class to extract objects(histograms and/or profiles) from
	22	a AliDielctronHF array of objects.
	23
	24
	25	How to use it:
	26
	27	AliDielectronHFhelper *hf = new AliDielectronHFhelper("path/to/the/output/file.root", "ConfigName");
	28	// print the structure
	29	hf->Print();
	30
	31	//apply some cuts and print them
	32	hf->SetRangeUser("cut1name",cutmin1,cutmax1);
	33	hf->SetRangeUser(AliDielectronVarManager::kPt,ptmin,ptmax);
	34	hf->PrintCuts();
	35
	36	// collect 1-,2- or 3-dim histograms or profiles with error option (default:"")
	37	TObjArray *arrHists = hf->CollectHistos(AliDielectronVarManager::kM);
	38	TObjArray *arrProfs = hf->CollectProfiles("",AliDielectronVarManager::kM,AliDielectronVarManager::kPt);
	39
	40	// then you are left with an array of histograms for all pair types or MC signals
	41
	42	*/
	43	// //
	44	///////////////////////////////////////////////////////////////////////////
	45
	46	#include <TObjArray.h>
	47	#include <TKey.h>
	48	#include <TList.h>
	49	#include <TClass.h>
	50	#include <TObject.h>
	51	#include <TFile.h>
	52	#include <TString.h>
	53	#include <TObjString.h>
	54	#include <TVectorD.h>
	55	#include <TMath.h>
	56	#include <TH1.h>
	57
	58	#include <AliLog.h>
	59
	60	#include "AliDielectron.h"
	61	#include "AliDielectronHFhelper.h"
	62	#include "AliDielectronHF.h"
	63
	64	ClassImp(AliDielectronHFhelper)
	65
	66	//const char* AliDielectronHFhelper::fCutVars[AliDielectronHFhelper::kMaxCuts] = {"","","","","","","","","",""};
	67
	68	//________________________________________________________________
	69	AliDielectronHFhelper::AliDielectronHFhelper(const char* filename, const char* container) :
	70	TNamed(),
	71	fMainArr(0x0),
	72	fCutVars(0x0),
	73	fCutLowLimits(0),
	74	fCutUpLimits(0)
	75	{
	76	//
	77	// get HF container(s) from file 'filename'
	78	//
	79	SetHFArray(filename, container);
	80
	81	}
	82
	83	//________________________________________________________________
	84	AliDielectronHFhelper::~AliDielectronHFhelper()
	85	{
	86	//
	87	// dtor
	88	//
	89	if(fMainArr) delete fMainArr;
	90	if(fCutVars) delete fCutVars;
	91
	92	}
	93
	94	//________________________________________________________________
	95	void AliDielectronHFhelper::SetHFArray(const char* filename, const char* container)
	96	{
	97	//
	98	// get HF container from file
	99	//
	100
	101	TFile *f = TFile::Open(filename);
	102
	103	TList *l=f->GetListOfKeys();
	104	TIter nextKey(l);
	105	TKey *k=0x0;
	106	while ( (k=static_cast<TKey*>(nextKey())) ){
	107
	108	TObject *o=k->ReadObj();
	109	if (o->IsA()==TList::Class()){
	110
	111	TList tlist=(TList)o;
	112
	113	TIter next(tlist);
	114	TObject *obj=0x0;
	115	while ((obj = next())) {
	116	TString objname(obj->GetName());
	117
	118	if( objname.Contains(Form("%s_HF",container)) && obj->IsA()==TObjArray::Class()) {
	119	fMainArr = new TObjArray( (dynamic_cast<TObjArray>(obj)) );
	120	fMainArr->SetOwner();
	121	// fMainArr->Print();
	122	return;
	123	}
	124	}
	125	}
	126	}
	127
	128	}
	129	//________________________________________________________________
	130	void AliDielectronHFhelper::SetRangeUser(const char *varname, Double_t min, Double_t max, Bool_t leg)
	131	{
	132	//
	133	// Set range from variable name
	134	//
	135
	136	Int_t size=fCutLowLimits.GetNrows();
	137
	138	// check if cut is already set
	139	for(Int_t icut=0; icut<size; icut++) {
	140	TString cutName = fCutVars->At(icut)->GetName();
	141	if(!cutName.CompareTo(Form("%s%s",(leg?"Leg":""),varname))) {
	142	UnsetRangeUser(varname,leg);
	143	SetRangeUser(varname, min, max, leg);
	144	return;
	145	}
	146	}
	147
	148	if(size>=kMaxCuts) return;
	149
	150	// arrays
	151	if(!fCutVars) {
	152	fCutVars = new TObjArray();
	153	fCutVars->SetOwner();
	154	}
	155	fCutLowLimits.ResizeTo(size+1);
	156	fCutUpLimits.ResizeTo(size+1);
	157
	158	// fill
	159	TObjString *str = new TObjString(Form("%s%s",(leg?"Leg":""),varname));
	160	fCutVars->Add(str);
	161	fCutLowLimits(size) = min;
	162	fCutUpLimits(size) = max;
	163	AliWarning(Form(" %s [%.2f,%.2f]",fCutVars->At(size)->GetName(),fCutLowLimits(size),fCutUpLimits(size)));
	164	}
	165
	166	//________________________________________________________________
	167	void AliDielectronHFhelper::SetRangeUser(AliDielectronVarManager::ValueTypes type, Double_t min, Double_t max, Bool_t leg)
	168	{
	169	//
	170	// Set range from AliDielectronVarManager
	171	//
	172	SetRangeUser(AliDielectronVarManager::GetValueName(type), min, max, leg);
	173	}
	174
	175	//________________________________________________________________
	176	void AliDielectronHFhelper::UnsetRangeUser(const char *varname, Bool_t leg)
	177	{
	178	//
	179	// unset range from variable name
	180	//
	181
	182	Int_t size=fCutLowLimits.GetNrows();
	183	TVectorD newlow;
	184	TVectorD newup;
	185
	186	// find cut and build new vectors w/o it
	187	Int_t ientries = 0;
	188	for(Int_t icut=0; icut<size; icut++) {
	189
	190	TString cutName = fCutVars->At(icut)->GetName();
	191	if(!cutName.CompareTo(Form("%s%s",(leg?"Leg":""),varname))) {
	192	fCutVars->AddAt(0x0,icut);
	193	continue;
	194	}
	195
	196	// fill new vectors
	197	newlow.ResizeTo(ientries+1);
	198	newup.ResizeTo(ientries+1);
	199	newlow(ientries) = fCutLowLimits(icut);
	200	newup(ientries) = fCutUpLimits(icut);
	201
	202	ientries++;
	203	}
	204
	205	// adapt new arrays/vectors
	206	fCutVars->Compress();
	207
	208	fCutLowLimits.ResizeTo(ientries);
	209	fCutUpLimits.ResizeTo(ientries);
	210	for(Int_t icut=0; icut<ientries; icut++) {
	211	fCutLowLimits(icut) = newlow(icut);
	212	fCutUpLimits(icut) = newup(icut);
	213	}
	214	// PrintCuts();
	215
	216	}
	217
	218	//________________________________________________________________
	219	void AliDielectronHFhelper::UnsetRangeUser(AliDielectronVarManager::ValueTypes type, Bool_t leg)
	220	{
	221	//
	222	// Unset range from AliDielectronVarManager
	223	//
	224	UnsetRangeUser(AliDielectronVarManager::GetValueName(type), leg);
	225	}
	226
	227	//________________________________________________________________
	228	TObjArray* AliDielectronHFhelper::CollectProfiles(TString option,
	229	AliDielectronVarManager::ValueTypes varx,
	230	AliDielectronVarManager::ValueTypes vary,
	231	AliDielectronVarManager::ValueTypes varz,
	232	AliDielectronVarManager::ValueTypes vart)
	233	{
	234	//
	235	// collect 1-3 dimensional TProfiles for all kind of pair types or sources
	236	//
	237
	238	// reconstruct the histogram/profile name resp. key
	239	Int_t dim = 0;
	240	if(varx < AliDielectronVarManager::kNMaxValues) dim++;
	241	if(vary < AliDielectronVarManager::kNMaxValues) dim++;
	242	if(varz < AliDielectronVarManager::kNMaxValues) dim++;
	243	if(vart < AliDielectronVarManager::kNMaxValues) dim++;
	244	Bool_t bPairClass=0;
	245	if( varx < AliDielectronVarManager::kPairMax \|\|
	246	vary < AliDielectronVarManager::kPairMax \|\|
	247	varz < AliDielectronVarManager::kPairMax \|\|
	248	vart < AliDielectronVarManager::kPairMax ) bPairClass=kTRUE;
	249	Bool_t bwght= (option.Contains("weight",TString::kIgnoreCase)); // weighted histogram
	250	Bool_t bprf = !(option.Contains("hist",TString::kIgnoreCase)); // tprofile
	251	Bool_t bStdOpt=kTRUE;
	252	if(option.Contains("S",TString::kIgnoreCase) \|\| option.Contains("rms",TString::kIgnoreCase)) bStdOpt=kFALSE;
	253
	254	TString key = "";
	255	if(bwght) dim--;
	256	if(bprf) dim--;
	257	switch(dim) {
	258	case 3:
	259	key+=Form("%s_",AliDielectronVarManager::GetValueName(varx));
	260	key+=Form("%s_",AliDielectronVarManager::GetValueName(vary));
	261	key+=Form("%s",AliDielectronVarManager::GetValueName(varz));
	262	if(bprf) key+=Form("-%s%s",AliDielectronVarManager::GetValueName(vart),(bStdOpt ? "avg" : "rms"));
	263	break;
	264	case 2:
	265	key+=Form("%s_",AliDielectronVarManager::GetValueName(varx));
	266	key+=Form("%s",AliDielectronVarManager::GetValueName(vary));
	267	if(bprf) key+=Form("-%s%s",AliDielectronVarManager::GetValueName(varz),(bStdOpt ? "avg" : "rms"));
	268	break;
	269	case 1:
	270	key+=Form("%s",AliDielectronVarManager::GetValueName(varx));
	271	if(bprf) key+=Form("-%s%s",AliDielectronVarManager::GetValueName(vary),(bStdOpt ? "avg" : "rms"));
	272	break;
	273	}
	274	// to differentiate btw. leg and pair histos
	275	if(bPairClass) key.Prepend("p");
	276	// prepend HF
	277	key.Prepend("HF_");
	278
	279	// add the weighted part to the key
	280	if(bwght) {
	281	if(bprf) dim++;
	282	switch(dim) {
	283	case 3: key+=Form("-wght%s",AliDielectronVarManager::GetValueName(vart)); break;
	284	case 2: key+=Form("-wght%s",AliDielectronVarManager::GetValueName(varz)); break;
	285	case 1: key+=Form("-wght%s",AliDielectronVarManager::GetValueName(vary)); break;
	286	case 4: AliError(Form(" NO weighted 3D profiles supported by the framework"));
	287	}
	288	}
	289
	290	//printf("--------> KEY: %s \n",key.Data());
	291	// init array of histograms of interest
	292	TObjArray *collection = new TObjArray(fMainArr->GetEntriesFast());
	293	collection->SetOwner(kTRUE);
	294
	295	TObjArray cloneArr = new TObjArray( (dynamic_cast<TObjArray*>(fMainArr)) );
	296	if(!cloneArr) return 0x0;
	297	cloneArr->SetOwner(kTRUE);
	298
	299	// loop over all pair types or sources
	300	for(Int_t i=0; i<cloneArr->GetEntriesFast(); i++) {
	301	if(!cloneArr->At(i)) continue;
	302	if(!((TObjArray*)cloneArr->At(i))->GetEntries()) continue;
	303
	304	// Get signal/source array with its histograms of interest
	305	AliDebug(1,Form(" Looking into step %s selected",cloneArr->At(i)->GetName()));
	306	TObjArray arr = FindObjects((TObjArray)cloneArr->At(i));
	307	if(arr) {
	308
	309	// find requested histogram
	310	collection->AddAt(arr->FindObject(key.Data()), i);
	311	if(!collection->At(i)) { AliError(Form("Object %s does not exist",key.Data())); return 0x0; }
	312
	313	// modify the signal/source name if needed (MConly)
	314	TString stepName = cloneArr->At(i)->GetName();
	315	stepName.ReplaceAll("(","");
	316	stepName.ReplaceAll(")","");
	317	stepName.ReplaceAll(": ","");
	318	stepName.ReplaceAll(" ","_");
	319	stepName.ReplaceAll("Signal","");
	320	((TH1*)collection->At(i))->SetName(Form("%s_%s",key.Data(),stepName.Data()));
	321	}
	322	else
	323	AliError(Form("Step %d not found",i));
	324
	325	}
	326
	327	//clean up
	328	//delete cloneArr;
	329	//cloneArr=0;
	330
	331	return collection;
	332	}
	333
	334	//________________________________________________________________
	335	TObjArray* AliDielectronHFhelper::FindObjects(TObjArray *stepArr)
	336	{
	337	// rename DoCuts, return values is a tobjarray
	338	// apply cuts and exclude objects from the array for merging (CUT selection)
	339	//
	340
	341	// debug
	342	// TString title = stepArr->At(0)->GetTitle();
	343	// TObjArray* vars = title.Tokenize(":");
	344	// AliDebug(1,Form(" number of cuts/vars: %d/%d",fCutLowLimits.GetNrows(),vars->GetEntriesFast()));
	345	if(!stepArr) { AliError("step is empty"); return 0x0;}
	346
	347	// check for missing cuts
	348	CheckCuts(stepArr);
	349
	350	// loop over all cuts
	351	for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++) {
	352
	353	Bool_t bFndBin = kFALSE; // bin with exact limits found
	354	const char *cutvar = fCutVars->At(icut)->GetName(); // cut variable
	355	Double_t min = fCutLowLimits(icut); // lower limit
	356	Double_t max = fCutUpLimits(icut); // upper limit
	357	AliDebug(5,Form(" Cut %d: %s [%.2f,%.2f]",icut,cutvar,min,max));
	358
	359	// loop over the full grid of the signalArray (pair,source)
	360	for(Int_t i=0; i<stepArr->GetEntriesFast(); i++) {
	361
	362	// continue if already empty
	363	if(!stepArr->At(i)) continue;
	364
	365	// collect bins from the name
	366	TString title = stepArr->At(i)->GetName();
	367	if(title.IsNull()) continue;
	368	AliDebug(5,Form(" %03d object name: %s",i,title.Data()));
	369
	370	// loop over all variables
	371	TObjArray *vars = title.Tokenize(":");
	372	for(Int_t ivar=0; ivar<vars->GetEntriesFast(); ivar++) {
	373	TString binvar = vars->At(ivar)->GetName();
	374	AliDebug(10,Form(" --> %d check bin var %s",ivar,binvar.Data()));
	375
	376	// check cuts set by the user, and compare to the bin limits
	377	if(binvar.Contains(cutvar)) {
	378	// bin limits
	379	TObjArray *limits = binvar.Tokenize("#");
	380	if(!limits) continue;
	381	Double_t binmin = atof(limits->At(1)->GetName()); // lower bin limit
	382	Double_t binmax = atof(limits->At(2)->GetName()); // upper bin limit
	383	AliDebug(10,Form(" bin %s var %s [%.2f,%.2f]",binvar.Data(),limits->At(0)->GetName(),binmin,binmax));
	384	delete limits;
	385
	386	// cut and remove objects from the array
	387	if(binmin < min \|\| binmax < min \|\| binmin > max \|\| binmax > max ) {
	388	AliDebug(10,Form(" removed, out of range! lower bin,cut: %.2f,%.2f or upper bin,cut: %.2f>%.2f",binmin,min,binmax,max));
	389	stepArr->AddAt(0x0,i);
	390	}
	391	if(bFndBin && !(binmin == min && binmax == max)) {
	392	stepArr->AddAt(0x0,i);
	393	AliDebug(10,Form(" removed, within range! lower bin,cut: %.2f,%.2f or upper bin,cut: %.2f,%.2f",binmin,min,binmax,max));
	394	}
	395
	396	// did we found a bin with exact the cut ranges
	397	// this can happen only once per variable
	398	if(binmin==min && binmax==max) bFndBin=kTRUE;
	399
	400	}
	401
	402	}
	403	// clean up
	404	if(vars) delete vars;
	405
	406	} //end loop: pair step
	407
	408	} //end: loop cuts
	409
	410	// compress the array by removing all empty entries
	411	stepArr->Compress();
	412	AliDebug(1,Form(" Compression: %d arrays left",stepArr->GetEntriesFast()));
	413
	414	// merge left objects
	415	TObjArray* hist = Merge(stepArr);
	416	if(hist) AliDebug(1,Form(" final array has %d objects",hist->GetEntries()));
	417	return hist;
	418	}
	419
	420	//________________________________________________________________
	421	TObjArray* AliDielectronHFhelper::Merge(TObjArray *arr)
	422	{
	423	//
	424	// merge left objects into a single one (LAST step)
	425	//
	426
	427	if(arr->GetEntriesFast()<1) { AliError(" No more objects left!"); return 0x0; }
	428
	429	TObjArray final=(TObjArray) arr->At(0)->Clone();
	430	if(!final) return 0x0;
	431	final->SetOwner();
	432
	433	TList listH;
	434	TString listHargs;
	435	listHargs.Form("((TCollection*)0x%lx)", (ULong_t)&listH);
	436	Int_t error = 0;
	437
	438	// final->Reset("CE");
	439	// final->SetTitle(""); //TODO: change in future
	440	for(Int_t i=1; i<arr->GetEntriesFast(); i++) {
	441	listH.Add(arr->At(i));
	442	// printf("%d: ent %.0f \n",i,((TH1)((TObjArray)arr->At(i))->At(0))->GetEntries());
	443	// final->Add((TH1*)arr->At(i));
	444	}
	445	// arr->Clear();
	446
	447	final->Execute("Merge", listHargs.Data(), &error);
	448	// returns a tobjarray of histograms/profiles
	449	return final;
	450	}
	451
	452	//________________________________________________________________
	453	void AliDielectronHFhelper::CheckCuts(TObjArray *arr)
	454	{
	455	//
	456	// Compare binning and cut variables. Add necessary cuts (full range, no exclusion)
	457	//
	458
	459	// build array with bin variable, minimum and maximum bin values
	460	TString titleFIRST = arr->First()->GetName();
	461	TString titleLAST = arr->Last()->GetName();
	462	TObjArray* binvarsF = titleFIRST.Tokenize(":#");
	463	TObjArray* binvarsL = titleLAST.Tokenize(":#");
	464	Double_t binmin[kMaxCuts]= {0.0};
	465	Double_t binmax[kMaxCuts]= {0.0};
	466	if(!binvarsF) { delete binvarsL; return; }
	467	if(!binvarsL) { delete binvarsF; return; }
	468	for(Int_t ivar=0; ivar<binvarsF->GetEntriesFast(); ivar++) {
	469
	470	TString elementF=binvarsF->At(ivar)->GetName();
	471	TString elementL=binvarsL->At(ivar)->GetName();
	472	AliDebug(1,Form(" binvar %d: %s,%s",ivar,elementF.Data(),elementL.Data()));
	473
	474	switch(ivar%3) {
	475	case 0: continue; break;
	476	case 1: binmin[(int)ivar/3]=atof(elementF.Data()); break;
	477	case 2: binmax[(int)ivar/3]=atof(elementL.Data()); break;
	478	}
	479
	480	binvarsF->AddAt(0x0,ivar);
	481	}
	482	binvarsF->Compress();
	483
	484	// loop over all vars and cuts, check for missing stuff
	485	for(Int_t ivar=0; ivar<binvarsF->GetEntriesFast(); ivar++) {
	486
	487	TString binvar=binvarsF->At(ivar)->GetName();
	488	Bool_t selected=kFALSE;
	489
	490	AliDebug(1,Form(" check cuts %d %s [%.2f,%.2f]",ivar,binvar.Data(),binmin[ivar],binmax[ivar]));
	491	// loop over all cuts and check for missing stuff
	492	for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++) {
	493	if(binvar.Contains(fCutVars->At(icut)->GetName())) { selected=kTRUE; break; }
	494	// else break;
	495	}
	496
	497	// add missing cut with max limits
	498	if(!selected) {
	499	AliWarning(Form(" Bin variable %s not covered. Add cut!",binvar.Data()));
	500	Bool_t leg = binvar.BeginsWith("Leg");
	501	if(leg) binvar.Remove(0,3);
	502	SetRangeUser(binvar.Data(),binmin[ivar],binmax[ivar], leg);
	503	}
	504
	505	}
	506
	507	// clean up
	508	delete binvarsF;
	509	delete binvarsL;
	510	}
	511
	512	//________________________________________________________________
	513	void AliDielectronHFhelper::Print(const Option_t* /option/) const
	514	{
	515
	516	//
	517	// Print out object contents
	518	//
	519	AliInfo(Form(" Container: %s",fMainArr->GetName()));
	520
	521	// pairtypes, steps and sources
	522	Int_t stepLast=0;
	523	AliInfo(Form(" Number of filled steps: %d",fMainArr->GetEntries()));
	524	for(Int_t istep=0; istep<fMainArr->GetEntriesFast(); istep++) {
	525	if(!fMainArr->At(istep)) continue;
	526	if(!((TObjArray*)fMainArr->At(istep))->GetEntries()) continue;
	527	AliInfo(Form(" step %d: %s",istep,fMainArr->At(istep)->GetName()));
	528	stepLast=istep;
	529	}
	530
	531	AliInfo(Form(" Number of objects: %d",
	532	((TObjArray) ((TObjArray)fMainArr->At(stepLast)) ->First())->GetEntriesFast()));
	533
	534	TString title = ((TObjArray*)fMainArr->At(stepLast))->First()->GetName();
	535	TObjArray* binvars = title.Tokenize(":");
	536	AliInfo(Form(" Number of variables: %d",binvars->GetEntriesFast()));
	537	delete binvars;
	538
	539	/* REACTIVATE
	540	// check for missing cuts
	541	CheckCuts(((TObjArray*)fMainArr->At(stepLast)));
	542	// loop over all cuts
	543	for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++) {
	544	const char *cutvar = fCutVars->At(icut)->GetName(); // cut variable
	545	Double_t min = fCutLowLimits(icut); // lower limit
	546	Double_t max = fCutUpLimits(icut); // upper limit
	547	AliInfo(Form(" variable %d: %s [%.2f,%.2f]",icut,cutvar,min,max));
	548	}
	549	*/
	550	TObjArray* binvars2 = title.Tokenize(":#");
	551	for(Int_t ivar=0; ivar<binvars2->GetEntriesFast(); ivar++) {
	552	if(ivar%3) continue;
	553	AliInfo(Form(" variable %.0f: %s",((Double_t)ivar)/3+1,binvars2->At(ivar)->GetName()));
	554	}
	555	delete binvars2;
	556
	557	}
	558
	559	//________________________________________________________________
	560	void AliDielectronHFhelper::PrintCuts()
	561	{
	562
	563	//
	564	// Print cuts
	565	//
	566
	567	// loop over all cuts
	568	AliInfo(" Selected cuts:");
	569	for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++)
	570	AliInfo(Form(" %d: %s [%.2f,%.2f]",icut,fCutVars->At(icut)->GetName(),fCutLowLimits(icut),fCutUpLimits(icut)));
	571
	572	}
	573