[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowAnalysisWithScalarProduct.cxx

/*************************************************************************
* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  * 
**************************************************************************/

#define AliFlowAnalysisWithScalarProduct_cxx
 
#include "Riostream.h"  //needed as include
#include "TFile.h"      //needed as include
#include "TList.h"
#include "TMath.h"
#include "TProfile.h"
#include "TVector2.h"

class TH1F;

#include "AliFlowCommonConstants.h"    //needed as include
#include "AliFlowEventSimple.h"
#include "AliFlowTrackSimple.h"
#include "AliFlowCommonHist.h"
#include "AliFlowCommonHistResults.h"
#include "AliFlowAnalysisWithScalarProduct.h"

class AliFlowVector;

// AliFlowAnalysisWithScalarProduct:
// Description: 
// Maker to analyze Flow with the Scalar product method.
//
// author: N. van der Kolk (kolk@nikhef.nl)

ClassImp(AliFlowAnalysisWithScalarProduct)

  //-----------------------------------------------------------------------
 
 AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
   fEventNumber(0),
   fDebug(kFALSE),
   fHistList(NULL),
   fHistProUQetaRP(NULL),
   fHistProUQetaPOI(NULL),
   fHistProUQPtRP(NULL),
   fHistProUQPtPOI(NULL),
   fHistProQaQb(NULL),
   fHistProM(NULL),
   fCommonHists(NULL),
   fCommonHistsRes(NULL)
{
  // Constructor.
  fHistList = new TList();
}
 //-----------------------------------------------------------------------


 AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct() 
 {
   //destructor
   delete fHistList;
 }
 

//-----------------------------------------------------------------------

void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
{
 //store the final results in output .root file

  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
  output->WriteObject(fHistList, "cobjSP","SingleKey");
  delete output;
}

//-----------------------------------------------------------------------

void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
{
 //store the final results in output .root file

  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
  output->WriteObject(fHistList, "cobjSP","SingleKey");
  delete output;
}

//-----------------------------------------------------------------------
void AliFlowAnalysisWithScalarProduct::Init() {

  //Define all histograms
  cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;

  Int_t iNbinsPt   = AliFlowCommonConstants::GetNbinsPt();
  Double_t dPtMin  = AliFlowCommonConstants::GetPtMin();	     
  Double_t dPtMax  = AliFlowCommonConstants::GetPtMax();
  Int_t iNbinsEta  = AliFlowCommonConstants::GetNbinsEta();
  Double_t dEtaMin = AliFlowCommonConstants::GetEtaMin();	     
  Double_t dEtaMax = AliFlowCommonConstants::GetEtaMax();

  fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
  fHistProUQetaRP->SetXTitle("{eta}");
  fHistProUQetaRP->SetYTitle("<uQ>");
  fHistList->Add(fHistProUQetaRP);

  fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
  fHistProUQetaPOI->SetXTitle("{eta}");
  fHistProUQetaPOI->SetYTitle("<uQ>");
  fHistList->Add(fHistProUQetaPOI);

  fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax);
  fHistProUQPtRP->SetXTitle("p_t (GeV)");
  fHistProUQPtRP->SetYTitle("<uQ>");
  fHistList->Add(fHistProUQPtRP);

  fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
  fHistProUQPtPOI->SetXTitle("p_t (GeV)");
  fHistProUQPtPOI->SetYTitle("<uQ>");
  fHistList->Add(fHistProUQPtPOI);

  fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, -0.5, 0.5);
  fHistProQaQb->SetYTitle("<QaQb>");
  fHistList->Add(fHistProQaQb);

  fHistProM = new TProfile("Flow_M_SP","Flow_M_SP",2,0.5, 2.5);
  fHistProM -> SetYTitle("<*>");
  fHistProM -> SetXTitle("<M-1>, <Ma*Mb>");
  fHistList->Add(fHistProM);

  fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
  fHistList->Add(fCommonHists);
  fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
  fHistList->Add(fCommonHistsRes);  

  fEventNumber = 0;  //set number of events to zero    
}

//-----------------------------------------------------------------------
 
void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {

  //Fill histogram
  if (anEvent) {

    //fill control histograms     
    fCommonHists->FillControlHistograms(anEvent);
    
    //get the Q vector from the FlowEvent
    AliFlowVector vQ = anEvent->GetQ();
    fHistProM -> Fill(1,vQ.GetMult()-1);
    //get Q vectors for the subevents
    AliFlowVector vQa = anEvent->GetQsub(-0.9,-0.1);  
    AliFlowVector vQb = anEvent->GetQsub(0.1,0.9);
    fHistProM -> Fill(2,vQa.GetMult()*vQb.GetMult());
    Double_t dQaQb = vQa*vQb; 
    fHistProQaQb -> Fill(0.,dQaQb);    
                
    //loop over the tracks of the event
    AliFlowTrackSimple*   pTrack = NULL; 
    Int_t iNumberOfTracks = anEvent->NumberOfTracks(); 
    for (Int_t i=0;i<iNumberOfTracks;i++) 
      {
	pTrack = anEvent->GetTrack(i) ; 
	if (pTrack){
	  Double_t dPhi = pTrack->Phi();
	  //calculate vU
	  TVector2 vU;
	  Double_t dUX = TMath::Cos(2*dPhi);
	  Double_t dUY = TMath::Sin(2*dPhi);
	  vU.Set(dUX,dUY);
	  Double_t dModulus = vU.Mod();
	  if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus);  // make length 1
	  else cerr<<"dModulus is zero!"<<endl;

	  TVector2 vQm = vQ;
	  //subtrackt particle from the flowvector if used to define it
	  if (pTrack->UseForIntegratedFlow()) {
	    Double_t dQmX = vQm.X() - dUX;
	    Double_t dQmY = vQm.Y() - dUY;
	    vQm.Set(dQmX,dQmY);
	  }

	  //dUQ = scalar product of vU and vQm
	  Double_t dUQ = vU * vQm;
	  Double_t dPt = pTrack->Pt();
	  Double_t dEta = pTrack->Eta();
	  //fill the profile histograms
	  if (pTrack->InRPSelection()) {
	    fHistProUQetaRP -> Fill(dEta,dUQ);
	    fHistProUQPtRP -> Fill(dPt,dUQ);
	  }
	  if (pTrack->InPOISelection()) {
	    fHistProUQetaPOI -> Fill(dEta,dUQ);
	    fHistProUQPtPOI -> Fill(dPt,dUQ);
	  }  
	}//track selected
      }//loop over tracks
	 
    fEventNumber++;
    //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
  }
}

  //--------------------------------------------------------------------    
void AliFlowAnalysisWithScalarProduct::Finish() {
   
  //calculate flow and fill the AliFlowCommonHistResults
  if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;

  cout<<"*************************************"<<endl;
  cout<<"*************************************"<<endl;
  cout<<"      Integrated flow from           "<<endl;
  cout<<"         Scalar product              "<<endl;
  cout<<endl;
  
  Int_t iNbinsPt  = AliFlowCommonConstants::GetNbinsPt();
  Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
  Double_t dMmin1 = fHistProM->GetBinContent(1);  //average over M-1
  Double_t dMaMb  = fHistProM->GetBinContent(2);  //average over Ma*Mb
  Double_t dQaQbAv  = fHistProQaQb->GetBinContent(1); //average over events
  Double_t dQaQbErr = fHistProQaQb->GetBinError(1);
  if (dQaQbAv <= 0.){
    //set v to -0
    fCommonHistsRes->FillIntegratedFlowRP(-0.,0.);
    fCommonHistsRes->FillIntegratedFlow(-0.,0.);
    cout<<"dV(RP) = -0. +- 0."<<endl;
    fCommonHistsRes->FillIntegratedFlowPOI(-0.,0.);
    cout<<"dV(POI) = -0. +- 0."<<endl;
  } else {
    Double_t dQaQbSqrt = TMath::Sqrt(dQaQbAv);
    if (dMaMb>0.) { dQaQbSqrt *= dMmin1/(TMath::Sqrt(dMaMb)); }
    else { dQaQbSqrt = 0.; }
    Double_t dQaQbSqrtErr = (dQaQbSqrt/2)*(dQaQbErr/dQaQbAv);
    Double_t dQaQbSqrtErrRel = 0.;
    if (dQaQbSqrt!=0.) {dQaQbSqrtErr/dQaQbSqrt; }
    Double_t dQaQbSqrtErrRel2 = dQaQbSqrtErrRel*dQaQbSqrtErrRel;

    //v as a function of eta for RP selection
    for(Int_t b=0;b<iNbinsEta;b++) {
      Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
      Double_t duQerr = fHistProUQetaRP->GetBinError(b); //copy error for now
      Double_t duQerrRel = 0.;
      if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
      Double_t duQerrRel2 = duQerrRel*duQerrRel;

      Double_t dv2pro     = 0.;
      if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
      Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
      Double_t dv2errRel  = 0.;
      if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
      Double_t dv2err     = dv2pro*dv2errRel; 
      //fill TH1D
      fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2err); 
    } //loop over bins b
    
    //v as a function of eta for POI selection
    for(Int_t b=0;b<iNbinsEta;b++) {
      Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
      Double_t duQerr = fHistProUQetaPOI->GetBinError(b); //copy error for now
      Double_t duQerrRel = 0.;
      if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
      Double_t duQerrRel2 = duQerrRel*duQerrRel;

      Double_t dv2pro     = 0.;
      if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
      Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
      Double_t dv2errRel  = 0.;
      if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
      Double_t dv2err     = dv2pro*dv2errRel; 
      //fill TH1D
      fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2err); 
    } //loop over bins b
    
    //v as a function of Pt for RP selection
    TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
    Double_t dVRP = 0.;
    Double_t dSum = 0.;
    Double_t dErrV =0.;

    for(Int_t b=0;b<iNbinsPt;b++) {
      Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
      Double_t duQerr = fHistProUQPtRP->GetBinError(b); //copy error for now
      Double_t duQerrRel = 0.;
      if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
      Double_t duQerrRel2 = duQerrRel*duQerrRel;

      Double_t dv2pro     = 0.;
      if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
      Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
      Double_t dv2errRel  = 0.;
      if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
      Double_t dv2err     = dv2pro*dv2errRel; 
      //fill TH1D
      fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2err);
      //calculate integrated flow for RP selection
      if (fHistPtRP){
	Double_t dYieldPt = fHistPtRP->GetBinContent(b);
	dVRP += dv2pro*dYieldPt;
	dSum +=dYieldPt;
	dErrV += dYieldPt*dYieldPt*dv2err*dv2err;
      } else { cout<<"fHistPtRP is NULL"<<endl; }
    } //loop over bins b

    if (dSum != 0.) {
      dVRP /= dSum; //the pt distribution should be normalised
      dErrV /= (dSum*dSum);
      dErrV = TMath::Sqrt(dErrV);
    }
    fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
    fCommonHistsRes->FillIntegratedFlow(dVRP,dErrV);

    cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
       
    //v as a function of Pt for POI selection 
    TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
    Double_t dVPOI = 0.;
    dSum = 0.;
    dErrV =0.;
  
    for(Int_t b=0;b<iNbinsPt;b++) {
      Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
      Double_t duQerr = fHistProUQPtPOI->GetBinError(b); //copy error for now
      Double_t duQerrRel = 0.;
      if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
      Double_t duQerrRel2 = duQerrRel*duQerrRel;

      Double_t dv2pro     = 0.;
      if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
      Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
      Double_t dv2errRel  = 0.;
      if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
      Double_t dv2err     = dv2pro*dv2errRel; 
      //fill TH1D
      fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2err); 

      //calculate integrated flow for POI selection
      if (fHistPtPOI){
	Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
	dVPOI += dv2pro*dYieldPt;
	dSum +=dYieldPt;
	dErrV += dYieldPt*dYieldPt*dv2err*dv2err;
      } else { cout<<"fHistPtPOI is NULL"<<endl; }
    } //loop over bins b

    if (dSum != 0.) {
      dVPOI /= dSum; //the pt distribution should be normalised
      dErrV /= (dSum*dSum);
      dErrV = TMath::Sqrt(dErrV);
    }
    fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);

    cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
  }
  cout<<endl;
  cout<<"*************************************"<<endl;
  cout<<"*************************************"<<endl;   	  

  //cout<<".....finished"<<endl;
 }
Commit	Line	Data
8d312f00	1	/*************************************************************************
	2	* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	#define AliFlowAnalysisWithScalarProduct_cxx
	17
	18	#include "Riostream.h" //needed as include
1dfa3c16	19	#include "TFile.h" //needed as include
d7eb18ec	20	#include "TList.h"
8d312f00	21	#include "TMath.h"
	22	#include "TProfile.h"
	23	#include "TVector2.h"
	24
	25	class TH1F;
	26
	27	#include "AliFlowCommonConstants.h" //needed as include
	28	#include "AliFlowEventSimple.h"
	29	#include "AliFlowTrackSimple.h"
	30	#include "AliFlowCommonHist.h"
395fadba	31	#include "AliFlowCommonHistResults.h"
8d312f00	32	#include "AliFlowAnalysisWithScalarProduct.h"
	33
	34	class AliFlowVector;
	35
	36	// AliFlowAnalysisWithScalarProduct:
	37	// Description:
	38	// Maker to analyze Flow with the Scalar product method.
	39	//
	40	// author: N. van der Kolk (kolk@nikhef.nl)
	41
	42	ClassImp(AliFlowAnalysisWithScalarProduct)
	43
	44	//-----------------------------------------------------------------------
	45
	46	AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
e35ddff0	47	fEventNumber(0),
	48	fDebug(kFALSE),
	49	fHistList(NULL),
395fadba	50	fHistProUQetaRP(NULL),
	51	fHistProUQetaPOI(NULL),
	52	fHistProUQPtRP(NULL),
	53	fHistProUQPtPOI(NULL),
e4cecfa0	54	fHistProQaQb(NULL),
92041674	55	fHistProM(NULL),
395fadba	56	fCommonHists(NULL),
395fadba	57	fCommonHistsRes(NULL)
8d312f00	58	{
8d312f00	59	// Constructor.
e35ddff0	60	fHistList = new TList();
8d312f00	61	}
	62	//-----------------------------------------------------------------------
	63
	64
	65	AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
	66	{
	67	//destructor
e35ddff0	68	delete fHistList;
8d312f00	69	}
	70
	71
1dfa3c16	72	//-----------------------------------------------------------------------
	73
	74	void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
	75	{
	76	//store the final results in output .root file
	77
	78	TFile *output = new TFile(outputFileName->Data(),"RECREATE");
	79	output->WriteObject(fHistList, "cobjSP","SingleKey");
	80	delete output;
	81	}
	82
b0fda271	83	//-----------------------------------------------------------------------
	84
	85	void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
	86	{
	87	//store the final results in output .root file
	88
	89	TFile *output = new TFile(outputFileName.Data(),"RECREATE");
	90	output->WriteObject(fHistList, "cobjSP","SingleKey");
	91	delete output;
	92	}
	93
8d312f00	94	//-----------------------------------------------------------------------
	95	void AliFlowAnalysisWithScalarProduct::Init() {
	96
	97	//Define all histograms
d7eb18ec	98	cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
8d312f00	99
395fadba	100	Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
	101	Double_t dPtMin = AliFlowCommonConstants::GetPtMin();
	102	Double_t dPtMax = AliFlowCommonConstants::GetPtMax();
	103	Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
	104	Double_t dEtaMin = AliFlowCommonConstants::GetEtaMin();
	105	Double_t dEtaMax = AliFlowCommonConstants::GetEtaMax();
8d312f00	106
395fadba	107	fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
	108	fHistProUQetaRP->SetXTitle("{eta}");
	109	fHistProUQetaRP->SetYTitle("<uQ>");
	110	fHistList->Add(fHistProUQetaRP);
	111
	112	fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
	113	fHistProUQetaPOI->SetXTitle("{eta}");
	114	fHistProUQetaPOI->SetYTitle("<uQ>");
	115	fHistList->Add(fHistProUQetaPOI);
	116
	117	fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax);
	118	fHistProUQPtRP->SetXTitle("p_t (GeV)");
	119	fHistProUQPtRP->SetYTitle("<uQ>");
	120	fHistList->Add(fHistProUQPtRP);
	121
	122	fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
	123	fHistProUQPtPOI->SetXTitle("p_t (GeV)");
	124	fHistProUQPtPOI->SetYTitle("<uQ>");
	125	fHistList->Add(fHistProUQPtPOI);
	126
e4cecfa0	127	fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, -0.5, 0.5);
	128	fHistProQaQb->SetYTitle("<QaQb>");
	129	fHistList->Add(fHistProQaQb);
8d312f00	130
92041674	131	fHistProM = new TProfile("Flow_M_SP","Flow_M_SP",2,0.5, 2.5);
	132	fHistProM -> SetYTitle("<*>");
	133	fHistProM -> SetXTitle("<M-1>, <Ma*Mb>");
	134	fHistList->Add(fHistProM);
	135
04f6283b	136	fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
099e1753	137	fHistList->Add(fCommonHists);
395fadba	138	fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
	139	fHistList->Add(fCommonHistsRes);
	140
e2d51347	141	fEventNumber = 0; //set number of events to zero
	142	}
	143
8d312f00	144	//-----------------------------------------------------------------------
8d312f00	145
8232a5ec	146	void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
8d312f00	147
8d312f00	148	//Fill histogram
8232a5ec	149	if (anEvent) {
8d312f00	150
8d312f00	151	//fill control histograms
8232a5ec	152	fCommonHists->FillControlHistograms(anEvent);
92041674	153
8d312f00	154	//get the Q vector from the FlowEvent
e35ddff0	155	AliFlowVector vQ = anEvent->GetQ();
92041674	156	fHistProM -> Fill(1,vQ.GetMult()-1);
395fadba	157	//get Q vectors for the subevents
92041674	158	AliFlowVector vQa = anEvent->GetQsub(-0.9,-0.1);
	159	AliFlowVector vQb = anEvent->GetQsub(0.1,0.9);
	160	fHistProM -> Fill(2,vQa.GetMult()*vQb.GetMult());
e4cecfa0	161	Double_t dQaQb = vQa*vQb;
e4cecfa0	162	fHistProQaQb -> Fill(0.,dQaQb);
e35ddff0	163
8d312f00	164	//loop over the tracks of the event
e35ddff0	165	AliFlowTrackSimple* pTrack = NULL;
	166	Int_t iNumberOfTracks = anEvent->NumberOfTracks();
	167	for (Int_t i=0;i<iNumberOfTracks;i++)
8d312f00	168	{
e35ddff0	169	pTrack = anEvent->GetTrack(i) ;
e35ddff0	170	if (pTrack){
e35ddff0	171	Double_t dPhi = pTrack->Phi();
e35ddff0	172	//calculate vU
	173	TVector2 vU;
	174	Double_t dUX = TMath::Cos(2*dPhi);
	175	Double_t dUY = TMath::Sin(2*dPhi);
	176	vU.Set(dUX,dUY);
	177	Double_t dModulus = vU.Mod();
	178	if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus); // make length 1
	179	else cerr<<"dModulus is zero!"<<endl;
	180
	181	TVector2 vQm = vQ;
8d312f00	182	//subtrackt particle from the flowvector if used to define it
e35ddff0	183	if (pTrack->UseForIntegratedFlow()) {
	184	Double_t dQmX = vQm.X() - dUX;
	185	Double_t dQmY = vQm.Y() - dUY;
	186	vQm.Set(dQmX,dQmY);
8d312f00	187	}
8d312f00	188
e35ddff0	189	//dUQ = scalar product of vU and vQm
	190	Double_t dUQ = vU * vQm;
	191	Double_t dPt = pTrack->Pt();
395fadba	192	Double_t dEta = pTrack->Eta();
395fadba	193	//fill the profile histograms
b7cb54d5	194	if (pTrack->InRPSelection()) {
395fadba	195	fHistProUQetaRP -> Fill(dEta,dUQ);
395fadba	196	fHistProUQPtRP -> Fill(dPt,dUQ);
395fadba	197	}
b7cb54d5	198	if (pTrack->InPOISelection()) {
395fadba	199	fHistProUQetaPOI -> Fill(dEta,dUQ);
395fadba	200	fHistProUQPtPOI -> Fill(dPt,dUQ);
8d312f00	201	}
	202	}//track selected
	203	}//loop over tracks
d7eb18ec	204
8d312f00	205	fEventNumber++;
deebd72f	206	// cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
8d312f00	207	}
	208	}
	209
	210	//--------------------------------------------------------------------
	211	void AliFlowAnalysisWithScalarProduct::Finish() {
	212
e4cecfa0	213	//calculate flow and fill the AliFlowCommonHistResults
395fadba	214	if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
	215
	216	cout<<"*************************************"<<endl;
	217	cout<<"*************************************"<<endl;
	218	cout<<" Integrated flow from "<<endl;
	219	cout<<" Scalar product "<<endl;
	220	cout<<endl;
e4cecfa0	221
395fadba	222	Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
395fadba	223	Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
92041674	224	Double_t dMmin1 = fHistProM->GetBinContent(1); //average over M-1
92041674	225	Double_t dMaMb = fHistProM->GetBinContent(2); //average over Ma*Mb
e4cecfa0	226	Double_t dQaQbAv = fHistProQaQb->GetBinContent(1); //average over events
	227	Double_t dQaQbErr = fHistProQaQb->GetBinError(1);
	228	if (dQaQbAv <= 0.){
	229	//set v to -0
	230	fCommonHistsRes->FillIntegratedFlowRP(-0.,0.);
	231	fCommonHistsRes->FillIntegratedFlow(-0.,0.);
	232	cout<<"dV(RP) = -0. +- 0."<<endl;
	233	fCommonHistsRes->FillIntegratedFlowPOI(-0.,0.);
	234	cout<<"dV(POI) = -0. +- 0."<<endl;
	235	} else {
92041674	236	Double_t dQaQbSqrt = TMath::Sqrt(dQaQbAv);
	237	if (dMaMb>0.) { dQaQbSqrt *= dMmin1/(TMath::Sqrt(dMaMb)); }
	238	else { dQaQbSqrt = 0.; }
e4cecfa0	239	Double_t dQaQbSqrtErr = (dQaQbSqrt/2)*(dQaQbErr/dQaQbAv);
92041674	240	Double_t dQaQbSqrtErrRel = 0.;
92041674	241	if (dQaQbSqrt!=0.) {dQaQbSqrtErr/dQaQbSqrt; }
e4cecfa0	242	Double_t dQaQbSqrtErrRel2 = dQaQbSqrtErrRel*dQaQbSqrtErrRel;
	243
	244	//v as a function of eta for RP selection
	245	for(Int_t b=0;b<iNbinsEta;b++) {
	246	Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
	247	Double_t duQerr = fHistProUQetaRP->GetBinError(b); //copy error for now
	248	Double_t duQerrRel = 0.;
	249	if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
	250	Double_t duQerrRel2 = duQerrRel*duQerrRel;
	251
92041674	252	Double_t dv2pro = 0.;
92041674	253	if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0	254	Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674	255	Double_t dv2errRel = 0.;
92041674	256	if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0	257	Double_t dv2err = dv2pro*dv2errRel;
	258	//fill TH1D
	259	fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2err);
	260	} //loop over bins b
395fadba	261
e4cecfa0	262	//v as a function of eta for POI selection
	263	for(Int_t b=0;b<iNbinsEta;b++) {
	264	Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
	265	Double_t duQerr = fHistProUQetaPOI->GetBinError(b); //copy error for now
	266	Double_t duQerrRel = 0.;
	267	if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
	268	Double_t duQerrRel2 = duQerrRel*duQerrRel;
	269
92041674	270	Double_t dv2pro = 0.;
92041674	271	if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0	272	Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674	273	Double_t dv2errRel = 0.;
92041674	274	if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0	275	Double_t dv2err = dv2pro*dv2errRel;
	276	//fill TH1D
	277	fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2err);
	278	} //loop over bins b
395fadba	279
e4cecfa0	280	//v as a function of Pt for RP selection
b7cb54d5	281	TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
e4cecfa0	282	Double_t dVRP = 0.;
	283	Double_t dSum = 0.;
	284	Double_t dErrV =0.;
	285
	286	for(Int_t b=0;b<iNbinsPt;b++) {
	287	Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
	288	Double_t duQerr = fHistProUQPtRP->GetBinError(b); //copy error for now
	289	Double_t duQerrRel = 0.;
	290	if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
	291	Double_t duQerrRel2 = duQerrRel*duQerrRel;
	292
92041674	293	Double_t dv2pro = 0.;
92041674	294	if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0	295	Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674	296	Double_t dv2errRel = 0.;
92041674	297	if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0	298	Double_t dv2err = dv2pro*dv2errRel;
	299	//fill TH1D
	300	fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2err);
	301	//calculate integrated flow for RP selection
b7cb54d5	302	if (fHistPtRP){
b7cb54d5	303	Double_t dYieldPt = fHistPtRP->GetBinContent(b);
e4cecfa0	304	dVRP += dv2pro*dYieldPt;
	305	dSum +=dYieldPt;
	306	dErrV += dYieldPtdYieldPtdv2err*dv2err;
b7cb54d5	307	} else { cout<<"fHistPtRP is NULL"<<endl; }
e4cecfa0	308	} //loop over bins b
	309
	310	if (dSum != 0.) {
	311	dVRP /= dSum; //the pt distribution should be normalised
	312	dErrV /= (dSum*dSum);
	313	dErrV = TMath::Sqrt(dErrV);
	314	}
	315	fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
	316	fCommonHistsRes->FillIntegratedFlow(dVRP,dErrV);
	317
	318	cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
395fadba	319
e4cecfa0	320	//v as a function of Pt for POI selection
b7cb54d5	321	TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
e4cecfa0	322	Double_t dVPOI = 0.;
	323	dSum = 0.;
	324	dErrV =0.;
395fadba	325
e4cecfa0	326	for(Int_t b=0;b<iNbinsPt;b++) {
	327	Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
	328	Double_t duQerr = fHistProUQPtPOI->GetBinError(b); //copy error for now
	329	Double_t duQerrRel = 0.;
	330	if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
	331	Double_t duQerrRel2 = duQerrRel*duQerrRel;
	332
92041674	333	Double_t dv2pro = 0.;
92041674	334	if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0	335	Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674	336	Double_t dv2errRel = 0.;
92041674	337	if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0	338	Double_t dv2err = dv2pro*dv2errRel;
	339	//fill TH1D
	340	fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2err);
	341
	342	//calculate integrated flow for POI selection
b7cb54d5	343	if (fHistPtPOI){
b7cb54d5	344	Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
e4cecfa0	345	dVPOI += dv2pro*dYieldPt;
	346	dSum +=dYieldPt;
	347	dErrV += dYieldPtdYieldPtdv2err*dv2err;
b7cb54d5	348	} else { cout<<"fHistPtPOI is NULL"<<endl; }
e4cecfa0	349	} //loop over bins b
	350
	351	if (dSum != 0.) {
	352	dVPOI /= dSum; //the pt distribution should be normalised
	353	dErrV /= (dSum*dSum);
	354	dErrV = TMath::Sqrt(dErrV);
	355	}
	356	fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
	357
	358	cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
395fadba	359	}
395fadba	360	cout<<endl;
	361	cout<<"*************************************"<<endl;
	362	cout<<"*************************************"<<endl;
8d312f00	363
395fadba	364	//cout<<".....finished"<<endl;
8d312f00	365	}
8232a5ec	366
8232a5ec	367