PWGDQ/dielectron/AliDielectronHFhelper.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////
  17 //       Dielectron Histogram framework helper                     //
  18 //                                                                       //
  19 /*
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29 */
  30 //                                                                       //
  31 ///////////////////////////////////////////////////////////////////////////
  32
  33 #include <TObjArray.h>
  34 #include <TKey.h>
  35 #include <TList.h>
  36 #include <TClass.h>
  37 #include <TObject.h>
  38 #include <TFile.h>
  39 #include <TString.h>
  40 #include <TObjString.h>
  41 #include <TVectorD.h>
  42 #include <TMath.h>
  43 #include <TH1.h>
  44
  45 #include <AliLog.h>
  46
  47 #include "AliDielectron.h"
  48 #include "AliDielectronHFhelper.h"
  49 #include "AliDielectronHF.h"
  50
  51 //ClassImp(AliDielectronHFhelper)
  52
  53 //const char* AliDielectronHFhelper::fCutVars[AliDielectronHFhelper::kMaxCuts] = {"","","","","","","","","",""};
  54
  55 //________________________________________________________________
  56 AliDielectronHFhelper::AliDielectronHFhelper(const char* filename, const char* container) :
  57   TNamed(),
  58   fArrPairType(0x0),
  59   fCutVars(0x0),
  60   fCutLowLimits(0),
  61   fCutUpLimits(0)
  62 {
  63   //
  64   // get HF container(s) from file 'filename'
  65   //
  66   SetHFArray(filename, container);
  67
  68 }
  69
  70 //________________________________________________________________
  71 AliDielectronHFhelper::~AliDielectronHFhelper()
  72 {
  73   //
  74   // dtor
  75   //
  76   if(fArrPairType) delete fArrPairType;
  77   if(fCutVars)     delete fCutVars;
  78
  79 }
  80
  81 //________________________________________________________________
  82 void AliDielectronHFhelper::SetHFArray(const char* filename, const char* container)
  83 {
  84   //
  85   // get HF containers from file
  86   //
  87
  88   TFile *f = TFile::Open(filename);
  89
  90   TList *l=f->GetListOfKeys();
  91   TIter nextKey(l);
  92   TKey *k=0x0;
  93   while ( (k=static_cast<TKey*>(nextKey())) ){
  94
  95     TObject *o=k->ReadObj();
  96     if (o->IsA()==TList::Class()){
  97
  98       TList *tlist=(TList*)o;
  99
 100       TIter next(tlist);
 101       TObject *obj=0x0;
 102       while ((obj = next())) {
 103         TString objname(obj->GetName());
 104
 105         if( objname.Contains(Form("%s_HF",container)) && obj->IsA()==TObjArray::Class()) {
 106           fArrPairType = new TObjArray( *(dynamic_cast<TObjArray*>(obj)) );
 107           //fArrPairType->Print();
 108           return;
 109         }
 110       }
 111     }
 112   }
 113
 114 }
 115 //________________________________________________________________
 116 void AliDielectronHFhelper::SetRangeUser(const char *varname, Double_t min, Double_t max, Bool_t leg)
 117 {
 118   //
 119   // Set range from variable name
 120   //
 121   //  Int_t size=sizeof(fCutVars)/sizeof(const char*);
 122
 123   Int_t size=fCutLowLimits.GetNrows();
 124
 125   // check if cut is already set
 126   for(Int_t icut=0; icut<size; icut++) {
 127     TString cutName = fCutVars->At(icut)->GetName();
 128     if(!cutName.CompareTo(Form("%s%s",(leg?"Leg":""),varname))) {
 129       UnsetRangeUser(varname,leg);
 130       SetRangeUser(varname, min, max, leg);
 131       return;
 132     }
 133   }
 134
 135   if(size>=kMaxCuts) return;
 136
 137   // arrays
 138   if(!fCutVars) {
 139     fCutVars = new TObjArray();
 140     fCutVars->SetOwner();
 141   }
 142   fCutLowLimits.ResizeTo(size+1);
 143   fCutUpLimits.ResizeTo(size+1);
 144
 145   // fill
 146   TObjString *str = new TObjString(Form("%s%s",(leg?"Leg":""),varname));
 147   fCutVars->Add(str);
 148   fCutLowLimits(size) = min;
 149   fCutUpLimits(size)  = max;
 150   AliWarning(Form(" %s [%.2f,%.2f]",fCutVars->At(size)->GetName(),fCutLowLimits(size),fCutUpLimits(size)));
 151 }
 152
 153 //________________________________________________________________
 154 void AliDielectronHFhelper::SetRangeUser(AliDielectronVarManager::ValueTypes type, Double_t min, Double_t max, Bool_t leg)
 155 {
 156   //
 157   // Set range from AliDielectronVarManager
 158   //
 159   SetRangeUser(AliDielectronVarManager::GetValueName(type), min, max, leg);
 160 }
 161
 162 //________________________________________________________________
 163 void AliDielectronHFhelper::UnsetRangeUser(const char *varname, Bool_t leg)
 164 {
 165   //
 166   // unset range from variable name
 167   //
 168   Int_t size=fCutLowLimits.GetNrows();
 169   //  PrintCuts();
 170   TVectorD newlow;
 171   TVectorD newup;
 172
 173   // find cut and build new vectors w/o it
 174   Int_t ientries = 0;
 175   for(Int_t icut=0; icut<size; icut++) {
 176
 177     TString cutName = fCutVars->At(icut)->GetName();
 178     if(!cutName.CompareTo(Form("%s%s",(leg?"Leg":""),varname))) {
 179       fCutVars->AddAt(0x0,icut);
 180       continue;
 181     }
 182
 183     // fill new vectors
 184     newlow.ResizeTo(ientries+1);
 185     newup.ResizeTo(ientries+1);
 186     newlow(ientries) = fCutLowLimits(icut);
 187     newup(ientries)  = fCutUpLimits(icut);
 188
 189     ientries++;
 190   }
 191
 192   // adapt new arrays/vectors
 193   fCutVars->Compress();
 194
 195   fCutLowLimits.ResizeTo(ientries);
 196   fCutUpLimits.ResizeTo(ientries);
 197   for(Int_t icut=0; icut<ientries; icut++) {
 198     fCutLowLimits(icut) = newlow(icut);
 199     fCutUpLimits(icut)  = newup(icut);
 200   }
 201   // PrintCuts();
 202
 203 }
 204
 205 //________________________________________________________________
 206 void AliDielectronHFhelper::UnsetRangeUser(AliDielectronVarManager::ValueTypes type, Bool_t leg)
 207 {
 208   //
 209   // Unset range from AliDielectronVarManager
 210   //
 211   UnsetRangeUser(AliDielectronVarManager::GetValueName(type), leg);
 212 }
 213
 214 //________________________________________________________________
 215 TObjArray* AliDielectronHFhelper::CollectHistos()
 216 {
 217   //
 218   // collect histograms for all kind of pair types or sources
 219   //
 220
 221   TObjArray *collection = new TObjArray(AliDielectron::kEv1PMRot+1);
 222
 223   TObjArray *histArr = (TObjArray*) fArrPairType->Clone("tmpArr");
 224   if(!histArr) return 0x0;
 225   histArr->SetOwner(kTRUE);
 226
 227   // loop over max. available pair types
 228   for(Int_t i=0; i<AliDielectron::kEv1PMRot+1; i++) {
 229
 230     collection->AddAt(GetHistogram(AliDielectron::PairClassName(i),histArr), i);
 231
 232   }
 233
 234   // clean up the clone
 235   if(histArr) {
 236     delete histArr;
 237     histArr=0;
 238   }
 239
 240   return collection;
 241 }
 242
 243 //________________________________________________________________
 244 TH1F* AliDielectronHFhelper::GetHistogram(const char *step, TObjArray *histArr)
 245 {
 246   //
 247   // main function to recieve a single histogram
 248   // TODO: check memory
 249
 250   AliDebug(1,Form(" Step %s selected",step));
 251
 252   TObjArray *histos= 0x0;
 253   TH1F *hist       = 0x0;
 254   if(!histArr) {
 255     histos = (TObjArray*) fArrPairType->FindObject(step)->Clone("tmpArr");
 256   }
 257   else {
 258     histos = (TObjArray*) histArr->FindObject(step);
 259   }
 260
 261   if(histos) hist   = FindHistograms(histos);
 262   return hist;
 263
 264 }
 265
 266 //________________________________________________________________
 267 TH1F* AliDielectronHFhelper::FindHistograms(TObjArray *histos)
 268 {
 269   //
 270   // exclude histograms
 271   //
 272
 273   // debug
 274   // TString title    = histos->At(0)->GetTitle();
 275   // TObjArray* vars  = title.Tokenize(":");
 276   // AliDebug(1,Form(" number of cuts/vars: %d/%d",fCutLowLimits.GetNrows(),vars->GetEntriesFast()));
 277
 278   // check for missing cuts
 279   CheckCuts(histos);
 280
 281   // loop over all cuts
 282   for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++) {
 283
 284     Bool_t bFndBin = kFALSE; // exact bin found
 285     const char *cutvar   = fCutVars->At(icut)->GetName();
 286     Double_t min   = fCutLowLimits(icut);
 287     Double_t max   = fCutUpLimits(icut);
 288     AliDebug(5,Form(" Cut %d: %s [%.2f,%.2f]",icut,cutvar,min,max));
 289
 290     // loop over all histograms
 291     for(Int_t i=0; i<histos->GetEntriesFast(); i++) {
 292
 293       // continue if already empty
 294       if(!histos->At(i)) continue;
 295
 296       // collect binning from histo title
 297       TString title    = histos->At(i)->GetTitle();
 298       if(title.IsNull()) continue;
 299       AliDebug(10,Form(" histo title: %s",title.Data()));
 300
 301       TObjArray *vars  = title.Tokenize(":");
 302       for(Int_t ivar=0; ivar<vars->GetEntriesFast(); ivar++) {
 303         TString binvar = vars->At(ivar)->GetName();
 304         AliDebug(10,Form(" Check ivar %d binvar %s",ivar,binvar.Data()));
 305
 306         // check for cuts and ranges by the user
 307         if(binvar.Contains(cutvar)) {
 308           TObjArray *limits = binvar.Tokenize("#");
 309
 310           Double_t binmin = atof(limits->At(1)->GetName());
 311           Double_t binmax = atof(limits->At(2)->GetName());
 312           AliDebug(10,Form(" bin %s var %s [%.2f,%.2f]",binvar.Data(),limits->At(0)->GetName(),binmin,binmax));
 313
 314           // remove histogram from array
 315           if(binmin < min || binmax < min || binmin > max || binmax > max ) {
 316             AliDebug(10,Form(" removed, out of range min %.2f,%.2f  max %.2f,%.2f",binmin,min,binmax,max));
 317             histos->AddAt(0x0,i);
 318           }
 319           if(bFndBin && !(binmin == min && binmax == max)) {
 320               histos->AddAt(0x0,i);
 321               AliDebug(10,Form(" removed, within range min %.2f,%.2f  max %.2f,%.2f",binmin,min,binmax,max));
 322             }
 323           // clean up
 324           if(limits) delete limits;
 325
 326           // do we have found an exact bin
 327           if(binmin==min && binmax==max) bFndBin=kTRUE;
 328
 329         }
 330
 331       }
 332       // clean up
 333       if(vars)   delete vars;
 334
 335     }
 336
 337   }
 338
 339   // compress the array by removing all empty histos
 340   histos->Compress();
 341   AliDebug(1,Form(" Compression: %d histograms left",histos->GetEntriesFast()));
 342
 343   // merge histograms
 344   TH1F* hist = MergeHistos(histos);
 345   if(hist) AliDebug(1,Form(" Merging: %e histogram entries",hist->GetEntries()));
 346   return hist;
 347 }
 348
 349 //________________________________________________________________
 350 TH1F* AliDielectronHFhelper::MergeHistos(TObjArray *arr)
 351 {
 352   //
 353   // merge histos to one single histogram
 354   //
 355
 356   if(arr->GetEntriesFast()<1) { AliError(" No more histosgrams left!"); return 0x0; }
 357
 358   TH1F *final=(TH1F*) arr->At(0)->Clone();
 359   if(!final) return 0x0;
 360
 361   final->Reset("CE");
 362   final->SetTitle(""); //TODO: change in future
 363   for(Int_t i=0; i<arr->GetEntriesFast(); i++) {
 364     final->Add((TH1F*)arr->At(i));
 365   }
 366   arr->Clear();
 367
 368   return final;
 369 }
 370
 371 //________________________________________________________________
 372 void AliDielectronHFhelper::CheckCuts(TObjArray *arr)
 373 {
 374   //
 375   // Compare histo binning and cut variables. Add necessary cuts (largest limits)
 376   //
 377
 378
 379   // build array with bin variable, minimum and maximum bin values
 380   TString titleFIRST       = arr->First()->GetTitle();
 381   TString titleLAST        = arr->Last()->GetTitle();
 382   TObjArray* binvarsF  = titleFIRST.Tokenize(":#");
 383   TObjArray* binvarsL  = titleLAST.Tokenize(":#");
 384   Double_t binmin[kMaxCuts]= {0.0};
 385   Double_t binmax[kMaxCuts]= {0.0};
 386   for(Int_t ivar=0; ivar<binvarsF->GetEntriesFast(); ivar++) {
 387
 388     TString elementF=binvarsF->At(ivar)->GetName();
 389     TString elementL=binvarsL->At(ivar)->GetName();
 390     AliDebug(1,Form(" binvar %d: %s,%s",ivar,elementF.Data(),elementL.Data()));
 391
 392     switch(ivar%3) {
 393     case 0: continue; break;
 394     case 1: binmin[(int)ivar/3]=atof(elementF.Data()); break;
 395     case 2: binmax[(int)ivar/3]=atof(elementL.Data()); break;
 396     }
 397
 398     binvarsF->AddAt(0x0,ivar);
 399   }
 400   binvarsF->Compress();
 401
 402   // loop over all vars and cuts, check for missing stuff
 403   for(Int_t ivar=0; ivar<binvarsF->GetEntriesFast(); ivar++) {
 404
 405     TString binvar=binvarsF->At(ivar)->GetName();
 406     Bool_t selected=kFALSE;
 407
 408     AliDebug(1,Form(" check cuts %d %s [%.2f,%.2f]",ivar,binvar.Data(),binmin[ivar],binmax[ivar]));
 409     // loop over all cuts and check for missing stuff
 410     for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++) {
 411       if(binvar.Contains(fCutVars->At(icut)->GetName())) { selected=kTRUE; break; }
 412       //      else break;
 413     }
 414
 415     // add missing cut with max limits
 416     if(!selected) {
 417       AliWarning(Form(" Bin variable %s not covered. Add cut!",binvar.Data()));
 418       Bool_t leg = binvar.BeginsWith("Leg");
 419       if(leg) binvar.Remove(0,3);
 420       SetRangeUser(binvar.Data(),binmin[ivar],binmax[ivar], leg);
 421     }
 422
 423   }
 424
 425   // clean up
 426   if(binvarsF) delete binvarsF;
 427   if(binvarsL) delete binvarsL;
 428 }
 429
 430 //________________________________________________________________
 431 void AliDielectronHFhelper::Print(const Option_t* /*option*/) const
 432 {
 433
 434   //
 435   // Print out object contents
 436   //
 437   AliInfo(Form(" Container:               %s",fArrPairType->GetName()));
 438
 439   // pairtypes, steps and sources
 440   AliInfo(Form(" Number of filled steps:  %d",fArrPairType->GetEntries()));
 441   for(Int_t istep=0; istep<fArrPairType->GetEntriesFast(); istep++) {
 442     if(fArrPairType->At(istep))
 443       AliInfo(Form(" step %d: %s",istep,fArrPairType->At(istep)->GetName()));
 444   }
 445
 446   AliInfo(Form(" Number of histograms:    %d",((TObjArray*)fArrPairType->At(0))->GetEntriesFast()));
 447
 448   TString title       = ((TObjArray*)fArrPairType->At(0))->First()->GetTitle();
 449   TObjArray* binvars  = title.Tokenize(":");
 450   AliInfo(Form(" Number of variables:     %d",binvars->GetEntriesFast()));
 451   delete binvars;
 452
 453   TObjArray* binvars2  = title.Tokenize(":#");
 454   for(Int_t ivar=0; ivar<binvars2->GetEntriesFast(); ivar++) {
 455     if(ivar%3) continue;
 456     AliInfo(Form(" variable %.0f: %s",((Double_t)ivar)/3+1,binvars2->At(ivar)->GetName()));
 457   }
 458   delete binvars2;
 459
 460 }
 461
 462 //________________________________________________________________
 463 void AliDielectronHFhelper::PrintCuts()
 464 {
 465
 466   //
 467   // Print out object contents
 468   //
 469
 470   // loop over all cuts
 471   AliInfo(" Selected cuts:");
 472   for(Int_t icut=0; icut<fCutLowLimits.GetNrows(); icut++)
 473     AliInfo(Form(" %d: %s [%.2f,%.2f]",icut,fCutVars->At(icut)->GetName(),fCutLowLimits(icut),fCutUpLimits(icut)));
 474
 475 }
 476