]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGDQ/dielectron/AliDielectronHFhelper.cxx
including switch to set on/off iso-track core removal, cleaning and bug fix
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronHFhelper.cxx
index 58faf81c93add70762657f585920be35e13aa67b..8e4e4602810f2affdfb6b5e5c26002a24771cc32 100644 (file)
 **************************************************************************/
 
 ///////////////////////////////////////////////////////////////////////////
-//       Dielectron Histogram framework helper                     //
+//       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 "AliDielectronHFhelper.h"
 #include "AliDielectronHF.h"
 
-//ClassImp(AliDielectronHFhelper)
+ClassImp(AliDielectronHFhelper)
 
 //const char* AliDielectronHFhelper::fCutVars[AliDielectronHFhelper::kMaxCuts] = {"","","","","","","","","",""};
 
 //________________________________________________________________
 AliDielectronHFhelper::AliDielectronHFhelper(const char* filename, const char* container) :
   TNamed(),
-  fArrPairType(0x0),
+  fMainArr(0x0),
   fCutVars(0x0),
   fCutLowLimits(0),
   fCutUpLimits(0)
@@ -73,7 +86,7 @@ AliDielectronHFhelper::~AliDielectronHFhelper()
   //
   // dtor
   //
-  if(fArrPairType) delete fArrPairType;
+  if(fMainArr) delete fMainArr;
   if(fCutVars)     delete fCutVars;
 
 }
@@ -82,7 +95,7 @@ AliDielectronHFhelper::~AliDielectronHFhelper()
 void AliDielectronHFhelper::SetHFArray(const char* filename, const char* container)
 {
   //
-  // get HF containers from file
+  // get HF container from file
   //
 
   TFile *f = TFile::Open(filename);
@@ -103,8 +116,9 @@ void AliDielectronHFhelper::SetHFArray(const char* filename, const char* contain
        TString objname(obj->GetName());
 
        if( objname.Contains(Form("%s_HF",container)) && obj->IsA()==TObjArray::Class()) {
-         fArrPairType = new TObjArray( *(dynamic_cast<TObjArray*>(obj)) );
-         //fArrPairType->Print();
+         fMainArr = new TObjArray( *(dynamic_cast<TObjArray*>(obj)) );
+         fMainArr->SetOwner();
+         //      fMainArr->Print();
          return;
        }
       }
@@ -118,7 +132,6 @@ void AliDielectronHFhelper::SetRangeUser(const char *varname, Double_t min, Doub
   //
   // Set range from variable name
   //
-  //  Int_t size=sizeof(fCutVars)/sizeof(const char*);
 
   Int_t size=fCutLowLimits.GetNrows();
 
@@ -165,8 +178,8 @@ void AliDielectronHFhelper::UnsetRangeUser(const char *varname, Bool_t leg)
   //
   // unset range from variable name
   //
+
   Int_t size=fCutLowLimits.GetNrows();
-  //  PrintCuts();
   TVectorD newlow;
   TVectorD newup;
 
@@ -212,117 +225,176 @@ void AliDielectronHFhelper::UnsetRangeUser(AliDielectronVarManager::ValueTypes t
 }
 
 //________________________________________________________________
-TObjArray* AliDielectronHFhelper::CollectHistos()
+TObjArray* AliDielectronHFhelper::CollectProfiles(TString option,
+                                                 AliDielectronVarManager::ValueTypes varx,
+                                                 AliDielectronVarManager::ValueTypes vary,
+                                                 AliDielectronVarManager::ValueTypes varz,
+                                                 AliDielectronVarManager::ValueTypes vart)
 {
   //
-  // collect histograms for all kind of pair types or sources
+  // collect 1-3 dimensional TProfiles for all kind of pair types or sources
   //
 
-  TObjArray *collection = new TObjArray(AliDielectron::kEv1PMRot+1);
-
-  TObjArray *histArr = (TObjArray*) fArrPairType->Clone("tmpArr");
-  histArr->SetOwner(kTRUE);
-
-  // loop over max. available pair types
-  for(Int_t i=0; i<AliDielectron::kEv1PMRot+1; i++) {
-
-    collection->AddAt(GetHistogram(AliDielectron::PairClassName(i),histArr), i);
-
+  // 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;
   }
-
-  // clean up the clone
-  if(histArr) {
-    delete histArr;
-    histArr=0;
+  // 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"));
+    }
   }
 
-  return collection;
-}
-
-//________________________________________________________________
-TH1F* AliDielectronHFhelper::GetHistogram(const char *step, TObjArray *histArr)
-{
-  //
-  // main function to recieve a single histogram
-  // TODO: check memory
-
-  AliDebug(1,Form(" Step %s selected",step));
+  //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));
 
-  TObjArray *histos= 0x0;
-  TH1F *hist       = 0x0;
-  if(!histArr) {
-    histos = (TObjArray*) fArrPairType->FindObject(step)->Clone("tmpArr");
-  }
-  else {
-    histos = (TObjArray*) histArr->FindObject(step);
   }
 
-  if(histos) hist   = FindHistograms(histos);
-  return hist;
+  //clean up
+  //delete cloneArr;
+  //cloneArr=0;
 
+  return collection;
 }
 
 //________________________________________________________________
-TH1F* AliDielectronHFhelper::FindHistograms(TObjArray *histos)
+TObjArray* AliDielectronHFhelper::FindObjects(TObjArray *stepArr)
 {
-  //
-  // exclude histograms
+  // rename DoCuts, return values is a tobjarray
+  // apply cuts and exclude objects from the array for merging (CUT selection)
   //
 
   // debug
-  // TString title    = histos->At(0)->GetTitle();
+  // 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(histos);
+  CheckCuts(stepArr);
 
   // loop over all cuts
   for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++) {
 
-    Bool_t bFndBin = kFALSE; // exact bin found
-    const char *cutvar   = fCutVars->At(icut)->GetName();
-    Double_t min   = fCutLowLimits(icut);
-    Double_t max   = fCutUpLimits(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 all histograms
-    for(Int_t i=0; i<histos->GetEntriesFast(); i++) {
+    // loop over the full grid of the signalArray (pair,source)
+    for(Int_t i=0; i<stepArr->GetEntriesFast(); i++) {
 
       // continue if already empty
-      if(!histos->At(i)) continue;
+      if(!stepArr->At(i)) continue;
 
-      // collect binning from histo title
-      TString title    = histos->At(i)->GetTitle();
+      // collect bins from the name
+      TString title    = stepArr->At(i)->GetName();
       if(title.IsNull()) continue;
-      AliDebug(10,Form(" histo title: %s",title.Data()));
+      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(" Check ivar %d binvar %s",ivar,binvar.Data()));
+       AliDebug(10,Form(" --> %d check bin var %s",ivar,binvar.Data()));
 
-       // check for cuts and ranges by the user
+       // check cuts set by the user, and compare to the bin limits
        if(binvar.Contains(cutvar)) {
+         // bin limits
          TObjArray *limits = binvar.Tokenize("#");
-
-         Double_t binmin = atof(limits->At(1)->GetName());
-         Double_t binmax = atof(limits->At(2)->GetName());
+         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;
 
-         // remove histogram from array
+         // cut and remove objects from the array
          if(binmin < min || binmax < min || binmin > max || binmax > max ) {
-           AliDebug(10,Form(" removed, out of range min %.2f,%.2f  max %.2f,%.2f",binmin,min,binmax,max));
-           histos->AddAt(0x0,i);
+           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)) {
-             histos->AddAt(0x0,i);
-             AliDebug(10,Form(" removed, within range min %.2f,%.2f  max %.2f,%.2f",binmin,min,binmax,max));
-           }
-         // clean up
-         if(limits) delete limits;
+           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));
+         }
 
-         // do we have found an exact bin
+         // did we found a bin with exact the cut ranges
+         // this can happen only once per variable
          if(binmin==min && binmax==max) bFndBin=kTRUE;
 
        }
@@ -331,39 +403,49 @@ TH1F* AliDielectronHFhelper::FindHistograms(TObjArray *histos)
       // clean up
       if(vars)   delete vars;
 
-    }
+    } //end loop: pair step
 
-  }
+  } //end: loop cuts
 
-  // compress the array by removing all empty histos
-  histos->Compress();
-  AliDebug(1,Form(" Compression: %d histograms left",histos->GetEntriesFast()));
+  // compress the array by removing all empty entries
+  stepArr->Compress();
+  AliDebug(1,Form(" Compression: %d arrays left",stepArr->GetEntriesFast()));
 
-  // merge histograms
-  TH1F* hist = MergeHistos(histos);
-  if(hist) AliDebug(1,Form(" Merging: %e histogram entries",hist->GetEntries()));
+  // merge left objects
+  TObjArray* hist = Merge(stepArr);
+  if(hist) AliDebug(1,Form(" final array has %d objects",hist->GetEntries()));
   return hist;
 }
 
 //________________________________________________________________
-TH1F* AliDielectronHFhelper::MergeHistos(TObjArray *arr)
+TObjArray* AliDielectronHFhelper::Merge(TObjArray *arr)
 {
   //
-  // merge histos to one single histogram
+  // merge left objects into a single one (LAST step)
   //
 
-  if(arr->GetEntriesFast()<1) { AliError(" No more histosgrams left!"); return 0x0; }
+  if(arr->GetEntriesFast()<1) { AliError(" No more objects left!"); return 0x0; }
 
-  TH1F *final=(TH1F*) arr->At(0)->Clone();
+  TObjArray *final=(TObjArray*) arr->At(0)->Clone();
   if(!final) return 0x0;
-
-  final->Reset("CE");
-  final->SetTitle(""); //TODO: change in future
-  for(Int_t i=0; i<arr->GetEntriesFast(); i++) {
-    final->Add((TH1F*)arr->At(i));
+  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();
+  //  arr->Clear();
 
+  final->Execute("Merge", listHargs.Data(), &error);
+  // returns a tobjarray of histograms/profiles
   return final;
 }
 
@@ -371,17 +453,18 @@ TH1F* AliDielectronHFhelper::MergeHistos(TObjArray *arr)
 void AliDielectronHFhelper::CheckCuts(TObjArray *arr)
 {
   //
-  // Compare histo binning and cut variables. Add necessary cuts (largest limits)
+  // 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()->GetTitle();
-  TString titleLAST        = arr->Last()->GetTitle();
+  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();
@@ -422,8 +505,8 @@ void AliDielectronHFhelper::CheckCuts(TObjArray *arr)
   }
 
   // clean up
-  if(binvarsF) delete binvarsF;
-  if(binvarsL) delete binvarsL;
+  delete binvarsF;
+  delete binvarsL;
 }
 
 //________________________________________________________________
@@ -433,26 +516,43 @@ void AliDielectronHFhelper::Print(const Option_t* /*option*/) const
   //
   // Print out object contents
   //
-  AliInfo(Form(" Container:               %s",fArrPairType->GetName()));
+  AliInfo(Form(" Container:               %s",fMainArr->GetName()));
 
   // pairtypes, steps and sources
-  AliInfo(Form(" Number of filled steps:  %d",fArrPairType->GetEntries()));
-  for(Int_t istep=0; istep<fArrPairType->GetEntriesFast(); istep++) {
-    if(fArrPairType->At(istep))
-      AliInfo(Form(" step %d: %s",istep,fArrPairType->At(istep)->GetName()));
+  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 histograms:    %d",((TObjArray*)fArrPairType->At(0))->GetEntriesFast()));
+  AliInfo(Form(" Number of objects:    %d",
+              ((TObjArray*) ((TObjArray*)fMainArr->At(stepLast)) ->First())->GetEntriesFast()));
 
-  TString title       = ((TObjArray*)fArrPairType->At(0))->First()->GetTitle();
+  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;
 
 }
 
@@ -461,7 +561,7 @@ void AliDielectronHFhelper::PrintCuts()
 {
 
   //
-  // Print out object contents
+  // Print cuts
   //
 
   // loop over all cuts