1 /*************************************************************************
2 * Copyright(c) 1998-2048, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2005
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 //=============================================================================
28 #include "AliUnicorEvent.h"
29 #include "AliUnicorHN.h"
30 #include "AliUnicorAnalCorrel.h"
32 ClassImp(AliUnicorAnalCorrel)
34 //=============================================================================
35 AliUnicorAnalCorrel::AliUnicorAnalCorrel(Char_t *nam, Double_t emi, Double_t ema,
36 Int_t pid0, Int_t pid1):
37 AliUnicorAnal(nam), fPid0(pid0), fPid1(pid1), fMass0(0), fMass1(0), fPa()
40 // emi and ema define the rapidity range for histogram
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 double pi = TMath::Pi();
48 // correlation function
50 double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.7,1.0};
52 for (int i=0;i<20;i++) qbins[i]=i*0.005;
53 for (int i=0;i<45;i++) qbins[20+i]=0.1+i*0.02;
56 ax[0] = new TAxis(3,-0.5,2.5);ax[0]->SetTitle("trumixrot");
57 ax[1] = new TAxis(5,0,1.0); ax[1]->SetTitle("centrality");
58 ax[2] = new TAxis(3,emi,ema); ax[2]->SetTitle("pair y");
59 //ax[3] = new TAxis(8,-pi,pi); ax[3]->SetTitle("pair phi"); // wrt event plane
60 ax[3] = new TAxis(1,-pi,pi); ax[3]->SetTitle("pair phi"); // wrt event plane
61 ax[4] = new TAxis(7,ptbins); ax[4]->SetTitle("(pair pt)/2 (GeV)");
62 ax[5] = new TAxis(8,0,pi); ax[5]->SetTitle("q-theta");
63 ax[6] = new TAxis(16,-pi,pi); ax[6]->SetTitle("q-phi");
64 //ax[7] = new TAxis(64,qbins); ax[7]->SetTitle("q (GeV/c)");
65 ax[7] = new TAxis(100,0,2.0); ax[7]->SetTitle("q (GeV/c)");
66 AliUnicorHN *pair = new AliUnicorHN("pair",8,ax);
67 for (int i=0; i<8; i++) delete ax[i];
70 // two-track resolution monitoring histogram
72 ax[0] = new TAxis(3,-0.5,2.5); ax[0]->SetTitle("trumixrot");
73 ax[1] = new TAxis(2,-0.5,1.5); ax[1]->SetTitle("cut applied");
74 ax[2] = new TAxis(7,ptbins); ax[2]->SetTitle("(pair pt)/2 (GeV)");
75 ax[3] = new TAxis(80,-0.02,0.02); ax[3]->SetTitle("dtheta");
76 ax[4] = new TAxis(80,-0.04,0.04); ax[4]->SetTitle("dphi");
77 AliUnicorHN *twot = new AliUnicorHN("twot",5,ax);
78 for (int i=0; i<5; i++) delete ax[i];
82 printf("%s object named %s created\n",ClassName(),GetName());
84 //=============================================================================
85 void AliUnicorAnalCorrel::Process(Int_t tmr, const AliUnicorEvent * const ev0, const AliUnicorEvent * const ev1, Double_t phirot)
87 // process pairs from one or two (if mixing) events
88 // tmr tells which histogram (bins) to fill: tru,mix,rot
90 // Could be possibly accelerated by checking the "good particle" only once
91 // and caching the result. (Maybe the optimizer does it already.)
94 AliUnicorHN *pair = (AliUnicorHN*) fHistos.At(0);
95 AliUnicorHN *twot = (AliUnicorHN*) fHistos.At(1);
97 // mixing-and-rotating-proof centrality and reaction plane angle
98 // (but not rotation-proof for rotation angles much different from 0 and 180)
99 // true and rotated pairs are within the triangle (j<i), mixed - all
100 // thus, proper rotation is either by 180, or by 170 AND 190, etc.
102 double cent = (ev0->Centrality()+ev1->Centrality())/2.0;
103 double q0x,q0y,q1x,q1y;
106 double rpphi = atan2(q0y+q1y,q0x+q1x);
110 for (int i=0; i<ev0->NParticles(); i++) {
111 if (!ev0->ParticleGood(i,fPid0)) continue;
112 for (int j=0; j<ev1->NParticles(); j++) {
113 if (ev0==ev1 && j<i && fPid0==fPid1 ) continue;
114 if (ev0==ev1 && j==i) continue; // beware, not even when rotated or non-identical
115 if (!ev1->ParticleGood(j,fPid1)) continue;
116 fPa.Set0(fMass0,ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i));
117 fPa.Set1(fMass1,ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot);
118 if (ev0==ev1 && fPid0==fPid1 && ran.Rndm()>=0.5) fPa.Swap();
119 twot->Fill((double) tmr, 0.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
120 if (!ev0->PairGood(ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i),
121 ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot)) continue;
122 twot->Fill((double) tmr, 1.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
125 if (fPa.QCM()==0) return; // should not be too frequent
126 double phi = TVector2::Phi_mpi_pi(fPa.Phi()-rpphi);
128 //if (tmr==0) weigth = 1.0+0.3*exp(-fPa.QCM()*fPa.QCM()*1*1/0.197/0.197);
129 pair->Fill((double) tmr, // 0 for tru, 1 for mix, 2 for rot
131 fPa.Rapidity(), // pair rapidity
132 phi, // pair phi wrt reaction plane
133 fPa.Pt()/2.0, // half of pair pt
134 fPa.QCMTheta(), // polar angle of Q
135 fPa.QCMPhiOut(), // azimuthal angle of Q w.r.t. out
136 fPa.QCM(), // |p2-p1| in c.m.s.
141 //=============================================================================