Coding violation corrected (G.Galazka)
[u/mrichter/AliRoot.git] / HBTAN / AliHBTWeightasCorrFctn.cxx
index 0430be7cf9021e444b8a14e78f4d903d07bf4505..54cc2aebf757e97762a6949b7fe003bb478b7d08 100644 (file)
@@ -1,6 +1,6 @@
 #include "AliHBTWeightasCorrFctn.h"
 #include <TH1.h>
-#include <TObjArray.h>
+#include <Riostream.h>
 
 ///////////////////////////////////////////////////////
 //                                                   //
@@ -129,11 +129,12 @@ void AliHBTWeightasCorrFctn::Init()
 
 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());
-     Double_t qout=GetValue(trackpair, partpair);
+     trackpair = CheckPair(trackpair);
+     
 
      int n = fNumberOfIntervals;
      
@@ -143,7 +144,8 @@ void AliHBTWeightasCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, Al
      TH1D &den = *((TH1D*)fDen->At(ntv));
      if ( trackpair && partpair)
      {
-          den.Fill(qout);
+     Double_t qout=GetValue(trackpair, partpair);      
+            den.Fill(qout);    
      }
 }
 
@@ -151,14 +153,16 @@ void AliHBTWeightasCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, Al
 /******************************************************************/
 
 
-void AliHBTWeightasCorrFctn::SetParams(Int_t nbins, Float_t maxXval, Float_t minXval){
+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;
  }
@@ -177,7 +181,7 @@ AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
 
 void AliHBTQOutWeightasCorrFctn::BuildHistos()
 {
-    
+//buils histograms (allocates memory etc)
      Int_t i;
      int n=GetNumberOfIntervals();
 
@@ -189,13 +193,13 @@ void AliHBTQOutWeightasCorrFctn::BuildHistos()
      
 
      
-     TH1D *Num;
-     TH1D *Den;
-     TH1D *Rat;
+     TH1D *num;
+     TH1D *den;
+     TH1D *rat;
      
-     TString nameNum = "NumOut";
-     TString nameDen = "DenOut";
-     TString nameRat = "RatOut";
+     TString nameNum = "WeightNumOut";
+     TString nameDen = "WeightDenOut";
+     TString nameRat = "WeightRatOut";
      
      for(i=0;i<n;i++){
          
@@ -207,25 +211,25 @@ void AliHBTQOutWeightasCorrFctn::BuildHistos()
          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 = 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->Sumw2();
+         den->Sumw2();
+         rat->Sumw2();
          
-         Num->Reset();
-         Den->Reset();
-         Rat->Reset();
+         num->Reset();
+         den->Reset();
+         rat->Reset();
          
-         fNum->Add(Num);
-         fDen->Add(Den);
-         fRat->Add(Rat);
+         fNum->Add(num);
+         fDen->Add(den);
+         fRat->Add(rat);
          
-         nameNum = TString("NumOut");
-         nameDen = TString("DenOut");
-         nameRat = TString("RatOut");
+         nameNum = TString("WeightNumOut");
+         nameDen = TString("WeightDenOut");
+         nameRat = TString("WeightRatOut");
          
      }
 
@@ -243,7 +247,7 @@ AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
 
 void AliHBTQSideWeightasCorrFctn::BuildHistos()
 {
-    
+//builds histograms
      Int_t i;
      int n=GetNumberOfIntervals();
      int nbins=Getnbins();
@@ -254,13 +258,14 @@ void AliHBTQSideWeightasCorrFctn::BuildHistos()
      
 
      
-     TH1D *Num;
-     TH1D *Den;
-     TH1D *Rat;
+     TH1D *num;
+     TH1D *den;
+     TH1D *rat;
      
-     TString nameNum = "NumSide";
-     TString nameDen = "DenSide";
-     TString nameRat = "RatSide";
+     
+     TString nameNum = "WeightNumSide";
+     TString nameDen = "WeightDenSide";
+     TString nameRat = "WeightRatSide";
      
      for(i=0;i<n;i++){
          
@@ -272,25 +277,25 @@ void AliHBTQSideWeightasCorrFctn::BuildHistos()
          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 = 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->Sumw2();
+         den->Sumw2();
+         rat->Sumw2();
          
-         Num->Reset();
-         Den->Reset();
-         Rat->Reset();
+         num->Reset();
+         den->Reset();
+         rat->Reset();
          
-         fNum->Add(Num);
-         fDen->Add(Den);
-         fRat->Add(Rat);
+         fNum->Add(num);
+         fDen->Add(den);
+         fRat->Add(rat);
          
-         nameNum = TString("NumSide");
-         nameDen = TString("DenSide");
-         nameRat = TString("RatSide");
+         nameNum = TString("WeightNumSide");
+         nameDen = TString("WeightDenSide");
+         nameRat = TString("WeightRatSide");
          
      }
 
@@ -308,7 +313,7 @@ AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
 
 void AliHBTQLongWeightasCorrFctn::BuildHistos()
 {
-    
+//builds histograms
      Int_t i;
      int n=GetNumberOfIntervals();
      int nbins=Getnbins();
@@ -319,13 +324,13 @@ void AliHBTQLongWeightasCorrFctn::BuildHistos()
      
 
      
-     TH1D *Num;
-     TH1D *Den;
-     TH1D *Rat;
+     TH1D *num;
+     TH1D *den;
+     TH1D *rat;
      
-     TString nameNum = "NumLong";
-     TString nameDen = "DenLong";
-     TString nameRat = "RatLong";
+     TString nameNum = "WeightNumLong";
+     TString nameDen = "WeightDenLong";
+     TString nameRat = "WeightRatLong";
      
      for(i=0;i<n;i++){
          
@@ -337,27 +342,327 @@ void AliHBTQLongWeightasCorrFctn::BuildHistos()
          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 = 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->Sumw2();
+         den->Sumw2();
+         rat->Sumw2();
          
-         Num->Reset();
-         Den->Reset();
-         Rat->Reset();
+         num->Reset();
+         den->Reset();
+         rat->Reset();
          
-         fNum->Add(Num);
-         fDen->Add(Den);
-         fRat->Add(Rat);
+         fNum->Add(num);
+         fDen->Add(den);
+         fRat->Add(rat);
          
-         nameNum = TString("NumLong");
-         nameDen = TString("DenLong");
-         nameRat = TString("RatLong");
+         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;
+}
+