UNICOR/AliDAnalCorrel.cxx

   1 // Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2005
   2
   3 //=============================================================================
   4 // two-particle correlation analyzer
   5 //=============================================================================
   6
   7 #include <TROOT.h>
   8 #include <TMath.h>
   9 #include <TRandom2.h>
  10 #include "AliDEvent.h"
  11 #include "AliDHN.h"
  12 #include "AliDAnalCorrel.h"
  13
  14 ClassImp(AliDAnalCorrel)
  15
  16 //=============================================================================
  17 AliDAnalCorrel::AliDAnalCorrel(Char_t *nam, Double_t emi, Double_t ema,
  18                          Int_t pid0, Int_t pid1):
  19   AliDAnal(nam), fPid0(pid0), fPid1(pid1), fMass0(0), fMass1(0), fPa()
  20 {
  21   // constructor
  22   // emi and ema define the rapidity range for histogram
  23
  24   TParticlePDG *part0 = AliDAnal::fgPDG.GetParticle(fPid0);
  25   TParticlePDG *part1 = AliDAnal::fgPDG.GetParticle(fPid1);
  26   fMass0 = part0? part0->Mass() : 0;
  27   fMass1 = part1? part1->Mass() : 0;
  28
  29   double pi = TMath::Pi();
  30   TAxis *ax[8];
  31   ax[0] = new TAxis(3,-0.5,2.5);ax[0]->SetTitle("trumixrot");
  32   ax[1] = new TAxis(5,0,0.5);   ax[1]->SetTitle("centrality");
  33   ax[2] = new TAxis(3,emi,ema); ax[2]->SetTitle("pair y");
  34   ax[3] = new TAxis(8,-pi,pi);  ax[3]->SetTitle("pair phi"); // wrt event plane
  35   double a0[]={0,0.1,0.2,0.3,0.4,0.5,0.7,1.0};
  36   ax[4] = new TAxis(7,a0);      ax[4]->SetTitle("(pair pt)/2 (GeV)");
  37   ax[5] = new TAxis(8,0,pi);    ax[5]->SetTitle("q-theta");
  38   ax[6] = new TAxis(16,-pi,pi); ax[6]->SetTitle("q-phi");
  39   double a1[100];
  40   for (int i=0;i<20;i++) a1[i]=i*0.005;
  41   for (int i=0;i<45;i++) a1[20+i]=0.1+i*0.02;
  42   ax[7] = new TAxis(64,a1);     ax[7]->SetTitle("q (GeV/c)");
  43   AliDHN *pair = new AliDHN("pair",8,ax);
  44   for (int i=0; i<8; i++) delete ax[i];
  45   fHistos.Add(pair);
  46   gROOT->cd();
  47   printf("%s object named %s created\n",ClassName(),GetName());
  48 }
  49 //=============================================================================
  50 void AliDAnalCorrel::Process(Int_t tmr, AliDEvent *ev0, AliDEvent *ev1, Double_t phirot)
  51 {
  52   // process pairs from one or two (if mixing) events
  53   // tmr tells which histogram (bins) to fill: tru,mix,rot
  54
  55   static TRandom2 ran;
  56   AliDHN *pair = (AliDHN*) fHistos.At(0);
  57
  58   // mixing-and-rotating-proof centrality and reaction plane angle
  59   // (but not rotation-proof for rotation angles much different from 0 and 180)
  60   // true and rotated pairs are within the triangle (j<i), mixed - all
  61   // thus, proper rotation is either by 180, or by 170 AND 190, etc.
  62
  63   double cent = (ev0->Centrality()+ev1->Centrality())/2.0;
  64   double q0x,q0y,q1x,q1y;
  65   ev0->RP(q0x,q0y);
  66   ev1->RP(q1x,q1y);
  67   double rpphi = atan2(q0y+q1y,q0x+q1x);
  68
  69   // loop over pairs
  70
  71   for (int i=0; i<ev0->NParticles(); i++) {
  72     if (!ev0->ParticleGood(i,fPid0)) continue;
  73     for (int j=0; j<ev1->NParticles(); j++) {
  74       if (ev0==ev1 && j<i && fPid0==fPid1 ) continue;
  75       if (ev0==ev1 && j==i) continue; // beware, not even when rotated or non-identical
  76       if (!ev1->ParticleGood(j,fPid1)) continue;
  77       if (!ev0->PairGood(ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i),
  78                          ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot)) continue;
  79       fPa.Set0(fMass0,ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i));
  80       fPa.Set1(fMass1,ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot);
  81       if (ev0==ev1 && fPid0==fPid1 && ran.Rndm()>=0.5) fPa.Swap();
  82       fPa.CalcLAB();
  83       fPa.CalcPairCM();
  84       if (fPa.QCM()==0) return; // should not be too frequent
  85       double phi = TVector2::Phi_mpi_pi(fPa.Phi()-rpphi);
  86       pair->Fill(tmr,                    // 0 for tru, 1 for mix, 2 for rot
  87                  cent,                   // centrality
  88                  fPa.Rapidity(),         // pair rapidity
  89                  phi,                    // pair phi wrt reaction plane
  90                  fPa.Pt()/2.0,           // half of pair pt
  91                  fPa.QCMTheta(),         // polar angle of Q
  92                  fPa.QCMPhiOut(),        // azimuthal angle of Q w.r.t. out
  93                  fPa.QCM(),              // |p2-p1| in c.m.s.
  94                  1.0);                   // weigth
  95     }
  96   }
  97 }
  98 //=============================================================================