PWGCF/FLOW/Attic/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<AliFlowConstants::kHars;j++)
 174   {
 175    fFlowSelect->SetEtaCut(0., 0.9, 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(-0.9,0.9);
 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.1) ;
 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 = "fof_flowWgts.root" ;
 224    TFile* wgtFile = new TFile(wgtFileName.Data(),"READ");
 225    if(!wgtFile || wgtFile->IsZombie()) { cout << " . NO phi Weights . " << endl ;}
 226    else { fFlowAnal->FillWgtArrays(wgtFile) ; 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() ;         // default kTRUE.
 234    fFlowAnal->SetSub(1) ; //eta  // ->SetChrSub() ; //charge  // ->SetRndSub() ; //random (default)
 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->SetCustomRespFunc();    // default kFALSE. uses the combined response function from the ESD
 244    fFlowAnal->Init() ;
 245   }
 246  }
 247
 248  cout << "                ... init is done " << endl ;
 249 }
 250
 251 //-----------------------------------------------------------------------
 252
 253 void AliSelectorFoF::SlaveBegin(TTree* tree)
 254 {
 255   // The SlaveBegin() function is called after the Begin() function.
 256   // When running with PROOF SlaveBegin() is called on each slave server.
 257   // The tree argument is deprecated (on PROOF 0 is passed).
 258
 259   CheckOptions();
 260
 261   AliDebug(AliLog::kDebug, "=======SLAVEBEGIN========");
 262   AliDebug(AliLog::kDebug, Form("Hostname: %s", gSystem->HostName()));
 263   AliDebug(AliLog::kDebug, Form("Time: %s", gSystem->Now().AsString()));
 264
 265   if (tree != 0) { Init(tree) ; }
 266 }
 267
 268 //-----------------------------------------------------------------------
 269
 270 void AliSelectorFoF::Init(TTree *tree)
 271 {
 272   // The Init() function is called when the selector needs to initialize
 273   // a new tree or chain. Typically here the branch addresses of the tree
 274   // will be set. It is normaly not necessary to make changes to the
 275   // generated code, but the routine can be extended by the user if needed.
 276   // Init() will be called many times when running with PROOF.
 277
 278   AliDebug(AliLog::kDebug, "=========Init==========");
 279
 280   fTree = tree;
 281
 282   if (fTree == 0)
 283   {
 284    AliDebug(AliLog::kError, "ERROR: tree argument is 0.");
 285    return;
 286   }
 287
 288   // Set branch address
 289   fTree->SetBranchAddress("ESD", &fESD);
 290   if (fESD != 0) { AliDebug(AliLog::kInfo, "INFO: Found ESD branch in chain.") ; }
 291 }
 292
 293 //-----------------------------------------------------------------------
 294
 295 Bool_t AliSelectorFoF::Notify()
 296 {
 297   // The Notify() function is called when a new file is opened. This
 298   // can be either for a new TTree in a TChain or when when a new TTree
 299   // is started when using PROOF. Typically here the branch pointers
 300   // will be retrieved. It is normaly not necessary to make changes
 301   // to the generated code, but the routine can be extended by the
 302   // user if needed.
 303
 304   AliDebug(AliLog::kDebug, "=========NOTIFY==========");
 305   AliDebug(AliLog::kDebug, Form("Hostname: %s", gSystem->HostName()));
 306   AliDebug(AliLog::kDebug, Form("Time: %s", TTimeStamp(time(0)).AsString()));
 307
 308   ++fCountFiles;
 309   if (fTree)
 310   {
 311     TFile *f = fTree->GetCurrentFile();
 312     AliDebug(AliLog::kInfo, Form("Processing %d. file %s", fCountFiles, f->GetName()));
 313   }
 314   else
 315   {
 316     AliDebug(AliLog::kError, "fTree not available");
 317   }
 318
 319   DeleteKinematicsFile();
 320
 321   return kTRUE;
 322 }
 323
 324 //-----------------------------------------------------------------------
 325
 326 Bool_t AliSelectorFoF::Process(Long64_t entry)
 327 {
 328   // The Process() function is called for each entry in the tree (or possibly
 329   // keyed object in the case of PROOF) to be processed. The entry argument
 330   // specifies which entry in the currently loaded tree is to be processed.
 331   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
 332   // to read either all or the required parts of the data. When processing
 333   // keyed objects with PROOF, the object is already loaded and is available
 334   // via the fObject pointer.
 335   //
 336   // This function should contain the "body" of the analysis. It can contain
 337   // simple or elaborate selection criteria, run algorithms on the data
 338   // of the event and typically fill histograms.
 339
 340   // WARNING when a selector is used with a TChain, you must use
 341   //  the pointer to the current TTree to call GetEntry(entry).
 342   //  The entry is always the local entry number in the current tree.
 343   //  Assuming that fTree is the pointer to the TChain being processed,
 344   //  use fTree->GetTree()->GetEntry(entry).
 345
 346   AliDebug(AliLog::kDebug, Form("=========PROCESS========== Entry %lld", entry));
 347
 348   if(!fTree)
 349   {
 350    AliDebug(AliLog::kError, "ERROR: fTree is 0.") ;
 351    return kFALSE ;
 352   }
 353
 354   fEventNumber = entry ;
 355   fTree->GetTree()->GetEntry(fEventNumber) ;
 356
 357   if(fESD) { AliDebug(AliLog::kDebug, Form("ESD: We have %d tracks.", fESD->GetNumberOfTracks())); }
 358   cout << " event !!! " << entry << endl ;
 359
 360   fRunID = fESD->GetRunNumber() ;
 361   //  fEventNumber = fESD->GetEventNumber() ;
 362   fEventNumber = -1 ;
 363   fNumberOfTracks = fESD->GetNumberOfTracks() ;
 364   fNumberOfV0s = fESD->GetNumberOfV0s() ;
 365
 366   cout << " *evt n. " << fEventNumber << " (run " << fRunID << ") " << endl ;
 367   cout << "  tracks: " << fNumberOfTracks << " ,   v0s " << fNumberOfV0s << endl ;
 368
 369  // Dummy Loop
 370   if(fDoNothing)  { cout << " ... doing nothing ... " << endl ; return kTRUE; }
 371
 372  // Instantiate a new AliFlowEvent
 373   cout << " filling the flow event :| " << endl ;
 374   fFlowEvent = fFlowMaker->FillFlowEvent(fESD) ;
 375   if(!fFlowEvent) { cout << "! something bad occurred !" << endl ; return kFALSE ; }
 376   else            { cout << "# event done :) " << entry << "                # ok ! #       " << endl ; }
 377
 378  // Saves the AliFlowEvent
 379   if(fSaveFlowEvents)
 380   {
 381    cout << " saving flow event :| " << endl ;
 382    TString strID = "" ; strID += entry ;
 383    fFlowfile->cd() ; fFlowEvent->Write(strID.Data()) ;
 384    cout << "# event saved :) " << strID.Data() << "                # ok ! #       " << endl ; cout << endl ;
 385   }
 386
 387  // Weights
 388   if(fOnFlyWeight)
 389   {
 390    cout << " weightening :| " << endl ;
 391    bool done = fFlowWeighter->WeightEvent(fFlowEvent) ;
 392    if(done) { cout << "# weights done :) " << entry << "                # ok ! #       " << endl ; }
 393   }
 394
 395  // On-Fly Analysis
 396   if(fOnFlyAnalysis)
 397   {
 398    cout << " doing analysis :| " << endl ;
 399    bool done = fFlowAnal->Analyze(fFlowEvent) ;
 400    fFlowAnal->PrintEventQuantities() ;
 401    if(done) { cout << "# analysis done :) " << entry << "                # ok ! #       " << endl ; }
 402   }
 403
 404   return kTRUE;
 405 }
 406
 407 //-----------------------------------------------------------------------
 408
 409 void AliSelectorFoF::SlaveTerminate()
 410 {
 411   // The SlaveTerminate() function is called after all entries or objects
 412   // have been processed. When running with PROOF SlaveTerminate() is called
 413   // on each slave server.
 414
 415   AliDebug(AliLog::kDebug, "=======SLAVETERMINATE=======");
 416
 417   DeleteKinematicsFile();
 418 }
 419
 420 //-----------------------------------------------------------------------
 421
 422 void AliSelectorFoF::Terminate()
 423 {
 424   // The Terminate() function is the last function to be called during
 425   // a query. It always runs on the client, it can be used to present
 426   // the results graphically or save the results to file.
 427
 428   AliDebug(AliLog::kDebug, "=========TERMINATE==========");
 429
 430   cout << " Finished ... " << endl ;
 431
 432   if(fDoNothing) { cout << "             ... & nothing was done ! " << endl ; cout << endl ; return ; }
 433
 434   cout << endl ;
 435   cout << "  nTracks:  " << fFlowMaker->GetNgoodTracks() << endl ;
 436   cout << "  nV0s:  " << fFlowMaker->GetNgoodV0s()  << endl ;
 437   cout << "  nTracks (|eta|<0.5):  " << fFlowMaker->GetNgoodTracksEta() << endl ;
 438   cout << "  nTracks+:  " << fFlowMaker->GetNposiTracks() << endl ;
 439   cout << "  nTracks-:  " << fFlowMaker->GetNnegaTracks() << endl ;
 440   cout << "  nTracks unconstrained:  " << fFlowMaker->GetNunconstrained() << endl ;
 441   cout << "  Bayesian (e,mu,pi,k,p) :  " ;
 442   for(int ii=0;ii<5;ii++) { cout << fFlowMaker->GetBayesianNorm(ii) << "   " ; }
 443   cout << " . " << endl ; cout << endl ;
 444
 445   if(fOnFlyWeight) { fFlowWeighter->Finish() ; delete fFlowWeighter ; }
 446
 447   if(fOnFlyAnalysis)
 448   {
 449    fFlowAnal->Resolution() ; cout << endl;
 450    fFlowAnal->Finish() ;  cout << endl;
 451    //cout << endl ; cout << "Particles normalized abundance:" << endl;  // shows the bayesian vector
 452    //fFlowAnal->PrintRunBayesian() ; cout << endl;
 453    delete fFlowAnal ;
 454   }
 455
 456   if(fSaveFlowEvents) { fFlowfile->Close() ; cout << " file closed . " << endl ; }
 457   delete fFlowMaker ;
 458
 459   delete fFlowSelect ;
 460
 461   cout << endl ;
 462   return ;
 463 }
 464
 465 //-----------------------------------------------------------------------
 466
 467 TTree* AliSelectorFoF::GetKinematics()
 468 {
 469   // Returns kinematics tree corresponding to current ESD active in fTree
 470   // Loads the kinematics from the kinematics file, the file is identified by replacing "AliESDs" to
 471   // "Kinematics" in the file path of the ESD file. This is a hack, to be changed!
 472
 473   if (!fKineFile)
 474   {
 475     if(!fTree->GetCurrentFile()) { return 0 ; }
 476
 477     TString fileName(fTree->GetCurrentFile()->GetName());
 478     fileName.ReplaceAll("AliESDs", "Kinematics");
 479
 480     // temporary workaround for PROOF bug #18505
 481     fileName.ReplaceAll("#Kinematics.root#Kinematics.root", "#Kinematics.root");
 482
 483     AliDebug(AliLog::kInfo, Form("Opening %s", fileName.Data()));
 484
 485     fKineFile = TFile::Open(fileName);
 486     if(!fKineFile) { return 0 ; }
 487   }
 488
 489   return dynamic_cast<TTree*> (fKineFile->Get(Form("Event%d/TreeK", fTree->GetTree()->GetReadEntry())));
 490 }
 491
 492 //-----------------------------------------------------------------------
 493
 494 void AliSelectorFoF::DeleteKinematicsFile()
 495 {
 496   //
 497   // Closes the kinematics file and deletes the pointer.
 498   //
 499
 500   if (fKineFile)
 501   {
 502     fKineFile->Close();
 503     delete fKineFile;
 504     fKineFile = 0;
 505   }
 506 }
 507
 508 //-----------------------------------------------------------------------
 509
 510