]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Fixes of older problems and adaptation to heavy ions.
authormisko <misko@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 10 Dec 2010 09:11:23 +0000 (09:11 +0000)
committermisko <misko@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 10 Dec 2010 09:11:23 +0000 (09:11 +0000)
- AliUnicorHN: Write replaced by Save, for reading use Retrieve rather
  than constructor, FindBin replaced by FindFixBin
- AliUnicorAnal: saving adapted correspondingly
- AliUnicorAnalCorrel: charges passed to two-track cut, bins modified,
  bimo has now matching bins to pair, qinv calculated and stored for the
  Coulomb correction (ugly, hopefully can be removed in future)
- AliUnicorCoulomb: added
- AliUnicorEventAliceESD: centrality using VZERO, cuts adjusted to ions
- AliTaskUnicor: list SetOwner
- AliTaskUnicorME: list SetOwner

15 files changed:
PWG2/CMakelibPWG2unicor.pkg
PWG2/UNICOR/AliAnalysisTaskUnicor.cxx
PWG2/UNICOR/AliAnalysisTaskUnicorME.cxx
PWG2/UNICOR/AliUnicorAnal.cxx
PWG2/UNICOR/AliUnicorAnalCorrel.cxx
PWG2/UNICOR/AliUnicorAnalCorrel.h
PWG2/UNICOR/AliUnicorAnalGlobal.cxx
PWG2/UNICOR/AliUnicorCoulomb.cxx [new file with mode: 0644]
PWG2/UNICOR/AliUnicorCoulomb.h [new file with mode: 0644]
PWG2/UNICOR/AliUnicorEvent.h
PWG2/UNICOR/AliUnicorEventAliceESD.cxx
PWG2/UNICOR/AliUnicorEventAliceESD.h
PWG2/UNICOR/AliUnicorHN.cxx
PWG2/UNICOR/AliUnicorHN.h
PWG2/libPWG2unicor.pkg

index c34ad0767bc74bc07d67c042daa8cfba94270c59..921ad679345c4a094b5dd06d7365047fe4c1821e 100644 (file)
@@ -25,7 +25,7 @@
 # SHLIBS - Shared Libraries and objects for linking (Executables only)           #
 #--------------------------------------------------------------------------------#
 
-set ( SRCS  UNICOR/AliAnalysisTaskUnicor.cxx UNICOR/AliUnicorAnalCorrel.cxx UNICOR/AliUnicorAnal.cxx UNICOR/AliUnicorAnalGlobal.cxx UNICOR/AliUnicorAnalHighpt.cxx UNICOR/AliUnicorAnalPtfluc.cxx UNICOR/AliUnicorAnalSingle.cxx UNICOR/AliUnicorEventAliceESD.cxx UNICOR/AliUnicorEvent.cxx UNICOR/AliUnicorHN.cxx UNICOR/AliUnicorPair.cxx)
+set ( SRCS  UNICOR/AliAnalysisTaskUnicor.cxx UNICOR/AliUnicorAnalCorrel.cxx UNICOR/AliUnicorAnal.cxx UNICOR/AliUnicorAnalGlobal.cxx UNICOR/AliUnicorAnalHighpt.cxx UNICOR/AliUnicorAnalPtfluc.cxx UNICOR/AliUnicorAnalSingle.cxx UNICOR/AliUnicorEventAliceESD.cxx UNICOR/AliUnicorEvent.cxx UNICOR/AliUnicorHN.cxx UNICOR/AliUnicorPair.cxx UNICOR/AliUnicorCoulomb.cxx)
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
 
index ea70fd256594f89eada690d8c3c9bab3b1cc5b01..162448627d4fd65d1ff03e7990e6b1b4d9c17840 100644 (file)
@@ -38,6 +38,7 @@ AliAnalysisTaskUnicor::AliAnalysisTaskUnicor(const char *name) :
   // constructor
 
   fEv0 = new AliUnicorEventAliceESD(); // needed for eta ranges only
+  fOutputList->SetOwner();
   DefineOutput(1, TList::Class());
 }
 //=============================================================================
@@ -50,7 +51,7 @@ void AliAnalysisTaskUnicor::UserCreateOutputObjects()
   fOutputList->Add(new AliUnicorAnalSingle("all",fEv0->Etamin(),fEv0->Etamax(),0));
   fOutputList->Add(new AliUnicorAnalSingle("pim",fEv0->Etamin(),fEv0->Etamax(),-211));
   fOutputList->Add(new AliUnicorAnalSingle("pip",fEv0->Etamin(),fEv0->Etamax(), 211));
-  int frame = AliUnicorAnalCorrel::kLCMS;
+  AliUnicorAnalCorrel::AnalysisFrame frame = AliUnicorAnalCorrel::kLCMS;
   fOutputList->Add(new AliUnicorAnalCorrel("cnn",fEv0->Etamin(),fEv0->Etamax(),-211,-211, frame));
   fOutputList->Add(new AliUnicorAnalCorrel("cpp",fEv0->Etamin(),fEv0->Etamax(), 211, 211, frame));
   fOutputList->Add(new AliUnicorAnalCorrel("cnp",fEv0->Etamin(),fEv0->Etamax(),-211, 211, frame));
index a074b25ea1c5d27bc880b2ca20ccac03c2f3057b..4d1617959125ee4f0771622f24a327442669e72f 100644 (file)
@@ -41,6 +41,7 @@ AliAnalysisTaskUnicorME::AliAnalysisTaskUnicorME(const char *name) :
 
   fEv0 = new AliUnicorEventAliceESD();
   fEv1 = new AliUnicorEventAliceESD();
+  fOutputList->SetOwner();
   DefineOutput(1, TList::Class());
 }
 //=============================================================================
@@ -53,7 +54,7 @@ void AliAnalysisTaskUnicorME::UserCreateOutputObjects()
   fOutputList->Add(new AliUnicorAnalSingle("all",fEv0->Etamin(),fEv0->Etamax(),0));
   fOutputList->Add(new AliUnicorAnalSingle("pim",fEv0->Etamin(),fEv0->Etamax(),-211));
   fOutputList->Add(new AliUnicorAnalSingle("pip",fEv0->Etamin(),fEv0->Etamax(), 211));
-  int frame = AliUnicorAnalCorrel::kLCMS;
+  AliUnicorAnalCorrel::AnalysisFrame frame = AliUnicorAnalCorrel::kLCMS;
   fOutputList->Add(new AliUnicorAnalCorrel("cnn",fEv0->Etamin(),fEv0->Etamax(),-211,-211, frame));
   fOutputList->Add(new AliUnicorAnalCorrel("cpp",fEv0->Etamin(),fEv0->Etamax(), 211, 211, frame));
   fOutputList->Add(new AliUnicorAnalCorrel("cnp",fEv0->Etamin(),fEv0->Etamax(),-211, 211, frame));
index 2184cbb9823ea1ab2e6490fec81d6a2fafd54566..e9f92fe2950044089747a1957c4fb6619d1c59b0 100644 (file)
@@ -26,6 +26,7 @@
 #include <TCollection.h>
 #include <TClass.h>
 #include <TH1.h>
+#include "AliUnicorHN.h"
 #include "AliUnicorAnal.h"
 
 ClassImp(AliUnicorAnal)
@@ -60,8 +61,8 @@ Long64_t AliUnicorAnal::Merge(const TCollection * const list) {
   return 0;
 }
 //=============================================================================
-void AliUnicorAnal::Save(const char *outfil, const char *mode) 
-{
+void AliUnicorAnal::Save(const char *outfil, const char *mode) {
+
   // store histograms on file in a directory named after the object
   // mode should be "update" (default) or "new"
 
@@ -69,7 +70,11 @@ void AliUnicorAnal::Save(const char *outfil, const char *mode)
   TFile * f = TFile::Open(outfil, mode);
   TDirectory *dest = f->mkdir(GetName());
   dest->cd();
-  fHistos.Write();
+  for (int i=0; i<fHistos.GetEntries(); i++) {
+    TObject *obj = fHistos.At(i);
+    if (obj->IsA()->InheritsFrom("AliUnicorHN")) ((AliUnicorHN*) obj)->Save();
+    else obj->Write();
+  }
   gROOT->cd();
   f->Close();
 }
index 75c2de3afb61c55f0e7c1c756e561779cb319f33..6c696f5b4a84241a9746cf8f5563c52cc251a121 100644 (file)
@@ -33,9 +33,9 @@ ClassImp(AliUnicorAnalCorrel)
  
 //=============================================================================
 AliUnicorAnalCorrel::AliUnicorAnalCorrel(const char *nam, Double_t emi, Double_t ema, 
-                        Int_t pid0, Int_t pid1, Int_t frame): 
-  AliUnicorAnal(nam), fPid0(pid0), fPid1(pid1), fMass0(0), fMass1(0), fFrame(frame), fPa() 
-{
+                        Int_t pid0, Int_t pid1, AnalysisFrame frame): 
+  AliUnicorAnal(nam), fPid0(pid0), fPid1(pid1), fMass0(0), fMass1(0), fZ0(0), fZ1(0), 
+  fFrame(frame), fPa() {
   // constructor
   // emi and ema define the rapidity range for histogram
 
@@ -43,45 +43,53 @@ AliUnicorAnalCorrel::AliUnicorAnalCorrel(const char *nam, Double_t emi, Double_t
   TParticlePDG *part1 = AliUnicorAnal::fgPDG.GetParticle(fPid1);
   fMass0 = part0? part0->Mass() : 0;
   fMass1 = part1? part1->Mass() : 0;
+  fZ0 = part0? part0->Charge()/3.0 : 0;
+  fZ1 = part1? part1->Charge()/3.0 : 0;
   double pi = TMath::Pi();
 
   // correlation function
 
+  int nce = 6; double cebins[]={0,0.05,0.1,0.2,0.4,0.6,1.0};  // centrality bins
+
   //int npt = 7; double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.7,1.0};
   //int npt = 6; double ptbins[]={0,0.1,0.25,0.35,0.55,1.0,2.0}; // like Adam, except last bin split
   //int npt = 7; double ptbins[]={0,0.1,0.25,0.40,0.55,0.7,1.0,2.0}; // like Adam in Mar-2010, + first+last
-  int npt = 8; double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.7,1.0,2.0};  // for Pb
+  int npt = 10; double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,1.0,2.0};  // for Pb
 
   double qbins[200] = {0};
-  //  for (int i=0;i<20;i++) qbins[i]=i*0.005;
+  //  for (int i=0;i<60;i++) qbins[i]=i*0.005;
   //  for (int i=0;i<45;i++) qbins[20+i]=0.1+i*0.02;
-  for (int i=0;i<30;i++) qbins[i]=i*0.010;
-  for (int i=0;i<35;i++) qbins[30+i]=0.3+i*0.02;
-  for (int i=0;i<20;i++) qbins[65+i]=1.0+i*0.05;
-  for (int i=0;i<11;i++) qbins[85+i]=2.0+i*0.20;
+  //  for (int i=0;i<30;i++) qbins[i]=i*0.010;
+  for (int i=0;i<60;i++) qbins[i]=i*0.005;
+  for (int i=0;i<35;i++) qbins[60+i]=0.3+i*0.02;
+  for (int i=0;i<20;i++) qbins[95+i]=1.0+i*0.05;
+  for (int i=0;i<11;i++) qbins[115+i]=2.0+i*0.20;
 
   TAxis *ax[8];
-  ax[0] = new TAxis(3,-0.5,2.5);ax[0]->SetTitle("trumixrot");
-  ax[1] = new TAxis(5,0,1.0);   ax[1]->SetTitle("centrality");
+  ax[0] = new TAxis(4,-0.5,3.5);ax[0]->SetTitle("trumixrot");
+  //  ax[1] = new TAxis(5,0,1.0);   ax[1]->SetTitle("centrality");
+  ax[1] = new TAxis(nce,cebins);ax[1]->SetTitle("centrality");
   ax[2] = new TAxis(3,emi,ema); ax[2]->SetTitle("y");          // pair y
   //  ax[3] = new TAxis(8,-pi,pi);  ax[3]->SetTitle("phi");      // wrt event plane
   ax[3] = new TAxis(1,-pi,pi);  ax[3]->SetTitle("phi");        // wrt event plane
   ax[4] = new TAxis(npt,ptbins);ax[4]->SetTitle("kt (GeV/c)"); // pair pt/2
   ax[5] = new TAxis(8,0,pi);    ax[5]->SetTitle("q-theta");
   ax[6] = new TAxis(16,-pi,pi); ax[6]->SetTitle("q-phi");
-  ax[7] = new TAxis(95,qbins);  ax[7]->SetTitle("q (GeV/c)");
+  ax[7] = new TAxis(125,qbins); ax[7]->SetTitle("q (GeV/c)");
   //  ax[7] = new TAxis(700,0,3.5); ax[7]->SetTitle("q (GeV/c)");
   AliUnicorHN *pair = new AliUnicorHN("pair",8,ax);
-  for (int i=0; i<8; i++) delete ax[i];
   fHistos.Add(pair);
 
   // correlation function bin monitor (needed to get <kt> etc.)
 
-  ax[0] = new TAxis(5,0,1);      ax[0]->SetTitle("centrality");
-  ax[1] = new TAxis(30,emi,ema); ax[1]->SetTitle("y");           // pair y
-  ax[2] = new TAxis(100,0,2);    ax[2]->SetTitle("kt (GeV/c)");  // pair pt/2
-  AliUnicorHN *bimo = new AliUnicorHN("bimo",3,ax);
-  for (int i=0; i<3; i++) delete ax[i];
+  TAxis *bx[3]={0};
+  //bx[0] = new TAxis(*(pair->GetAxis(1))); // wait until root bug (516-00..527-06) fixed. For now, do: 
+  bx[0] = new TAxis(ax[1]->GetNbins(), ax[1]->GetXmin(), ax[1]->GetXmax()); bx[0]->SetTitle("centrality"); 
+  bx[1] = new TAxis(10*ax[2]->GetNbins(),emi,ema); bx[1]->SetTitle("y");         // pair y
+  bx[2] = new TAxis(100,0,2); bx[2]->SetTitle("kt (GeV/c)");                     // pair pt/2
+  AliUnicorHN *bimo = new AliUnicorHN("bimo",3,bx);
+  for (int i=0; i<8; i++) delete ax[i];
+  for (int i=0; i<3; i++) delete bx[i];
   fHistos.Add(bimo);
 
   // two-track resolution monitoring histogram
@@ -89,8 +97,8 @@ AliUnicorAnalCorrel::AliUnicorAnalCorrel(const char *nam, Double_t emi, Double_t
   ax[0] = new TAxis(3,-0.5,2.5);    ax[0]->SetTitle("trumixrot");
   ax[1] = new TAxis(2,-0.5,1.5);    ax[1]->SetTitle("cut applied");
   ax[2] = new TAxis(npt,ptbins);    ax[2]->SetTitle("(pair pt)/2 (GeV)");
-  ax[3] = new TAxis(80,-0.02,0.02); ax[3]->SetTitle("dtheta");
-  ax[4] = new TAxis(80,-0.04,0.04); ax[4]->SetTitle("dphi");
+  ax[3] = new TAxis(80,-0.08,0.08); ax[3]->SetTitle("dtheta");
+  ax[4] = new TAxis(80,-0.20,0.20); ax[4]->SetTitle("dphi");
   AliUnicorHN *twot = new AliUnicorHN("twot",5,ax);
   for (int i=0; i<5; i++) delete ax[i];
   fHistos.Add(twot);
@@ -135,17 +143,34 @@ void AliUnicorAnalCorrel::Process(Int_t tmr, const AliUnicorEvent * const ev0, c
       fPa.Set1(fMass1,ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot);
       if (ev0==ev1 && fPid0==fPid1 && ran.Rndm()>=0.5) fPa.Swap();
       twot->Fill((double) tmr, 0.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
-      if (!ev0->PairGood(ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i),
-                        ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot)) continue;
+      if (!ev0->PairGood(ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i),fZ0,
+                        ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot,fZ1)) continue;
       twot->Fill((double) tmr, 1.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
+      fPa.CalcLAB(); // this could be organized better. AliUnicorPair should give k*?
+      fPa.CalcPairCM();
+      double qcm = fPa.QCM(); // momdif in pair cm - argument for Coulomb correction
+
       fPa.CalcLAB();
       if (fFrame == kPairFrame) fPa.CalcPairCM();
       if (fFrame == kLCMS)      fPa.CalcLcmsCM();
       if (fPa.QCM()==0) {printf("AliUnicorAnalCorrel: Q=0\n"); return;} // should not be too frequent
       double phi = TVector2::Phi_mpi_pi(fPa.Phi()-rpphi);
       double weigth = 1.0;
-      //      if (tmr==0 && fPid0==fPid1) weigth = 1.0+0.6*exp(-pow(fPa.QCM()*5/0.197,2)); 
-      //      if (tmr==0 && fPid0==fPid1) weigth = 1.0+0.6*exp(-fPa.QCM()*2/0.197); 
+      /*
+      static TH2D *coul = 0; 
+      if (!coul) {
+       TFile::Open("coulomb.root","read");
+       coul = (TH2D*) gFile->Get("co");
+       coul->SetDirectory(gROOT);
+       gFile->Close();
+      }
+      if (tmr==0 && fPid0==fPid1) {
+       double co = 0;
+       if (qcm>0.999) co = 1;
+       else if (qcm>0.001) co = coul->Interpolate(7,qcm);
+       weigth = 1.0-0.5+0.5*co*(1+exp(-pow(fPa.QCM()*7/0.197,2))); 
+      }
+      */
       pair->Fill((double) tmr,           // 0 for tru, 1 for mix, 2 for rot
                 cent,                   // centrality
                 fPa.Rapidity(),         // pair rapidity
@@ -155,9 +180,16 @@ void AliUnicorAnalCorrel::Process(Int_t tmr, const AliUnicorEvent * const ev0, c
                 fPa.QCMPhiOut(),        // azimuthal angle of Q w.r.t. out
                 fPa.QCM(),              // |p2-p1| in c.m.s.
                 weigth);                // weigth
-      if (tmr) continue;
-      if (fPa.QCM()>0.2) continue;
-      bimo->Fill(cent, fPa.Rapidity(), fPa.Pt()/2.0, 1.0);
+      if (tmr==0 && fPa.QCM()<0.2) bimo->Fill(cent, fPa.Rapidity(), fPa.Pt()/2.0, weigth);
+      if (tmr==0) pair->Fill((double) 3,             // this is for Coulomb correction, maybe not necessary
+                            cent,                   // centrality
+                            fPa.Rapidity(),         // pair rapidity
+                            phi,                    // pair phi wrt reaction plane
+                            fPa.Pt()/2.0,           // half of pair pt
+                            fPa.QCMTheta(),         // polar angle of Q
+                            fPa.QCMPhiOut(),        // azimuthal angle of Q w.r.t. out
+                            fPa.QCM(),              // |p2-p1| in c.m.s.
+                            weigth*qcm);            // weigth*Q
     }
   }
 }
index b81e513e173ea13574132183ce328d15bedb6de8..10457a9b7d5c9a7cc086710e748d2c0a2ed27910 100644 (file)
@@ -19,10 +19,10 @@ class AliUnicorEvent;
 class AliUnicorAnalCorrel : public AliUnicorAnal {
    
  public:
-  enum AnalysisFrames {kPairFrame, kLCMS, kLAB};
+  enum AnalysisFrame {kPairFrame, kLCMS, kLAB};
   AliUnicorAnalCorrel(const char *nam="correl", Double_t emi=-1, Double_t ema=1, 
-             Int_t pid0=0, Int_t pid1=0, Int_t frame=0); // constructor
-  virtual ~AliUnicorAnalCorrel(){}                                // destructor
+             Int_t pid0=0, Int_t pid1=0, AnalysisFrame frame=kPairFrame); 
+  virtual ~AliUnicorAnalCorrel(){}
   // process one (tru) or two (mix) events
   void Process(Int_t tmr, const AliUnicorEvent * const ev0, const AliUnicorEvent * const ev1, Double_t phirot);
 
@@ -31,6 +31,8 @@ class AliUnicorAnalCorrel : public AliUnicorAnal {
   Int_t    fPid1;                       // particle species 1
   Double_t fMass0;                      // mass 0
   Double_t fMass1;                      // mass 1
+  double   fZ0;                         // charge 0 in units of |e|
+  double   fZ1;                         // charge 1 in units of |e|
   Int_t    fFrame;                      // analysis frame
   AliUnicorPair    fPa;                         // pair buffer for calculations
 
index 8b175220af91d7843d9445774d44dab2f5005e86..4b116e117974c4903146b53b462f0d7da00263a2 100644 (file)
@@ -27,7 +27,6 @@
 #include "AliUnicorEvent.h"
 #include "AliUnicorAnalGlobal.h"
 
-
 class TH2;
 ClassImp(AliUnicorAnalGlobal)
   
@@ -66,9 +65,10 @@ void AliUnicorAnalGlobal::Process(AliUnicorEvent *ev) const
   TH2D *dire = (TH2D*) fHistos.At(3);
   TH1D *zver = (TH1D*) fHistos.At(4);
 
-  mult->Fill(ev->NGoodParticles(),1.0);
+  double n = ev->NGoodParticles();
+  mult->Fill(n,1.0);
   cent->Fill(ev->Centrality(),1.0);
-  cemu->Fill(ev->Centrality(),ev->NGoodParticles());
+  cemu->Fill(ev->Centrality(),n);
   Double_t qx=0,qy=0;
   ev->RP(qx,qy);
   dire->Fill(qx,qy,1.0);
diff --git a/PWG2/UNICOR/AliUnicorCoulomb.cxx b/PWG2/UNICOR/AliUnicorCoulomb.cxx
new file mode 100644 (file)
index 0000000..6c40938
--- /dev/null
@@ -0,0 +1,168 @@
+/************************************************************************* 
+* Copyright(c) 1998-2048, 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.                  * 
+**************************************************************************/
+
+//Author: Dariusz Miskowiec
+//Date:   2010
+
+//=============================================================================
+// Coulomb correlation function
+//=============================================================================
+
+#include <TRandom2.h>
+#include <TComplex.h>
+#include <TFile.h>
+#include <TH2.h>
+#include "AliUnicorCoulomb.h"
+
+ClassImp(AliUnicorCoulomb)
+
+//=============================================================================
+AliUnicorCoulomb::AliUnicorCoulomb(int zz, double mass, double R) : TGraph(1001) {
+
+  // constructor
+  // zz = +1 for pi+pi+, -1 for pi-pi-, etc. 
+  // mass is the reduced mass: m1*m2/(m1+m2)
+
+  SetMarkerStyle(20);
+  int nstep = 10000;
+
+  double rx = R;
+  double ry = R;
+  double rz = R;
+  double r0 = 0;
+  TRandom2 ran;
+
+  if (zz>0) SetPoint(0,0,0);
+  else SetPoint(0,0,1e9);
+  for (int i=1; i<1000; i++) {
+    double qinv = i/1000.0;
+    double k = qinv/2.0;
+    double sum = 0;
+    for (int j=0; j<nstep; j++) {
+      double x0 = ran.Gaus(0.0,rx);
+      double y0 = ran.Gaus(0.0,ry);
+      double z0 = ran.Gaus(0.0,rz);
+      double t0 = ran.Gaus(0.0,r0);
+      double x1 = ran.Gaus(0.0,rx);
+      double y1 = ran.Gaus(0.0,ry);
+      double z1 = ran.Gaus(0.0,rz);
+      double t1 = ran.Gaus(0.0,r0);
+      double dx = x1-x0;
+      double dy = y1-y0;
+      double dz = z1-z0;
+      double dt = t1-t0;
+      dz += dt*k/mass;
+      sum += WaveFunction2(zz,mass,k,dx,dy,dz);
+    }
+    SetPoint(i,qinv,sum/nstep);
+  }
+  SetPoint(1000,1000,1);
+}
+//=============================================================================
+double AliUnicorCoulomb::WaveFunction2(int zz, double mass, double k, double x, double y, double z) {
+
+  // Square of normalized Coulomb wave function. 
+  // Hypergeometrical function diverges for numerical reasons for large Qinv, z and x. 
+  // Out of the convergence limit I will return 1.0
+
+  double qinv12 = 1.2 * 2.0 * k;
+  double p0 = -9.1569/qinv12;
+  double p2 = 2.735e-2*qinv12;
+  double wf2 = 1.0;
+  int converg = (z>p0+p2*(x*x+y*y));
+  if (converg) {
+    TComplex co = WaveFunction(zz,mass,k,x,y,z);
+    wf2 =  co.Rho2()*Gamow(zz,mass,k);
+  }
+  return wf2;
+}
+//=============================================================================
+TComplex AliUnicorCoulomb::WaveFunction(int zz, double mass, double k, double x, double y, double z) {
+
+  // Coulomb wave function in parabolic coordinates
+  // Merzbacher, Quantum Mechanics, p.248
+  // Landau-Lifshitz, Quantum Mechanics (Non-Relativistic Theory), p.518
+
+  double r = sqrt(x*x+y*y+z*z);
+  //  double ksi = r+z;
+  double eta = r-z;
+  TComplex jjj(0,1);
+  TComplex a = -zz*mass/k/137.036*jjj;
+  TComplex b(1.0,0.0);
+  TComplex c = jjj*k*eta/0.197327;
+  TComplex wf = TComplex::Exp(jjj*k*z/0.197327)*F1(a,b,c);
+  //printf("%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f \n",a.Re(),a.Im(),b.Re(),b.Im(),c.Re(),c.Im(),F1(a,b,c).Rho());
+  return wf;
+}
+//=============================================================================
+double AliUnicorCoulomb::Gamow(int zz, double m, double k) {
+
+  // Squared Coulomb wave function in infinity = Gamow
+  // zz - product of charges
+  // m - reduced mass
+  // k = Qinv/2 for equal masses
+
+  double n = -zz*m/k/137.036;
+  double a = 2*TMath::Pi()*n;
+  double ga = a/(1.0-exp(-a));
+  return ga;
+}
+//=============================================================================
+TComplex AliUnicorCoulomb::F1(TComplex alpha, TComplex gamma, TComplex z) {
+
+  // confluent hypergeometric function 1F1
+
+  double prec = 0.0000001;
+  TComplex term(1.0,0.0);
+  TComplex sum(1.0,0.0);
+  for (int n=1; n<1000; n++) {
+    double u = (double) n;
+    TComplex v = (alpha+u-1.0)/(gamma+u-1.0);
+    TComplex w = z/u;
+    term *= v;
+    term *= w;
+    sum += term;
+    //    printf("%10d %10.3f %10.3f %10.3f %10.3f\n", n, v.Rho(), w.Rho(), term.Rho(), sum.Rho());
+    if (fabs(term/sum) < prec) return sum;
+  }
+  printf("F1 Maximum number of iterations reached\n");
+  return sum;
+}
+//=============================================================================
+ void AliUnicorCoulomb::Makehist(int zz, double m, const char *outfil) {
+
+  // Make a two-dim histogram coulomb(R,Q) and store it on outfil.
+  // Later to be used via TH2::Interpolate(). 
+  // zz and m are the product of charges and the reduced mass. 
+
+  TH2D *hi = new TH2D("co","coulomb",11,-0.5,10.5,999,0.0005,0.9995);
+  hi->SetXTitle("R (fm)");
+  hi->SetYTitle("Q (GeV/c)");
+  AliUnicorCoulomb *co = 0;
+  for (int i=1; i<=hi->GetNbinsX(); i++) {
+    double R = hi->GetXaxis()->GetBinCenter(i);
+    printf("R = %.1f fm\n",R); 
+    co = new AliUnicorCoulomb(zz, m, R);
+    for (int j=1; j<=hi->GetNbinsY(); j++) {
+      double Q = hi->GetYaxis()->GetBinCenter(j);
+      hi->SetBinContent(i,j,co->Eval(Q));
+    }
+    delete co;
+  }
+  TFile::Open(outfil,"recreate");
+  hi->Write();
+  gFile->Close();
+}
+//=============================================================================
diff --git a/PWG2/UNICOR/AliUnicorCoulomb.h b/PWG2/UNICOR/AliUnicorCoulomb.h
new file mode 100644 (file)
index 0000000..890dfc2
--- /dev/null
@@ -0,0 +1,36 @@
+#ifndef ALIUNICORCOULOMB_H
+#define ALIUNICORCOULOMB_H
+
+/* Copyright(c) 1998-2048, ALICE Experiment at CERN, All rights reserved. *
+* See cxx source for full Copyright notice */
+/* $Id$ */
+
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2010
+
+//=============================================================================
+// Coulomb correlation function
+//=============================================================================
+
+#include <TComplex.h>
+#include <TGraph.h>
+
+//=============================================================================
+class AliUnicorCoulomb : public TGraph {
+
+ public: 
+  AliUnicorCoulomb(TRootIOCtor *) : TGraph() {}             // default constructor
+  AliUnicorCoulomb(int sign, double mass, double R);        // constructor
+  virtual ~AliUnicorCoulomb() {}                            // destructor
+  double Cf(double qinv) const {return Eval(qinv);} // value of the correlation function
+  static double Gamow(int zz, double m, double k);  // poin source case - Gamow function
+  static void Makehist(int zz, double m, const char *outfil); // make TH2(R,Q)
+
+ protected:
+  static double   WaveFunction2(int zz, double mass, double k, double x, double y, double z);
+  static TComplex WaveFunction(int zz, double mass, double k, double x, double y, double z);
+  static TComplex F1(TComplex alpha, TComplex gamma, TComplex z);
+
+  ClassDef(AliUnicorCoulomb,1)
+};
+//=============================================================================
+#endif
index fc2c7cff0416309a24fca19ff12a2b9d695805b8..c61c0e399746f223dc6fb5fb055f54f935101e82 100644 (file)
@@ -39,8 +39,8 @@ class AliUnicorEvent : public TObject {
   virtual Double_t    ParticleTheta(Int_t i) const = 0;
   virtual Double_t    ParticlePhi(Int_t i) const = 0;
   virtual Double_t    ParticleDedx(Int_t i) const = 0;
-  virtual Bool_t      PairGood(Double_t p0, Double_t the0, Double_t phi0, 
-                              Double_t p1, Double_t the1, Double_t phi1) const = 0;
+  virtual Bool_t      PairGood(double p0, double the0, double phi0, double z0, 
+                              double p1, double the1, double phi1, double z1) const = 0;
 
   // toolkit part
 
index a4c6b1a91c827c843e533967cca67467e1a0ee33..6fa6c8c462ed53218435b2dc674c3054bb9c4833 100644 (file)
 //=============================================================================
 
 #include <cmath>
+#include "TFile.h"
+#include "TH1.h"
+#include "TGraph.h"
+#include "AliESDVZERO.h"
 #include "AliESDEvent.h"
+#include "AliMultiplicity.h"
 #include "AliUnicorEventAliceESD.h"
 
 ClassImp(AliUnicorEventAliceESD)
 
 //=============================================================================
-  AliUnicorEventAliceESD::AliUnicorEventAliceESD(AliESDEvent *esd) : AliUnicorEvent(), fESD(esd)//, fPhysicsSelection(0) 
+AliUnicorEventAliceESD::AliUnicorEventAliceESD(AliESDEvent *esd) : AliUnicorEvent(), fViper(0), fESD(esd)//, fPhysicsSelection(0) 
 {
   // constructor 
 
   //  printf("%s object created\n",ClassName());
   if (!fESD) fESD = new AliESDEvent();
   //  fPhysicsSelection = new AliPhysicsSelection();
+
+  // V0 percentile graph for centrality determination
+
+  //  TFile::Open("/lustre/alice/misko/AliCentralityBy1D_137161_v4.root","read");
+  //  TFile::Open("/lustre/alice/misko/AliCentralityBy1D_137161_v5.root","read"); not checked
+  TFile::Open("/lustre/alice/misko/AliCentralityBy1D_137366_v2.root","read");
+  //TFile::Open("/lustre/alice/misko/AliCentralityBy1D_137366_v4.root","read");
+  const TH1F *hi = (const TH1F*) gFile->Get("hmultV0_percentile");
+  //const TH1D *hi = (const TH1D*) gFile->Get("hNtracks_percentile");
+  //const TH1F *hi = (const TH1F*) gFile->Get("hNclusters1_percentile");
+  fViper = new TGraph((const TH1*) hi);
+  gFile->Close();
 }
 //=============================================================================
-AliUnicorEventAliceESD::~AliUnicorEventAliceESD()
-{
+AliUnicorEventAliceESD::~AliUnicorEventAliceESD() {
+
   // destructor
 
 }
@@ -49,19 +66,36 @@ Bool_t AliUnicorEventAliceESD::Good() const
   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
   if (!vtx->GetStatus()) return kFALSE;
   if (fabs(Zver())>1) return kFALSE;
+  if (NGoodParticles()<9) return kFALSE;
+  return kTRUE;
+}
+//=============================================================================
+Double_t AliUnicorEventAliceESD::Centrality() const {
 
-  /*
-  int hard=0;
-  for (int i=0; i<NParticles(); i++) {
-    if (!ParticleGood(i)) continue;
-    if (ParticlePt(i)<1.0) continue;
-    hard = 1;
-    break;
-  }
-  return hard;
-  */
+  // centrality between 0 (central) and 1 (very peripheral)
 
-  return kTRUE;
+  // V0 multiplicity, not linearized
+  
+  AliESDVZERO *v0 = fESD->GetVZEROData();
+  float multa = v0->GetMTotV0A();
+  float multc = v0->GetMTotV0C();
+  double cent = fViper->Eval(multa+multc)/100.0;
+
+  // number of tracks
+
+  // double cent = fViper->Eval(((AliVEvent*) fESD)->GetNumberOfTracks())/100.0;  
+
+  // number of clusters in second layer of SPD
+
+  //  double nspd2 = fESD->GetMultiplicity()->GetNumberOfITSClusters(1);
+  //  double zv = fESD->GetPrimaryVertexSPD()->GetZ();
+  //  const double pars[] = {8.10030e-01,-2.80364e-03,-7.19504e-04};
+  //  zv -= pars[0];
+  //  float corr = 1 + zv*(pars[1] + zv*pars[2]);
+  //  nspd2 = corr>0? nspd2/corr : -1;
+  //  double cent = fViper->Eval(nspd2)/100.0;  
+
+  return cent;
 }
 //=============================================================================
 Bool_t AliUnicorEventAliceESD::ParticleGood(Int_t i, Int_t pidi) const 
@@ -74,7 +108,7 @@ Bool_t AliUnicorEventAliceESD::ParticleGood(Int_t i, Int_t pidi) const
   AliESDtrack *track = fESD->GetTrack(i);
   if (!track->IsOn(AliESDtrack::kTPCrefit)) return 0;        // TPC refit
   if (!track->IsOn(AliESDtrack::kITSrefit)) return 0;        // ITS refit
-  if (track->GetTPCNcls() < 70) return 0;                    // number of TPC clusters
+  if (track->GetTPCNcls() < 90) return 0;                    // number of TPC clusters
   if (track->GetKinkIndex(0) > 0) return 0;                  // no kink daughters
   const AliExternalTrackParam *tp = GetTrackParam(i);
   if (!tp) return 0;                                         // track param
@@ -89,8 +123,9 @@ Bool_t AliUnicorEventAliceESD::ParticleGood(Int_t i, Int_t pidi) const
 
   Float_t r,z;
   track->GetImpactParametersTPC(r,z);
-  if (fabs(z)>3.2) return 0;                          // impact parameter in z
-  if (fabs(r)>2.4) return 0;                          // impact parameter in xy
+  if (fabs(z)>1.0) return 0;                          // impact parameter in z
+  if (fabs(r)>1.0) return 0;                          // impact parameter in xy
+  if (r==0) return 0;
 
   //TBits shared = track->GetTPCSharedMap();
   //if (shared.CountBits()) return 0;                 // no shared clusters; pragmatic but dangerous
@@ -115,16 +150,26 @@ Bool_t AliUnicorEventAliceESD::ParticleGood(Int_t i, Int_t pidi) const
   else return 0;
 }
 //=============================================================================
-Bool_t AliUnicorEventAliceESD::PairGood(Double_t /*p0*/, Double_t /*the0*/, Double_t /*phi0*/, 
-                               Double_t /*p1*/, Double_t /*the1*/, Double_t /*phi1*/) const {
+Bool_t AliUnicorEventAliceESD::PairGood(double p0, double the0, double phi0, double z0,  
+                               double p1, double the1, double phi1, double z1) const {
 
   // two-track separation cut
 
-  // double dthe = the1-the0;
-  // double dphi = TVector2::Phi_mpi_pi(phi1-phi0);
-  // double dpt = p1*sin(the1) - p0*sin(the0);
-  // return (fabs(dthe)>0.010 || fabs(dphi)>0.060);
-  // return (dpt*dphi<0);
-  return 1;
+  double dthe = the1-the0;
+  if (fabs(dthe)>0.010) return kTRUE;
+  double B = -0.5; // magnetic field in T, needed for helix; should be gotten in the proper way
+  double pt0 = p0*sin(the0);
+  double pt1 = p1*sin(the1);
+  double r = 1.2; // we will calculate the distance between the two tracks at r=1.2 m
+  // for (double r=0.8; r<2.5; r+=0.2) {
+  double si0 = -0.3*B*z0*r/2/pt0;
+  double si1 = -0.3*B*z1*r/2/pt1;
+  if (fabs(si0)>=1.0) return kTRUE; // could be done better
+  if (fabs(si1)>=1.0) return kTRUE;
+  double dphi = phi1 - phi0 + asin(si1) - asin(si0);
+  dphi = TVector2::Phi_mpi_pi(dphi);
+  if (fabs(dphi)<0.020) return kFALSE;
+  // }
+  return kTRUE;
 }
 //=============================================================================
index 79984792cf1786e02e1ba6185b0769299415eaf1..d08e6b622e526429b1d7d9a62f71b771d063653e 100644 (file)
 #include "AliUnicorEvent.h"
 
 //class AliPhysicsSelection;
+class TGraph;
 
 //=============================================================================
 class AliUnicorEventAliceESD : public AliUnicorEvent {
 
  public:
               AliUnicorEventAliceESD(AliESDEvent *esd=0);
-             //AliUnicorEventAliceESD(const AliUnicorEventAliceESD &ev): AliUnicorEvent(ev), fESD(ev.fESD), fPhysicsSelection(ev.fPhysicsSelection){}
-              AliUnicorEventAliceESD(const AliUnicorEventAliceESD &ev): AliUnicorEvent(ev), fESD(ev.fESD) {}
+             //AliUnicorEventAliceESD(const AliUnicorEventAliceESD &ev): AliUnicorEvent(ev), fViper(ev.fViper), fESD(ev.fESD), fPhysicsSelection(ev.fPhysicsSelection){}
+             AliUnicorEventAliceESD(const AliUnicorEventAliceESD &ev): AliUnicorEvent(ev), fViper(ev.fViper), fESD(ev.fESD) {}
   virtual     ~AliUnicorEventAliceESD();
-  AliUnicorEventAliceESD &operator=(const AliUnicorEventAliceESD &source) {fESD = source.fESD; return *this;}
+  AliUnicorEventAliceESD &operator=(const AliUnicorEventAliceESD &source) {fViper=source.fViper; fESD=source.fESD; return *this;}
   Double_t    Etamin() const {return -0.75;}
   Double_t    Etamax() const {return  0.75;}
   void        AttachTree(TTree *tr) {fESD->ReadFromTree(tr);}
   Bool_t      Good() const;
-  Double_t    Centrality() const {return 0.9999*exp(-NGoodParticles()/2000.0);} // 7 for pp 900 GeV, 12 for pp 7 TeV, 20 extr, 2000  PbPb
+  //  Double_t    Centrality() const {return 0.9999*exp(-NGoodParticles()/1000.0);} // 7 for pp 900 GeV, 12 for pp 7 TeV, 1000  PbPb
+  Double_t    Centrality() const; 
   void        RP(Double_t &qx, Double_t &qy) const {AliUnicorEvent::RP(qx,qy,2);}
   Double_t    RPphi() const {Double_t qx,qy; RP(qx,qy); return atan2(qy,qx);}
-  Double_t    Zver() const {return fESD->GetPrimaryVertex()->GetZv()/10.0;}
+  Double_t    Zver() const {return fESD->GetPrimaryVertex()->GetZv()/12.0;}
   Int_t       NParticles() const {return fESD->GetNumberOfTracks();}
-  //  Int_t       NGoodParticles() const {int n=0; for (int i=0; i<fESD->GetMultiplicity()->GetNumberOfTracklets(); i++) if (fabs(fESD->GetMultiplicity()->GetEta(i))<0.8) n++; return n;}
 
   Bool_t      ParticleGood(Int_t i, Int_t pidi=0) const;
   Double_t    ParticleP(Int_t i)     const {return GetTrackParam(i)->P();}
   Double_t    ParticleTheta(Int_t i) const {return GetTrackParam(i)->Theta();}
   Double_t    ParticlePhi(Int_t i)   const {return TVector2::Phi_mpi_pi(GetTrackParam(i)->Phi());}
   Double_t    ParticleDedx(Int_t i)  const {return fESD->GetTrack(i)->GetTPCsignal()/50.0;}
-  Bool_t      PairGood(Double_t p0, Double_t the0, Double_t phi0, 
-                      Double_t p1, Double_t the1, Double_t phi1) const;
+  Bool_t      PairGood(double p0, double the0, double phi0, double z0, 
+                      double p1, double the1, double phi1, double z1) const;
   void        SetESD(AliESDEvent * const esd) {fESD = esd;}
   AliESDEvent *GetESD() const {return fESD;}
   //const AliExternalTrackParam *GetTrackParam(Int_t i) const {return fESD->GetTrack(i);}
@@ -51,6 +52,7 @@ class AliUnicorEventAliceESD : public AliUnicorEvent {
   const AliExternalTrackParam *GetTrackParam(Int_t i) const {return fESD->GetTrack(i)->GetTPCInnerParam();}
 
  protected:
+  TGraph      *fViper;                    // V0 centrality percentile
   AliESDEvent *fESD;                      //! pointer to the actual source of data
   //  AliPhysicsSelection *fPhysicsSelection; //! interaction event filter
 
index 7d041a4d4a65a24f625d856f12c03c83ab5f883d..38fb824624853fee829ab6cac78e86c9452788c6 100644 (file)
 #include <stdlib.h>
 #include <TFile.h>
 #include <TDirectory.h>
+#include <TROOT.h>
 #include <TAxis.h>
 #include <TH2.h>
+#include <TObjArray.h>
+#include <TObjString.h>
+#include <TString.h>
 #include "AliUnicorHN.h"
 
 ClassImp(AliUnicorHN)
@@ -36,64 +40,64 @@ ClassImp(AliUnicorHN)
 //=============================================================================
 AliUnicorHN::AliUnicorHN(const char *nam, Int_t ndim, TAxis **ax) 
   : TH1D(nam, nam, Albins(ndim,ax), 0.5, Albins(ndim,ax)+0.5), 
-    fNdim(ndim)
-{
-  // Constructor for building from scratch.
+    fNdim(ndim) {
+
+  // constructor
  
   // Above, we just have managed to create a single one-dimensional histogram 
   // with number of bins equal to the product of the numbers of bins in all 
   // dimensions. For easy inspection the histogram range was set to -0.5,n-0.5.
 
-  for (int i=0; i<fNdim; i++) ax[i]->SetName(Form("axis%d",i));
   for (int i=0; i<fNdim; i++) ax[i]->Copy(fAxis[i]); 
-
+  for (int i=0; i<fNdim; i++) fAxis[i].SetName(Form("axis%d",i));
+  
   // for speed reasons, number of bins of each axis is stored on fNbins
   // and the running product of the latter on fMbins
   // so fMbins = {...,fNbins[fNdim-2]*fNbins[fNdim-1],fNbins[fNdim-1],1}
-
+  
   for (int i=0; i<fgkMaxNdim; i++) fNbins[i] = fAxis[i].GetNbins();
   for (int i=0; i<fgkMaxNdim; i++) fMbins[i] = 1;
   for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i];
   printf("   %d-dimensional histogram %s with %d bins created\n",fNdim,nam,GetNbinsX());
 }
 //=============================================================================
-AliUnicorHN::AliUnicorHN(const char *filnam, const char *nam) : 
-  TH1D(*(TFile::Open(filnam,"read")?
-        TFile::Open(filnam,"read")->GetDirectory(nam)?
-        (TH1D*) TFile::Open(filnam,"read")->GetDirectory(nam)->Get("histo"):new TH1D():new TH1D()
-        )), 
-  fNdim(0) 
-{
-  // Constructor for reading from file.
+AliUnicorHN* AliUnicorHN::Retrieve(const char *filnam, const char *nam) { 
 
-  TFile *f = TFile::Open(filnam,"read");
-  if (f) if (f->cd(nam)) {
-    TAxis *ax[fgkMaxNdim];
-    for (fNdim=0; fNdim<fgkMaxNdim; fNdim++) {
-      ax[fNdim] = (TAxis *) gDirectory->Get(Form("axis%d",fNdim));
-      if (ax[fNdim]) ax[fNdim]->Copy(fAxis[fNdim]); 
-      else break;
-    }
-    f->Close();
-
-    for (int i=0; i<fgkMaxNdim; i++) fNbins[i] = fAxis[i].GetNbins();
-    for (int i=0; i<fgkMaxNdim; i++) fMbins[i] = 1;
-    for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i];
-    
-    if (GetNbinsX()!=Albins(fNdim,ax)) {
-      printf("number of bins of histo %d differs from product of nbins of axes %d\n",
-            GetNbinsX(),Albins(fNdim,ax));
-      printf("bombing\n");
-      exit(-1);
-    }
+  // retrieve a multidimensional histogram from file
 
-    printf("%d-dimensional histogram %s with %d bins read from file %s\n",
-          fNdim,nam,GetNbinsX(),filnam);
+  TFile *f = TFile::Open(filnam,"read");
+  if (!f) return 0;
+  if (!f->cd(nam)) {f->Close(); return 0;}
+  TH1D *hist = (TH1D*) gDirectory->Get("histo");
+  if (!hist) {printf("No histogram histo on file %s in directory %s\n",filnam,nam); return 0;}
+  hist->SetDirectory(gROOT);
+  TAxis *ax[fgkMaxNdim]={0};
+  int n=0;
+  while ((ax[n] = (TAxis *) gDirectory->Get(Form("axis%d",n)))) n++; 
+  f->Close();
+  
+  if (hist->GetNbinsX()!=Albins(n,ax)) {
+    printf("number of bins of histo %d differs from product of nbins of axes %d\n",
+          hist->GetNbinsX(),Albins(n,ax));
+    return 0;
   }
+
+  // derive the name from the path (ccc from aaa/bbb/ccc)
+  TString path=nam;
+  TObjArray *ar = path.Tokenize("/");
+  TObjString *os = (TObjString*) ar->Last();
+  const char *lastnam = os->GetString().Data();
+
+  AliUnicorHN *hi = new AliUnicorHN(lastnam,n,ax);
+  //  *((TH1D*) hi) = *hist;
+  hist->Copy(*((TH1D*)hi));
+  hi->SetName(lastnam);
+  delete hist;
+  return hi;
 }
 //=============================================================================
-Int_t AliUnicorHN::Albins(Int_t n, TAxis **ax) 
-{
+Int_t AliUnicorHN::Albins(Int_t n, TAxis **ax) {
+
   // Calculate product of nbins of ax[0]...ax[n-1]
   // (= total number of bins of the multidimensional histogram to be made). 
 
@@ -103,8 +107,8 @@ Int_t AliUnicorHN::Albins(Int_t n, TAxis **ax)
   return nbins;
 }
 //=============================================================================
-Int_t AliUnicorHN::MulToOne(const Int_t * const k) const 
-{
+Int_t AliUnicorHN::MulToOne(const Int_t * const k) const {
+
   // Calculate the 1-dim index n from n-dim indices k[fNdim].
   // Valid k[i] should be between 0 and fNbins[i]-1.
   // Valid result will be between 0 and GetNbinsX()-1. 
@@ -119,19 +123,19 @@ Int_t AliUnicorHN::MulToOne(const Int_t * const k) const
   return n;
 }
 //=============================================================================
-Int_t AliUnicorHN::MulToOne(Double_t *x) 
-{
+Int_t AliUnicorHN::MulToOne(Double_t *x) {
+
   // Calculate the 1-dim index n from n-dim vector x, representing the 
   // abscissa of the n-dim histogram. The result will be between 0 and 
   // GetNbinsX()-1. 
 
   Int_t k[fgkMaxNdim] = {0};
-  for (int i=0; i<fNdim; i++) k[i] = fAxis[i].FindBin(x[i])-1;
+  for (int i=0; i<fNdim; i++) k[i] = fAxis[i].FindFixBin(x[i])-1;
   return MulToOne(k);
 }
 //=============================================================================
-void AliUnicorHN::OneToMul(Int_t n, Int_t *k) const 
-{
+void AliUnicorHN::OneToMul(Int_t n, Int_t *k) const {
+
   // Calculate the n-dim indices k[fNdim] from 1-dim index n.
   // Valid n should be between 0 and GetNbinsX()-1. 
   // Valid results will be between 0 and fNbins[i]-1.
@@ -144,8 +148,8 @@ void AliUnicorHN::OneToMul(Int_t n, Int_t *k) const
   }
 }
 //=============================================================================
-Int_t AliUnicorHN::Fill(Double_t *xx, Double_t w)  
-{
+Int_t AliUnicorHN::Fill(Double_t *xx, Double_t w) {
+
   // Fill the histogram. The array xx holds the abscissa information, w is the 
   // weigth. The 1-dim histogram is filled using the standard Fill method 
   // (direct access to the arrays was tried and was not faster). 
@@ -155,21 +159,21 @@ Int_t AliUnicorHN::Fill(Double_t *xx, Double_t w)
   return TH1D::Fill(nbin+1,w); 
 }
 //=============================================================================
-Int_t AliUnicorHN::Fill(Double_t x0, Double_t x1, ...) 
-{
+Int_t AliUnicorHN::Fill(Double_t x0, Double_t x1, ...) {
+
   // Fill the histogram. Arguments are passed as doubles rather than array. 
+
   va_list ap;
   Double_t xx[fgkMaxNdim] = {x0, x1};
   va_start(ap,x1);
   for (int i=2; i<fNdim; i++) xx[i] = va_arg(ap,Double_t);
   Double_t weigth = va_arg(ap,Double_t);
   va_end(ap);
-  //for (int i=0; i<fgkMaxNdim; i++) printf("%8.2f ",xx[i]); printf("%6.2f \n",weigth);
   return Fill(xx,weigth);
 }
 //=============================================================================
-Int_t AliUnicorHN::Write() const  
-{
+Int_t AliUnicorHN::Save() const {
+
   // Save the 1-dim histo and the axes in a subdirectory on file. This might 
   // not be the most elegant way but it is very simple and backward and forward 
   // compatible. 
@@ -189,8 +193,8 @@ Int_t AliUnicorHN::Write() const
   return nbytes;
 }
 //=============================================================================
-AliUnicorHN *AliUnicorHN::ProjectAlong(const char *nam, Int_t dim, Int_t first, Int_t last) 
-{
+AliUnicorHN *AliUnicorHN::ProjectAlong(const char *nam, Int_t dim, Int_t first, Int_t last) {
+
   // Reduce dimension dim by summing up its bins between first and last. 
   // Use root convention: bin=1 is the first bin, bin=nbins is the last. 
   // last=0 means till the last bin
@@ -241,8 +245,9 @@ AliUnicorHN *AliUnicorHN::ProjectAlong(const char *nam, Int_t dim, Int_t first,
   return his;
 }
 //=============================================================================
-TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Int_t * const first, const Int_t * const last) const 
-{
+TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Int_t * const first, 
+                    const Int_t * const last) const {
+
   // Project on dimension dim. Use only bins between first[i] and last[i]. 
   // Use root convention: bin=1 is the first bin, bin=nbins is the last. 
   // first[i]=0 means from the first bin
@@ -275,7 +280,6 @@ TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Int_t * const fir
 
   TH1D *his;
   const char *name = strlen(nam)? nam : Form("%s_proj%d",GetName(),dim);
-  //  if (gDirectory->Get(name)) gDirectory->Get(name)->Delete(); // for some reason leads to troubles
   if (fAxis[dim].IsVariableBinSize()) 
     his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXbins()->GetArray());
   else 
@@ -302,22 +306,24 @@ TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Int_t * const fir
   return his;
 }
 //=============================================================================
-TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Double_t * const first, const Double_t * const last) 
-{
+TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Double_t * const first, 
+                    const Double_t * const last) const {
+
   // Project on dimension dim. Use only bins between first[i] and last[i]. 
 
   if (dim<0 || dim>fNdim-1) return 0;
   Int_t kfirst[fgkMaxNdim] = {0};
   Int_t klast[fgkMaxNdim] = {0};
   for (int i=0; i<fNdim; i++) {
-    kfirst[i] = fAxis[i].FindBin(first[i]);
-    klast[i] = fAxis[i].FindBin(last[i]);
+    kfirst[i] = fAxis[i].FindFixBin(first[i]);
+    klast[i] = fAxis[i].FindFixBin(last[i]);
   }
   return ProjectOn(nam,dim,kfirst,klast);
 }
 //=============================================================================
-TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Int_t * const first, const Int_t * const last) const
-{
+TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Int_t * 
+                    const first, const Int_t * const last) const {
+
   // Project on dim1 vs dim0. Use only bins between first[i] and last[i]. 
   // Use root convention: bin=1 is the first bin, bin=nbins is the last. 
   // first[i]=0 means from the first bin
@@ -396,3 +402,19 @@ TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Int_
   return his;
 }
 //=============================================================================
+TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Double_t * 
+                    const first, const Double_t * const last) const {
+
+  // Project on dim1 vs dim0. Use only bins between first[i] and last[i]. 
+
+  if (dim0<0 || dim0>fNdim-1) return 0;
+  if (dim1<0 || dim1>fNdim-1) return 0;
+  Int_t kfirst[fgkMaxNdim] = {0};
+  Int_t klast[fgkMaxNdim] = {0};
+  for (int i=0; i<fNdim; i++) {
+    kfirst[i] = fAxis[i].FindFixBin(first[i]);
+    klast[i] = fAxis[i].FindFixBin(last[i]);
+  }
+  return ProjectOn(nam,dim0,dim1,kfirst,klast);
+}
+//=============================================================================
index a86f99e2b1122cbe57451fc4cd10fe833081e030..ee0b97f9762b557397888634ec27eb11dc159a03 100644 (file)
@@ -19,9 +19,11 @@ class TAxis;
 class AliUnicorHN : public TH1D {
 
  public:
-  AliUnicorHN(const char *nam="muhi", Int_t ndim=0, TAxis **ax=0); // constructor from scratch
-  AliUnicorHN(const char *filename, const char *name);             // constructor from file
-  virtual ~AliUnicorHN() {}                                        // destructor
+  AliUnicorHN(const char *nam="muhi", Int_t ndim=0, TAxis **ax=0);     // constructor
+  AliUnicorHN(TRootIOCtor *) : TH1D(), fNdim(0) {}                     // default constructor
+  virtual ~AliUnicorHN() {}                                            // destructor
+  static AliUnicorHN* Retrieve(const char *filnam, const char *nam);   // read from file
+
   Int_t GetNdim() const                     {return fNdim;}
   TAxis *GetAxis(Int_t i) const             {return (TAxis*) &fAxis[i];}
 
@@ -31,20 +33,17 @@ class AliUnicorHN : public TH1D {
   Int_t Fill(Double_t x0, Double_t x1, ...);// 2 or more dim histo fill
   Int_t Fill(const char*, Double_t)         {return -1;} // overload TH1
 
-  Int_t Write() const;                      // save histo and axis on file 
-  Int_t Write()                             {return ((const AliUnicorHN*)this)->Write();}
-  Int_t Write(const char *, Int_t, Int_t)   {return Write();} // overload TObject
-  Int_t Write(const char *, Int_t, Int_t) const {return Write();} 
+  Int_t Save() const;                      // save histo and axis on file 
 
   // project along (integrate over) one axis
   AliUnicorHN  *ProjectAlong(const char *nam, Int_t dim, Int_t first=0, Int_t last=0);
   // project on 1-dim histogram
   TH1D *ProjectOn(const char *nam, Int_t dim, const Int_t * const first=0, const Int_t * const last=0) const;
-  // project on 1-dim histogram
-  TH1D *ProjectOn(const char *nam, Int_t dim, const Double_t * const first, const Double_t * const last);
+  TH1D *ProjectOn(const char *nam, Int_t dim, const Double_t * const first, const Double_t * const last) const;
   // project on 2-dim histogram
   TH2D *ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Int_t * const first=0, const Int_t * const last=0) const;
-  void  OneToMul(Int_t n, Int_t *k) const;      // calc n-dim indices from 1-dim index
+  TH2D *ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Double_t * const first, const Double_t * const last) const;
+  void OneToMul(Int_t n, Int_t *k) const;      // calc n-dim indices from 1-dim index
 
  protected:
 
index 8e702e0727e9ba3ca47db1010393e3093a00d0d5..a3ccac8b66638b4945ebbbaf7ec939157ac481de 100644 (file)
@@ -10,7 +10,8 @@ SRCS= UNICOR/AliAnalysisTaskUnicor.cxx \
       UNICOR/AliUnicorEventAliceESD.cxx \
       UNICOR/AliUnicorEvent.cxx \
       UNICOR/AliUnicorHN.cxx \
-      UNICOR/AliUnicorPair.cxx 
+      UNICOR/AliUnicorPair.cxx \
+      UNICOR/AliUnicorCoulomb.cxx 
 
 HDRS= $(SRCS:.cxx=.h)