PairAnalysisHF.cxx

Go to the documentation of this file.

/*************************************************************************
* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
**************************************************************************/
 
//                PairAnalysis HF                                  //
//                                                                       //
//                                                                       //
/*
Detailed description
 
 
*/
//                                                                       //
 
#include <TAxis.h>
#include <TFile.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2.h>
#include <TH3.h>
#include <THnSparse.h>
#include <TKey.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TProfile3D.h>
#include <TString.h>
#include <TVectorD.h>
 
#include "CbmMCTrack.h"
 
#include "PairAnalysis.h"
#include "PairAnalysisHelper.h"
#include "PairAnalysisMC.h"
#include "PairAnalysisPair.h"
#include "PairAnalysisSignalMC.h"
#include "PairAnalysisStyler.h"
 
#include "PairAnalysisHF.h"
#include "PairAnalysisHistos.h"
 
#include "PairAnalysisTrack.h"
 
ClassImp(PairAnalysisHF)
 
  PairAnalysisHF::PairAnalysisHF()
  :  //TNamed(),
  PairAnalysisHistos()
  , fUsedVars(new TBits(PairAnalysisVarManager::kNMaxValuesMC))
  , fArrDielHistos()
  , fSignalsMC(0x0)
  , fVarCutType(new TBits(kMaxCuts))
  , fAxes(kMaxCuts) {
  //
  // Default Constructor
  //
  for (Int_t i = 0; i < kMaxCuts; ++i) {
    fVarCuts[i] = 0;
    //    fVarCutType[i]=0;
  }
  fAxes.SetOwner(kTRUE);
  fArrDielHistos.SetOwner(kTRUE);
}
 
//______________________________________________
PairAnalysisHF::PairAnalysisHF(const char* name, const char* title)
  :  //  TNamed(name, title),
  PairAnalysisHistos(name, title)
  , fUsedVars(new TBits(PairAnalysisVarManager::kNMaxValuesMC))
  , fArrDielHistos()
  , fSignalsMC(0x0)
  , fVarCutType(new TBits(kMaxCuts))
  , fAxes(kMaxCuts) {
  //
  // Named Constructor
  //
  for (Int_t i = 0; i < kMaxCuts; ++i) {
    fVarCuts[i] = 0;
    //    fVarCutType[i]=0;
  }
  fAxes.SetOwner(kTRUE);
  fArrDielHistos.SetOwner(kTRUE);
}
 
//______________________________________________
PairAnalysisHF::~PairAnalysisHF() {
  //
  // Default Destructor
  //
  if (fUsedVars) delete fUsedVars;
  if (fVarCutType) delete fVarCutType;
  fAxes.Delete();
  fArrDielHistos.Delete();  //TODO: better Clear?
}
 
//________________________________________________________________
void PairAnalysisHF::AddCutVariable(PairAnalysisVarManager::ValueTypes type,
                                    TVectorD* binLimits,
                                    Bool_t leg) {
  //
  // Add a variable to the histogram array 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;
  fVarCutType->SetBitNumber(size, leg);
  fAxes.Add(binLimits);
  fUsedVars->SetBitNumber(type, kTRUE);
}
 
//_____________________________________________________________________________
void PairAnalysisHF::FillClass(const char* histClass, const Double_t* values) {
  //
  // Fill the histograms if cuts are passed
  //
 
  // find cell described by values
  Int_t cell = FindCell(values);
  //  printf("cell: %d \n",cell);
  if (cell < 0) return;
  //  printf("  --> %s \n",fArrDielHistos.UncheckedAt(cell)->GetName());
 
  // do NOT set the ownership otherwise you will delete all histos!!
  SetHistogramList(*static_cast<THashList*>(fArrDielHistos.UncheckedAt(cell)),
                   kFALSE);
  PairAnalysisHistos::FillClass(histClass, values);
}
 
//______________________________________________
void PairAnalysisHF::ReadFromFile(const char* file,
                                  const char* task,
                                  const char* config) {
  //
  // Read histos from file
  //
  // TODO: to be tested!
  TFile f(file);
  TIter nextKey(f.GetListOfKeys());
  TKey* key = 0;
  while ((key = (TKey*) nextKey())) {
    TString name = key->GetName();
    if (!name.Contains("PairAnalysisHistos")) continue;
    if (!strlen(task) && !name.Contains(task)) continue;
    TObject* o  = f.Get(key->GetName());
    TList* list = dynamic_cast<TList*>(o);
    if (!list) continue;
    TObjArray* listCfg =
      dynamic_cast<TObjArray*>(list->FindObject(Form("%s_HF", config)));
    if (!listCfg) continue;
    break;
  }
  f.Close();
  // load style
  PairAnalysisStyler::LoadStyle();
}
 
//______________________________________________
void PairAnalysisHF::Fill(Int_t /*label1*/,
                          Int_t /*label2*/,
                          Int_t /*nSignal*/) {
  return;  //TODO: OBSOLETE?
  /*
  //
  // 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 || fEventArray) return;
 
  CbmMCTrack* part1 = PairAnalysisMC::Instance()->GetMCTrackFromMCEvent(label1);
  CbmMCTrack* part2 = PairAnalysisMC::Instance()->GetMCTrackFromMCEvent(label2);
  if(!part1 || !part2) return;
 
  PairAnalysisMC* papaMC = PairAnalysisMC::Instance();
 
  Int_t mLabel1 = papaMC->GetMothersLabel(label1);
  Int_t mLabel2 = papaMC->GetMothersLabel(label2);
 
  // check the same mother option
  PairAnalysisSignalMC* sigMC = (PairAnalysisSignalMC*)fSignalsMC->At(nSignal);
  if(sigMC->GetMothersRelation()==PairAnalysisSignalMC::kSame      && mLabel1!=mLabel2) return;
  if(sigMC->GetMothersRelation()==PairAnalysisSignalMC::kDifferent && mLabel1==mLabel2) return;
 
  PairAnalysisVarManager::SetFillMap(fUsedVars);
  // fill the leg variables
  Double_t valuesLeg1[PairAnalysisVarManager::kNMaxValuesMC];
  Double_t valuesLeg2[PairAnalysisVarManager::kNMaxValuesMC];
  PairAnalysisVarManager::Fill(part1,valuesLeg1);
  PairAnalysisVarManager::Fill(part2,valuesLeg2);
 
  // fill the pair and event variables
  Double_t valuesPair[PairAnalysisVarManager::kNMaxValuesMC];
  //TODO PairAnalysisVarManager::Fill(papaMC->GetMCEvent(), valuesPair);
  PairAnalysisVarManager::FillVarMCParticle(part1,part2,valuesPair);
 
  // if pair types are filled, fill mc sources at the end
  Int_t istep=0;
  if(fPairType!=kMConly) istep=PairAnalysis::kSEPMRot+1;
 
  // only OS at the moment
  if(part1->GetCharge()*part2->GetCharge()<0) {
    Fill(istep+nSignal+fSignalsMC->GetEntries(), valuesPair,  valuesLeg1, valuesLeg2);
  }
 
  return;
  */
}
//______________________________________________
void PairAnalysisHF::Fill(Int_t /*pairIndex*/,
                          const PairAnalysisPair* /*particle*/) {
  //
  // fill histograms for event, pair and daughter cuts and pair types
  //
  return;  //TODO: OBSOLETE?
  /*
  // only OS pairs in case of MC
 
  // only selected pair types in case of data
  if(!IsPairTypeSelected(pairIndex) || fEventArray) return;
 
  // get event and pair variables
  Double_t valuesPair[PairAnalysisVarManager::kNMaxValuesMC];
  PairAnalysisVarManager::SetFillMap(fUsedVars);
  PairAnalysisVarManager::Fill(particle,valuesPair);
 
  // get leg variables (TODO: do not fill for the moment since leg cuts are not opened)
  Double_t valuesLeg1[PairAnalysisVarManager::kNMaxValuesMC]={0};
  if(fVarCutType->CountBits())  PairAnalysisVarManager::Fill(particle->GetFirstDaughter(),valuesLeg1);
  Double_t valuesLeg2[PairAnalysisVarManager::kNMaxValuesMC]={0};
  if(fVarCutType->CountBits())  PairAnalysisVarManager::Fill(particle->GetSecondDaughter(),valuesLeg2);
 
  // fill
 
  // if pair types are filled, fill mc sources at the end
  Int_t istep = 0;
  if(fPairType!=kMConly) istep=PairAnalysis::kSEPMRot+1;
 
  // mc source steps (only OS SE pairs)
  if(fHasMC && fSignalsMC && pairIndex==PairAnalysis::kSEPM) {
    for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
        Fill(istep+i, valuesPair,  valuesLeg1, valuesLeg2);
    }
  }
 
  // all pair types w/o use of mc information
  if(fPairType==kMConly) return;
 
  // remove comments
  Fill(pairIndex, valuesPair,  valuesLeg1, valuesLeg2); 
 
  return;
  */
}
 
//______________________________________________
void PairAnalysisHF::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
  //
  return;  //TODO: OBSOLETE?
  /*
  TObjArray *histArr = static_cast<TObjArray*>(fArrDielHistos.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->TestBitNumber(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 object with Pair and event values
    TObjArray *tmp = (TObjArray*) histArr->At(ihist);
    TString title = tmp->GetName();
    for(Int_t i=0; i<tmp->GetEntriesFast(); i++) {
      PairAnalysisHistos::FillValues(tmp->At(i), valuesPair);
    }
    //    Debug(10,Form("Fill var %d %s value %f in %s \n",fVar,PairAnalysisVarManager::GetValueName(fVar),valuesPair[fVar],tmp->GetName()));
  } //end of hist loop
  */
}
 
//______________________________________________
void PairAnalysisHF::Init() {
  //
  // initialise the HF array
  //
 
  // create an TObjArray of 'size' with PairAnalysisHistos objects
  Int_t size = GetNumberOfBins();
 
  fArrDielHistos.SetName(Form("%s_HF", GetName()));
  fArrDielHistos.Expand(size);
 
  // add 'PairAnalysisHistos'-list to each cell
  for (Int_t icell = 0; icell < size; icell++) {
    fArrDielHistos.AddAt(fHistoList.Clone(""), icell);
  }
 
  // loop over all cut variables and do the naming according to its bin cell
  Int_t sizeAdd = 1;  // counter for processed cells
  Int_t nvars   = fAxes.GetEntriesFast();
  for (Int_t ivar = 0; ivar < nvars; ivar++) {
 
    // get bin limits for variabvle ivar
    TVectorD* bins = static_cast<TVectorD*>(fAxes.At(ivar));
    Int_t nbins    = bins->GetNrows() - 1;
 
    // loop over all cells
    // Add 'variable' name and current 'bin limits' to
    // the title of the 'PairAnalysisHistos'-list title
    // which is unique
    for (Int_t icell = 0; icell < size; icell++) {
 
      // get the lower limit for current ivar bin
      Int_t ibin       = (icell / sizeAdd) % nbins;
      Double_t lowEdge = (*bins)[ibin];
      Double_t upEdge  = (*bins)[ibin + 1];
 
      // modify title of hashlist
      TString title = fArrDielHistos.UncheckedAt(icell)->GetName();
      if (!ivar)
        title = "";
      else
        title += ":";  // delimiter for new variable
      if (fVarCutType->TestBitNumber(ivar))
        title += "Leg";  // for leg variable Identification
      title += PairAnalysisVarManager::GetValueName(fVarCuts[ivar]);
      title += Form("#%.2f#%.2f", lowEdge, upEdge);  // TODO: precision enough?
      static_cast<THashList*>(fArrDielHistos.UncheckedAt(icell))
        ->SetName(title.Data());
 
    }  // end: array of cells
 
    sizeAdd *= nbins;
 
  }  //end: loop variables
 
  // clean up
  // if(fArrDielHistos) {
  //   fArrDielHistos.Delete();
  //   fArrDielHistos=0;
  // }
}
 
//______________________________________________
Int_t PairAnalysisHF::GetNumberOfBins() const {
  //
  // return the number of bins this histogram grid 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;
}
 
//______________________________________________
Int_t PairAnalysisHF::FindCell(const Double_t* values) {
  //
  // cell described by 'values'
  // if the values are outside the binning range -1 is returned
  //
 
  if (fAxes.GetEntriesFast() == 0) return 0;
 
  Int_t sizeAdd = 1;
  Int_t cell    = 0;
  for (Int_t i = 0; i < fAxes.GetEntriesFast(); ++i) {
    Double_t val   = values[fVarCuts[i]];
    TVectorD* bins = static_cast<TVectorD*>(fAxes.At(i));
    Int_t nRows    = bins->GetNrows();
    if ((val < (*bins)[0]) || (val > (*bins)[nRows - 1])) { return -1; }
 
    Int_t pos = TMath::BinarySearch(nRows, bins->GetMatrixArray(), val);
    cell += sizeAdd * pos;
    // printf(" \t %s: %.2f-%.2f %d \n",
    //     PairAnalysisVarManager::GetValueName(fVarCuts[i]), (*bins)[pos], (*bins)[pos+1], pos);
    sizeAdd *= (nRows - 1);
  }
 
  return cell;
}