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