<commit_before>/*************************************************************************
* 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"
#include "TFile.h"      
#include "TList.h"
#include "TMath.h"
#include "TProfile.h"
#include "TVector2.h"
#include "TH1D.h"
#include "TH2D.h"

#include "AliFlowCommonConstants.h"
#include "AliFlowEventSimple.h"
#include "AliFlowVector.h"
#include "AliFlowTrackSimple.h"
#include "AliFlowCommonHist.h"
#include "AliFlowCommonHistResults.h"
#include "AliFlowAnalysisWithScalarProduct.h"

//////////////////////////////////////////////////////////////////////////////
// AliFlowAnalysisWithScalarProduct:
// Description: 
// Maker to analyze Flow with the Scalar product method.
//
// authors: <NAME> (<EMAIL>), <NAME> (<EMAIL>)
//////////////////////////////////////////////////////////////////////////////

ClassImp(AliFlowAnalysisWithScalarProduct)

  //-----------------------------------------------------------------------
  AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
   fEventNumber(0),
   fDebug(kFALSE),
   fApplyCorrectionForNUA(kFALSE),
   fHarmonic(2),
   fTotalQvector(NULL),
   fRelDiffMsub(1.),
   fWeightsList(NULL),
   fUsePhiWeights(kFALSE),
   fPhiWeightsSub0(NULL),
   fPhiWeightsSub1(NULL),
   fHistProFlags(NULL),
   fHistProUQetaRP(NULL),
   fHistProUQetaPOI(NULL),
   fHistProUQetaAllEventsPOI(NULL),
   fHistProUQPtRP(NULL),
   fHistProUQPtPOI(NULL),
   fHistProUQPtAllEventsPOI(NULL),
   fHistProQNorm(NULL),
   fHistProQaQb(NULL),
   fHistProQaQbNorm(NULL),
   fHistProQaQbReImNorm(NULL),
   fHistProNonIsotropicTermsQ(NULL),
   fHistSumOfLinearWeights(NULL),
   fHistSumOfQuadraticWeights(NULL),
   fHistProUQQaQbPtRP(NULL),
   fHistProUQQaQbEtaRP(NULL),
   fHistProUQQaQbPtPOI(NULL),
   fHistProUQQaQbEtaPOI(NULL),
   fCommonHistsSP(NULL),
   fCommonHistsResSP(NULL),
   fCommonHistsmuQ(NULL),
   fHistQNorm(NULL),
   fHistQaQb(NULL),
   fHistQaQbNorm(NULL),
   fHistQNormvsQaQbNorm(NULL),
   fHistQaQbCos(NULL),
   fHistResolution(NULL),
   fHistQaNorm(NULL),
   fHistQaNormvsMa(NULL),
   fHistQbNorm(NULL),
   fHistQbNormvsMb(NULL),
   fHistMavsMb(NULL),
   fHistList(NULL)
{
  // Constructor.
  fWeightsList = new TList();
  fHistList = new TList();
  fHistList->SetOwner(kTRUE);
  
  // Total Q-vector is: "QaQb" (means Qa+Qb), "Qa"  or "Qb"
  fTotalQvector = new TString("QaQb");
  
  // Initialize arrays:
  for(Int_t i=0;i<3;i++)
  {
   fHistSumOfWeightsPtRP[i] = NULL;
   fHistSumOfWeightsEtaRP[i] = NULL;
   fHistSumOfWeightsPtPOI[i] = NULL;
   fHistSumOfWeightsEtaPOI[i] = NULL;
  }
  for(Int_t rp=0;rp<2;rp++)
  {
   for(Int_t pe=0;pe<2;pe++)
   {
    for(Int_t sc=0;sc<2;sc++)
    {
     fHistProNonIsotropicTermsU[rp][pe][sc] = NULL;
    }
   } 
  }
}
 //-----------------------------------------------------------------------
 AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct() 
 {
   //destructor
   delete fWeightsList;
   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::WriteHistograms(TDirectoryFile *outputFileName)
{
 //store the final results in output .root file
 fHistList->SetName("cobjSP");
 fHistList->SetOwner(kTRUE);
 outputFileName->Add(fHistList);
 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
}

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

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

  //save old value and prevent histograms from being added to directory
  //to avoid name clashes in case multiple analaysis objects are used
  //in an analysis

  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);
 
  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();

  fHistProFlags = new TProfile("FlowPro_Flags_SP","Flow_Flags_SP",1,0,1,"s");
  fHistProFlags->GetXaxis()->SetBinLabel(1,"fApplyCorrectionForNUA");
  fHistList->Add(fHistProFlags);
  
  fHistProUQetaRP = new TProfile("FlowPro_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
  fHistProUQetaRP->SetXTitle("{eta}");
  fHistProUQetaRP->SetYTitle("<uQ>");
  fHistList->Add(fHistProUQetaRP);

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

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

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

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

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

  fHistProQNorm = new TProfile("FlowPro_QNorm_SP","FlowPro_QNorm_SP", 1, 0.5, 1.5,"s");
  fHistProQNorm ->SetYTitle("<|Qa+Qb|>");
  fHistList->Add(fHistProQNorm); 

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

  fHistProQaQbNorm = new TProfile("FlowPro_QaQbNorm_SP","FlowPro_QaQbNorm_SP", 1, 0.5, 1.5,"s");
  fHistProQaQbNorm->SetYTitle("<QaQb/MaMb>");
  fHistList->Add(fHistProQaQbNorm);
  
  fHistProQaQbReImNorm = new TProfile("FlowPro_QaQbReImNorm_SP","FlowPro_QaQbReImNorm_SP", 4, 0.5, 4.5,"s");
  fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(1,"#LT#LTsin(#phi_{a})#GT#GT");
  fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(2,"#LT#LTcos(#phi_{a})#GT#GT");
  fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(3,"#LT#LTsin(#phi_{b})#GT#GT");
  fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(4,"#LT#LTcos(#phi_{b})#GT#GT");
  fHistList->Add(fHistProQaQbReImNorm); 
  
  fHistProNonIsotropicTermsQ = new TProfile("FlowPro_NonIsotropicTermsQ_SP","FlowPro_NonIsotropicTermsQ_SP", 2, 0.5, 2.5,"s");
  fHistProNonIsotropicTermsQ->GetXaxis()->SetBinLabel(1,"#LT#LTsin(#phi_{a+b})#GT#GT");
  fHistProNonIsotropicTermsQ->GetXaxis()->SetBinLabel(2,"#LT#LTcos(#phi_{a+b})#GT#GT");
  fHistList->Add(fHistProNonIsotropicTermsQ); 
  
  TString rpPoi[2] = {"RP","POI"};
  TString ptEta[2] = {"Pt","Eta"};
  TString sinCos[2] = {"sin","cos"};
  Int_t nBinsPtEta[2] = {iNbinsPt,iNbinsEta};
  Double_t minPtEta[2] = {dPtMin,dEtaMin};
  Double_t maxPtEta[2] = {dPtMax,dEtaMax};
  for(Int_t rp=0;rp<2;rp++)
  {
   for(Int_t pe=0;pe<2;pe++)
   {
    for(Int_t sc=0;sc<2;sc++)
    {  
     fHistProNonIsotropicTermsU[rp][pe][sc] = new TProfile(Form("FlowPro_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data()),Form("FlowPro_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data()),nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]);
     fHistList->Add(fHistProNonIsotropicTermsU[rp][pe][sc]);
    } 
   }
  } 
   
  fHistSumOfLinearWeights = new TH1D("Flow_SumOfLinearWeights_SP","Flow_SumOfLinearWeights_SP",1,-0.5, 0.5);
  fHistSumOfLinearWeights -> SetYTitle("sum (*)");
  fHistSumOfLinearWeights -> SetXTitle("sum (Ma*Mb)");
  fHistList->Add(fHistSumOfLinearWeights);
  
  fHistSumOfQuadraticWeights = new TH1D("Flow_SumOfQuadraticWeights_SP","Flow_SumOfQuadraticWeights_SP",1,-0.5, 0.5);
  fHistSumOfQuadraticWeights -> SetYTitle("sum (*)");
  fHistSumOfQuadraticWeights -> SetXTitle("sum (Ma*Mb)^2");
  fHistList->Add(fHistSumOfQuadraticWeights);
  
  fHistProUQQaQbPtRP = new TProfile("FlowPro_UQQaQbPtRP_SP","FlowPro_UQQaQbPtRP_SP",iNbinsPt,dPtMin,dPtMax);
  fHistProUQQaQbPtRP -> SetYTitle("<*>");
  fHistProUQQaQbPtRP -> SetXTitle("<Qu QaQb>");
  fHistList->Add(fHistProUQQaQbPtRP);
  
  fHistProUQQaQbEtaRP = new TProfile("FlowPro_UQQaQbEtaRP_SP","FlowPro_UQQaQbEtaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
  fHistProUQQaQbEtaRP -> SetYTitle("<*>");
  fHistProUQQaQbEtaRP -> SetXTitle("<Qu QaQb>");
  fHistList->Add(fHistProUQQaQbEtaRP);
  
  fHistProUQQaQbPtPOI = new TProfile("FlowPro_UQQaQbPtPOI_SP","FlowPro_UQQaQbPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
  fHistProUQQaQbPtPOI -> SetYTitle("<*>");
  fHistProUQQaQbPtPOI -> SetXTitle("<Qu QaQb>");
  fHistList->Add(fHistProUQQaQbPtPOI);
  
  fHistProUQQaQbEtaPOI = new TProfile("FlowPro_UQQaQbEtaPOI_SP","FlowPro_UQQaQbEtaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
  fHistProUQQaQbEtaPOI -> SetYTitle("<*>");
  fHistProUQQaQbEtaPOI -> SetXTitle("<Qu QaQb>");
  fHistList->Add(fHistProUQQaQbEtaPOI);
   
  TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"}; 
  for(Int_t i=0;i<3;i++)
  {
   fHistSumOfWeightsPtRP[i] = new TH1D(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),
                              Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
   fHistSumOfWeightsPtRP[i] -> SetYTitle("sum (*)");
   fHistSumOfWeightsPtRP[i] -> SetXTitle("p_{T}");
   fHistList->Add(fHistSumOfWeightsPtRP[i]);
 
   fHistSumOfWeightsEtaRP[i] = new TH1D(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),
                               Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
   fHistSumOfWeightsEtaRP[i] -> SetYTitle("sum (*)");
   fHistSumOfWeightsEtaRP[i] -> SetXTitle("#eta");
   fHistList->Add(fHistSumOfWeightsEtaRP[i]);
  
   fHistSumOfWeightsPtPOI[i] = new TH1D(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),
                               Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
   fHistSumOfWeightsPtPOI[i] -> SetYTitle("sum (*)");
   fHistSumOfWeightsPtPOI[i] -> SetXTitle("p_{T}");
   fHistList->Add(fHistSumOfWeightsPtPOI[i]);
 
   fHistSumOfWeightsEtaPOI[i] = new TH1D(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),
                                Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
   fHistSumOfWeightsEtaPOI[i] -> SetYTitle("sum (*)");
   fHistSumOfWeightsEtaPOI[i] -> SetXTitle("#eta");
   fHistList->Add(fHistSumOfWeightsEtaPOI[i]);
  }
      
  fCommonHistsSP = new AliFlowCommonHist("AliFlowCommonHistSP");
  fHistList->Add(fCommonHistsSP);
  fCommonHistsResSP = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
  fHistList->Add(fCommonHistsResSP);  
  fCommonHistsmuQ = new AliFlowCommonHist("AliFlowCommonHistmuQ");
  fHistList->Add(fCommonHistsmuQ);

  (fCommonHistsSP->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic 
  (fCommonHistsmuQ->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic 

  fHistQNorm = new TH1D("Flow_QNorm_SP","Flow_QNorm_SP",110,0.,1.1);
  fHistQNorm -> SetYTitle("dN/d(|(Qa+Qb)/(Ma+Mb)|)");
  fHistQNorm -> SetXTitle("|(Qa+Qb)/(Ma+Mb)|");
  fHistList->Add(fHistQNorm);

  fHistQaQb = new TH1D("Flow_QaQb_SP","Flow_QaQb_SP",200,-100.,100.);
  fHistQaQb -> SetYTitle("dN/dQaQb");
  fHistQaQb -> SetXTitle("QaQb");
  fHistList->Add(fHistQaQb);

  fHistQaQbNorm = new TH1D("Flow_QaQbNorm_SP","Flow_QaQbNorm_SP",44,-1.1,1.1);
  fHistQaQbNorm -> SetYTitle("dN/d(QaQb/MaMb)");
  fHistQaQbNorm -> SetXTitle("QaQb/MaMb");
  fHistList->Add(fHistQaQbNorm);

  fHistQNormvsQaQbNorm = new TH2D("Flow_QNormvsQaQbNorm_SP","Flow_QNormvsQaQbNorm_SP",88,-1.1,1.1,22,0.,1.1);
  fHistQNormvsQaQbNorm -> SetYTitle("|Q/Mq|");
  fHistQNormvsQaQbNorm -> SetXTitle("QaQb/MaMb");
  fHistList->Add(fHistQNormvsQaQbNorm);

  fHistQaQbCos = new TH1D("Flow_QaQbCos_SP","Flow_QaQbCos_SP",63,0.,TMath::Pi());
  fHistQaQbCos -> SetYTitle("dN/d(#phi)");
  fHistQaQbCos -> SetXTitle("#phi");
  fHistList->Add(fHistQaQbCos);

  fHistResolution = new TH1D("Flow_resolution_SP","Flow_resolution_SP",100,-1.0,1.0);
  fHistResolution -> SetYTitle("dN/d(cos(2(#phi_a - #phi_b))");
  fHistResolution -> SetXTitle("cos(2*(#phi_a - #phi_b))");
  fHistList->Add(fHistResolution);

  fHistQaNorm = new TH1D("Flow_QaNorm_SP","Flow_QaNorm_SP",22,0.,1.1);
  fHistQaNorm -> SetYTitle("dN/d(|Qa/Ma|)");
  fHistQaNorm -> SetXTitle("|Qa/Ma|");
  fHistList->Add(fHistQaNorm);

  fHistQaNormvsMa = new TH2D("Flow_QaNormvsMa_SP","Flow_QaNormvsMa_SP",100,0.,100.,22,0.,1.1);
  fHistQaNormvsMa -> SetYTitle("|Qa/Ma|");
  fHistQaNormvsMa -> SetXTitle("Ma");
  fHistList->Add(fHistQaNormvsMa);

  fHistQbNorm = new TH1D("Flow_QbNorm_SP","Flow_QbNorm_SP",22,0.,1.1);
  fHistQbNorm -> SetYTitle("dN/d(|Qb/Mb|)");
  fHistQbNorm -> SetXTitle("|Qb/Mb|");
  fHistList->Add(fHistQbNorm);

  fHistQbNormvsMb = new TH2D("Flow_QbNormvsMb_SP","Flow_QbNormvsMb_SP",100,0.,100.,22,0.,1.1);
  fHistQbNormvsMb -> SetYTitle("|Qb/Mb|");
  fHistQbNormvsMb -> SetXTitle("|Mb|");
  fHistList->Add(fHistQbNormvsMb);

  fHistMavsMb = new TH2D("Flow_MavsMb_SP","Flow_MavsMb_SP",100,0.,100.,100,0.,100.);
  fHistMavsMb -> SetYTitle("Ma");
  fHistMavsMb -> SetXTitle("Mb");
  fHistList->Add(fHistMavsMb);


  //weights
  if(fUsePhiWeights) {
    if(!fWeightsList) {
      cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
      exit(0);  
    }
    if(fWeightsList->FindObject("phi_weights_sub0"))  {
      fPhiWeightsSub0 = dynamic_cast<TH1F*>
	(fWeightsList->FindObject("phi_weights_sub0"));
      fHistList->Add(fPhiWeightsSub0);
    } else {
      cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
      exit(0);
    }
    if(fWeightsList->FindObject("phi_weights_sub1"))  {
      fPhiWeightsSub1 = dynamic_cast<TH1F*>
	(fWeightsList->FindObject("phi_weights_sub1"));
      fHistList->Add(fPhiWeightsSub1);
    } else {
      cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
      exit(0);
    }

  } // end of if(fUsePhiWeights)

  fEventNumber = 0;  //set number of events to zero 
  
  //store all boolean flags needed in Finish():
  this->StoreFlags();   

  TH1::AddDirectory(oldHistAddStatus);
}

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


  if (anEvent) {

  //Calculate muQ (for comparing pp and PbPb)
  FillmuQ(anEvent);

  //Calculate flow based on  <QaQb/MaMb> = <v^2>
  FillSP(anEvent);

  }
}

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

  //Calculate flow based on  <QaQb/MaMb> = <v^2>

  //Fill histograms
  if (anEvent) {

    //get Q vectors for the eta-subevents
    AliFlowVector* vQarray = new AliFlowVector[2];
    if (fUsePhiWeights) {
      anEvent->Get2Qsub(vQarray,fHarmonic,fWeightsList,kTRUE);
    } else {
      anEvent->Get2Qsub(vQarray,fHarmonic);
    }
    //subevent a
    AliFlowVector vQa = vQarray[0];
    //subevent b
    AliFlowVector vQb = vQarray[1];

    //For calculating v2 only events should be taken where both subevents are not empty
    //check that the subevents are not empty:
    Double_t dMa = vQa.GetMult();
    Double_t dMb = vQb.GetMult();
    if (dMa > 0. && dMb > 0.) {
      
      //request that the subevent multiplicities are not too different
      //fRelDiffMsub can be set from the configuration macro
      Double_t dRelDiff = TMath::Abs((dMa - dMb)/(dMa + dMb));
      if (dRelDiff < fRelDiffMsub) {

	//fill control histograms 	   
	if (fUsePhiWeights) {
	  fCommonHistsSP->FillControlHistograms(anEvent,fWeightsList,kTRUE);
	} else {
	  fCommonHistsSP->FillControlHistograms(anEvent);
	}

	//fill some SP control histograms
	fHistProQaQb -> Fill(1.,vQa*vQb,1.); //Fill with weight 1 -> Weight with MaMb????
	fHistQaQbCos ->Fill(TMath::ACos((vQa/vQa.Mod())*(vQb/vQb.Mod())));  //Acos(Qa*Qb) = angle
	fHistResolution -> Fill(TMath::Cos( vQa.Phi()- vQb.Phi() ));  //vQa.Phi() returns 2*phi
	fHistQaQb -> Fill(vQa*vQb);
	fHistMavsMb -> Fill(dMb,dMa);

	//get total Q vector = the sum of subevent a and subevent b
	AliFlowVector vQ;
	if(!strcmp(fTotalQvector->Data(),"QaQb"))
	{
	 vQ = vQa + vQb;
	} else if(!strcmp(fTotalQvector->Data(),"Qa"))
	  {
	   vQ = vQa; 
	  } else if(!strcmp(fTotalQvector->Data(),"Qb"))
	    {
	     vQ = vQb; 
	    }

	//needed to correct for non-uniform acceptance:
	fHistProNonIsotropicTermsQ->Fill(1.,vQ.Y()/(dMa+dMb),dMa+dMb);
	fHistProNonIsotropicTermsQ->Fill(2.,vQ.X()/(dMa+dMb),dMa+dMb);

	//weight the Q vectors for the subevents by the multiplicity
	//Note: Weight Q only in the particle loop when it is clear 
	//if it should be (m-1) or M
	Double_t dQXa = vQa.X()/dMa; 
	Double_t dQYa = vQa.Y()/dMa;
	vQa.Set(dQXa,dQYa);
	
	Double_t dQXb = vQb.X()/dMb; 
	Double_t dQYb = vQb.Y()/dMb;
	vQb.Set(dQXb,dQYb);
        
	//scalar product of the two subevents
	Double_t dQaQb = (vQa*vQb);
	fHistProQaQbNorm -> Fill(1.,dQaQb,dMa*dMb);  //Fill (QaQb/MaMb) with weight (MaMb). 
	//needed for the error calculation:
	fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
	fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
	//needed for correcting non-uniform acceptance: 
	fHistProQaQbReImNorm->Fill(1.,dQYa,dMa); // to get <<sin(phi_a)>>
	fHistProQaQbReImNorm->Fill(2.,dQXa,dMa); // to get <<cos(phi_a)>>
	fHistProQaQbReImNorm->Fill(3.,dQYb,dMb); // to get <<sin(phi_b)>>
	fHistProQaQbReImNorm->Fill(4.,dQXb,dMb); // to get <<cos(phi_b)>>
	
	//fill some SP control histograms
	fHistQaQbNorm ->Fill(vQa*vQb);
	fHistQaNorm ->Fill(vQa.Mod());
	fHistQaNormvsMa->Fill(dMa,vQa.Mod());
	fHistQbNorm ->Fill(vQb.Mod());
	fHistQbNormvsMb->Fill(dMb,vQb.Mod());
	
	//loop over the tracks of the event
	AliFlowTrackSimple*   pTrack = NULL; 
	Int_t iNumberOfTracks = anEvent->NumberOfTracks(); 
	Double_t dMq =  vQ.GetMult();
	
	for (Int_t i=0;i<iNumberOfTracks;i++) 
	  {
	    pTrack = anEvent->GetTrack(i) ; 
	    if (pTrack){
	      Double_t dPhi = pTrack->Phi();
	      Double_t dWeightUQ = 1.; // weight for u*Q 	    
	      //calculate vU
	      TVector2 vU;
	      //do not need to use weight for v as the length will be made 1
	      Double_t dUX = TMath::Cos(fHarmonic*dPhi);
	      Double_t dUY = TMath::Sin(fHarmonic*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;
	    
	      //redefine the Q vector and devide by its multiplicity
	      TVector2 vQm;
	      Double_t dQmX = 0.;
	      Double_t dQmY = 0.;
	      //subtract particle from the flowvector if used to define it
	      if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) { 
		//set default phi weight to 1
		Double_t dW = 1.; 
		//if phi weights are used
		if(fUsePhiWeights && fPhiWeightsSub0 && fPhiWeightsSub1) 
		{
		  if(strcmp(fTotalQvector->Data(),"QaQb"))
		  {
		   printf("\n WARNING (SP): If you use phi-weights total Q-vector has to be Qa+Qb in the current implementation!!!! \n");
		   exit(0);
		  }
		  //value of the center of the phi bin
		  Double_t dPhiCenter = 0.;  
		  if (pTrack->InSubevent(0) ) {
		    Int_t iNBinsPhiSub0 = fPhiWeightsSub0->GetNbinsX();
		    Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub0/TMath::TwoPi()));
		    dW = fPhiWeightsSub0->GetBinContent(phiBin); 
		    dPhiCenter = fPhiWeightsSub0->GetBinCenter(phiBin);
		    dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter) )/(dMq-1);
		    dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter) )/(dMq-1);
		    
		    vQm.Set(dQmX,dQmY);
		  }

		  else if ( pTrack->InSubevent(1)) { 
		    Int_t iNBinsPhiSub1 = fPhiWeightsSub1->GetNbinsX();
		    Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub1/TMath::TwoPi()));
		    dW = fPhiWeightsSub1->GetBinContent(phiBin);
		    dPhiCenter = fPhiWeightsSub1->GetBinCenter(phiBin);
		    dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter) )/(dMq-1);
		    dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter) )/(dMq-1);
		    
		    vQm.Set(dQmX,dQmY);
		  }
		  //bin = 1 + value*nbins/range
		  //TMath::Floor rounds to the lower integer
		}     
		// if no phi weights are used
		else 
		{
		 if(!strcmp(fTotalQvector->Data(),"QaQb"))
		 {
		  dQmX = (vQ.X() - (pTrack->Weight())*dUX)/(dMq-1);
		  dQmY = (vQ.Y() - (pTrack->Weight())*dUY)/(dMq-1);
		  dWeightUQ = dMq-1;
		  vQm.Set(dQmX,dQmY);
		 } else if((!strcmp(fTotalQvector->Data(),"Qa") && pTrack->InSubevent(0)) ||
		           (!strcmp(fTotalQvector->Data(),"Qb") && pTrack->InSubevent(1)))
		   {
		    dQmX = (vQ.X() - (pTrack->Weight())*dUX)/(dMq-1);
		    dQmY = (vQ.Y() - (pTrack->Weight())*dUY)/(dMq-1);
		    dWeightUQ = dMq-1;
		    vQm.Set(dQmX,dQmY);
		   } else if((!strcmp(fTotalQvector->Data(),"Qa") && pTrack->InSubevent(1)) ||
		             (!strcmp(fTotalQvector->Data(),"Qb") && pTrack->InSubevent(0)))
		     {
   		      dQmX = vQ.X()/dMq;
		      dQmY = vQ.Y()/dMq;
		      dWeightUQ = dMq;
		      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,dWeightUQ); //Fill (Qu/(Mq-1)) with weight (Mq-1) 
		  //needed for the error calculation:
		  fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dWeightUQ*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb	    
		  fHistProUQPtRP -> Fill(dPt,dUQ,dWeightUQ);                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
		  fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dWeightUQ*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb	
		  <commit_msg>fix for using particle weights<commit_after>/*************************************************************************
* 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"
#include "TFile.h"      
#include "TList.h"
#include "TMath.h"
#include "TProfile.h"
#include "TVector2.h"
#include "TH1D.h"
#include "TH2D.h"

#include "AliFlowCommonConstants.h"
#include "AliFlowEventSimple.h"
#include "AliFlowVector.h"
#include "AliFlowTrackSimple.h"
#include "AliFlowCommonHist.h"
#include "AliFlowCommonHistResults.h"
#include "AliFlowAnalysisWithScalarProduct.h"

//////////////////////////////////////////////////////////////////////////////
// AliFlowAnalysisWithScalarProduct:
// Description: 
// Maker to analyze Flow with the Scalar product method.
//
// authors: <NAME> (<EMAIL>), <NAME> (<EMAIL>)
//////////////////////////////////////////////////////////////////////////////

ClassImp(AliFlowAnalysisWithScalarProduct)

  //-----------------------------------------------------------------------
  AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
   fEventNumber(0),
   fDebug(kFALSE),
   fApplyCorrectionForNUA(kFALSE),
   fHarmonic(2),
   fTotalQvector(NULL),
   fRelDiffMsub(1.),
   fWeightsList(NULL),
   fUsePhiWeights(kFALSE),
   fPhiWeightsSub0(NULL),
   fPhiWeightsSub1(NULL),
   fHistProFlags(NULL),
   fHistProUQetaRP(NULL),
   fHistProUQetaPOI(NULL),
   fHistProUQetaAllEventsPOI(NULL),
   fHistProUQPtRP(NULL),
   fHistProUQPtPOI(NULL),
   fHistProUQPtAllEventsPOI(NULL),
   fHistProQNorm(NULL),
   fHistProQaQb(NULL),
   fHistProQaQbNorm(NULL),
   fHistProQaQbReImNorm(NULL),
   fHistProNonIsotropicTermsQ(NULL),
   fHistSumOfLinearWeights(NULL),
   fHistSumOfQuadraticWeights(NULL),
   fHistProUQQaQbPtRP(NULL),
   fHistProUQQaQbEtaRP(NULL),
   fHistProUQQaQbPtPOI(NULL),
   fHistProUQQaQbEtaPOI(NULL),
   fCommonHistsSP(NULL),
   fCommonHistsResSP(NULL),
   fCommonHistsmuQ(NULL),
   fHistQNorm(NULL),
   fHistQaQb(NULL),
   fHistQaQbNorm(NULL),
   fHistQNormvsQaQbNorm(NULL),
   fHistQaQbCos(NULL),
   fHistResolution(NULL),
   fHistQaNorm(NULL),
   fHistQaNormvsMa(NULL),
   fHistQbNorm(NULL),
   fHistQbNormvsMb(NULL),
   fHistMavsMb(NULL),
   fHistList(NULL)
{
  // Constructor.
  fWeightsList = new TList();
  fHistList = new TList();
  fHistList->SetOwner(kTRUE);
  
  // Total Q-vector is: "QaQb" (means Qa+Qb), "Qa"  or "Qb"
  fTotalQvector = new TString("QaQb");
  
  // Initialize arrays:
  for(Int_t i=0;i<3;i++)
  {
   fHistSumOfWeightsPtRP[i] = NULL;
   fHistSumOfWeightsEtaRP[i] = NULL;
   fHistSumOfWeightsPtPOI[i] = NULL;
   fHistSumOfWeightsEtaPOI[i] = NULL;
  }
  for(Int_t rp=0;rp<2;rp++)
  {
   for(Int_t pe=0;pe<2;pe++)
   {
    for(Int_t sc=0;sc<2;sc++)
    {
     fHistProNonIsotropicTermsU[rp][pe][sc] = NULL;
    }
   } 
  }
}
 //-----------------------------------------------------------------------
 AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct() 
 {
   //destructor
   delete fWeightsList;
   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::WriteHistograms(TDirectoryFile *outputFileName)
{
 //store the final results in output .root file
 fHistList->SetName("cobjSP");
 fHistList->SetOwner(kTRUE);
 outputFileName->Add(fHistList);
 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
}

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

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

  //save old value and prevent histograms from being added to directory
  //to avoid name clashes in case multiple analaysis objects are used
  //in an analysis

  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);
 
  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();

  fHistProFlags = new TProfile("FlowPro_Flags_SP","Flow_Flags_SP",1,0,1,"s");
  fHistProFlags->GetXaxis()->SetBinLabel(1,"fApplyCorrectionForNUA");
  fHistList->Add(fHistProFlags);
  
  fHistProUQetaRP = new TProfile("FlowPro_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
  fHistProUQetaRP->SetXTitle("{eta}");
  fHistProUQetaRP->SetYTitle("<uQ>");
  fHistList->Add(fHistProUQetaRP);

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

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

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

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

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

  fHistProQNorm = new TProfile("FlowPro_QNorm_SP","FlowPro_QNorm_SP", 1, 0.5, 1.5,"s");
  fHistProQNorm ->SetYTitle("<|Qa+Qb|>");
  fHistList->Add(fHistProQNorm); 

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

  fHistProQaQbNorm = new TProfile("FlowPro_QaQbNorm_SP","FlowPro_QaQbNorm_SP", 1, 0.5, 1.5,"s");
  fHistProQaQbNorm->SetYTitle("<QaQb/MaMb>");
  fHistList->Add(fHistProQaQbNorm);
  
  fHistProQaQbReImNorm = new TProfile("FlowPro_QaQbReImNorm_SP","FlowPro_QaQbReImNorm_SP", 4, 0.5, 4.5,"s");
  fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(1,"#LT#LTsin(#phi_{a})#GT#GT");
  fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(2,"#LT#LTcos(#phi_{a})#GT#GT");
  fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(3,"#LT#LTsin(#phi_{b})#GT#GT");
  fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(4,"#LT#LTcos(#phi_{b})#GT#GT");
  fHistList->Add(fHistProQaQbReImNorm); 
  
  fHistProNonIsotropicTermsQ = new TProfile("FlowPro_NonIsotropicTermsQ_SP","FlowPro_NonIsotropicTermsQ_SP", 2, 0.5, 2.5,"s");
  fHistProNonIsotropicTermsQ->GetXaxis()->SetBinLabel(1,"#LT#LTsin(#phi_{a+b})#GT#GT");
  fHistProNonIsotropicTermsQ->GetXaxis()->SetBinLabel(2,"#LT#LTcos(#phi_{a+b})#GT#GT");
  fHistList->Add(fHistProNonIsotropicTermsQ); 
  
  TString rpPoi[2] = {"RP","POI"};
  TString ptEta[2] = {"Pt","Eta"};
  TString sinCos[2] = {"sin","cos"};
  Int_t nBinsPtEta[2] = {iNbinsPt,iNbinsEta};
  Double_t minPtEta[2] = {dPtMin,dEtaMin};
  Double_t maxPtEta[2] = {dPtMax,dEtaMax};
  for(Int_t rp=0;rp<2;rp++)
  {
   for(Int_t pe=0;pe<2;pe++)
   {
    for(Int_t sc=0;sc<2;sc++)
    {  
     fHistProNonIsotropicTermsU[rp][pe][sc] = new TProfile(Form("FlowPro_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data()),Form("FlowPro_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data()),nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]);
     fHistList->Add(fHistProNonIsotropicTermsU[rp][pe][sc]);
    } 
   }
  } 
   
  fHistSumOfLinearWeights = new TH1D("Flow_SumOfLinearWeights_SP","Flow_SumOfLinearWeights_SP",1,-0.5, 0.5);
  fHistSumOfLinearWeights -> SetYTitle("sum (*)");
  fHistSumOfLinearWeights -> SetXTitle("sum (Ma*Mb)");
  fHistList->Add(fHistSumOfLinearWeights);
  
  fHistSumOfQuadraticWeights = new TH1D("Flow_SumOfQuadraticWeights_SP","Flow_SumOfQuadraticWeights_SP",1,-0.5, 0.5);
  fHistSumOfQuadraticWeights -> SetYTitle("sum (*)");
  fHistSumOfQuadraticWeights -> SetXTitle("sum (Ma*Mb)^2");
  fHistList->Add(fHistSumOfQuadraticWeights);
  
  fHistProUQQaQbPtRP = new TProfile("FlowPro_UQQaQbPtRP_SP","FlowPro_UQQaQbPtRP_SP",iNbinsPt,dPtMin,dPtMax);
  fHistProUQQaQbPtRP -> SetYTitle("<*>");
  fHistProUQQaQbPtRP -> SetXTitle("<Qu QaQb>");
  fHistList->Add(fHistProUQQaQbPtRP);
  
  fHistProUQQaQbEtaRP = new TProfile("FlowPro_UQQaQbEtaRP_SP","FlowPro_UQQaQbEtaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
  fHistProUQQaQbEtaRP -> SetYTitle("<*>");
  fHistProUQQaQbEtaRP -> SetXTitle("<Qu QaQb>");
  fHistList->Add(fHistProUQQaQbEtaRP);
  
  fHistProUQQaQbPtPOI = new TProfile("FlowPro_UQQaQbPtPOI_SP","FlowPro_UQQaQbPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
  fHistProUQQaQbPtPOI -> SetYTitle("<*>");
  fHistProUQQaQbPtPOI -> SetXTitle("<Qu QaQb>");
  fHistList->Add(fHistProUQQaQbPtPOI);
  
  fHistProUQQaQbEtaPOI = new TProfile("FlowPro_UQQaQbEtaPOI_SP","FlowPro_UQQaQbEtaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
  fHistProUQQaQbEtaPOI -> SetYTitle("<*>");
  fHistProUQQaQbEtaPOI -> SetXTitle("<Qu QaQb>");
  fHistList->Add(fHistProUQQaQbEtaPOI);
   
  TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"}; 
  for(Int_t i=0;i<3;i++)
  {
   fHistSumOfWeightsPtRP[i] = new TH1D(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),
                              Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
   fHistSumOfWeightsPtRP[i] -> SetYTitle("sum (*)");
   fHistSumOfWeightsPtRP[i] -> SetXTitle("p_{T}");
   fHistList->Add(fHistSumOfWeightsPtRP[i]);
 
   fHistSumOfWeightsEtaRP[i] = new TH1D(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),
                               Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
   fHistSumOfWeightsEtaRP[i] -> SetYTitle("sum (*)");
   fHistSumOfWeightsEtaRP[i] -> SetXTitle("#eta");
   fHistList->Add(fHistSumOfWeightsEtaRP[i]);
  
   fHistSumOfWeightsPtPOI[i] = new TH1D(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),
                               Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
   fHistSumOfWeightsPtPOI[i] -> SetYTitle("sum (*)");
   fHistSumOfWeightsPtPOI[i] -> SetXTitle("p_{T}");
   fHistList->Add(fHistSumOfWeightsPtPOI[i]);
 
   fHistSumOfWeightsEtaPOI[i] = new TH1D(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),
                                Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
   fHistSumOfWeightsEtaPOI[i] -> SetYTitle("sum (*)");
   fHistSumOfWeightsEtaPOI[i] -> SetXTitle("#eta");
   fHistList->Add(fHistSumOfWeightsEtaPOI[i]);
  }
      
  fCommonHistsSP = new AliFlowCommonHist("AliFlowCommonHistSP");
  fHistList->Add(fCommonHistsSP);
  fCommonHistsResSP = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
  fHistList->Add(fCommonHistsResSP);  
  fCommonHistsmuQ = new AliFlowCommonHist("AliFlowCommonHistmuQ");
  fHistList->Add(fCommonHistsmuQ);

  (fCommonHistsSP->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic 
  (fCommonHistsmuQ->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic 

  fHistQNorm = new TH1D("Flow_QNorm_SP","Flow_QNorm_SP",110,0.,1.1);
  fHistQNorm -> SetYTitle("dN/d(|(Qa+Qb)/(Ma+Mb)|)");
  fHistQNorm -> SetXTitle("|(Qa+Qb)/(Ma+Mb)|");
  fHistList->Add(fHistQNorm);

  fHistQaQb = new TH1D("Flow_QaQb_SP","Flow_QaQb_SP",200,-100.,100.);
  fHistQaQb -> SetYTitle("dN/dQaQb");
  fHistQaQb -> SetXTitle("QaQb");
  fHistList->Add(fHistQaQb);

  fHistQaQbNorm = new TH1D("Flow_QaQbNorm_SP","Flow_QaQbNorm_SP",44,-1.1,1.1);
  fHistQaQbNorm -> SetYTitle("dN/d(QaQb/MaMb)");
  fHistQaQbNorm -> SetXTitle("QaQb/MaMb");
  fHistList->Add(fHistQaQbNorm);

  fHistQNormvsQaQbNorm = new TH2D("Flow_QNormvsQaQbNorm_SP","Flow_QNormvsQaQbNorm_SP",88,-1.1,1.1,22,0.,1.1);
  fHistQNormvsQaQbNorm -> SetYTitle("|Q/Mq|");
  fHistQNormvsQaQbNorm -> SetXTitle("QaQb/MaMb");
  fHistList->Add(fHistQNormvsQaQbNorm);

  fHistQaQbCos = new TH1D("Flow_QaQbCos_SP","Flow_QaQbCos_SP",63,0.,TMath::Pi());
  fHistQaQbCos -> SetYTitle("dN/d(#phi)");
  fHistQaQbCos -> SetXTitle("#phi");
  fHistList->Add(fHistQaQbCos);

  fHistResolution = new TH1D("Flow_resolution_SP","Flow_resolution_SP",100,-1.0,1.0);
  fHistResolution -> SetYTitle("dN/d(cos(2(#phi_a - #phi_b))");
  fHistResolution -> SetXTitle("cos(2*(#phi_a - #phi_b))");
  fHistList->Add(fHistResolution);

  fHistQaNorm = new TH1D("Flow_QaNorm_SP","Flow_QaNorm_SP",22,0.,1.1);
  fHistQaNorm -> SetYTitle("dN/d(|Qa/Ma|)");
  fHistQaNorm -> SetXTitle("|Qa/Ma|");
  fHistList->Add(fHistQaNorm);

  fHistQaNormvsMa = new TH2D("Flow_QaNormvsMa_SP","Flow_QaNormvsMa_SP",100,0.,100.,22,0.,1.1);
  fHistQaNormvsMa -> SetYTitle("|Qa/Ma|");
  fHistQaNormvsMa -> SetXTitle("Ma");
  fHistList->Add(fHistQaNormvsMa);

  fHistQbNorm = new TH1D("Flow_QbNorm_SP","Flow_QbNorm_SP",22,0.,1.1);
  fHistQbNorm -> SetYTitle("dN/d(|Qb/Mb|)");
  fHistQbNorm -> SetXTitle("|Qb/Mb|");
  fHistList->Add(fHistQbNorm);

  fHistQbNormvsMb = new TH2D("Flow_QbNormvsMb_SP","Flow_QbNormvsMb_SP",100,0.,100.,22,0.,1.1);
  fHistQbNormvsMb -> SetYTitle("|Qb/Mb|");
  fHistQbNormvsMb -> SetXTitle("|Mb|");
  fHistList->Add(fHistQbNormvsMb);

  fHistMavsMb = new TH2D("Flow_MavsMb_SP","Flow_MavsMb_SP",100,0.,100.,100,0.,100.);
  fHistMavsMb -> SetYTitle("Ma");
  fHistMavsMb -> SetXTitle("Mb");
  fHistList->Add(fHistMavsMb);


  //weights
  if(fUsePhiWeights) {
    if(!fWeightsList) {
      cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
      exit(0);  
    }
    if(fWeightsList->FindObject("phi_weights_sub0"))  {
      fPhiWeightsSub0 = dynamic_cast<TH1F*>
	(fWeightsList->FindObject("phi_weights_sub0"));
      fHistList->Add(fPhiWeightsSub0);
    } else {
      cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
      exit(0);
    }
    if(fWeightsList->FindObject("phi_weights_sub1"))  {
      fPhiWeightsSub1 = dynamic_cast<TH1F*>
	(fWeightsList->FindObject("phi_weights_sub1"));
      fHistList->Add(fPhiWeightsSub1);
    } else {
      cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
      exit(0);
    }

  } // end of if(fUsePhiWeights)

  fEventNumber = 0;  //set number of events to zero 
  
  //store all boolean flags needed in Finish():
  this->StoreFlags();   

  TH1::AddDirectory(oldHistAddStatus);
}

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


  if (anEvent) {

  //Calculate muQ (for comparing pp and PbPb)
  FillmuQ(anEvent);

  //Calculate flow based on  <QaQb/MaMb> = <v^2>
  FillSP(anEvent);

  }
}

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

  //Calculate flow based on  <QaQb/MaMb> = <v^2>

  //Fill histograms
  if (anEvent) {

    //get Q vectors for the eta-subevents
    AliFlowVector* vQarray = new AliFlowVector[2];
    if (fUsePhiWeights) {
      anEvent->Get2Qsub(vQarray,fHarmonic,fWeightsList,kTRUE);
    } else {
      anEvent->Get2Qsub(vQarray,fHarmonic);
    }
    //subevent a
    AliFlowVector vQa = vQarray[0];
    //subevent b
    AliFlowVector vQb = vQarray[1];

    //For calculating v2 only events should be taken where both subevents are not empty
    //check that the subevents are not empty:
    Double_t dMa = vQa.GetMult();
    Double_t dMb = vQb.GetMult();
    if (dMa > 0. && dMb > 0.) {
      
      //request that the subevent multiplicities are not too different
      //fRelDiffMsub can be set from the configuration macro
      Double_t dRelDiff = TMath::Abs((dMa - dMb)/(dMa + dMb));
      if (dRelDiff < fRelDiffMsub) {

	//fill control histograms 	   
	if (fUsePhiWeights) {
	  fCommonHistsSP->FillControlHistograms(anEvent,fWeightsList,kTRUE);
	} else {
	  fCommonHistsSP->FillControlHistograms(anEvent);
	}

	//fill some SP control histograms
	fHistProQaQb -> Fill(1.,vQa*vQb,1.); //Fill with weight 1 -> Weight with MaMb????
	fHistQaQbCos ->Fill(TMath::ACos((vQa/vQa.Mod())*(vQb/vQb.Mod())));  //Acos(Qa*Qb) = angle
	fHistResolution -> Fill(TMath::Cos( vQa.Phi()- vQb.Phi() ));  //vQa.Phi() returns 2*phi
	fHistQaQb -> Fill(vQa*vQb);
	fHistMavsMb -> Fill(dMb,dMa);

	//get total Q vector = the sum of subevent a and subevent b
	AliFlowVector vQ;
	if(!strcmp(fTotalQvector->Data(),"QaQb"))
	{
	 vQ = vQa + vQb;
	} else if(!strcmp(fTotalQvector->Data(),"Qa"))
	  {
	   vQ = vQa; 
	  } else if(!strcmp(fTotalQvector->Data(),"Qb"))
	    {
	     vQ = vQb; 
	    }

	//needed to correct for non-uniform acceptance:
	fHistProNonIsotropicTermsQ->Fill(1.,vQ.Y()/(dMa+dMb),dMa+dMb);
	fHistProNonIsotropicTermsQ->Fill(2.,vQ.X()/(dMa+dMb),dMa+dMb);

	//weight the Q vectors for the subevents by the multiplicity
	//Note: Weight Q only in the particle loop when it is clear 
	//if it should be (m-1) or M
	Double_t dQXa = vQa.X()/dMa; 
	Double_t dQYa = vQa.Y()/dMa;
	vQa.Set(dQXa,dQYa);
	
	Double_t dQXb = vQb.X()/dMb; 
	Double_t dQYb = vQb.Y()/dMb;
	vQb.Set(dQXb,dQYb);
        
	//scalar product of the two subevents
	Double_t dQaQb = (vQa*vQb);
	fHistProQaQbNorm -> Fill(1.,dQaQb,dMa*dMb);  //Fill (QaQb/MaMb) with weight (MaMb). 
	//needed for the error calculation:
	fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
	fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
	//needed for correcting non-uniform acceptance: 
	fHistProQaQbReImNorm->Fill(1.,dQYa,dMa); // to get <<sin(phi_a)>>
	fHistProQaQbReImNorm->Fill(2.,dQXa,dMa); // to get <<cos(phi_a)>>
	fHistProQaQbReImNorm->Fill(3.,dQYb,dMb); // to get <<sin(phi_b)>>
	fHistProQaQbReImNorm->Fill(4.,dQXb,dMb); // to get <<cos(phi_b)>>
	
	//fill some SP control histograms
	fHistQaQbNorm ->Fill(vQa*vQb);
	fHistQaNorm ->Fill(vQa.Mod());
	fHistQaNormvsMa->Fill(dMa,vQa.Mod());
	fHistQbNorm ->Fill(vQb.Mod());
	fHistQbNormvsMb->Fill(dMb,vQb.Mod());
	
	//loop over the tracks of the event
	AliFlowTrackSimple*   pTrack = NULL; 
	Int_t iNumberOfTracks = anEvent->NumberOfTracks(); 
	Double_t dMq =  vQ.GetMult();
	
	for (Int_t i=0;i<iNumberOfTracks;i++) 
	  {
	    pTrack = anEvent->GetTrack(i) ; 
	    if (pTrack){
	      Double_t dPhi = pTrack->Phi();
	      Double_t dWeightUQ = 1.; // weight for u*Q 	    
	      //calculate vU
	      TVector2 vU;
	      //do not need to use weight for v as the length will be made 1
	      Double_t dUX = TMath::Cos(fHarmonic*dPhi);
	      Double_t dUY = TMath::Sin(fHarmonic*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;
	    
	      //redefine the Q vector and devide by its multiplicity
	      TVector2 vQm;
	      Double_t dQmX = 0.;
	      Double_t dQmY = 0.;
	      //subtract particle from the flowvector if used to define it
	      if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) { 
		//set default phi weight to 1
		Double_t dW = 1.; 
		//if phi weights are used
		if(fUsePhiWeights && fPhiWeightsSub0 && fPhiWeightsSub1) 
		{
		  if(strcmp(fTotalQvector->Data(),"QaQb"))
		  {
		   printf("\n WARNING (SP): If you use phi-weights total Q-vector has to be Qa+Qb in the current implementation!!!! \n");
		   exit(0);
		  }
		  //value of the center of the phi bin
		  Double_t dPhiCenter = 0.;  
		  if (pTrack->InSubevent(0) ) {
		    Int_t iNBinsPhiSub0 = fPhiWeightsSub0->GetNbinsX();
		    Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub0/TMath::TwoPi()));
		    dW = fPhiWeightsSub0->GetBinContent(phiBin); 
		    dPhiCenter = fPhiWeightsSub0->GetBinCenter(phiBin);
		    dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter) )/(dMq-dW*pTrack->Weight());
		    dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter) )/(dMq-dW*pTrack->Weight());
		    
		    vQm.Set(dQmX,dQmY);
		  }

		  else if ( pTrack->InSubevent(1)) { 
		    Int_t iNBinsPhiSub1 = fPhiWeightsSub1->GetNbinsX();
		    Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub1/TMath::TwoPi()));
		    dW = fPhiWeightsSub1->GetBinContent(phiBin);
		    dPhiCenter = fPhiWeightsSub1->GetBinCenter(phiBin);
		    dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter) )/(dMq-dW*pTrack->Weight());
		    dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter) )/(dMq-dW*pTrack->Weight());
		    
		    vQm.Set(dQmX,dQmY);
		  }
		  //bin = 1 + value*nbins/range
		  //TMath::Floor rounds to the lower integer
		}     
		// if no phi weights are used
		else 
		{
		 if(!strcmp(fTotalQvector->Data(),"QaQb"))
		 {
		  dQmX = (vQ.X() - (pTrack->Weight())*dUX)/(dMq-pTrack->Weight());
		  dQmY = (vQ.Y() - (pTrack->Weight())*dUY)/(dMq-pTrack->Weight());
		  dWeightUQ = dMq-pTrack->Weight();
		  vQm.Set(dQmX,dQmY);
		 } else if((!strcmp(fTotalQvector->Data(),"Qa") && pTrack->InSubevent(0)) ||
		           (!strcmp(fTotalQvector->Data(),"Qb") && pTrack->InSubevent(1)))
		   {
		    dQmX = (vQ.X() - (pTrack->Weight())*dUX)/(dMq-pTrack->Weight());
		    dQmY = (vQ.Y() - (pTrack->Weight())*dUY)/(dMq-pTrack->Weight());
		    dWeightUQ = dMq-pTrack->Weight();
		    vQm.Set(dQmX,dQmY);
		   } else if((!strcmp(fTotalQvector->Data(),"Qa") && pTrack->InSubevent(1)) ||
		             (!strcmp(fTotalQvector->Data(),"Qb") && pTrack->InSubevent(0)))
		     {
   		      dQmX = vQ.X()/dMq;
		      dQmY = vQ.Y()/dMq;
		      dWeightUQ = dMq;
		      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,dWeightUQ); //Fill (Qu/(Mq-1)) with weight (Mq-1) 
		  //needed for the error calculation:
		  fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dWeightUQ*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb	    
		  fHistProUQPtRP -> Fill(dPt,dUQ,dWeightUQ);                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
		  fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dWeightUQ*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb	
		  
{'fraction_non_alphanumeric': 0.18960895455196547, 'fraction_numerical': 0.011868464670854733, 'mean_word_length': 3.815252734259533, 'pattern_counts': {'":': 0, '<': 133, '<?xml version=': 0, '>': 503, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4890799', 'n_tokens_mistral': 20701, 'n_tokens_neox': 19614, 'n_words': 2678}