[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowWeighter.cxx

//////////////////////////////////////////////////////////////////////
//
// $Id$
//
// Author: Emanuele Simili
//
//////////////////////////////////////////////////////////////////////
//_____________________________________________________________
//
// Description: 
//         the AliFlowWeighter class generates the phi-weights which 
// are later used to correct for azimuthal anisotropy in the reconstruction 
// efficiency of the ALICE TPC. It also fills an histogram of normalised 
// particle abundancies, which can be used as bayesian weights for particle Id. 
//
// - The method Init() generates the histograms and opens a new fPhiWgt file.
// - The method WeightEvent(AliFlowEvent*) fills phi and PId histograms. 
//   It must be inserted in a loop over the event sample.
// - The method Finish() calculates the weights, saves the histograms and 
//   closes the file. The AliFlowSelection object is saved as well.
//
//////////////////////////////////////////////////////////////////////

#ifndef ALIFLOWWEIGHTER_CXX
#define ALIFLOWWEIGHTER_CXX

// ROOT things
#include <TROOT.h>
#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include <TObject.h>
#include <TObjArray.h>
#include <TOrdCollection.h>
#include <TH1.h>
#include <TVector2.h>

// Flow things
#include "AliFlowEvent.h"
#include "AliFlowTrack.h"
#include "AliFlowV0.h"
#include "AliFlowConstants.h"
#include "AliFlowSelection.h"
#include "AliFlowWeighter.h"

// ANSI things
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
using namespace std; //required for resolving the 'cout' symbol

ClassImp(AliFlowWeighter) 
//-----------------------------------------------------------------------
AliFlowWeighter::AliFlowWeighter(const AliFlowSelection* flowSelect)
{
 // default constructor (selection given or default selection)

 if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect) ; }
 else 		{ fFlowSelect = new AliFlowSelection() ; }
 
 // output file (histograms)
 fWgtFileName = "flowPhiWgt.hist.root" ;
 fWgtFile     = 0 ;
}
//-----------------------------------------------------------------------
AliFlowWeighter::~AliFlowWeighter()
{
 // default distructor (no actions) 
}
//-----------------------------------------------------------------------
Bool_t AliFlowWeighter::Init() 
{
// sets some defaults for the analysis
 
 cout << "* FlowWeighter *  -  Init()" << endl ; cout << endl ; 
 
 // Open output files (->plots)
 fWgtFile = new TFile(fWgtFileName.Data(), "RECREATE");
 fWgtFile->cd() ;

 // counters and pointers to 0
 fEventNumber = 0 ;
 fTrackNumber = 0 ; 	
 fNumberOfTracks = 0 ;	
 fNumberOfV0s    = 0 ;	
 fFlowEvent  = 0 ;
 fFlowTrack  = 0 ;
 fFlowTracks = 0 ;
 //for(Int_t ii=0;ii<3;ii++) { fVertex[ii] = 0 ; }

 // Histogram settings
 fPhiBins = AliFlowConstants::kPhiBins ;     
 fPhiMin  = 0.;
 fPhiMax  = 2*TMath::Pi() ; 

 TString* histTitle ;
 for(int k = 0; k < AliFlowConstants::kSels; k++)
 {
  histTitle = new TString("Flow_BayPidMult_Sel");
  *histTitle += k+1;
  fHistFull[k].fHistBayPidMult = new TH1F(histTitle->Data(), histTitle->Data(),AliFlowConstants::kPid,-0.5,((Float_t)AliFlowConstants::kPid-0.5));
  fHistFull[k].fHistBayPidMult->Sumw2() ;
  fHistFull[k].fHistBayPidMult->SetXTitle("e+/-  ,  mu+/-  ,  pi+/-  ,  K+/-  ,  p+/-  ,  d+/- ");
  fHistFull[k].fHistBayPidMult->SetYTitle("Counts");
  delete histTitle;

  for(int j = 0; j < AliFlowConstants::kHars; j++) 
  {
   // Tpc - Phi lab
   histTitle = new TString("Flow_Phi_TPC_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhi = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhi->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhi->SetYTitle("Counts");
   delete histTitle;
   // Tpc (plus)
   histTitle = new TString("Flow_Phi_TPCplus_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetYTitle("Counts");
   delete histTitle;
   // Tpc (minus)
   histTitle = new TString("Flow_Phi_TPCminus_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetYTitle("Counts");
   delete histTitle;
   // Tpc (cross)
   histTitle = new TString("Flow_Phi_TPCcross_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiAll->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiAll->SetYTitle("Counts");
   delete histTitle;

   // Tpc - Phi lab flattened
   histTitle = new TString("Flow_Phi_Flat_TPC_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiFlat = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetYTitle("Counts");
   delete histTitle;
   // Tpc Plus
   histTitle = new TString("Flow_Phi_Flat_TPCplus_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetYTitle("Counts");
   delete histTitle;
   // Tpc Minus
   histTitle = new TString("Flow_Phi_Flat_TPCminus_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetYTitle("Counts");
   delete histTitle;
   // Tpc cross
   histTitle = new TString("Flow_Phi_Flat_TPCcross_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiFlatAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetYTitle("Counts");
   delete histTitle;
  }
 }
 
 return kTRUE ;
} 
//-----------------------------------------------------------------------

//-----------------------------------------------------------------------
Bool_t AliFlowWeighter::Finish() 
{
 // Calls the method to fill wgt histograms, then saves them 
 // on the fWgtFile and closes the file .
 
 cout << "* FlowWeighter *  -  Finish()" << endl ; cout << endl ;

 Weightening() ;
 fWgtFile->cd() ; 

 // Write PhiWgt histograms
 fPhiWgtHistList->Write();
 delete fPhiWgtHistList ;
 
 // Write Bayesian Weights for P.Id.
 for(int k=0;k<AliFlowConstants::kSels;k++) { fHistFull[k].fHistBayPidMult->Write() ; }

 // Write the AliFlowSelection object
 fFlowSelect->Write();
 delete fFlowSelect ;

 fWgtFile->Close();

 cout << "    Finish()  -  Wgt file closed : " << fWgtFileName.Data() << endl ; cout << endl ; 

 return kTRUE ;
}
//-----------------------------------------------------------------------
// ###
//-----------------------------------------------------------------------
Bool_t AliFlowWeighter::WeightEvent(AliFlowEvent* fFlowEvent)         
{
 // Reads the AliFlowEvent (* fFlowEvent) and loops over the 
 // AliFlowTraks to fill phi histograms .
 
 cout << " AliFlowWeighter::WeightEvent(" << fFlowEvent << " )   -   " << fEventNumber << endl ;
 if(!fFlowEvent) { return kFALSE ; }

 if(fFlowSelect->Select(fFlowEvent))	 // event selected - here below the ANALYSIS FLAGS are setted -
 {
  fFlowTracks = fFlowEvent->TrackCollection() ; 
  fNumberOfTracks = fFlowTracks->GetEntries() ;
  cout << "       event ID = " << fFlowEvent->EventID() << " :  found " << fNumberOfTracks << " AliFlowTracks . " << endl ;  
  fFlowEvent->SetNoWgt() ;
  fFlowEvent->SetSelections(fFlowSelect) ; // does the selection of tracks for r.p. calculation (sets flags in AliFlowTrack)
  
  TracksLoop(fFlowTracks) ;  
 }
 else 
 {
  cout << "       * " << fEventNumber << " (event ID = " << fFlowEvent->EventID() << ") discarded . " << endl ; 
  delete fFlowEvent  ; fFlowEvent = 0 ; 
  return kFALSE ;
 }
 fEventNumber++ ;
 
 return kTRUE ;
}
//-----------------------------------------------------------------------
// ###
//-----------------------------------------------------------------------
void AliFlowWeighter::TracksLoop(TObjArray* fFlowTracks) 
{
 // fills phi and PId histograms

 cout << " Tracks Loop . " << endl ; 
  
 Float_t phi, eta, zFirstPoint ; // , zLastPoint ;
 Int_t  fPidId ;
 Char_t pid[10] = "0" ;
 for(fTrackNumber=0;fTrackNumber<fNumberOfTracks;fTrackNumber++) 
 {
  fFlowTrack = (AliFlowTrack*)fFlowTracks->At(fTrackNumber) ;
  // cout << "Track n. " << fTrackNumber << endl ; fFlowTrack->Dump() ; 
 
  phi = fFlowTrack->Phi() ;		
  eta = fFlowTrack->Eta() ;		
  zFirstPoint = fFlowTrack->ZFirstPoint() ; 
  //zLastPoint = fFlowTrack->ZLastPoint() ;
  strcpy(pid,fFlowTrack->Pid()) ; 

  fPidId = -1 ;
  if(strstr(pid,"e"))  	    { fPidId = 0 ; }
  else if(strstr(pid,"mu")) { fPidId = 1 ; }
  else if(strstr(pid,"pi")) { fPidId = 2 ; }
  else if(strstr(pid,"k"))  { fPidId = 3 ; }
  else if(strstr(pid,"pr")) { fPidId = 4 ; }
  else if(strstr(pid,"d"))  { fPidId = 5 ; }

  // Looping over Selections and Harmonics
  for (int k = 0; k < AliFlowConstants::kSels; k++) 
  {
   fFlowSelect->SetSelection(k) ;
   for (int j = 0; j < AliFlowConstants::kHars; j++) 
   {
    fFlowSelect->SetHarmonic(j);
    if(fFlowSelect->Select(fFlowTrack))
    {
     Bool_t kTpcPlus  = kFALSE ;
     Bool_t kTpcMinus = kFALSE ;
     Bool_t kTpcAll   = kFALSE ;

     // Set Tpc (+ and -)
     if(fFlowTrack->FitPtsTPC())
     {
      if(zFirstPoint >= 0. && eta > 0.)	     { kTpcPlus  = kTRUE ; } 
      else if(zFirstPoint <= 0. && eta < 0.) { kTpcMinus = kTRUE ; }
      else				     { kTpcAll   = kTRUE ; }	  
     }
    
     // Phi distribution (particle for R.P.)
     Float_t wt = 1. ; // TMath::Abs(fFlowEvent->Weight(k, j, fFlowTrack)) ; 
     if(kTpcPlus)       { fHistFull[k].fHistFullHar[j].fHistPhiPlus->Fill(phi,wt) ;  } 
     else if(kTpcMinus) { fHistFull[k].fHistFullHar[j].fHistPhiMinus->Fill(phi,wt) ; } 
     else if(kTpcAll)   { fHistFull[k].fHistFullHar[j].fHistPhiAll->Fill(phi,wt) ;   } 
     fHistFull[k].fHistFullHar[j].fHistPhi->Fill(phi,wt) ;

     // PID Multiplicities (particle for R.P.) - just once for each selection
     if(j==0) 		{ fHistFull[k].fHistBayPidMult->Fill(fPidId) ; }
    }
   }
  }
 }
  
 return ;
}
//-----------------------------------------------------------------------
Bool_t AliFlowWeighter::Weightening()
{
 // Calculates weights, and fills PhiWgt histograms .
 // This is called at the end of the event loop .
 
 cout << " AliFlowWeighter::Weightening() " << endl ; cout << endl ;
 
 // PhiWgt histogram collection
 fPhiWgtHistList = new TOrdCollection(4*AliFlowConstants::kSels*AliFlowConstants::kHars) ;
 
 // Creates PhiWgt Histograms
 TString* histTitle ;
 for(Int_t k = 0; k < AliFlowConstants::kSels; k++)
 {
  for(Int_t j = 0; j < AliFlowConstants::kHars; j++) 
  {
   // Tpc plus
   histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Sumw2();
   fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetYTitle("PhiWgt");
   delete histTitle;
   // Tpc minus
   histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Sumw2();
   fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetYTitle("PhiWgt");
   delete histTitle;
   // Tpc cross
   histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiWgtAll = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Sumw2();
   fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetYTitle("PhiWgt");
   delete histTitle;
   // Tpc
   histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
   *histTitle += k+1;
   histTitle->Append("_Har");
   *histTitle += j+1;
   fHistFull[k].fHistFullHar[j].fHistPhiWgt = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
   fHistFull[k].fHistFullHar[j].fHistPhiWgt->Sumw2();
   fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetXTitle("Azimuthal Angles (rad)");
   fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetYTitle("PhiWgt");
   delete histTitle;

   // Calculate PhiWgt
   Double_t meanPlus  = fHistFull[k].fHistFullHar[j].fHistPhiPlus->Integral() / (Double_t)fPhiBins ;
   Double_t meanMinus = fHistFull[k].fHistFullHar[j].fHistPhiMinus->Integral() / (Double_t)fPhiBins ;
   Double_t meanCross = fHistFull[k].fHistFullHar[j].fHistPhiAll->Integral() / (Double_t)fPhiBins ;
   Double_t meanTPC = fHistFull[k].fHistFullHar[j].fHistPhi->Integral() / (Double_t)fPhiBins ;

   // Tpc
   for(Int_t i=0;i<fPhiBins;i++) 
   {
    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinContent(i+1,meanPlus);
    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinError(i+1, 0.);
    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinContent(i+1,meanMinus);
    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinError(i+1, 0.);
    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinContent(i+1,meanCross);
    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinError(i+1, 0.);
    fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinContent(i+1,meanTPC);
    fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinError(i+1, 0.);
   }
   
   if(meanTPC==0) { cout << " Sel." << k << " , Har." << j << " :  empty histogram ! " << endl ; }
   else 
   {
    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiPlus);
    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiMinus);
    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Divide(fHistFull[k].fHistFullHar[j].fHistPhiAll);
    fHistFull[k].fHistFullHar[j].fHistPhiWgt->Divide(fHistFull[k].fHistFullHar[j].fHistPhi);
   }

   fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus);
   fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus);
   fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtAll);
   fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgt);
  }
 }

 return kTRUE ;
}
//-----------------------------------------------------------------------
void AliFlowWeighter::PrintBayesian(Int_t selN) const
{
 // Prints the normalized particle abundance of all events (selection selN).
 // Call this at the end of the loop, just before Finish() .

 if(selN>AliFlowConstants::kSels) { selN = 0 ; } 
 Char_t* names[AliFlowConstants::kPid] = {"e","mu","pi","k","p","d"} ;
 Double_t bayes = 0. ;
 Double_t totCount = (fHistFull[selN].fHistBayPidMult)->GetSumOfWeights() ;
 if(totCount) 
 { 
  cout << " Sel." << selN << " particles normalized abundance (tot. " << totCount << " tracks) : " ;
  for(Int_t ii=0;ii<AliFlowConstants::kPid;ii++)
  {
   bayes = (fHistFull[selN].fHistBayPidMult->GetBinContent(ii+1) / totCount) ; 
   cout << bayes << "_" << names[ii] << " ; " ;
  }
 }
 else {  cout << " Sel." << selN << " :  empty P.Id. histogram ! " << endl ; }
 cout << endl ;
 
 return ;
}
//-----------------------------------------------------------------------


#endif
Commit	Line	Data
327288af	1	//////////////////////////////////////////////////////////////////////
9777bfcb	2	//
	3	// $Id$
	4	//
	5	// Author: Emanuele Simili
	6	//
	7	//////////////////////////////////////////////////////////////////////
	8	//_____________________________________________________________
	9	//
	10	// Description:
	11	// the AliFlowWeighter class generates the phi-weights which
	12	// are later used to correct for azimuthal anisotropy in the reconstruction
	13	// efficiency of the ALICE TPC. It also fills an histogram of normalised
	14	// particle abundancies, which can be used as bayesian weights for particle Id.
	15	//
	16	// - The method Init() generates the histograms and opens a new fPhiWgt file.
	17	// - The method WeightEvent(AliFlowEvent*) fills phi and PId histograms.
	18	// It must be inserted in a loop over the event sample.
	19	// - The method Finish() calculates the weights, saves the histograms and
	20	// closes the file. The AliFlowSelection object is saved as well.
	21	//
	22	//////////////////////////////////////////////////////////////////////
	23
327288af	24	#ifndef ALIFLOWWEIGHTER_CXX
327288af	25	#define ALIFLOWWEIGHTER_CXX
9777bfcb	26
	27	// ROOT things
	28	#include <TROOT.h>
	29	#include <TMath.h>
	30	#include <TFile.h>
	31	#include <TTree.h>
	32	#include <TString.h>
	33	#include <TObject.h>
	34	#include <TObjArray.h>
	35	#include <TOrdCollection.h>
	36	#include <TH1.h>
	37	#include <TVector2.h>
	38
	39	// Flow things
	40	#include "AliFlowEvent.h"
	41	#include "AliFlowTrack.h"
	42	#include "AliFlowV0.h"
	43	#include "AliFlowConstants.h"
	44	#include "AliFlowSelection.h"
	45	#include "AliFlowWeighter.h"
	46
	47	// ANSI things
51a9b7d8	48	#include <math.h>
9777bfcb	49	#include <stdio.h>
	50	#include <stdlib.h>
	51	#include <iostream>
	52	using namespace std; //required for resolving the 'cout' symbol
	53
da4d7f17	54	ClassImp(AliFlowWeighter)
9777bfcb	55	//-----------------------------------------------------------------------
	56	AliFlowWeighter::AliFlowWeighter(const AliFlowSelection* flowSelect)
	57	{
	58	// default constructor (selection given or default selection)
	59
	60	if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect) ; }
	61	else { fFlowSelect = new AliFlowSelection() ; }
	62
	63	// output file (histograms)
	64	fWgtFileName = "flowPhiWgt.hist.root" ;
	65	fWgtFile = 0 ;
	66	}
	67	//-----------------------------------------------------------------------
	68	AliFlowWeighter::~AliFlowWeighter()
	69	{
	70	// default distructor (no actions)
	71	}
	72	//-----------------------------------------------------------------------
	73	Bool_t AliFlowWeighter::Init()
	74	{
	75	// sets some defaults for the analysis
	76
	77	cout << "* FlowWeighter * - Init()" << endl ; cout << endl ;
	78
	79	// Open output files (->plots)
	80	fWgtFile = new TFile(fWgtFileName.Data(), "RECREATE");
	81	fWgtFile->cd() ;
	82
	83	// counters and pointers to 0
	84	fEventNumber = 0 ;
	85	fTrackNumber = 0 ;
	86	fNumberOfTracks = 0 ;
	87	fNumberOfV0s = 0 ;
	88	fFlowEvent = 0 ;
	89	fFlowTrack = 0 ;
	90	fFlowTracks = 0 ;
	91	//for(Int_t ii=0;ii<3;ii++) { fVertex[ii] = 0 ; }
	92
	93	// Histogram settings
327288af	94	fPhiBins = AliFlowConstants::kPhiBins ;
9777bfcb	95	fPhiMin = 0.;
	96	fPhiMax = 2*TMath::Pi() ;
	97
	98	TString* histTitle ;
327288af	99	for(int k = 0; k < AliFlowConstants::kSels; k++)
9777bfcb	100	{
	101	histTitle = new TString("Flow_BayPidMult_Sel");
	102	*histTitle += k+1;
327288af	103	fHistFull[k].fHistBayPidMult = new TH1F(histTitle->Data(), histTitle->Data(),AliFlowConstants::kPid,-0.5,((Float_t)AliFlowConstants::kPid-0.5));
9777bfcb	104	fHistFull[k].fHistBayPidMult->Sumw2() ;
	105	fHistFull[k].fHistBayPidMult->SetXTitle("e+/- , mu+/- , pi+/- , K+/- , p+/- , d+/- ");
	106	fHistFull[k].fHistBayPidMult->SetYTitle("Counts");
	107	delete histTitle;
	108
327288af	109	for(int j = 0; j < AliFlowConstants::kHars; j++)
9777bfcb	110	{
	111	// Tpc - Phi lab
	112	histTitle = new TString("Flow_Phi_TPC_Sel");
	113	*histTitle += k+1;
	114	histTitle->Append("_Har");
	115	*histTitle += j+1;
	116	fHistFull[k].fHistFullHar[j].fHistPhi = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	117	fHistFull[k].fHistFullHar[j].fHistPhi->SetXTitle("Azimuthal Angles (rad)");
	118	fHistFull[k].fHistFullHar[j].fHistPhi->SetYTitle("Counts");
	119	delete histTitle;
	120	// Tpc (plus)
	121	histTitle = new TString("Flow_Phi_TPCplus_Sel");
	122	*histTitle += k+1;
	123	histTitle->Append("_Har");
	124	*histTitle += j+1;
	125	fHistFull[k].fHistFullHar[j].fHistPhiPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	126	fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetXTitle("Azimuthal Angles (rad)");
	127	fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetYTitle("Counts");
	128	delete histTitle;
	129	// Tpc (minus)
	130	histTitle = new TString("Flow_Phi_TPCminus_Sel");
	131	*histTitle += k+1;
	132	histTitle->Append("_Har");
	133	*histTitle += j+1;
	134	fHistFull[k].fHistFullHar[j].fHistPhiMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	135	fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetXTitle("Azimuthal Angles (rad)");
	136	fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetYTitle("Counts");
	137	delete histTitle;
	138	// Tpc (cross)
	139	histTitle = new TString("Flow_Phi_TPCcross_Sel");
	140	*histTitle += k+1;
	141	histTitle->Append("_Har");
	142	*histTitle += j+1;
	143	fHistFull[k].fHistFullHar[j].fHistPhiAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	144	fHistFull[k].fHistFullHar[j].fHistPhiAll->SetXTitle("Azimuthal Angles (rad)");
	145	fHistFull[k].fHistFullHar[j].fHistPhiAll->SetYTitle("Counts");
	146	delete histTitle;
	147
	148	// Tpc - Phi lab flattened
	149	histTitle = new TString("Flow_Phi_Flat_TPC_Sel");
	150	*histTitle += k+1;
	151	histTitle->Append("_Har");
	152	*histTitle += j+1;
	153	fHistFull[k].fHistFullHar[j].fHistPhiFlat = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	154	fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetXTitle("Azimuthal Angles (rad)");
	155	fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetYTitle("Counts");
	156	delete histTitle;
	157	// Tpc Plus
	158	histTitle = new TString("Flow_Phi_Flat_TPCplus_Sel");
	159	*histTitle += k+1;
	160	histTitle->Append("_Har");
	161	*histTitle += j+1;
	162	fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	163	fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetXTitle("Azimuthal Angles (rad)");
	164	fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetYTitle("Counts");
	165	delete histTitle;
	166	// Tpc Minus
	167	histTitle = new TString("Flow_Phi_Flat_TPCminus_Sel");
	168	*histTitle += k+1;
	169	histTitle->Append("_Har");
	170	*histTitle += j+1;
	171	fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	172	fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetXTitle("Azimuthal Angles (rad)");
	173	fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetYTitle("Counts");
174	delete histTitle;
175	// Tpc cross
176	histTitle = new TString("Flow_Phi_Flat_TPCcross_Sel");
177	*histTitle += k+1;
178	histTitle->Append("_Har");
179	*histTitle += j+1;
180	fHistFull[k].fHistFullHar[j].fHistPhiFlatAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
181	fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetXTitle("Azimuthal Angles (rad)");
182	fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetYTitle("Counts");
183	delete histTitle;
184	}
185	}
186
187	return kTRUE ;
188	}
189	//-----------------------------------------------------------------------
190
191	//-----------------------------------------------------------------------
192	Bool_t AliFlowWeighter::Finish()
193	{
194	// Calls the method to fill wgt histograms, then saves them
195	// on the fWgtFile and closes the file .
196
197	cout << "* FlowWeighter * - Finish()" << endl ; cout << endl ;
198
199	Weightening() ;
200	fWgtFile->cd() ;
201
202	// Write PhiWgt histograms
203	fPhiWgtHistList->Write();
204	delete fPhiWgtHistList ;
205
206	// Write Bayesian Weights for P.Id.
327288af	207	for(int k=0;k<AliFlowConstants::kSels;k++) { fHistFull[k].fHistBayPidMult->Write() ; }
9777bfcb	208
	209	// Write the AliFlowSelection object
	210	fFlowSelect->Write();
	211	delete fFlowSelect ;
	212
	213	fWgtFile->Close();
	214
	215	cout << " Finish() - Wgt file closed : " << fWgtFileName.Data() << endl ; cout << endl ;
	216
	217	return kTRUE ;
	218	}
	219	//-----------------------------------------------------------------------
	220	// ###
	221	//-----------------------------------------------------------------------
	222	Bool_t AliFlowWeighter::WeightEvent(AliFlowEvent* fFlowEvent)
	223	{
	224	// Reads the AliFlowEvent (* fFlowEvent) and loops over the
	225	// AliFlowTraks to fill phi histograms .
	226
	227	cout << " AliFlowWeighter::WeightEvent(" << fFlowEvent << " ) - " << fEventNumber << endl ;
	228	if(!fFlowEvent) { return kFALSE ; }
	229
	230	if(fFlowSelect->Select(fFlowEvent)) // event selected - here below the ANALYSIS FLAGS are setted -
	231	{
	232	fFlowTracks = fFlowEvent->TrackCollection() ;
	233	fNumberOfTracks = fFlowTracks->GetEntries() ;
	234	cout << " event ID = " << fFlowEvent->EventID() << " : found " << fNumberOfTracks << " AliFlowTracks . " << endl ;
	235	fFlowEvent->SetNoWgt() ;
	236	fFlowEvent->SetSelections(fFlowSelect) ; // does the selection of tracks for r.p. calculation (sets flags in AliFlowTrack)
	237
	238	TracksLoop(fFlowTracks) ;
	239	}
	240	else
	241	{
	242	cout << " * " << fEventNumber << " (event ID = " << fFlowEvent->EventID() << ") discarded . " << endl ;
	243	delete fFlowEvent ; fFlowEvent = 0 ;
	244	return kFALSE ;
	245	}
	246	fEventNumber++ ;
	247
	248	return kTRUE ;
	249	}
	250	//-----------------------------------------------------------------------
	251	// ###
	252	//-----------------------------------------------------------------------
	253	void AliFlowWeighter::TracksLoop(TObjArray* fFlowTracks)
	254	{
	255	// fills phi and PId histograms
	256
	257	cout << " Tracks Loop . " << endl ;
	258
	259	Float_t phi, eta, zFirstPoint ; // , zLastPoint ;
	260	Int_t fPidId ;
	261	Char_t pid[10] = "0" ;
	262	for(fTrackNumber=0;fTrackNumber<fNumberOfTracks;fTrackNumber++)
	263	{
	264	fFlowTrack = (AliFlowTrack*)fFlowTracks->At(fTrackNumber) ;
	265	// cout << "Track n. " << fTrackNumber << endl ; fFlowTrack->Dump() ;
	266
	267	phi = fFlowTrack->Phi() ;
	268	eta = fFlowTrack->Eta() ;
	269	zFirstPoint = fFlowTrack->ZFirstPoint() ;
	270	//zLastPoint = fFlowTrack->ZLastPoint() ;
	271	strcpy(pid,fFlowTrack->Pid()) ;
272
273	fPidId = -1 ;
274	if(strstr(pid,"e")) { fPidId = 0 ; }
275	else if(strstr(pid,"mu")) { fPidId = 1 ; }
276	else if(strstr(pid,"pi")) { fPidId = 2 ; }
277	else if(strstr(pid,"k")) { fPidId = 3 ; }
278	else if(strstr(pid,"pr")) { fPidId = 4 ; }
279	else if(strstr(pid,"d")) { fPidId = 5 ; }
280
281	// Looping over Selections and Harmonics
327288af	282	for (int k = 0; k < AliFlowConstants::kSels; k++)
9777bfcb	283	{
9777bfcb	284	fFlowSelect->SetSelection(k) ;
327288af	285	for (int j = 0; j < AliFlowConstants::kHars; j++)
9777bfcb	286	{
	287	fFlowSelect->SetHarmonic(j);
	288	if(fFlowSelect->Select(fFlowTrack))
	289	{
	290	Bool_t kTpcPlus = kFALSE ;
	291	Bool_t kTpcMinus = kFALSE ;
	292	Bool_t kTpcAll = kFALSE ;
	293
	294	// Set Tpc (+ and -)
	295	if(fFlowTrack->FitPtsTPC())
	296	{
	297	if(zFirstPoint >= 0. && eta > 0.) { kTpcPlus = kTRUE ; }
	298	else if(zFirstPoint <= 0. && eta < 0.) { kTpcMinus = kTRUE ; }
	299	else { kTpcAll = kTRUE ; }
	300	}
	301
	302	// Phi distribution (particle for R.P.)
	303	Float_t wt = 1. ; // TMath::Abs(fFlowEvent->Weight(k, j, fFlowTrack)) ;
	304	if(kTpcPlus) { fHistFull[k].fHistFullHar[j].fHistPhiPlus->Fill(phi,wt) ; }
	305	else if(kTpcMinus) { fHistFull[k].fHistFullHar[j].fHistPhiMinus->Fill(phi,wt) ; }
	306	else if(kTpcAll) { fHistFull[k].fHistFullHar[j].fHistPhiAll->Fill(phi,wt) ; }
	307	fHistFull[k].fHistFullHar[j].fHistPhi->Fill(phi,wt) ;
	308
	309	// PID Multiplicities (particle for R.P.) - just once for each selection
	310	if(j==0) { fHistFull[k].fHistBayPidMult->Fill(fPidId) ; }
	311	}
	312	}
	313	}
	314	}
	315
	316	return ;
	317	}
	318	//-----------------------------------------------------------------------
	319	Bool_t AliFlowWeighter::Weightening()
	320	{
	321	// Calculates weights, and fills PhiWgt histograms .
	322	// This is called at the end of the event loop .
	323
	324	cout << " AliFlowWeighter::Weightening() " << endl ; cout << endl ;
	325
	326	// PhiWgt histogram collection
327288af	327	fPhiWgtHistList = new TOrdCollection(4AliFlowConstants::kSelsAliFlowConstants::kHars) ;
9777bfcb	328
	329	// Creates PhiWgt Histograms
	330	TString* histTitle ;
327288af	331	for(Int_t k = 0; k < AliFlowConstants::kSels; k++)
9777bfcb	332	{
327288af	333	for(Int_t j = 0; j < AliFlowConstants::kHars; j++)
9777bfcb	334	{
	335	// Tpc plus
	336	histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
	337	*histTitle += k+1;
	338	histTitle->Append("_Har");
	339	*histTitle += j+1;
	340	fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	341	fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Sumw2();
	342	fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetXTitle("Azimuthal Angles (rad)");
	343	fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetYTitle("PhiWgt");
	344	delete histTitle;
	345	// Tpc minus
	346	histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
	347	*histTitle += k+1;
	348	histTitle->Append("_Har");
	349	*histTitle += j+1;
	350	fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	351	fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Sumw2();
	352	fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetXTitle("Azimuthal Angles (rad)");
	353	fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetYTitle("PhiWgt");
	354	delete histTitle;
	355	// Tpc cross
	356	histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
	357	*histTitle += k+1;
	358	histTitle->Append("_Har");
	359	*histTitle += j+1;
	360	fHistFull[k].fHistFullHar[j].fHistPhiWgtAll = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	361	fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Sumw2();
	362	fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetXTitle("Azimuthal Angles (rad)");
	363	fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetYTitle("PhiWgt");
	364	delete histTitle;
	365	// Tpc
	366	histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
	367	*histTitle += k+1;
	368	histTitle->Append("_Har");
	369	*histTitle += j+1;
	370	fHistFull[k].fHistFullHar[j].fHistPhiWgt = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
	371	fHistFull[k].fHistFullHar[j].fHistPhiWgt->Sumw2();
	372	fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetXTitle("Azimuthal Angles (rad)");
	373	fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetYTitle("PhiWgt");
	374	delete histTitle;
	375
	376	// Calculate PhiWgt
	377	Double_t meanPlus = fHistFull[k].fHistFullHar[j].fHistPhiPlus->Integral() / (Double_t)fPhiBins ;
	378	Double_t meanMinus = fHistFull[k].fHistFullHar[j].fHistPhiMinus->Integral() / (Double_t)fPhiBins ;
	379	Double_t meanCross = fHistFull[k].fHistFullHar[j].fHistPhiAll->Integral() / (Double_t)fPhiBins ;
	380	Double_t meanTPC = fHistFull[k].fHistFullHar[j].fHistPhi->Integral() / (Double_t)fPhiBins ;
	381
	382	// Tpc
	383	for(Int_t i=0;i<fPhiBins;i++)
	384	{
	385	fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinContent(i+1,meanPlus);
	386	fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinError(i+1, 0.);
	387	fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinContent(i+1,meanMinus);
	388	fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinError(i+1, 0.);
	389	fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinContent(i+1,meanCross);
	390	fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinError(i+1, 0.);
	391	fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinContent(i+1,meanTPC);
	392	fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinError(i+1, 0.);
	393	}
	394
	395	if(meanTPC==0) { cout << " Sel." << k << " , Har." << j << " : empty histogram ! " << endl ; }
	396	else
	397	{
398	fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiPlus);
399	fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiMinus);
400	fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Divide(fHistFull[k].fHistFullHar[j].fHistPhiAll);
401	fHistFull[k].fHistFullHar[j].fHistPhiWgt->Divide(fHistFull[k].fHistFullHar[j].fHistPhi);
402	}
403
404	fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus);
405	fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus);
406	fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtAll);
407	fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgt);
408	}
409	}
410
411	return kTRUE ;
412	}
413	//-----------------------------------------------------------------------
327288af	414	void AliFlowWeighter::PrintBayesian(Int_t selN) const
9777bfcb	415	{
	416	// Prints the normalized particle abundance of all events (selection selN).
	417	// Call this at the end of the loop, just before Finish() .
	418
327288af	419	if(selN>AliFlowConstants::kSels) { selN = 0 ; }
327288af	420	Char_t* names[AliFlowConstants::kPid] = {"e","mu","pi","k","p","d"} ;
9777bfcb	421	Double_t bayes = 0. ;
	422	Double_t totCount = (fHistFull[selN].fHistBayPidMult)->GetSumOfWeights() ;
	423	if(totCount)
	424	{
	425	cout << " Sel." << selN << " particles normalized abundance (tot. " << totCount << " tracks) : " ;
327288af	426	for(Int_t ii=0;ii<AliFlowConstants::kPid;ii++)
9777bfcb	427	{
	428	bayes = (fHistFull[selN].fHistBayPidMult->GetBinContent(ii+1) / totCount) ;
	429	cout << bayes << "_" << names[ii] << " ; " ;
	430	}
	431	}
	432	else { cout << " Sel." << selN << " : empty P.Id. histogram ! " << endl ; }
	433	cout << endl ;
	434
	435	return ;
	436	}
	437	//-----------------------------------------------------------------------
	438
	439
	440	#endif