[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");
  fHistList->SetName("cobjSP");
  fHistList->SetOwner(kTRUE);
  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
  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");
  fHistList->SetName("cobjSP");
  fHistList->SetOwner(kTRUE);
  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
  delete output;
}

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

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

  Int_t iNbinsPt   = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
  Double_t dPtMin  = AliFlowCommonConstants::GetMaster()->GetPtMin();	     
  Double_t dPtMax  = AliFlowCommonConstants::GetMaster()->GetPtMax();
  Int_t iNbinsEta  = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
  Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();	     
  Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->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 Q vectors for the eta-subevents
    AliFlowVector* vQarray = new AliFlowVector[2];
    anEvent->GetQsub(vQarray);
    AliFlowVector vQa = vQarray[0];
    AliFlowVector vQb = vQarray[1];
    //get total Q vector
    AliFlowVector vQ = vQa + vQb;
    
    //fill the multiplicity histograms for the prefactor
    fHistProM -> Fill(1,vQ.GetMult()-1);                //<M-1>
    fHistProM -> Fill(2,vQa.GetMult()*vQb.GetMult());   //<Ma*Mb>
    //scalar product of the two subevents
    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;
	  //subtract particle from the flowvector if used to define it
	  if (pTrack->InRPSelection()) {
	    if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) { 
	      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;
    delete [] vQarray;
  }
}

  //--------------------------------------------------------------------  
void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
  
  //get pointers to all output histograms (called before Finish())
  if (outputListHistos) {
  //Get the common histograms from the output list
    AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*> 
      (outputListHistos->FindObject("AliFlowCommonHistSP"));
    AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*> 
      (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
    TProfile* pHistProQaQb     = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
    TProfile* pHistProM        = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_M_SP"));
    TProfile* pHistProUQetaRP  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
    TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
    TProfile* pHistProUQPtRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
    TProfile* pHistProUQPtPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
    if (pCommonHist && pCommonHistResults && pHistProQaQb && pHistProM &&
	pHistProUQetaRP	&& pHistProUQetaPOI && pHistProUQPtRP && pHistProUQPtPOI) {
      this -> SetCommonHists(pCommonHist);
      this -> SetCommonHistsRes(pCommonHistResults);
      this -> SetHistProQaQb(pHistProQaQb);
      this -> SetHistProM(pHistProM);
      this -> SetHistProUQetaRP(pHistProUQetaRP);
      this -> SetHistProUQetaPOI(pHistProUQetaPOI);
      this -> SetHistProUQPtRP(pHistProUQPtRP);
      this -> SetHistProUQPtPOI(pHistProUQPtPOI);
      }  
   }
}            

//--------------------------------------------------------------------            
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::GetMaster()->GetNbinsPt();
  Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();

  Double_t dMmin1    = fHistProM->GetBinContent(1);  //average over M-1
  Double_t dMmin1Err = fHistProM->GetBinError(1);    //error on average over M-1
  Double_t dMaMb     = fHistProM->GetBinContent(2);  //average over Ma*Mb
  Double_t dMaMbErr  = fHistProM->GetBinError(2);    //error on average over Ma*Mb

  Double_t dMcorrection = 0.;     //correction factor for Ma != Mb
  Double_t dMcorrectionErr = 0.;  
  Double_t dMcorrectionErrRel = 0.; 
  Double_t dMcorrectionErrRel2 = 0.; 

  if (dMaMb != 0. && dMmin1 != 0.) {
    dMcorrection    = dMmin1/(TMath::Sqrt(dMaMb)); 
    dMcorrectionErr = dMcorrection*(dMmin1Err/dMmin1 + dMaMbErr/(2*dMaMb));
    dMcorrectionErrRel = dMcorrectionErr/dMcorrection;
    dMcorrectionErrRel2 = dMcorrectionErrRel*dMcorrectionErrRel;
  }

  Double_t dQaQbAv  = TMath::Abs(fHistProQaQb->GetBinContent(1)); //average over events //TEST TAKE ABS
  Double_t dQaQbErr = fHistProQaQb->GetBinError(1);
  Double_t dQaQbErrRel = 0.;
  if (dQaQbAv != 0.) {
    dQaQbErrRel = dQaQbErr/dQaQbAv; }
  Double_t dQaQbErrRel2 = dQaQbErrRel*dQaQbErrRel;

  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);  //DOES NOT WORK IF dQaQbAv IS NEGATIVE
    if (dMaMb>0.) { dQaQbSqrt *= dMcorrection; }
    else { dQaQbSqrt = 0.; }
    Double_t dQaQbSqrtErrRel2 = dMcorrectionErrRel2 + (1/4)*dQaQbErrRel2;
    
    //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");
0fe80f88	79	//output->WriteObject(fHistList, "cobjSP","SingleKey");
0fe80f88	80	fHistList->SetName("cobjSP");
9455e15e	81	fHistList->SetOwner(kTRUE);
0fe80f88	82	fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
1dfa3c16	83	delete output;
	84	}
	85
b0fda271	86	//-----------------------------------------------------------------------
	87
	88	void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
	89	{
	90	//store the final results in output .root file
	91
	92	TFile *output = new TFile(outputFileName.Data(),"RECREATE");
0fe80f88	93	//output->WriteObject(fHistList, "cobjSP","SingleKey");
0fe80f88	94	fHistList->SetName("cobjSP");
9455e15e	95	fHistList->SetOwner(kTRUE);
0fe80f88	96	fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
b0fda271	97	delete output;
	98	}
	99
8d312f00	100	//-----------------------------------------------------------------------
	101	void AliFlowAnalysisWithScalarProduct::Init() {
	102
	103	//Define all histograms
d7eb18ec	104	cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
8d312f00	105
bee2efdc	106	Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
	107	Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
	108	Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
	109	Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
	110	Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
	111	Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
8d312f00	112
395fadba	113	fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
	114	fHistProUQetaRP->SetXTitle("{eta}");
	115	fHistProUQetaRP->SetYTitle("<uQ>");
	116	fHistList->Add(fHistProUQetaRP);
	117
	118	fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
	119	fHistProUQetaPOI->SetXTitle("{eta}");
	120	fHistProUQetaPOI->SetYTitle("<uQ>");
	121	fHistList->Add(fHistProUQetaPOI);
	122
	123	fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax);
	124	fHistProUQPtRP->SetXTitle("p_t (GeV)");
	125	fHistProUQPtRP->SetYTitle("<uQ>");
	126	fHistList->Add(fHistProUQPtRP);
	127
	128	fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
	129	fHistProUQPtPOI->SetXTitle("p_t (GeV)");
	130	fHistProUQPtPOI->SetYTitle("<uQ>");
	131	fHistList->Add(fHistProUQPtPOI);
	132
e4cecfa0	133	fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, -0.5, 0.5);
	134	fHistProQaQb->SetYTitle("<QaQb>");
	135	fHistList->Add(fHistProQaQb);
8d312f00	136
92041674	137	fHistProM = new TProfile("Flow_M_SP","Flow_M_SP",2,0.5, 2.5);
	138	fHistProM -> SetYTitle("<*>");
	139	fHistProM -> SetXTitle("<M-1>, <Ma*Mb>");
	140	fHistList->Add(fHistProM);
	141
04f6283b	142	fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
099e1753	143	fHistList->Add(fCommonHists);
395fadba	144	fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
	145	fHistList->Add(fCommonHistsRes);
	146
e2d51347	147	fEventNumber = 0; //set number of events to zero
	148	}
	149
8d312f00	150	//-----------------------------------------------------------------------
8d312f00	151
8232a5ec	152	void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
8d312f00	153
8d312f00	154	//Fill histogram
8232a5ec	155	if (anEvent) {
8d312f00	156
8d312f00	157	//fill control histograms
8232a5ec	158	fCommonHists->FillControlHistograms(anEvent);
cbbaf54a	159
cbbaf54a	160	//get Q vectors for the eta-subevents
b125a454	161	AliFlowVector* vQarray = new AliFlowVector[2];
	162	anEvent->GetQsub(vQarray);
	163	AliFlowVector vQa = vQarray[0];
	164	AliFlowVector vQb = vQarray[1];
cbbaf54a	165	//get total Q vector
	166	AliFlowVector vQ = vQa + vQb;
	167
	168	//fill the multiplicity histograms for the prefactor
	169	fHistProM -> Fill(1,vQ.GetMult()-1); //<M-1>
	170	fHistProM -> Fill(2,vQa.GetMult()vQb.GetMult()); //<MaMb>
	171	//scalar product of the two subevents
e4cecfa0	172	Double_t dQaQb = vQa*vQb;
e4cecfa0	173	fHistProQaQb -> Fill(0.,dQaQb);
e35ddff0	174
8d312f00	175	//loop over the tracks of the event
e35ddff0	176	AliFlowTrackSimple* pTrack = NULL;
	177	Int_t iNumberOfTracks = anEvent->NumberOfTracks();
	178	for (Int_t i=0;i<iNumberOfTracks;i++)
8d312f00	179	{
e35ddff0	180	pTrack = anEvent->GetTrack(i) ;
e35ddff0	181	if (pTrack){
e35ddff0	182	Double_t dPhi = pTrack->Phi();
e35ddff0	183	//calculate vU
	184	TVector2 vU;
	185	Double_t dUX = TMath::Cos(2*dPhi);
	186	Double_t dUY = TMath::Sin(2*dPhi);
	187	vU.Set(dUX,dUY);
	188	Double_t dModulus = vU.Mod();
	189	if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus); // make length 1
	190	else cerr<<"dModulus is zero!"<<endl;
	191
	192	TVector2 vQm = vQ;
cbbaf54a	193	//subtract particle from the flowvector if used to define it
6a1a854d	194	if (pTrack->InRPSelection()) {
b125a454	195	if (pTrack->InSubevent(0) \|\| pTrack->InSubevent(1)) {
	196	Double_t dQmX = vQm.X() - dUX;
	197	Double_t dQmY = vQm.Y() - dUY;
	198	vQm.Set(dQmX,dQmY);
	199	}
8d312f00	200	}
8d312f00	201
e35ddff0	202	//dUQ = scalar product of vU and vQm
	203	Double_t dUQ = vU * vQm;
	204	Double_t dPt = pTrack->Pt();
395fadba	205	Double_t dEta = pTrack->Eta();
395fadba	206	//fill the profile histograms
b7cb54d5	207	if (pTrack->InRPSelection()) {
395fadba	208	fHistProUQetaRP -> Fill(dEta,dUQ);
395fadba	209	fHistProUQPtRP -> Fill(dPt,dUQ);
395fadba	210	}
b7cb54d5	211	if (pTrack->InPOISelection()) {
395fadba	212	fHistProUQetaPOI -> Fill(dEta,dUQ);
395fadba	213	fHistProUQPtPOI -> Fill(dPt,dUQ);
8d312f00	214	}
	215	}//track selected
	216	}//loop over tracks
d7eb18ec	217
8d312f00	218	fEventNumber++;
deebd72f	219	// cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
b125a454	220	delete [] vQarray;
8d312f00	221	}
	222	}
	223
fd46c3dd	224	//--------------------------------------------------------------------
	225	void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
	226
	227	//get pointers to all output histograms (called before Finish())
	228	if (outputListHistos) {
	229	//Get the common histograms from the output list
	230	AliFlowCommonHist pCommonHist = dynamic_cast<AliFlowCommonHist>
	231	(outputListHistos->FindObject("AliFlowCommonHistSP"));
	232	AliFlowCommonHistResults pCommonHistResults = dynamic_cast<AliFlowCommonHistResults>
	233	(outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
	234	TProfile* pHistProQaQb = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
	235	TProfile* pHistProM = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_M_SP"));
	236	TProfile* pHistProUQetaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
	237	TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
	238	TProfile* pHistProUQPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
	239	TProfile* pHistProUQPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
	240	if (pCommonHist && pCommonHistResults && pHistProQaQb && pHistProM &&
	241	pHistProUQetaRP && pHistProUQetaPOI && pHistProUQPtRP && pHistProUQPtPOI) {
	242	this -> SetCommonHists(pCommonHist);
	243	this -> SetCommonHistsRes(pCommonHistResults);
	244	this -> SetHistProQaQb(pHistProQaQb);
	245	this -> SetHistProM(pHistProM);
	246	this -> SetHistProUQetaRP(pHistProUQetaRP);
	247	this -> SetHistProUQetaPOI(pHistProUQetaPOI);
	248	this -> SetHistProUQPtRP(pHistProUQPtRP);
	249	this -> SetHistProUQPtPOI(pHistProUQPtPOI);
	250	}
	251	}
	252	}
	253
	254	//--------------------------------------------------------------------
8d312f00	255	void AliFlowAnalysisWithScalarProduct::Finish() {
8d312f00	256
e4cecfa0	257	//calculate flow and fill the AliFlowCommonHistResults
395fadba	258	if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
	259
	260	cout<<"*************************************"<<endl;
	261	cout<<"*************************************"<<endl;
	262	cout<<" Integrated flow from "<<endl;
	263	cout<<" Scalar product "<<endl;
	264	cout<<endl;
e4cecfa0	265
bee2efdc	266	Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
bee2efdc	267	Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
80f72ed6	268
	269	Double_t dMmin1 = fHistProM->GetBinContent(1); //average over M-1
	270	Double_t dMmin1Err = fHistProM->GetBinError(1); //error on average over M-1
c09d90fd	271	Double_t dMaMb = fHistProM->GetBinContent(2); //average over Ma*Mb
c09d90fd	272	Double_t dMaMbErr = fHistProM->GetBinError(2); //error on average over Ma*Mb
80f72ed6	273
	274	Double_t dMcorrection = 0.; //correction factor for Ma != Mb
	275	Double_t dMcorrectionErr = 0.;
	276	Double_t dMcorrectionErrRel = 0.;
	277	Double_t dMcorrectionErrRel2 = 0.;
	278
	279	if (dMaMb != 0. && dMmin1 != 0.) {
	280	dMcorrection = dMmin1/(TMath::Sqrt(dMaMb));
	281	dMcorrectionErr = dMcorrection(dMmin1Err/dMmin1 + dMaMbErr/(2dMaMb));
	282	dMcorrectionErrRel = dMcorrectionErr/dMcorrection;
	283	dMcorrectionErrRel2 = dMcorrectionErrRel*dMcorrectionErrRel;
	284	}
	285
c09d90fd	286	Double_t dQaQbAv = TMath::Abs(fHistProQaQb->GetBinContent(1)); //average over events //TEST TAKE ABS
80f72ed6	287	Double_t dQaQbErr = fHistProQaQb->GetBinError(1);
	288	Double_t dQaQbErrRel = 0.;
	289	if (dQaQbAv != 0.) {
	290	dQaQbErrRel = dQaQbErr/dQaQbAv; }
	291	Double_t dQaQbErrRel2 = dQaQbErrRel*dQaQbErrRel;
	292
e4cecfa0	293	if (dQaQbAv <= 0.){
	294	//set v to -0
	295	fCommonHistsRes->FillIntegratedFlowRP(-0.,0.);
	296	fCommonHistsRes->FillIntegratedFlow(-0.,0.);
	297	cout<<"dV(RP) = -0. +- 0."<<endl;
	298	fCommonHistsRes->FillIntegratedFlowPOI(-0.,0.);
	299	cout<<"dV(POI) = -0. +- 0."<<endl;
	300	} else {
c09d90fd	301	Double_t dQaQbSqrt = TMath::Sqrt(dQaQbAv); //DOES NOT WORK IF dQaQbAv IS NEGATIVE
80f72ed6	302	if (dMaMb>0.) { dQaQbSqrt *= dMcorrection; }
92041674	303	else { dQaQbSqrt = 0.; }
80f72ed6	304	Double_t dQaQbSqrtErrRel2 = dMcorrectionErrRel2 + (1/4)*dQaQbErrRel2;
80f72ed6	305
e4cecfa0	306	//v as a function of eta for RP selection
	307	for(Int_t b=0;b<iNbinsEta;b++) {
	308	Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
	309	Double_t duQerr = fHistProUQetaRP->GetBinError(b); //copy error for now
	310	Double_t duQerrRel = 0.;
	311	if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
	312	Double_t duQerrRel2 = duQerrRel*duQerrRel;
	313
92041674	314	Double_t dv2pro = 0.;
92041674	315	if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0	316	Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674	317	Double_t dv2errRel = 0.;
92041674	318	if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0	319	Double_t dv2err = dv2pro*dv2errRel;
	320	//fill TH1D
	321	fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2err);
	322	} //loop over bins b
395fadba	323
e4cecfa0	324	//v as a function of eta for POI selection
	325	for(Int_t b=0;b<iNbinsEta;b++) {
	326	Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
	327	Double_t duQerr = fHistProUQetaPOI->GetBinError(b); //copy error for now
	328	Double_t duQerrRel = 0.;
	329	if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
	330	Double_t duQerrRel2 = duQerrRel*duQerrRel;
	331
92041674	332	Double_t dv2pro = 0.;
92041674	333	if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0	334	Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674	335	Double_t dv2errRel = 0.;
92041674	336	if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0	337	Double_t dv2err = dv2pro*dv2errRel;
	338	//fill TH1D
	339	fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2err);
	340	} //loop over bins b
395fadba	341
e4cecfa0	342	//v as a function of Pt for RP selection
b7cb54d5	343	TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
e4cecfa0	344	Double_t dVRP = 0.;
	345	Double_t dSum = 0.;
	346	Double_t dErrV =0.;
	347
	348	for(Int_t b=0;b<iNbinsPt;b++) {
	349	Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
	350	Double_t duQerr = fHistProUQPtRP->GetBinError(b); //copy error for now
	351	Double_t duQerrRel = 0.;
	352	if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
	353	Double_t duQerrRel2 = duQerrRel*duQerrRel;
	354
92041674	355	Double_t dv2pro = 0.;
92041674	356	if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0	357	Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674	358	Double_t dv2errRel = 0.;
92041674	359	if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0	360	Double_t dv2err = dv2pro*dv2errRel;
	361	//fill TH1D
	362	fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2err);
	363	//calculate integrated flow for RP selection
b7cb54d5	364	if (fHistPtRP){
b7cb54d5	365	Double_t dYieldPt = fHistPtRP->GetBinContent(b);
e4cecfa0	366	dVRP += dv2pro*dYieldPt;
	367	dSum +=dYieldPt;
	368	dErrV += dYieldPtdYieldPtdv2err*dv2err;
b7cb54d5	369	} else { cout<<"fHistPtRP is NULL"<<endl; }
e4cecfa0	370	} //loop over bins b
	371
	372	if (dSum != 0.) {
	373	dVRP /= dSum; //the pt distribution should be normalised
	374	dErrV /= (dSum*dSum);
	375	dErrV = TMath::Sqrt(dErrV);
	376	}
	377	fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
	378	fCommonHistsRes->FillIntegratedFlow(dVRP,dErrV);
	379
	380	cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
395fadba	381
e4cecfa0	382	//v as a function of Pt for POI selection
b7cb54d5	383	TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
e4cecfa0	384	Double_t dVPOI = 0.;
	385	dSum = 0.;
	386	dErrV =0.;
395fadba	387
e4cecfa0	388	for(Int_t b=0;b<iNbinsPt;b++) {
	389	Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
	390	Double_t duQerr = fHistProUQPtPOI->GetBinError(b); //copy error for now
	391	Double_t duQerrRel = 0.;
	392	if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
	393	Double_t duQerrRel2 = duQerrRel*duQerrRel;
	394
92041674	395	Double_t dv2pro = 0.;
92041674	396	if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0	397	Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674	398	Double_t dv2errRel = 0.;
92041674	399	if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0	400	Double_t dv2err = dv2pro*dv2errRel;
	401	//fill TH1D
	402	fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2err);
	403
	404	//calculate integrated flow for POI selection
b7cb54d5	405	if (fHistPtPOI){
b7cb54d5	406	Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
e4cecfa0	407	dVPOI += dv2pro*dYieldPt;
	408	dSum +=dYieldPt;
	409	dErrV += dYieldPtdYieldPtdv2err*dv2err;
b7cb54d5	410	} else { cout<<"fHistPtPOI is NULL"<<endl; }
e4cecfa0	411	} //loop over bins b
	412
	413	if (dSum != 0.) {
	414	dVPOI /= dSum; //the pt distribution should be normalised
	415	dErrV /= (dSum*dSum);
	416	dErrV = TMath::Sqrt(dErrV);
	417	}
	418	fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
	419
	420	cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
c09d90fd	421	}
395fadba	422	cout<<endl;
	423	cout<<"*************************************"<<endl;
	424	cout<<"*************************************"<<endl;
8d312f00	425
395fadba	426	//cout<<".....finished"<<endl;
8d312f00	427	}
8232a5ec	428
8232a5ec	429