[u/mrichter/AliRoot.git] / HBTAN / AliHBTWeightasCorrFctn.cxx

#include "AliHBTWeightasCorrFctn.h"
#include <TH1.h>
#include <Riostream.h>

///////////////////////////////////////////////////////
//                                                   //
// AliHBTWeightasCorrFctn.h                             //
//                                                   //
// Class for calculating 3D Weightas correlation        //
// functions                                         //
//                                                   //
///////////////////////////////////////////////////////

ClassImp(AliHBTWeightasCorrFctn)

     AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const char* name, const char* title):
 AliHBTTwoPairFctn1D(name,title),

     fNum(0x0),
     fDen(0x0),
     fRat(0x0)

{
//ctor
}
     
/******************************************************************/
AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
 AliHBTTwoPairFctn1D(name,title,nbins,maxXval,minXval),


     fNum(new TObjArray()),
     fDen(new TObjArray()),
     fRat(new TObjArray())
{
     SetParams(nbins,maxXval, minXval);
}

/******************************************************************/
AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const AliHBTWeightasCorrFctn& in):
 AliHBTTwoPairFctn1D(in),


     fNum((in.fNum)?(TObjArray*)in.fNum->Clone():0x0),
     fDen((in.fDen)?(TObjArray*)in.fDen->Clone():0x0),
     fRat((in.fRat)?(TObjArray*)in.fRat->Clone():0x0)
 {
//ctor
}

/******************************************************************/

AliHBTWeightasCorrFctn::~AliHBTWeightasCorrFctn()
{
 //dtor

     delete fNum;
     delete fDen;
     delete fRat;
     
}

/******************************************************************/
void AliHBTWeightasCorrFctn::Write()
{
//out    
     Int_t i;
//     Int_t n=GetNumberOfIntervals();
     Double_t scale;

     for(i=0;i<fNumberOfIntervals;i++){
	  TH1D *num = ((TH1D*)fNum->At(i));
	  TH1D *den = ((TH1D*)fDen->At(i));
	  TH1D &rat = *((TH1D*)fRat->At(i));
	  scale = Scale(num,den);
	  Info("Write():","Scale in interval %d = %lf",i,scale);
	  rat.Divide(num,den,scale);
	  
	  
	  num->Write();
	  den->Write();
	  rat.Write();
     }

}

//-------------------------------------
void AliHBTWeightasCorrFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
{
    //Fills the numerator using pair from the same event
     trackpair = CheckPair(trackpair);
     if(partpair == 0x0) return;
     if(trackpair == 0x0) return;
     
     Double_t weight = 1.0;

     int n = fNumberOfIntervals;
     Double_t rplane=0.;   //reaction plane angle - 2 B determined
     Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka

     phi=phi*360/(2*TMath::Pi());
     Double_t q=GetValue(trackpair, partpair);
     Int_t ntv;
     ntv =  (int)(phi*n/(360.));
     
     TH1D *num = ((TH1D*)fNum->At(ntv));
     
        
     if ( trackpair && partpair)
     {
	  if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) &&
	       ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode())    )
	  {
	       weight=partpair->GetWeight();
//	       Info("ProcessSameEvent","weight=%lf",weight);
	  }
	  num->Fill(q,weight);
     }
}

/****************************************************************/
void AliHBTWeightasCorrFctn::Init()
{
     BuildHistos();
}
/****************************************************************/


void AliHBTWeightasCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
{
     //process particles from different events 
     Double_t rplane=0.;   //reaction plane angle - 2 B determined
     Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
     phi=phi*360/(2*TMath::Pi());
     trackpair = CheckPair(trackpair);
     

     int n = fNumberOfIntervals;
     
     Int_t ntv;
     ntv =  (int)(phi*n/(360.));

     TH1D &den = *((TH1D*)fDen->At(ntv));
     if ( trackpair && partpair)
     {
     Double_t qout=GetValue(trackpair, partpair);	
	     den.Fill(qout);	
     }
}


/******************************************************************/


void AliHBTWeightasCorrFctn::SetParams(Int_t nbins, Float_t maxXval, Float_t minXval)
{
//set histogram parameters
     fnbins=nbins;
     fmaxXval= maxXval;
     fminXval=minXval;
}
TH1* AliHBTWeightasCorrFctn::GetResult()
{
//just for compatibility       
     TH1D *den = ((TH1D*)fDen->UncheckedAt(1));
     return den;
 }


ClassImp(AliHBTQOutWeightasCorrFctn)
     
     AliHBTQOutWeightasCorrFctn::AliHBTQOutWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)

{
//ct0r
}

void AliHBTQOutWeightasCorrFctn::BuildHistos()
{
//buils histograms (allocates memory etc)
     Int_t i;
     int n=GetNumberOfIntervals();

     int nbins=Getnbins();

     double max = GetmaxXval();
     double min = GetminXval();
     char buff[10];
     

     TH1D *num;
     TH1D *den;
     TH1D *rat;
     
     TString nameNum = "WeightNumOut";
     TString nameDen = "WeightDenOut";
     TString nameRat = "WeightRatOut";
     
     for(i=0;i<n;i++){
	  
	  sprintf(buff,"%d",i);

	  nameNum +=TString(buff);

	  nameDen +=TString(buff);
	  nameRat +=TString(buff);
	  
	  
	  num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
	  den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
	  rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
	  
	  num->Sumw2();
	  den->Sumw2();
	  rat->Sumw2();
	  
	  num->Reset();
	  den->Reset();
	  rat->Reset();
	  
	  fNum->Add(num);
	  fDen->Add(den);
	  fRat->Add(rat);
	  
	  nameNum = TString("WeightNumOut");
	  nameDen = TString("WeightDenOut");
	  nameRat = TString("WeightRatOut");
	  
     }

     
 }

ClassImp(AliHBTQSideWeightasCorrFctn)
     
     AliHBTQSideWeightasCorrFctn::AliHBTQSideWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)

{
//ct0r
}

void AliHBTQSideWeightasCorrFctn::BuildHistos()
{
//builds histograms
     Int_t i;
     int n=GetNumberOfIntervals();
     int nbins=Getnbins();

     double max = GetmaxXval();
     double min = GetminXval();
     char buff[10];
     

     TH1D *num;
     TH1D *den;
     TH1D *rat;
     
     
     TString nameNum = "WeightNumSide";
     TString nameDen = "WeightDenSide";
     TString nameRat = "WeightRatSide";
     
     for(i=0;i<n;i++){
	  
	  sprintf(buff,"%d",i);

	  nameNum +=TString(buff);

	  nameDen +=TString(buff);
	  nameRat +=TString(buff);
	  
	  
	  num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
	  den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
	  rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
	  
	  num->Sumw2();
	  den->Sumw2();
	  rat->Sumw2();
	  
	  num->Reset();
	  den->Reset();
	  rat->Reset();
	  
	  fNum->Add(num);
	  fDen->Add(den);
	  fRat->Add(rat);
	  
	  nameNum = TString("WeightNumSide");
	  nameDen = TString("WeightDenSide");
	  nameRat = TString("WeightRatSide");
	  
     }

     
 }

ClassImp(AliHBTQLongWeightasCorrFctn)
     
     AliHBTQLongWeightasCorrFctn::AliHBTQLongWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)

{
//ct0r
}

void AliHBTQLongWeightasCorrFctn::BuildHistos()
{
//builds histograms
     Int_t i;
     int n=GetNumberOfIntervals();
     int nbins=Getnbins();

     double max = GetmaxXval();
     double min = GetminXval();
     char buff[10];
     

     TH1D *num;
     TH1D *den;
     TH1D *rat;
     
     TString nameNum = "WeightNumLong";
     TString nameDen = "WeightDenLong";
     TString nameRat = "WeightRatLong";
     
     for(i=0;i<n;i++){
	  
	  sprintf(buff,"%d",i);

	  nameNum +=TString(buff);

	  nameDen +=TString(buff);
	  nameRat +=TString(buff);
	  
	  
	  num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
	  den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
	  rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
	  
	  num->Sumw2();
	  den->Sumw2();
	  rat->Sumw2();
	  
	  num->Reset();
	  den->Reset();
	  rat->Reset();
	  
	  fNum->Add(num);
	  fDen->Add(den);
	  fRat->Add(rat);
	  
	  nameNum = TString("WeightNumLong");
	  nameDen = TString("WeightDenLong");
	  nameRat = TString("WeightRatLong");
	  
     }

     
 }

/********************************************************************/
ClassImp(AliHBTasWeightQOSLCorrFctn)
     AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const char* name, const char* title):
 AliHBTTwoPairFctn3D(name,title),

    fNum(0x0),
     fDen(0x0),
     fRat(0x0)
//	fNumberOfIntervals(1)
    

{
//ctor
}
     
/******************************************************************/
AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const char* name, const char* title, Int_t nXbins, Double_t maxXval, Double_t minXval, Int_t nYbins, Double_t maxYval, Double_t minYval,Int_t nZbins, Double_t maxZval, Double_t minZval):
 AliHBTTwoPairFctn3D(name,title,nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval),


     fNum(new TObjArray()),
     fDen(new TObjArray()),
     fRat(new TObjArray())
//	fNumberOfIntervals(1)
     
{
     SetParams(nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval);
}

/******************************************************************/
AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const AliHBTasWeightQOSLCorrFctn& in):
 AliHBTTwoPairFctn3D(in),


     fNum((in.fNum)?(TObjArray*)in.fNum->Clone():0x0),
     fDen((in.fDen)?(TObjArray*)in.fDen->Clone():0x0),
     fRat((in.fRat)?(TObjArray*)in.fRat->Clone():0x0)
 {
//ctor
}

/******************************************************************/

AliHBTasWeightQOSLCorrFctn::~AliHBTasWeightQOSLCorrFctn()
{
 //dtor

     delete fNum;
     delete fDen;
     delete fRat;
     
}

/******************************************************************/
void AliHBTasWeightQOSLCorrFctn::Write()
{
//out    
     Int_t i;
//     Int_t n=GetNumberOfIntervals();
    // Double_t scale;

     for(i=0;i<fNumberOfIntervals;i++){
Info("Write()","%d",i);
	     TH3D *num = ((TH3D*)fNum->UncheckedAt(i));
	  TH3D *den = ((TH3D*)fDen->UncheckedAt(i));
	  TH3D &rat = *((TH3D*)fRat->At(i));
	 // scale = Scale(num,den);
//	  Info("Write():","Scale in interval %d = %lf",i,scale);
	  rat.Divide(num,den);//,scale);
		  

	  num->Write();
	  den->Write();
	  rat.Write();
     }

}

//-------------------------------------
void AliHBTasWeightQOSLCorrFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
{
    //Fills the numerator using pair from the same event
     trackpair = CheckPair(trackpair);
     if(partpair == 0x0) return;
     if(trackpair == 0x0) return;
     
     Double_t weight = 1.0;

     int n = fNumberOfIntervals;
     Double_t rplane=0.;   //reaction plane angle - 2 B determined
     Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka

     phi=phi*360/(2*TMath::Pi());
     Double_t out = TMath::Abs(trackpair->GetQOutLCMS());
     Double_t side = TMath::Abs(trackpair->GetQSideLCMS());
     Double_t lon = TMath::Abs(trackpair->GetQLongLCMS());
	       
     
     Int_t ntv;
     ntv =  (int)(phi*n/(360.));
     
     TH3D *num = ((TH3D*)fNum->At(ntv));
     
        
     if ( trackpair && partpair)
     {
	  if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) &&
	       ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode())    )
	  {
	       weight=partpair->GetWeight();
//	       Info("ProcessSameEvent","weight=%lf",weight);
	  }
	  num->Fill(out,side,lon,weight);
     }
}

/****************************************************************/
void AliHBTasWeightQOSLCorrFctn::Init()
{
     BuildHistos();
}
/****************************************************************/


void AliHBTasWeightQOSLCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
{
//fill denominator     
     Double_t rplane=0.;   //reaction plane angle - 2 B determined
     Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
     phi=phi*360/(2*TMath::Pi());

     Double_t out = TMath::Abs(trackpair->GetQOutLCMS());
     Double_t side = TMath::Abs(trackpair->GetQSideLCMS());
     Double_t lon = TMath::Abs(trackpair->GetQLongLCMS());

     int n = fNumberOfIntervals;
     
     Int_t ntv;
     ntv =  (int)(phi*n/(360.));

     TH3D &den = *((TH3D*)fDen->At(ntv));
     if ( trackpair && partpair)
     {
          den.Fill(out,side,lon);
     }
}


/******************************************************************/


void AliHBTasWeightQOSLCorrFctn::SetParams(Int_t nXbins, Float_t maxXval, Float_t minXval,Int_t nYbins, Float_t maxYval, Float_t minYval,Int_t nZbins, Float_t maxZval, Float_t minZval)
{
//set histogram parametera
     fnXbins=nXbins;
     fmaxXval= maxXval;
     fminXval=minXval;
     fnYbins=nYbins;
     fmaxYval= maxYval;
     fminYval=minYval;

     fnZbins=nZbins;
     fmaxZval= maxZval;
     fminZval=minZval;


}
TH1* AliHBTasWeightQOSLCorrFctn::GetResult()
{
 //just for compatibility      
     TH1D *den = ((TH1D*)fDen->UncheckedAt(1));
     return den;
 }

void AliHBTasWeightQOSLCorrFctn::BuildHistos()
{
    //builds histograms
     Int_t i;
     int n=GetNumberOfIntervals();
     int nXbins=GetnXbins();
     double maxX = GetmaxXval();
     double minX = GetminXval();
     
     int nYbins=GetnYbins();
     double maxY = GetmaxXval();
     double minY = GetminXval();
     int nZbins=GetnZbins();
     double maxZ = GetmaxZval();
     double minZ = GetminZval();
     char buff[10];
        
     TH3D *num;
     TH3D *den;
     TH3D *rat;
     
     TString nameNum = "OSLWeightNum";
     TString nameDen = "OSLWeightDen";
     TString nameRat = "OSLWeightRat";
     
     for(i=0;i<n;i++){
	  
	  sprintf(buff,"%d",i);

	  nameNum +=TString(buff);

	  nameDen +=TString(buff);
	  nameRat +=TString(buff);
	  
	  
	  num = new TH3D(nameNum.Data(),nameNum.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
	  den = new TH3D(nameDen.Data(),nameDen.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
	  rat = new TH3D(nameRat.Data(),nameRat.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
	  
	  num->Sumw2();
	  den->Sumw2();
	  rat->Sumw2();
	  
	  num->Reset();
	  den->Reset();
	  rat->Reset();
	  
	  fNum->Add(num);
	  fDen->Add(den);
	  fRat->Add(rat);
	  
	  nameNum = TString("OSLWeightNumLong");
	  nameDen = TString("OSLWeightDenLong");
	  nameRat = TString("OSLWeightRatLong");
	  
     }
}
Double_t AliHBTasWeightQOSLCorrFctn::Scale(TH3D *num, TH3D *den)
{
//scale histograms method, the default one didn't like me so I've done this one :)
if (AliVAODParticle::GetDebug()>0) Info("Scale","Enetered Scale()");
  if(!num)
     {
         Error("Scale","No numerator");
         return 0.0;
     }
  if(!den)
     {
         Error("Scale","No denominator");
         return 0.0;
     }
                                       
   if( (fNBinsToScaleX < 1) || (fNBinsToScaleY < 1) || (fNBinsToScaleZ < 1))
    {
         Error("Scale","Number of bins for scaling is smaller thnan 1");
         return 0.0;
					    }
    UInt_t nbinsX = num->GetNbinsX();
  if (fNBinsToScaleX > nbinsX)
    {
        Error("Scale","Number of X bins for scaling is bigger thnan number of bins in histograms");
     return 0.0;
    }
		     
    UInt_t nbinsY = num->GetNbinsX();
   if (fNBinsToScaleY > nbinsY)
        {
         Error("Scale","Number of Y bins for scaling is bigger thnan number of bins in histograms");
        return 0.0;
        }

   UInt_t nbinsZ = num->GetNbinsZ();
  if (fNBinsToScaleZ > nbinsZ)
    {
      Error("Scale","Number of Z bins for scaling is bigger thnan number of bins in histograms");
   return 0.0;
     }
		    
  if (AliVAODParticle::GetDebug()>0) Info("Scale","No errors detected");
    Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
    Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
    Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
    Double_t densum = 0.0;
    Double_t numsum = 0.0;

   for (UInt_t k = offsetZ; k<nbinsZ; k++)
     for (UInt_t j = offsetY; j<nbinsY; j++)
       for (UInt_t i = offsetX; i<nbinsX; i++)
         {
            if ( num->GetBinContent(i,j,k) > 0.0 )
              {
              densum += den->GetBinContent(i,j,k);
              numsum += num->GetBinContent(i,j,k);
               }
         }
     if(AliVAODParticle::GetDebug() > 0)
     Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
              numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ);
   if (numsum == 0) return 0.0;
     Double_t ret = densum/numsum;
  if(AliVAODParticle::GetDebug() > 0) Info("Scale","returning %f",ret);
  return ret;
}
Commit	Line	Data
94a709e1	1	#include "AliHBTWeightasCorrFctn.h"
94a709e1	2	#include <TH1.h>
265f04f7	3	#include <Riostream.h>
94a709e1	4
	5	///////////////////////////////////////////////////////
	6	// //
	7	// AliHBTWeightasCorrFctn.h //
	8	// //
	9	// Class for calculating 3D Weightas correlation //
	10	// functions //
	11	// //
	12	///////////////////////////////////////////////////////
	13
	14	ClassImp(AliHBTWeightasCorrFctn)
	15
	16	AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const char* name, const char* title):
	17	AliHBTTwoPairFctn1D(name,title),
	18
	19	fNum(0x0),
	20	fDen(0x0),
	21	fRat(0x0)
	22
	23	{
	24	//ctor
	25	}
	26
	27	/******************************************************************/
	28	AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
	29	AliHBTTwoPairFctn1D(name,title,nbins,maxXval,minXval),
	30
	31
	32	fNum(new TObjArray()),
	33	fDen(new TObjArray()),
	34	fRat(new TObjArray())
	35	{
	36	SetParams(nbins,maxXval, minXval);
	37	}
	38
	39	/******************************************************************/
	40	AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const AliHBTWeightasCorrFctn& in):
	41	AliHBTTwoPairFctn1D(in),
	42
	43
	44
	45	fNum((in.fNum)?(TObjArray*)in.fNum->Clone():0x0),
	46	fDen((in.fDen)?(TObjArray*)in.fDen->Clone():0x0),
	47	fRat((in.fRat)?(TObjArray*)in.fRat->Clone():0x0)
	48	{
	49	//ctor
	50	}
	51
	52	/******************************************************************/
	53
	54	AliHBTWeightasCorrFctn::~AliHBTWeightasCorrFctn()
	55	{
	56	//dtor
	57
	58	delete fNum;
	59	delete fDen;
	60	delete fRat;
	61
	62	}
	63
	64	/******************************************************************/
	65	void AliHBTWeightasCorrFctn::Write()
	66	{
	67	//out
68	Int_t i;
69	// Int_t n=GetNumberOfIntervals();
70	Double_t scale;
71
72	for(i=0;i<fNumberOfIntervals;i++){
73	TH1D num = ((TH1D)fNum->At(i));
74	TH1D den = ((TH1D)fDen->At(i));
75	TH1D &rat = ((TH1D)fRat->At(i));
76	scale = Scale(num,den);
77	Info("Write():","Scale in interval %d = %lf",i,scale);
78	rat.Divide(num,den,scale);
79
80
81	num->Write();
82	den->Write();
83	rat.Write();
84	}
85
86	}
87
88	//-------------------------------------
89	void AliHBTWeightasCorrFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
90	{
91	//Fills the numerator using pair from the same event
92	trackpair = CheckPair(trackpair);
93	if(partpair == 0x0) return;
94	if(trackpair == 0x0) return;
95
96	Double_t weight = 1.0;
97
98	int n = fNumberOfIntervals;
99	Double_t rplane=0.; //reaction plane angle - 2 B determined
100	Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
101
102	phi=phi360/(2TMath::Pi());
103	Double_t q=GetValue(trackpair, partpair);
104	Int_t ntv;
105	ntv = (int)(phi*n/(360.));
106
107	TH1D num = ((TH1D)fNum->At(ntv));
108
109
110	if ( trackpair && partpair)
111	{
112	if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) &&
113	( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) )
114	{
115	weight=partpair->GetWeight();
116	// Info("ProcessSameEvent","weight=%lf",weight);
117	}
118	num->Fill(q,weight);
119	}
120	}
121
122	/****************************************************************/
123	void AliHBTWeightasCorrFctn::Init()
124	{
125	BuildHistos();
126	}
127	/****************************************************************/
128
129
130	void AliHBTWeightasCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
131	{
265f04f7	132	//process particles from different events
94a709e1	133	Double_t rplane=0.; //reaction plane angle - 2 B determined
	134	Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
	135	phi=phi360/(2TMath::Pi());
265f04f7	136	trackpair = CheckPair(trackpair);
265f04f7	137
94a709e1	138
	139	int n = fNumberOfIntervals;
	140
	141	Int_t ntv;
	142	ntv = (int)(phi*n/(360.));
	143
	144	TH1D &den = ((TH1D)fDen->At(ntv));
	145	if ( trackpair && partpair)
	146	{
265f04f7	147	Double_t qout=GetValue(trackpair, partpair);
265f04f7	148	den.Fill(qout);
94a709e1	149	}
	150	}
	151
	152
	153	/******************************************************************/
	154
	155
265f04f7	156	void AliHBTWeightasCorrFctn::SetParams(Int_t nbins, Float_t maxXval, Float_t minXval)
	157	{
	158	//set histogram parameters
94a709e1	159	fnbins=nbins;
	160	fmaxXval= maxXval;
	161	fminXval=minXval;
	162	}
	163	TH1* AliHBTWeightasCorrFctn::GetResult()
	164	{
265f04f7	165	//just for compatibility
94a709e1	166	TH1D den = ((TH1D)fDen->UncheckedAt(1));
	167	return den;
	168	}
	169
	170
	171
	172
	173	ClassImp(AliHBTQOutWeightasCorrFctn)
	174
	175	AliHBTQOutWeightasCorrFctn::AliHBTQOutWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
	176	AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
	177
	178	{
	179	//ct0r
	180	}
	181
	182	void AliHBTQOutWeightasCorrFctn::BuildHistos()
	183	{
265f04f7	184	//buils histograms (allocates memory etc)
94a709e1	185	Int_t i;
	186	int n=GetNumberOfIntervals();
	187
	188	int nbins=Getnbins();
	189
	190	double max = GetmaxXval();
	191	double min = GetminXval();
	192	char buff[10];
	193
	194
	195
265f04f7	196	TH1D *num;
	197	TH1D *den;
	198	TH1D *rat;
94a709e1	199
265f04f7	200	TString nameNum = "WeightNumOut";
	201	TString nameDen = "WeightDenOut";
	202	TString nameRat = "WeightRatOut";
94a709e1	203
	204	for(i=0;i<n;i++){
	205
	206	sprintf(buff,"%d",i);
	207
	208	nameNum +=TString(buff);
	209
	210	nameDen +=TString(buff);
	211	nameRat +=TString(buff);
	212
	213
265f04f7	214	num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
	215	den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
	216	rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
94a709e1	217
265f04f7	218	num->Sumw2();
	219	den->Sumw2();
	220	rat->Sumw2();
94a709e1	221
265f04f7	222	num->Reset();
	223	den->Reset();
	224	rat->Reset();
94a709e1	225
265f04f7	226	fNum->Add(num);
	227	fDen->Add(den);
	228	fRat->Add(rat);
94a709e1	229
265f04f7	230	nameNum = TString("WeightNumOut");
	231	nameDen = TString("WeightDenOut");
	232	nameRat = TString("WeightRatOut");
94a709e1	233
	234	}
	235
	236
	237	}
	238
	239	ClassImp(AliHBTQSideWeightasCorrFctn)
	240
	241	AliHBTQSideWeightasCorrFctn::AliHBTQSideWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
	242	AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
	243
	244	{
	245	//ct0r
	246	}
	247
	248	void AliHBTQSideWeightasCorrFctn::BuildHistos()
	249	{
265f04f7	250	//builds histograms
94a709e1	251	Int_t i;
	252	int n=GetNumberOfIntervals();
	253	int nbins=Getnbins();
	254
	255	double max = GetmaxXval();
	256	double min = GetminXval();
	257	char buff[10];
	258
	259
	260
265f04f7	261	TH1D *num;
	262	TH1D *den;
	263	TH1D *rat;
94a709e1	264
265f04f7	265
	266	TString nameNum = "WeightNumSide";
	267	TString nameDen = "WeightDenSide";
	268	TString nameRat = "WeightRatSide";
94a709e1	269
	270	for(i=0;i<n;i++){
	271
	272	sprintf(buff,"%d",i);
	273
	274	nameNum +=TString(buff);
	275
	276	nameDen +=TString(buff);
	277	nameRat +=TString(buff);
	278
	279
265f04f7	280	num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
	281	den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
	282	rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
94a709e1	283
265f04f7	284	num->Sumw2();
	285	den->Sumw2();
	286	rat->Sumw2();
94a709e1	287
265f04f7	288	num->Reset();
	289	den->Reset();
	290	rat->Reset();
94a709e1	291
265f04f7	292	fNum->Add(num);
	293	fDen->Add(den);
	294	fRat->Add(rat);
94a709e1	295
265f04f7	296	nameNum = TString("WeightNumSide");
	297	nameDen = TString("WeightDenSide");
	298	nameRat = TString("WeightRatSide");
94a709e1	299
	300	}
	301
	302
	303	}
	304
	305	ClassImp(AliHBTQLongWeightasCorrFctn)
	306
	307	AliHBTQLongWeightasCorrFctn::AliHBTQLongWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
	308	AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
	309
	310	{
	311	//ct0r
	312	}
	313
	314	void AliHBTQLongWeightasCorrFctn::BuildHistos()
	315	{
265f04f7	316	//builds histograms
94a709e1	317	Int_t i;
	318	int n=GetNumberOfIntervals();
	319	int nbins=Getnbins();
	320
	321	double max = GetmaxXval();
	322	double min = GetminXval();
	323	char buff[10];
	324
	325
	326
265f04f7	327	TH1D *num;
	328	TH1D *den;
	329	TH1D *rat;
94a709e1	330
265f04f7	331	TString nameNum = "WeightNumLong";
	332	TString nameDen = "WeightDenLong";
	333	TString nameRat = "WeightRatLong";
94a709e1	334
	335	for(i=0;i<n;i++){
	336
	337	sprintf(buff,"%d",i);
	338
	339	nameNum +=TString(buff);
	340
	341	nameDen +=TString(buff);
	342	nameRat +=TString(buff);
	343
	344
265f04f7	345	num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
	346	den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
	347	rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
94a709e1	348
265f04f7	349	num->Sumw2();
	350	den->Sumw2();
	351	rat->Sumw2();
94a709e1	352
265f04f7	353	num->Reset();
	354	den->Reset();
	355	rat->Reset();
94a709e1	356
265f04f7	357	fNum->Add(num);
	358	fDen->Add(den);
	359	fRat->Add(rat);
94a709e1	360
265f04f7	361	nameNum = TString("WeightNumLong");
	362	nameDen = TString("WeightDenLong");
	363	nameRat = TString("WeightRatLong");
94a709e1	364
	365	}
	366
	367
	368	}
265f04f7	369
	370	/********************************************************************/
	371	ClassImp(AliHBTasWeightQOSLCorrFctn)
	372	AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const char* name, const char* title):
	373	AliHBTTwoPairFctn3D(name,title),
	374
	375	fNum(0x0),
	376	fDen(0x0),
	377	fRat(0x0)
	378	// fNumberOfIntervals(1)
	379
	380
	381	{
	382	//ctor
	383	}
	384
	385	/******************************************************************/
	386	AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const char* name, const char* title, Int_t nXbins, Double_t maxXval, Double_t minXval, Int_t nYbins, Double_t maxYval, Double_t minYval,Int_t nZbins, Double_t maxZval, Double_t minZval):
	387	AliHBTTwoPairFctn3D(name,title,nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval),
	388
	389
	390	fNum(new TObjArray()),
	391	fDen(new TObjArray()),
	392	fRat(new TObjArray())
	393	// fNumberOfIntervals(1)
	394
	395	{
	396	SetParams(nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval);
	397	}
	398
	399	/******************************************************************/
	400	AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const AliHBTasWeightQOSLCorrFctn& in):
	401	AliHBTTwoPairFctn3D(in),
	402
	403
	404
	405	fNum((in.fNum)?(TObjArray*)in.fNum->Clone():0x0),
	406	fDen((in.fDen)?(TObjArray*)in.fDen->Clone():0x0),
	407	fRat((in.fRat)?(TObjArray*)in.fRat->Clone():0x0)
	408	{
	409	//ctor
	410	}
	411
	412	/******************************************************************/
	413
	414	AliHBTasWeightQOSLCorrFctn::~AliHBTasWeightQOSLCorrFctn()
	415	{
	416	//dtor
	417
	418	delete fNum;
	419	delete fDen;
	420	delete fRat;
	421
	422	}
	423
	424	/******************************************************************/
	425	void AliHBTasWeightQOSLCorrFctn::Write()
	426	{
	427	//out
	428	Int_t i;
	429	// Int_t n=GetNumberOfIntervals();
	430	// Double_t scale;
	431
	432	for(i=0;i<fNumberOfIntervals;i++){
433	Info("Write()","%d",i);
434	TH3D num = ((TH3D)fNum->UncheckedAt(i));
435	TH3D den = ((TH3D)fDen->UncheckedAt(i));
436	TH3D &rat = ((TH3D)fRat->At(i));
437	// scale = Scale(num,den);
438	// Info("Write():","Scale in interval %d = %lf",i,scale);
439	rat.Divide(num,den);//,scale);
440
441
442	num->Write();
443	den->Write();
444	rat.Write();
445	}
446
447	}
448
449	//-------------------------------------
450	void AliHBTasWeightQOSLCorrFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
451	{
452	//Fills the numerator using pair from the same event
453	trackpair = CheckPair(trackpair);
454	if(partpair == 0x0) return;
455	if(trackpair == 0x0) return;
456
457	Double_t weight = 1.0;
458
459	int n = fNumberOfIntervals;
460	Double_t rplane=0.; //reaction plane angle - 2 B determined
461	Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
462
463	phi=phi360/(2TMath::Pi());
464	Double_t out = TMath::Abs(trackpair->GetQOutLCMS());
465	Double_t side = TMath::Abs(trackpair->GetQSideLCMS());
466	Double_t lon = TMath::Abs(trackpair->GetQLongLCMS());
467
468
469	Int_t ntv;
470	ntv = (int)(phi*n/(360.));
471
472	TH3D num = ((TH3D)fNum->At(ntv));
473
474
475	if ( trackpair && partpair)
476	{
477	if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) &&
478	( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) )
479	{
480	weight=partpair->GetWeight();
481	// Info("ProcessSameEvent","weight=%lf",weight);
482	}
483	num->Fill(out,side,lon,weight);
484	}
485	}
486
487	/****************************************************************/
488	void AliHBTasWeightQOSLCorrFctn::Init()
489	{
490	BuildHistos();
491	}
492	/****************************************************************/
493
494
495	void AliHBTasWeightQOSLCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
496	{
497	//fill denominator
498	Double_t rplane=0.; //reaction plane angle - 2 B determined
499	Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
500	phi=phi360/(2TMath::Pi());
501
502	Double_t out = TMath::Abs(trackpair->GetQOutLCMS());
503	Double_t side = TMath::Abs(trackpair->GetQSideLCMS());
504	Double_t lon = TMath::Abs(trackpair->GetQLongLCMS());
505
506	int n = fNumberOfIntervals;
507
508	Int_t ntv;
509	ntv = (int)(phi*n/(360.));
510
511	TH3D &den = ((TH3D)fDen->At(ntv));
512	if ( trackpair && partpair)
513	{
514	den.Fill(out,side,lon);
515	}
516	}
517
518
519	/******************************************************************/
520
521
522	void AliHBTasWeightQOSLCorrFctn::SetParams(Int_t nXbins, Float_t maxXval, Float_t minXval,Int_t nYbins, Float_t maxYval, Float_t minYval,Int_t nZbins, Float_t maxZval, Float_t minZval)
523	{
524	//set histogram parametera
525	fnXbins=nXbins;
526	fmaxXval= maxXval;
527	fminXval=minXval;
528	fnYbins=nYbins;
529	fmaxYval= maxYval;
530	fminYval=minYval;
531
532	fnZbins=nZbins;
533	fmaxZval= maxZval;
534	fminZval=minZval;
535
536
537	}
538	TH1* AliHBTasWeightQOSLCorrFctn::GetResult()
539	{
540	//just for compatibility
541	TH1D den = ((TH1D)fDen->UncheckedAt(1));
542	return den;
543	}
544
545	void AliHBTasWeightQOSLCorrFctn::BuildHistos()
546	{
547	//builds histograms
548	Int_t i;
549	int n=GetNumberOfIntervals();
550	int nXbins=GetnXbins();
551	double maxX = GetmaxXval();
552	double minX = GetminXval();
553
554	int nYbins=GetnYbins();
555	double maxY = GetmaxXval();
556	double minY = GetminXval();
557	int nZbins=GetnZbins();
558	double maxZ = GetmaxZval();
559	double minZ = GetminZval();
560	char buff[10];
561
562	TH3D *num;
563	TH3D *den;
564	TH3D *rat;
565
566	TString nameNum = "OSLWeightNum";
567	TString nameDen = "OSLWeightDen";
568	TString nameRat = "OSLWeightRat";
569
570	for(i=0;i<n;i++){
571
572	sprintf(buff,"%d",i);
573
574	nameNum +=TString(buff);
575
576	nameDen +=TString(buff);
577	nameRat +=TString(buff);
578
579
580	num = new TH3D(nameNum.Data(),nameNum.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
581	den = new TH3D(nameDen.Data(),nameDen.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
582	rat = new TH3D(nameRat.Data(),nameRat.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
583
584	num->Sumw2();
585	den->Sumw2();
586	rat->Sumw2();
587
588	num->Reset();
589	den->Reset();
590	rat->Reset();
591
592	fNum->Add(num);
593	fDen->Add(den);
594	fRat->Add(rat);
595
596	nameNum = TString("OSLWeightNumLong");
597	nameDen = TString("OSLWeightDenLong");
598	nameRat = TString("OSLWeightRatLong");
599
600	}
601	}
602	Double_t AliHBTasWeightQOSLCorrFctn::Scale(TH3D num, TH3D den)
603	{
604	//scale histograms method, the default one didn't like me so I've done this one :)
605	if (AliVAODParticle::GetDebug()>0) Info("Scale","Enetered Scale()");
606	if(!num)
607	{
608	Error("Scale","No numerator");
609	return 0.0;
610	}
611	if(!den)
612	{
613	Error("Scale","No denominator");
614	return 0.0;
615	}
616
617	if( (fNBinsToScaleX < 1) \|\| (fNBinsToScaleY < 1) \|\| (fNBinsToScaleZ < 1))
618	{
619	Error("Scale","Number of bins for scaling is smaller thnan 1");
620	return 0.0;
621	}
622	UInt_t nbinsX = num->GetNbinsX();
623	if (fNBinsToScaleX > nbinsX)
624	{
625	Error("Scale","Number of X bins for scaling is bigger thnan number of bins in histograms");
626	return 0.0;
627	}
628
629	UInt_t nbinsY = num->GetNbinsX();
630	if (fNBinsToScaleY > nbinsY)
631	{
632	Error("Scale","Number of Y bins for scaling is bigger thnan number of bins in histograms");
633	return 0.0;
634	}
635
636	UInt_t nbinsZ = num->GetNbinsZ();
637	if (fNBinsToScaleZ > nbinsZ)
638	{
639	Error("Scale","Number of Z bins for scaling is bigger thnan number of bins in histograms");
640	return 0.0;
641	}
642
643	if (AliVAODParticle::GetDebug()>0) Info("Scale","No errors detected");
644	Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
645	Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
646	Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
647	Double_t densum = 0.0;
648	Double_t numsum = 0.0;
649
650	for (UInt_t k = offsetZ; k<nbinsZ; k++)
651	for (UInt_t j = offsetY; j<nbinsY; j++)
652	for (UInt_t i = offsetX; i<nbinsX; i++)
653	{
654	if ( num->GetBinContent(i,j,k) > 0.0 )
655	{
656	densum += den->GetBinContent(i,j,k);
657	numsum += num->GetBinContent(i,j,k);
658	}
659	}
660	if(AliVAODParticle::GetDebug() > 0)
661	Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
662	numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ);
663	if (numsum == 0) return 0.0;
664	Double_t ret = densum/numsum;
665	if(AliVAODParticle::GetDebug() > 0) Info("Scale","returning %f",ret);
666	return ret;
667	}
668