PWG2/FLOW/AliSelectorFoF.cxx

   1 /* AliSelectorFoF.cxx, v1.1 01/02/2007 esimili Exp */
   2 /* derived from AliSelector.cxx,v 1.17 2006/08/31 jgrosseo Exp */
   3
   4 // The class definition in esdV0.h has been generated automatically
   5 // by the ROOT utility TTree::MakeSelector(). This class is derived
   6 // from the ROOT class TSelector. For more information on the TSelector
   7 // framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
   8
   9 // The following methods are defined in this file:
  10 //    Begin():        called everytime a loop on the tree starts,
  11 //                    a convenient place to create your histograms.
  12 //    SlaveBegin():   called after Begin(), when on PROOF called only on the
  13 //                    slave servers.
  14 //    Process():      called for each event, in this function you decide what
  15 //                    to read and fill your histograms.
  16 //    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
  17 //                    called only on the slave servers.
  18 //    Terminate():    called at the end of the loop on the tree,
  19 //                    a convenient place to draw/fit your histograms.
  20 //
  21 // To use this file, try the following session on your Tree T:
  22 //
  23 // Root > T->Process("AliSelector.C")
  24 // Root > T->Process("AliSelector.C","some options")
  25 // Root > T->Process("AliSelector.C+")
  26 //
  27
  28 #include "AliSelectorFoF.h"
  29
  30 #include <TStyle.h>
  31 #include <TSystem.h>
  32 #include <TCanvas.h>
  33 #include <TRegexp.h>
  34 #include <TTime.h>
  35 #include <TFriendElement.h>
  36 #include <TTree.h>
  37 #include <TChain.h>
  38 #include <TFile.h>
  39 #include <TTimeStamp.h>
  40 #include <TMath.h>
  41 #include <TObject.h>
  42 #include <TObjArray.h>
  43 #include <TSelector.h>
  44
  45 #include "AliLog.h"
  46 #include "AliESD.h"
  47 #include "AliESDtrack.h"
  48 #include "AliESDv0.h"
  49
  50 #include <stdio.h>
  51 #include <stdlib.h>
  52 #include <iostream>
  53 using namespace std; //required for resolving the 'cout' symbol
  54
  55 #include "AliFlowEvent.h"
  56 #include "AliFlowTrack.h"
  57 #include "AliFlowV0.h"
  58 #include "AliFlowConstants.h"
  59 #include "AliFlowSelection.h"
  60 #include "AliFlowMaker.h"
  61 #include "AliFlowAnalyser.h"
  62 #include "AliFlowWeighter.h"
  63
  64 ClassImp(AliSelectorFoF)
  65
  66 //-----------------------------------------------------------------------
  67
  68 AliSelectorFoF::AliSelectorFoF() :
  69   TSelector(),
  70   fTree(0),
  71   fESD(0),
  72   fCountFiles(0),
  73   fKineFile(0)
  74 {
  75   //
  76   // Constructor. Initialization of pointers
  77   //
  78
  79   fFlowEventFileName    = "fof_flowEvts.root" ;
  80   fFlowWgtFileName      = "fof_flowWgts.root" ;
  81   fFlowAnalysisFileName = "fof_flowAnal.root" ;
  82 }
  83
  84 //-----------------------------------------------------------------------
  85
  86 AliSelectorFoF::~AliSelectorFoF()
  87 {
  88   //
  89   // Destructor
  90   //
  91
  92  if (fTree) { fTree->ResetBranchAddresses() ; }
  93
  94  if (fESD)
  95  {
  96    delete fESD;
  97    fESD = 0;
  98  }
  99 }
 100
 101 //-----------------------------------------------------------------------
 102
 103 void AliSelectorFoF::CheckOptions()
 104 {
 105   // checks the option string for the debug flag
 106
 107   AliLog::SetClassDebugLevel(ClassName(), AliLog::kInfo);
 108
 109   TString option = GetOption();
 110
 111   if (option.Contains("moredebug"))
 112   {
 113     printf("Enabling verbose debug mode for %s\n", ClassName());
 114     AliLog::SetClassDebugLevel(ClassName(), AliLog::kDebug+1);
 115     AliInfo(Form("Called with option %s.", option.Data()));
 116   }
 117   else if (option.Contains("debug"))
 118   {
 119     printf("Enabling debug mode for %s\n", ClassName());
 120     AliLog::SetClassDebugLevel(ClassName(), AliLog::kDebug);
 121     AliInfo(Form("Called with option %s.", option.Data()));
 122   }
 123 }
 124
 125 //-----------------------------------------------------------------------
 126
 127 void AliSelectorFoF::Begin(TTree*)
 128 {
 129   // The Begin() function is called at the start of the query.
 130   // When running with PROOF Begin() is only called on the client.
 131   // The tree argument is deprecated (on PROOF 0 is passed).
 132
 133   cout << " HERE I begin !!! " << endl ; cout << endl ;
 134
 135   CheckOptions();
 136
 137   AliDebug(AliLog::kDebug, "============BEGIN===========");
 138
 139  // flags
 140  fDoNothing     = kFALSE ;  // kFALSE ;
 141  fOnFlyAnalysis = kTRUE ;  // kTRUE ;
 142  fOnFlyWeight   = kFALSE ;  // kTRUE ;
 143  fSaveFlowEvents= kTRUE ;   // kTRUE ;
 144
 145 // Maker part :
 146  fFlowMaker = new AliFlowMaker() ;
 147  // ESD Cuts
 148  fFlowMaker->SetNHitsCut(1) ;
 149  fFlowMaker->SetECut(0.01,100.) ;
 150  fFlowMaker->PrintCutList() ;
 151
 152 // Opens flow event (output) file
 153  if(fSaveFlowEvents)
 154  {
 155  // flow events file (output)
 156   fFlowfile = new TFile(fFlowEventFileName.Data(),"RECREATE") ;
 157   fFlowfile->cd() ;
 158   cout << " . Writing AliFlowEvents on  : " << fFlowEventFileName.Data()   << "  . " << endl ;
 159  }
 160
 161 // Analysis part :
 162  if(fOnFlyWeight || fOnFlyAnalysis)
 163  {
 164   cout << " . Here the flow selection ... " << endl ;
 165   cout << endl ;
 166
 167  // AliFlowSelection...
 168   fFlowSelect = new AliFlowSelection() ;
 169   // Event Cuts
 170   fFlowSelect->SetCentralityCut(-1) ;
 171   fFlowSelect->SetRunIdCut(-1) ;
 172   // R.P. calculation cuts
 173   for(int j=0;j<Flow::nHars;j++)
 174   {
 175    fFlowSelect->SetEtaCut(0., 1.1, j, 1) ;
 176    fFlowSelect->SetPtCut(0.1, 10. , j, 1);
 177   }
 178   fFlowSelect->SetConstrainCut(kTRUE) ;
 179   fFlowSelect->SetDcaGlobalCut(0.,0.1);
 180   // Correlation analysis cuts (not all of them)
 181   fFlowSelect->SetEtaPart(-1.1,1.1);
 182   fFlowSelect->SetPtPart(0.1,10.);
 183   fFlowSelect->SetConstrainablePart(kTRUE);
 184   fFlowSelect->SetDcaGlobalPart(0.,0.1);
 185   // V0 analysis cuts (not all of them ... they are useless anyway)
 186   fFlowSelect->SetV0Mass(0.4875,0.5078) ;        // Mk0 = 0.49765
 187   fFlowSelect->SetV0SideBands(0.08) ;
 188   fFlowSelect->SetV0Pt(0.1,10.) ;
 189   fFlowSelect->SetV0Eta(-2.1,2.1) ;
 190   // print list :
 191   cout << " . Selection for R.P. calculation: " << endl ;
 192   fFlowSelect->PrintSelectionList() ;
 193   cout << " . Selection for correlation analysis: " << endl ;
 194   fFlowSelect->PrintList() ;
 195   cout << " . Selection for V0 analysis: " << endl ;
 196   fFlowSelect->PrintV0List() ;
 197
 198 // Weight part :
 199   if(fOnFlyWeight)
 200   {
 201    cout << " . Here the flow weighter ... " << endl ;
 202    cout << endl ;
 203
 204   // AliFlowWeighter...
 205    fFlowWeighter = new AliFlowWeighter(fFlowSelect) ;
 206
 207   // Wgt file (create)
 208    fFlowWgtFile = new TFile(fFlowWgtFileName.Data(),"RECREATE") ;
 209    fFlowWeighter->SetWgtFile(fFlowWgtFile) ;
 210    cout << " . Writing Wgt Histograms on  : " << fFlowWeighter->GetWgtFileName()  << "  . " << endl ;
 211    fFlowWeighter->Init() ;
 212   }
 213
 214   if(fOnFlyAnalysis)
 215   {
 216    cout << " . Here the flow analysis ... " << endl ;
 217    cout << endl ;
 218
 219   // AliFlowAnalyser...
 220    fFlowAnal = new AliFlowAnalyser(fFlowSelect) ;
 221
 222   // Wgt file (read)
 223    //TString wgtFileName = "fFlowPhiWgt.root" ;
 224    //TFile* wgtFile = new TFile(wgtFileName.Data(),"READ");
 225    //fFlowAnal->FillWgtArrays(wgtFile) ; // fix this !!!
 226    //cout << " . Weights from  : " << fFlowAnal->GetWgtFileName() << "  . " << endl ;
 227
 228   // analysis file (output)
 229    fFlowAnal->SetHistFileName(fFlowAnalysisFileName.Data()) ;
 230    cout << " . Writing Analysis Histograms on  : " << fFlowAnal->GetHistFileName()   << "  . " << endl ;
 231
 232   // Analysis settings
 233    fFlowAnal->SetFlowForV0(kFALSE) ;         // default kTRUE.
 234    fFlowAnal->SetEtaSub() ;            // default kFALSE
 235    //fFlowAnal->SetV1Ep1Ep2() ;        // default kFALSE.
 236    fFlowAnal->SetShuffle() ;           // default kFALSE. shuffles track array
 237    //fFlowAnal->SetRedoWgt();             // default kFALSE. recalculates phiWgt (even if phiWgt file is already there)
 238    fFlowAnal->SetUsePhiWgt(kFALSE) ;   // default kTRUE if phiWgt file is there, kFALSE if not (& phiWgt file is created)
 239    fFlowAnal->SetUseOnePhiWgt() ; // or // fFlowAnal->SetUseFirstLastPhiWgt() ; // uses 1 or 3 wgt histograms (default is 1)
 240    fFlowAnal->SetUseBayWgt(kFALSE) ;   // default kFALSE. uses bayesian weights in P.id.
 241    //fFlowAnal->SetUsePtWgt();  // default kFALSE. uses pT as a weight for RP determination
 242    //fFlowAnal->SetUseEtaWgt(); // default kFALSE. uses eta as a weight for RP determination
 243    fFlowAnal->Init() ;
 244   }
 245  }
 246
 247  cout << "                ... init is done " << endl ;
 248 }
 249
 250 //-----------------------------------------------------------------------
 251
 252 void AliSelectorFoF::SlaveBegin(TTree* tree)
 253 {
 254   // The SlaveBegin() function is called after the Begin() function.
 255   // When running with PROOF SlaveBegin() is called on each slave server.
 256   // The tree argument is deprecated (on PROOF 0 is passed).
 257
 258   CheckOptions();
 259
 260   AliDebug(AliLog::kDebug, "=======SLAVEBEGIN========");
 261   AliDebug(AliLog::kDebug, Form("Hostname: %s", gSystem->HostName()));
 262   AliDebug(AliLog::kDebug, Form("Time: %s", gSystem->Now().AsString()));
 263
 264   if (tree != 0) { Init(tree) ; }
 265 }
 266
 267 //-----------------------------------------------------------------------
 268
 269 void AliSelectorFoF::Init(TTree *tree)
 270 {
 271   // The Init() function is called when the selector needs to initialize
 272   // a new tree or chain. Typically here the branch addresses of the tree
 273   // will be set. It is normaly not necessary to make changes to the
 274   // generated code, but the routine can be extended by the user if needed.
 275   // Init() will be called many times when running with PROOF.
 276
 277   AliDebug(AliLog::kDebug, "=========Init==========");
 278
 279   fTree = tree;
 280
 281   if (fTree == 0)
 282   {
 283    AliDebug(AliLog::kError, "ERROR: tree argument is 0.");
 284    return;
 285   }
 286
 287   // Set branch address
 288   fTree->SetBranchAddress("ESD", &fESD);
 289   if (fESD != 0) { AliDebug(AliLog::kInfo, "INFO: Found ESD branch in chain.") ; }
 290 }
 291
 292 //-----------------------------------------------------------------------
 293
 294 Bool_t AliSelectorFoF::Notify()
 295 {
 296   // The Notify() function is called when a new file is opened. This
 297   // can be either for a new TTree in a TChain or when when a new TTree
 298   // is started when using PROOF. Typically here the branch pointers
 299   // will be retrieved. It is normaly not necessary to make changes
 300   // to the generated code, but the routine can be extended by the
 301   // user if needed.
 302
 303   AliDebug(AliLog::kDebug, "=========NOTIFY==========");
 304   AliDebug(AliLog::kDebug, Form("Hostname: %s", gSystem->HostName()));
 305   AliDebug(AliLog::kDebug, Form("Time: %s", TTimeStamp(time(0)).AsString()));
 306
 307   ++fCountFiles;
 308   if (fTree)
 309   {
 310     TFile *f = fTree->GetCurrentFile();
 311     AliDebug(AliLog::kInfo, Form("Processing %d. file %s", fCountFiles, f->GetName()));
 312   }
 313   else
 314   {
 315     AliDebug(AliLog::kError, "fTree not available");
 316   }
 317
 318   DeleteKinematicsFile();
 319
 320   return kTRUE;
 321 }
 322
 323 //-----------------------------------------------------------------------
 324
 325 Bool_t AliSelectorFoF::Process(Long64_t entry)
 326 {
 327   // The Process() function is called for each entry in the tree (or possibly
 328   // keyed object in the case of PROOF) to be processed. The entry argument
 329   // specifies which entry in the currently loaded tree is to be processed.
 330   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
 331   // to read either all or the required parts of the data. When processing
 332   // keyed objects with PROOF, the object is already loaded and is available
 333   // via the fObject pointer.
 334   //
 335   // This function should contain the "body" of the analysis. It can contain
 336   // simple or elaborate selection criteria, run algorithms on the data
 337   // of the event and typically fill histograms.
 338
 339   // WARNING when a selector is used with a TChain, you must use
 340   //  the pointer to the current TTree to call GetEntry(entry).
 341   //  The entry is always the local entry number in the current tree.
 342   //  Assuming that fTree is the pointer to the TChain being processed,
 343   //  use fTree->GetTree()->GetEntry(entry).
 344
 345   AliDebug(AliLog::kDebug, Form("=========PROCESS========== Entry %lld", entry));
 346
 347   if(!fTree)
 348   {
 349    AliDebug(AliLog::kError, "ERROR: fTree is 0.") ;
 350    return kFALSE ;
 351   }
 352
 353   fEventNumber = entry ;
 354   fTree->GetTree()->GetEntry(fEventNumber) ;
 355
 356   if(fESD) { AliDebug(AliLog::kDebug, Form("ESD: We have %d tracks.", fESD->GetNumberOfTracks())); }
 357   cout << " event !!! " << entry << endl ;
 358
 359   fRunID = fESD->GetRunNumber() ;
 360   fEventNumber = fESD->GetEventNumber() ;
 361   fNumberOfTracks = fESD->GetNumberOfTracks() ;
 362   fNumberOfV0s = fESD->GetNumberOfV0s() ;
 363
 364   cout << " *evt n. " << fEventNumber << " (run " << fRunID << ") " << endl ;
 365   cout << "  tracks: " << fNumberOfTracks << " ,   v0s " << fNumberOfV0s << endl ;
 366
 367  // Dummy Loop
 368   if(fDoNothing)  { cout << " ... doing nothing ... " << endl ; return kTRUE; }
 369
 370  // Instantiate a new AliFlowEvent
 371   cout << " filling the flow event :| " << endl ;
 372   fFlowEvent = fFlowMaker->FillFlowEvent(fESD) ;
 373   if(!fFlowEvent) { cout << "!something bad occurred" << endl ; return kFALSE ; }
 374   else            { cout << " event filled :) " << entry << "                # t.G. ! #       "<< endl ; }
 375
 376  // Saves the AliFlowEvent
 377   if(fSaveFlowEvents)
 378   {
 379    cout << " saving flow event :| " << endl ;
 380    TString strID = "" ; strID += entry ;
 381    fFlowfile->cd() ; fFlowEvent->Write(strID.Data()) ;
 382    cout << " saved :) " << strID.Data() << "                # t.g. ! #       " << endl ; cout << endl ;
 383   }
 384
 385  // Weights
 386   if(fOnFlyWeight)
 387   {
 388    cout << " weightening :| " << endl ;
 389    bool done = fFlowWeighter->WeightEvent(fFlowEvent) ;
 390    if(done) { cout << " weights done :) " << entry << "                # t.h.v. ! #       " << endl ; }
 391   }
 392
 393  // On-Fly Analysis
 394   if(fOnFlyAnalysis)
 395   {
 396    cout << " doing analysis :| " << endl ;
 397    bool done = fFlowAnal->Analyse(fFlowEvent) ;
 398    if(done) { cout << " analysis done :) " << entry << "                # t.h.v. ! #       " << endl ; }
 399   }
 400
 401   return kTRUE;
 402 }
 403
 404 //-----------------------------------------------------------------------
 405
 406 void AliSelectorFoF::SlaveTerminate()
 407 {
 408   // The SlaveTerminate() function is called after all entries or objects
 409   // have been processed. When running with PROOF SlaveTerminate() is called
 410   // on each slave server.
 411
 412   AliDebug(AliLog::kDebug, "=======SLAVETERMINATE=======");
 413
 414   DeleteKinematicsFile();
 415 }
 416
 417 //-----------------------------------------------------------------------
 418
 419 void AliSelectorFoF::Terminate()
 420 {
 421   // The Terminate() function is the last function to be called during
 422   // a query. It always runs on the client, it can be used to present
 423   // the results graphically or save the results to file.
 424
 425   AliDebug(AliLog::kDebug, "=========TERMINATE==========");
 426
 427   cout << " Finished ... " << endl ;
 428
 429   if(fDoNothing) { cout << "             ... & nothing was done ! " << endl ; cout << endl ; return ; }
 430
 431   cout << "  nTracks:  " << fFlowMaker->GetNgoodTracks() << endl ;
 432   cout << "  nV0s:  " << fFlowMaker->GetNgoodV0s()  << endl ;
 433   cout << "  nTracks (|eta|<0.5):  " << fFlowMaker->GetNgoodTracksEta() << endl ;
 434   cout << "  nTracks+:  " << fFlowMaker->GetNposiTracks() << endl ;
 435   cout << "  nTracks-:  " << fFlowMaker->GetNnegaTracks() << endl ;
 436   cout << "  nTracks unconstrained:  " << fFlowMaker->GetNunconstrained() << endl ;
 437   cout << "  Bayesian :  " ;
 438   for(int ii=0;ii<5;ii++) { cout << fFlowMaker->GetBayesianNorm(ii) << "   " ; }
 439   cout << " . " << endl ;
 440
 441   if(fOnFlyWeight) { fFlowWeighter->Finish() ; delete fFlowWeighter ; }
 442
 443   if(fOnFlyAnalysis)
 444   {
 445    cout << endl ; cout << "Particles normalized abundance:" << endl;  // shows the bayesian vector
 446    fFlowAnal->Resolution() ; cout << endl;
 447    fFlowAnal->PrintRunBayesian() ; cout << endl;
 448    fFlowAnal->Finish() ;  cout << endl;
 449    delete fFlowAnal ;
 450   }
 451
 452   if(fSaveFlowEvents) { fFlowfile->Close() ; cout << " file closed . " << endl ; }
 453   delete fFlowMaker ;
 454
 455   delete fFlowSelect ;
 456
 457   cout << endl ;
 458   return ;
 459 }
 460
 461 //-----------------------------------------------------------------------
 462
 463 TTree* AliSelectorFoF::GetKinematics()
 464 {
 465   // Returns kinematics tree corresponding to current ESD active in fTree
 466   // Loads the kinematics from the kinematics file, the file is identified by replacing "AliESDs" to
 467   // "Kinematics" in the file path of the ESD file. This is a hack, to be changed!
 468
 469   if (!fKineFile)
 470   {
 471     if(!fTree->GetCurrentFile()) { return 0 ; }
 472
 473     TString fileName(fTree->GetCurrentFile()->GetName());
 474     fileName.ReplaceAll("AliESDs", "Kinematics");
 475
 476     // temporary workaround for PROOF bug #18505
 477     fileName.ReplaceAll("#Kinematics.root#Kinematics.root", "#Kinematics.root");
 478
 479     AliDebug(AliLog::kInfo, Form("Opening %s", fileName.Data()));
 480
 481     fKineFile = TFile::Open(fileName);
 482     if(!fKineFile) { return 0 ; }
 483   }
 484
 485   return dynamic_cast<TTree*> (fKineFile->Get(Form("Event%d/TreeK", fTree->GetTree()->GetReadEntry())));
 486 }
 487
 488 //-----------------------------------------------------------------------
 489
 490 void AliSelectorFoF::DeleteKinematicsFile()
 491 {
 492   //
 493   // Closes the kinematics file and deletes the pointer.
 494   //
 495
 496   if (fKineFile)
 497   {
 498     fKineFile->Close();
 499     delete fKineFile;
 500     fKineFile = 0;
 501   }
 502 }
 503
 504 //-----------------------------------------------------------------------
 505
 506