PWG2/UNICOR/AliUnicorAnalCorrel.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2048, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 // Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2005
  17
  18 //=============================================================================
  19 // two-particle correlation analyzer
  20 // Loop over pairs and fill pair histograms. The first particle is always
  21 // taken from ev0 and the second from ev1 so for ev0=ev1 one is getting true
  22 // pairs, otherwise mixed ones.
  23 //=============================================================================
  24
  25 #include <TROOT.h>
  26 #include <TMath.h>
  27 #include <TRandom2.h>
  28 #include "AliUnicorEvent.h"
  29 #include "AliUnicorHN.h"
  30 #include "AliUnicorAnalCorrel.h"
  31
  32 ClassImp(AliUnicorAnalCorrel)
  33
  34 //=============================================================================
  35 AliUnicorAnalCorrel::AliUnicorAnalCorrel(const char *nam, Double_t emi, Double_t ema,
  36                          Int_t pid0, Int_t pid1, AnalysisFrame frame):
  37   AliUnicorAnal(nam), fPid0(pid0), fPid1(pid1), fMass0(0), fMass1(0), fZ0(0), fZ1(0),
  38   fFrame(frame), fPa() {
  39   // constructor
  40   // emi and ema define the rapidity range for histogram
  41
  42   TParticlePDG *part0 = AliUnicorAnal::fgPDG.GetParticle(fPid0);
  43   TParticlePDG *part1 = AliUnicorAnal::fgPDG.GetParticle(fPid1);
  44   fMass0 = part0? part0->Mass() : 0;
  45   fMass1 = part1? part1->Mass() : 0;
  46   fZ0 = part0? part0->Charge()/3.0 : 0;
  47   fZ1 = part1? part1->Charge()/3.0 : 0;
  48   double pi = TMath::Pi();
  49
  50   // correlation function
  51
  52   int nce = 6; double cebins[]={0,0.05,0.1,0.2,0.4,0.6,1.0};  // centrality bins
  53
  54   //int npt = 7; double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.7,1.0};
  55   //int npt = 6; double ptbins[]={0,0.1,0.25,0.35,0.55,1.0,2.0}; // like Adam, except last bin split
  56   //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
  57   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
  58
  59   double qbins[200] = {0};
  60   //  for (int i=0;i<60;i++) qbins[i]=i*0.005;
  61   //  for (int i=0;i<45;i++) qbins[20+i]=0.1+i*0.02;
  62   //  for (int i=0;i<30;i++) qbins[i]=i*0.010;
  63   for (int i=0;i<60;i++) qbins[i]=i*0.005;
  64   for (int i=0;i<35;i++) qbins[60+i]=0.3+i*0.02;
  65   for (int i=0;i<20;i++) qbins[95+i]=1.0+i*0.05;
  66   for (int i=0;i<11;i++) qbins[115+i]=2.0+i*0.20;
  67
  68   TAxis *ax[8];
  69   ax[0] = new TAxis(4,-0.5,3.5);ax[0]->SetTitle("trumixrot");
  70   //  ax[1] = new TAxis(5,0,1.0);   ax[1]->SetTitle("centrality");
  71   ax[1] = new TAxis(nce,cebins);ax[1]->SetTitle("centrality");
  72   ax[2] = new TAxis(3,emi,ema); ax[2]->SetTitle("y");          // pair y
  73   //  ax[3] = new TAxis(8,-pi,pi);  ax[3]->SetTitle("phi");      // wrt event plane
  74   ax[3] = new TAxis(1,-pi,pi);  ax[3]->SetTitle("phi");        // wrt event plane
  75   ax[4] = new TAxis(npt,ptbins);ax[4]->SetTitle("kt (GeV/c)"); // pair pt/2
  76   ax[5] = new TAxis(8,0,pi);    ax[5]->SetTitle("q-theta");
  77   ax[6] = new TAxis(16,-pi,pi); ax[6]->SetTitle("q-phi");
  78   ax[7] = new TAxis(125,qbins); ax[7]->SetTitle("q (GeV/c)");
  79   //  ax[7] = new TAxis(700,0,3.5); ax[7]->SetTitle("q (GeV/c)");
  80   AliUnicorHN *pair = new AliUnicorHN("pair",8,ax);
  81   fHistos.Add(pair);
  82
  83   // correlation function bin monitor (needed to get <kt> etc.)
  84
  85   TAxis *bx[3]={0};
  86   //bx[0] = new TAxis(*(pair->GetAxis(1))); // wait until root bug (516-00..527-06) fixed. For now, do:
  87   bx[0] = new TAxis(ax[1]->GetNbins(), ax[1]->GetXmin(), ax[1]->GetXmax()); bx[0]->SetTitle("centrality");
  88   bx[1] = new TAxis(10*ax[2]->GetNbins(),emi,ema); bx[1]->SetTitle("y");         // pair y
  89   bx[2] = new TAxis(100,0,2); bx[2]->SetTitle("kt (GeV/c)");                     // pair pt/2
  90   AliUnicorHN *bimo = new AliUnicorHN("bimo",3,bx);
  91   for (int i=0; i<8; i++) delete ax[i];
  92   for (int i=0; i<3; i++) delete bx[i];
  93   fHistos.Add(bimo);
  94
  95   // two-track resolution monitoring histogram
  96
  97   ax[0] = new TAxis(3,-0.5,2.5);    ax[0]->SetTitle("trumixrot");
  98   ax[1] = new TAxis(2,-0.5,1.5);    ax[1]->SetTitle("cut applied");
  99   ax[2] = new TAxis(npt,ptbins);    ax[2]->SetTitle("(pair pt)/2 (GeV)");
 100   ax[3] = new TAxis(80,-0.08,0.08); ax[3]->SetTitle("dtheta");
 101   ax[4] = new TAxis(80,-0.20,0.20); ax[4]->SetTitle("dphi");
 102   AliUnicorHN *twot = new AliUnicorHN("twot",5,ax);
 103   for (int i=0; i<5; i++) delete ax[i];
 104   fHistos.Add(twot);
 105
 106   gROOT->cd();
 107   printf("%s object named %s created\n",ClassName(),GetName());
 108 }
 109 //=============================================================================
 110 void AliUnicorAnalCorrel::Process(Int_t tmr, const AliUnicorEvent * const ev0, const AliUnicorEvent * const ev1, Double_t phirot)
 111 {
 112   // process pairs from one or two (if mixing) events
 113   // tmr tells which histogram (bins) to fill: tru,mix,rot
 114
 115   // Could be possibly accelerated by checking the "good particle" only once
 116   // and caching the result. (Maybe the optimizer does it already.)
 117
 118   static TRandom2 ran;
 119   AliUnicorHN *pair = (AliUnicorHN*) fHistos.At(0);
 120   AliUnicorHN *bimo = (AliUnicorHN*) fHistos.At(1);
 121   AliUnicorHN *twot = (AliUnicorHN*) fHistos.At(2);
 122
 123   // mixing-and-rotating-proof centrality and reaction plane angle
 124   // (but not rotation-proof for rotation angles much different from 0 and 180)
 125   // true and rotated pairs are within the triangle (j<i), mixed - all
 126   // thus, proper rotation is either by 180, or by 170 AND 190, etc.
 127
 128   double cent = (ev0->Centrality()+ev1->Centrality())/2.0;
 129   double q0x=0,q0y=0,q1x=0,q1y=0;
 130   ev0->RP(q0x,q0y);
 131   ev1->RP(q1x,q1y);
 132   double rpphi = atan2(q0y+q1y,q0x+q1x);
 133
 134   // loop over pairs
 135
 136   for (int i=0; i<ev0->NParticles(); i++) {
 137     if (!ev0->ParticleGood(i,fPid0)) continue;
 138     for (int j=0; j<ev1->NParticles(); j++) {
 139       if (ev0==ev1 && j<i && fPid0==fPid1 ) continue;
 140       if (ev0==ev1 && j==i) continue; // beware, not even when rotated or non-identical
 141       if (!ev1->ParticleGood(j,fPid1)) continue;
 142       fPa.Set0(fMass0,ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i));
 143       fPa.Set1(fMass1,ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot);
 144       if (ev0==ev1 && fPid0==fPid1 && ran.Rndm()>=0.5) fPa.Swap();
 145       twot->Fill((double) tmr, 0.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
 146       if (!ev0->PairGood(ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i),fZ0,
 147                          ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot,fZ1)) continue;
 148       twot->Fill((double) tmr, 1.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
 149       fPa.CalcLAB(); // this could be organized better. AliUnicorPair should give k*?
 150       fPa.CalcPairCM();
 151       double qcm = fPa.QCM(); // momdif in pair cm - argument for Coulomb correction
 152
 153       fPa.CalcLAB();
 154       if (fFrame == kPairFrame) fPa.CalcPairCM();
 155       if (fFrame == kLCMS)      fPa.CalcLcmsCM();
 156       if (fPa.QCM()==0) {printf("AliUnicorAnalCorrel: Q=0\n"); return;} // should not be too frequent
 157       double phi = TVector2::Phi_mpi_pi(fPa.Phi()-rpphi);
 158       double weigth = 1.0;
 159       /*
 160       static TH2D *coul = 0;
 161       if (!coul) {
 162         TFile::Open("coulomb.root","read");
 163         coul = (TH2D*) gFile->Get("co");
 164         coul->SetDirectory(gROOT);
 165         gFile->Close();
 166       }
 167       if (tmr==0 && fPid0==fPid1) {
 168         double co = 0;
 169         if (qcm>0.999) co = 1;
 170         else if (qcm>0.001) co = coul->Interpolate(7,qcm);
 171         weigth = 1.0-0.5+0.5*co*(1+exp(-pow(fPa.QCM()*7/0.197,2)));
 172       }
 173       */
 174       pair->Fill((double) tmr,           // 0 for tru, 1 for mix, 2 for rot
 175                  cent,                   // centrality
 176                  fPa.Rapidity(),         // pair rapidity
 177                  phi,                    // pair phi wrt reaction plane
 178                  fPa.Pt()/2.0,           // half of pair pt
 179                  fPa.QCMTheta(),         // polar angle of Q
 180                  fPa.QCMPhiOut(),        // azimuthal angle of Q w.r.t. out
 181                  fPa.QCM(),              // |p2-p1| in c.m.s.
 182                  weigth);                // weigth
 183       if (tmr==0 && fPa.QCM()<0.2) bimo->Fill(cent, fPa.Rapidity(), fPa.Pt()/2.0, weigth);
 184       if (tmr==0) pair->Fill((double) 3,             // this is for Coulomb correction, maybe not necessary
 185                              cent,                   // centrality
 186                              fPa.Rapidity(),         // pair rapidity
 187                              phi,                    // pair phi wrt reaction plane
 188                              fPa.Pt()/2.0,           // half of pair pt
 189                              fPa.QCMTheta(),         // polar angle of Q
 190                              fPa.QCMPhiOut(),        // azimuthal angle of Q w.r.t. out
 191                              fPa.QCM(),              // |p2-p1| in c.m.s.
 192                              weigth*qcm);            // weigth*Q
 193     }
 194   }
 195 }
 196 //=============================================================================