PWGHF/correlationHF/AliDxHFEParticleSelection.cxx

   1 // $Id$
   2
   3 //**************************************************************************
   4 //* This file is property of and copyright by the ALICE Project            *
   5 //* ALICE Experiment at CERN, All rights reserved.                         *
   6 //*                                                                        *
   7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
   8 //*                  Sedat Altinpinar <Sedat.Altinpinar@cern.ch>           *
   9 //*                  Hege Erdal       <hege.erdal@gmail.com>               *
  10 //*                                                                        *
  11 //* Permission to use, copy, modify and distribute this software and its   *
  12 //* documentation strictly for non-commercial purposes is hereby granted   *
  13 //* without fee, provided that the above copyright notice appears in all   *
  14 //* copies and that both the copyright notice and this permission notice   *
  15 //* appear in the supporting documentation. The authors make no claims     *
  16 //* about the suitability of this software for any purpose. It is          *
  17 //* provided "as is" without express or implied warranty.                  *
  18 //**************************************************************************
  19
  20 /// @file   AliDxHFEParticleSelection.cxx
  21 /// @author Sedat Altinpinar, Hege Erdal, Matthias Richter
  22 /// @date   2012-03-19
  23 /// @brief  Base class for particle selection
  24 ///
  25
  26 #include "AliDxHFEParticleSelection.h"
  27 #include "AliLog.h"
  28 #include "AliVEvent.h"
  29 #include "AliVParticle.h"
  30 #include "TObjArray.h"
  31 #include "TList.h"
  32 #include "TMath.h"
  33 #include "TH1D.h"
  34 #include "TH2D.h"
  35 #include "THnSparse.h"
  36 #include "AliReducedParticle.h"
  37 #include "TFile.h"
  38 #include "TPRegexp.h" // for TStringToken
  39 #include <iostream>
  40 #include <cerrno>
  41 #include <memory>
  42
  43 using namespace std;
  44
  45 /// ROOT macro for the implementation of ROOT specific class methods
  46 ClassImp(AliDxHFEParticleSelection)
  47
  48 AliDxHFEParticleSelection::AliDxHFEParticleSelection(const char* name, const char* opt)
  49   : TNamed(name?name:"AliDxHFEParticleSelection", name?name:"AliDxHFEParticleSelection")
  50   , fOption(opt)
  51   , fSelectedTracks(NULL)
  52   , fControlObjects(NULL)
  53   , fhEventControl(NULL)
  54   , fhTrackControl(NULL)
  55   , fUseMC(0)
  56   , fVerbosity(0)
  57   , fDimThn(-1)
  58   , fParticleProperties(NULL)
  59   , fSystem(0)
  60 {
  61   // constructor
  62   //
  63   //
  64   //
  65   //
  66 }
  67
  68 const char* AliDxHFEParticleSelection::fgkEventControlBinNames[]={
  69   "nEventsAll",
  70   "nEventsSelected",
  71   "nEventsWParticle"
  72 };
  73
  74 const char* AliDxHFEParticleSelection::fgkTrackControlBinNames[]={
  75   "nTrackAll",
  76   "nTrackSelected",
  77 };
  78
  79 AliDxHFEParticleSelection::~AliDxHFEParticleSelection()
  80 {
  81   // destructor
  82   if (fParticleProperties) delete[] fParticleProperties;
  83   fParticleProperties=NULL;
  84   if (fSelectedTracks) delete fSelectedTracks;
  85   fSelectedTracks=NULL;
  86   if (fControlObjects) delete fControlObjects;
  87   fControlObjects=NULL;
  88   fhEventControl=NULL;
  89   fhTrackControl=NULL;
  90 }
  91
  92 int AliDxHFEParticleSelection::Init()
  93 {
  94   //
  95   // Init part sel. Calls InitControlObjects()
  96   //
  97
  98   InitControlObjects();
  99   return 0;
 100 }
 101
 102
 103 int AliDxHFEParticleSelection::InitControlObjects()
 104 {
 105   //
 106   // init control objects
 107   // TODO: Change to private now that have Init()?
 108   //
 109   if (fVerbosity>0) {
 110     AliInfo("Setting up control objects");
 111   }
 112
 113   /// init the control objects, can be overloaded by childs which should
 114   /// call AliDxHFEParticleSelection::InitControlObjects() explicitly
 115   std::auto_ptr<TH1D> hEventControl(new TH1D("hEventControl", "hEventControl", 10, 0, 10));
 116   std::auto_ptr<TH1D> hTrackControl(new TH1D("hTrackControl", "hTrackControl", 10, 0, 10));
 117
 118   fhEventControl=hEventControl.release();
 119   for (int iLabel=0; iLabel<kNEventPropertyLabels; iLabel++)
 120     fhEventControl->GetXaxis()->SetBinLabel(iLabel+1, fgkEventControlBinNames[iLabel]);
 121   AddControlObject(fhEventControl);
 122   fhTrackControl=hTrackControl.release();
 123   for (int iLabel=0; iLabel<kNTrackPropertyLabels; iLabel++)
 124     fhTrackControl->GetXaxis()->SetBinLabel(iLabel+1, fgkTrackControlBinNames[iLabel]);
 125   AddControlObject(fhTrackControl);
 126
 127   return 0;
 128 }
 129
 130
 131 TH1* AliDxHFEParticleSelection::CreateControlHistogram(const char* name,
 132                                                        const char* title,
 133                                                        int nBins,
 134                                                        double min,
 135                                                        double max,
 136                                                        const char** binLabels) const
 137 {
 138   /// create control histogram
 139   std::auto_ptr<TH1> h(new TH1D(name, title, nBins, min, max));
 140   if (!h.get()) return NULL;
 141   if (binLabels) {
 142   for (int iLabel=0; iLabel<nBins; iLabel++) {
 143     h->GetXaxis()->SetBinLabel(iLabel+1, binLabels[iLabel]);
 144   }
 145   }
 146
 147   return h.release();
 148 }
 149
 150
 151 TH2* AliDxHFEParticleSelection::CreateControl2DHistogram(const char* name,
 152                                                          const char* title,
 153                                                          double* nBins,
 154                                                          const char* xaxis,
 155                                                          const char* yaxis
 156                                                          ) const
 157 {
 158   /// create control 2D histogram. Requires as input:
 159   // name = name of histogram
 160   // title = title of histogram
 161   // nBins (array with 6 elements) containing apropriate binning and range for x and y axis
 162   // xaxis = title of x axis
 163   // yaxis = title of y axis
 164
 165   std::auto_ptr<TH2> h(new TH2D(name, title, (Int_t)nBins[0], nBins[1], nBins[2], (Int_t)nBins[3], nBins[4],nBins[5]));
 166   if (!h.get()) return NULL;
 167   h->GetXaxis()->SetTitle(xaxis);
 168   h->GetYaxis()->SetTitle(yaxis);
 169
 170   return h.release();
 171 }
 172
 173 THnSparse* AliDxHFEParticleSelection::CreateControlTHnSparse(const char* name,
 174                                                              int thnSize,
 175                                                              int* thnBins,
 176                                                              double* thnMin,
 177                                                              double* thnMax,
 178                                                              const char** binLabels) const
 179 {
 180   //
 181   // Creates THnSparse.
 182   //
 183
 184   AliInfo("Setting up THnSparse");
 185
 186   std::auto_ptr<THnSparseF> th(new THnSparseF(name, name, thnSize, thnBins, thnMin, thnMax));
 187   if (th.get()==NULL) {
 188     return NULL;
 189   }
 190   for (int iLabel=0; iLabel<thnSize; iLabel++) {
 191     th->GetAxis(iLabel)->SetTitle(binLabels[iLabel]);
 192
 193   }
 194  return th.release();
 195
 196 }
 197
 198 THnSparse* AliDxHFEParticleSelection::DefineTHnSparse()
 199 {
 200   //
 201   // Defines the THnSparse. For now, only calls CreatControlTHnSparse
 202
 203   //TODO: Should make it more general. Or maybe one can use this here, and skip in PartSelEl?
 204
 205   // here is the only place to change the dimension
 206   const int thnSize = 3;
 207   InitTHnSparseArray(thnSize);
 208
 209   const double Pi=TMath::Pi();
 210   TString name;
 211   //                                0    1       2
 212   //                                Pt   Phi    Eta
 213   int         thnBins [thnSize] = { 1000,  200, 500};
 214   double      thnMin  [thnSize] = {    0,    0, -1.};
 215   double      thnMax  [thnSize] = {  100, 2*Pi,  1.};
 216   const char* thnNames[thnSize] = { "Pt","Phi","Eta"};
 217
 218   name.Form("%s info", GetName());
 219
 220   return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
 221 }
 222
 223 int AliDxHFEParticleSelection::AddControlObject(TObject* pObj)
 224 {
 225   /// add control object to list, the base class becomes owner of the object
 226   if (!pObj) return -EINVAL;
 227   if (!fControlObjects) {
 228     fControlObjects=new TList;
 229     if (!fControlObjects) return -ENOMEM;
 230     fControlObjects->SetOwner();
 231   }
 232   if (fControlObjects->FindObject(pObj->GetName())) {
 233     AliError(Form("ignoring duplicate object '%s' of type %s", pObj->GetName(), pObj->ClassName()));
 234     return -EEXIST;
 235   }
 236   if (GetVerbosity()>0) {
 237     AliInfo(Form("Adding object '%s' of type %s",pObj->GetName(),pObj->ClassName()));
 238   }
 239   fControlObjects->Add(pObj);
 240   return 0;
 241 }
 242
 243 int AliDxHFEParticleSelection::HistogramEventProperties(int bin)
 244 {
 245   /// histogram event properties
 246   if (!fControlObjects) return 0;
 247
 248   // TODO: use enums for the bins of the control histogram
 249   // for now: 0=all, 1=events with D0s, 2=events with correlated D0s
 250   fhEventControl->Fill(bin);
 251   return 0;
 252 }
 253
 254 int AliDxHFEParticleSelection::HistogramParticleProperties(AliVParticle* p, int selected)
 255 {
 256   /// histogram particle properties
 257
 258   if (!p) return -EINVAL;
 259   if (!fControlObjects) return 0;
 260
 261   // TODO: use enums for the bins of the control histogram
 262   fhTrackControl->Fill(kTrackAll);
 263   if (selected) fhTrackControl->Fill(kTrackSel);
 264   return 0;
 265 }
 266
 267 int AliDxHFEParticleSelection::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
 268 {
 269   // fill the data array from the particle data
 270   if (!data) return -EINVAL;
 271   int i=0;
 272   if (dimension!=GetDimTHnSparse()) {
 273     // TODO: think about filling only the available data and throwing a warning
 274     return -ENOSPC;
 275   }
 276   data[i++]=p->Pt();
 277   data[i++]=p->Phi();
 278   data[i++]=p->Eta();
 279   return i;
 280 }
 281
 282 TObjArray* AliDxHFEParticleSelection::Select(const AliVEvent* pEvent)
 283 {
 284   /// create selection from 'Tracks' member of the event,
 285   /// array contains only pointers but does not own the objects
 286   /// object array needs to be deleted by caller
 287   if (!pEvent) return NULL;
 288   TObjArray* selectedTracks=new TObjArray;
 289   if (!selectedTracks) return NULL;
 290   selectedTracks->SetOwner(kFALSE); // creating new track objects below
 291   int nofTracks=pEvent->GetNumberOfTracks();
 292   for (int itrack=0; itrack<nofTracks; itrack++) {
 293     AliVParticle* track=pEvent->GetTrack(itrack);
 294     int selectionCode=IsSelected(track,pEvent);
 295     HistogramParticleProperties(track, selectionCode);
 296     if (selectionCode==0) continue;
 297     selectedTracks->Add(CreateParticle(track));
 298   }
 299   return selectedTracks;
 300 }
 301
 302 TObjArray* AliDxHFEParticleSelection::Select(TObjArray* pParticles, const AliVEvent* pEvent)
 303 {
 304   /// create selection from the array of particles,
 305   /// array contains only pointers but does not own the objects
 306   /// object array needs to be deleted by caller
 307   if (!pParticles) return NULL;
 308   TObjArray* selectedTracks=new TObjArray;
 309   if (!selectedTracks) return NULL;
 310   selectedTracks->SetOwner(kFALSE); // creating new track objects below
 311   TIter next(pParticles);
 312   TObject* pObj=NULL;
 313   while ((pObj=next())) {
 314     AliVParticle* track=dynamic_cast<AliVParticle*>(pObj);
 315     if (!track) continue;
 316     int selectionCode=IsSelected(track, pEvent);
 317     HistogramParticleProperties(track, selectionCode);
 318     if (selectionCode ==0) continue;
 319     selectedTracks->Add(CreateParticle(track));
 320   }
 321   return selectedTracks;
 322 }
 323
 324 int AliDxHFEParticleSelection::CheckAndAdd(AliVParticle* /*p*/)
 325 {
 326   /// check and add track to internal array
 327   /// TODO: check if needed
 328   return -ENOSYS;
 329 }
 330
 331 int AliDxHFEParticleSelection::IsSelected(AliVParticle* /*p*/, const AliVEvent* /*e*/)
 332 {
 333   /// check particle if it passes the selection criteria
 334   /// childs can overload, by default all tracks are selected
 335   return 1;
 336 }
 337
 338 void AliDxHFEParticleSelection::AliDxHFEParticleSelection::Clear(Option_t * /*option*/)
 339 {
 340   /// inherited from TObject: cleanup
 341 }
 342
 343 void AliDxHFEParticleSelection::Print(Option_t */*option*/) const
 344 {
 345   /// inherited from TObject: print info
 346   cout << "====================================================================" << endl;
 347   TNamed::Print();
 348   if (fControlObjects) fControlObjects->Print();
 349 }
 350
 351 void AliDxHFEParticleSelection::SaveAs(const char* filename, Option_t */*option*/) const
 352 {
 353   /// inherited from TObject: save selection criteria
 354   TString fileoption;
 355   // TODO: options recreate
 356   fileoption="RECREATE";
 357   //else fileoption="UPDATE";
 358
 359   std::auto_ptr<TFile> output(TFile::Open(filename,fileoption));
 360   if (!output.get() || output->IsZombie()) {
 361     AliError(Form("can not open file %s from writing", filename));
 362     return;
 363   }
 364   output->cd();
 365   if (fControlObjects) fControlObjects->Write();
 366   output->Close();
 367 }
 368
 369 void AliDxHFEParticleSelection::Draw(Option_t* /*option*/)
 370 {
 371   /// inherited from TObject: draw content
 372
 373   // TODO: implement drawing code
 374   // - create canvas objects
 375   // - plot internal objects
 376   // - optionally save canvases to file
 377   //
 378   // It might be appropriate to have another Draw function taking a
 379   // TList as argument and implementing the actual drawing. If this
 380   // function is 'static', it can be used stand-alone also from macros
 381 }
 382
 383 TObject* AliDxHFEParticleSelection::FindObject(const char* name) const
 384 {
 385   /// inherited from TObject: find object by name
 386
 387   if (fControlObjects) {
 388     return fControlObjects->FindObject(name);
 389   }
 390   return NULL;
 391 }
 392
 393 TObject* AliDxHFEParticleSelection::FindObject(const TObject* obj) const
 394 {
 395   /// inherited from TObject: find object by pointer
 396   if (fControlObjects) {
 397     return fControlObjects->FindObject(obj);
 398   }
 399   return NULL;
 400 }
 401
 402 AliVParticle *AliDxHFEParticleSelection::CreateParticle(AliVParticle* track)
 403 {
 404   // Creating object with reduced particle properties
 405   AliReducedParticle *part = new AliReducedParticle(track->Eta(), track->Phi(), track->Pt(), track->Charge(), 0);
 406
 407   return part;
 408
 409 }
 410
 411 int AliDxHFEParticleSelection::ParseArguments(const char* arguments)
 412 {
 413   // parse arguments and set internal flags
 414   const char* key=NULL;
 415   const TString delimiter(" ");
 416   TStringToken token(arguments, delimiter);
 417   while (token.NextToken()) {
 418     TString argument=token;
 419     key="system=";
 420     if (argument.BeginsWith(key)) {
 421       argument.ReplaceAll(key, "");
 422       if (argument.CompareTo("pp")==0) fSystem=0;
 423       else if (argument.CompareTo("Pb-Pb")==0) fSystem=1;
 424       else if (argument.CompareTo("p-Pb")==0) fSystem=2;
 425       else {
 426         AliWarning(Form("can not set collision system, unknown parameter '%s'", argument.Data()));
 427         // TODO: check what makes sense
 428         fSystem=0;
 429       }
 430       continue;
 431     }
 432
 433     AliWarning(Form("unknown argument '%s'", argument.Data()));
 434   }
 435
 436   return 0;
 437 }