1 #include "AliSPDMCTrackDensity.h"
2 #include <AliMCEvent.h>
3 #include <AliTrackReference.h>
13 #include "AliGenHijingEventHeader.h"
17 //____________________________________________________________________
18 AliSPDMCTrackDensity::AliSPDMCTrackDensity()
20 fUseOnlyPrimary(false),
23 fMinZ(-15), // -14.1),
24 fMaxZ(+15), // +14.1),
30 // Default constructor
33 //____________________________________________________________________
34 AliSPDMCTrackDensity::AliSPDMCTrackDensity(const char*)
35 : TNamed("spdMCTrackDensity","spdMCTrackDensity"),
36 fUseOnlyPrimary(false),
46 // Normal constructor constructor
49 //____________________________________________________________________
50 AliSPDMCTrackDensity::AliSPDMCTrackDensity(const AliSPDMCTrackDensity& o)
52 fUseOnlyPrimary(o.fUseOnlyPrimary),
62 // Normal constructor constructor
65 //____________________________________________________________________
67 AliSPDMCTrackDensity::operator=(const AliSPDMCTrackDensity& o)
69 // Assignment operator
70 if (&o == this) return *this;
72 fUseOnlyPrimary = o.fUseOnlyPrimary;
80 fBinFlow = o.fBinFlow;
84 //____________________________________________________________________
86 AliSPDMCTrackDensity::DefineOutput(TList* l)
88 fRZ = new TH2D("rz", "(r,z) of used track references",
89 30, -15, 15, 50, 0, 5);
90 fRZ->SetXTitle("z [cm]");
91 fRZ->SetYTitle("r [cm]");
95 fXYZ = new TH3D("xyz", "(x,y,z) of used track references",
96 100, -5., 5., 100, -5., 5., 30, -15., 15.);
97 fXYZ->SetXTitle("x [cm]");
98 fXYZ->SetYTitle("y [cm]");
99 fXYZ->SetZTitle("z [cm]");
100 fXYZ->SetDirectory(0);
103 fNRefs = new TH1D("nrefs", "Number of references used per track",
105 fNRefs->SetXTitle("# of references per track");
106 fNRefs->SetDirectory(0);
109 fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow",
110 120, -3, 3, 40, -180, 180);
111 fBinFlow->SetXTitle("#Delta#eta");
112 fBinFlow->SetYTitle("#Delta#varphi");
113 fBinFlow->SetDirectory(0);
117 //____________________________________________________________________
119 AliSPDMCTrackDensity::StoreParticle(AliMCParticle* particle,
120 const AliMCParticle* mother,
127 if (refNo < 0) return;
129 AliTrackReference* ref = particle->GetTrackReference(refNo);
132 Double_t r = ref->R();
133 Double_t x = ref->X();
134 Double_t y = ref->Y();
135 Double_t z = ref->Z();
138 Double_t th = TMath::ATan2(r,zr);
139 if (th < 0) th += 2*TMath::Pi();
140 Double_t et = -TMath::Log(TMath::Tan(th/2));
141 Double_t ph = TMath::ATan2(y,x);
142 if (ph < 0) ph += 2*TMath::Pi();
143 output.Fill(et,ph,w);
145 const AliMCParticle* mp = (mother ? mother : particle);
146 Double_t dEta = mp->Eta() - et;
147 Double_t dPhi = (mp->Phi() - ph) * 180 / TMath::Pi();
148 if (dPhi > 180) dPhi -= 360;
149 if (dPhi < -180) dPhi += 360;
150 fBinFlow->Fill(dEta, dPhi,w);
154 //____________________________________________________________________
156 AliSPDMCTrackDensity::GetMother(Int_t iTr,
157 const AliMCEvent& event) const
160 // Track down primary mother
164 const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
165 if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
173 // #define USE_FLOW_WEIGHTS
174 //____________________________________________________________________
176 AliSPDMCTrackDensity::Calculate(const AliMCEvent& event,
182 // Filter the input kinematics and track references, using
183 // some of the ESD information
186 // input Input ESD event
187 // event Input MC event
188 // vz Vertex position
189 // output Output ESD-like object
190 // primary Per-event histogram of primaries
193 // True on succes, false otherwise
196 #ifdef USE_FLOW_WEIGHTS
197 AliGenHijingEventHeader* hd = dynamic_cast<AliGenHijingEventHeader*>
198 (event.GenEventHeader());
199 Double_t rp = (hd ? hd->ReactionPlaneAngle() : 0.);
200 Double_t b = (hd ? hd->ImpactParameter() : -1 );
203 AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
204 Int_t nTracks = stack->GetNtrack();
205 Int_t nPrim = stack->GetNtrack();
206 for (Int_t iTr = 0; iTr < nTracks; iTr++) {
207 AliMCParticle* particle =
208 static_cast<AliMCParticle*>(event.GetTrack(iTr));
210 // Check the returned particle
211 if (!particle) continue;
213 // Check if this charged and a primary
214 Bool_t isCharged = particle->Charge() != 0;
215 if (!isCharged) continue;
217 Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
219 // Fill 'dn/deta' histogram
220 if (isPrimary && iTr < nPrim) {
222 Double_t eta = particle->Eta();
223 Double_t phi = particle->Phi();
224 primary->Fill(eta, phi);
228 // Bail out if we're only processing primaries - perhaps we should
229 // track back to the original primary?
230 if (fUseOnlyPrimary && !isPrimary) continue;
232 Int_t nTrRef = particle->GetNumberOfTrackReferences();
234 for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
235 AliTrackReference* ref = particle->GetTrackReference(iTrRef);
240 // Check that we hit an FMD element
241 if (ref->DetectorId() != AliTrackReference::kITS)
244 // Get radius and z where the track reference was made
245 Double_t r = ref->R();
246 Double_t x = ref->X();
247 Double_t y = ref->Y();
248 Double_t z = ref->Z();
249 if (r > fMaxR || r < fMinR) continue;
250 if (z > fMaxZ || z < fMinZ) continue;
255 // Only fill first reference
257 const AliMCParticle* mother = GetMother(iTr, event);
258 #ifdef USE_FLOW_WEIGHTS
259 Double_t phi = (mother ? mother->Phi() : particle->Phi());
260 Double_t eta = (mother ? mother->Eta() : particle->Eta());
261 Double_t pt = (mother ? mother->Pt() : particle->Pt());
262 Int_t id = (mother ? mother->PdgCode() : 2212);
263 Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id);
265 Double_t weight = 1.;
267 StoreParticle(particle, mother, iTrRef, vz, output, weight);
275 //____________________________________________________________________
277 AliSPDMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt, Double_t b,
278 Double_t phi, Double_t rp, Int_t id) const
280 static TF1 gaus = TF1("gaus", "gaus", -6, 6);
281 gaus.SetParameters(0.1, 0., 9);
282 // gaus.SetParameters(0.1, 0., 3);
283 // gaus.SetParameters(0.1, 0., 15);
285 const Double_t xCumulant2nd4050ALICE[] = {0.00, 0.25, 0.350,
292 const Double_t yCumulant2nd4050ALICE[] = {0.00000, 0.043400,
302 const Int_t nPointsCumulant2nd4050ALICE =
303 sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);
304 static TGraph alicePointsPt2(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
306 const Double_t xCumulant4th3040ALICE[] = {0.00,0.250,0.35,
314 const Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,
326 const Double_t xCumulant4th4050ALICE[] = {0.00,0.25,0.350,
333 const Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,
343 const Int_t nPointsCumulant4th4050ALICE =
344 sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);
345 static TGraph alicePointsPt4(nPointsCumulant4th4050ALICE,
346 xCumulant4th4050ALICE,
347 yCumulant4th4050ALICE);
349 const Double_t xCumulant4thTPCrefMultTPConlyAll[] = {1.75,
357 const Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,
361 const Int_t nPointsCumulant4thTPCrefMultTPConlyAll =
362 sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
363 TGraph aliceCent(nPointsCumulant4thTPCrefMultTPConlyAll,
364 xCumulant4thTPCrefMultTPConlyAll,
365 yCumulant4thTPCrefMultTPConlyAll);
368 Double_t weight = (20. * gaus.Eval(eta) * (alicePointsPt2.Eval(pt) * 0.5 +
369 alicePointsPt4.Eval(pt) * 0.5)
370 * (aliceCent.Eval(b) / aliceCent.Eval(10.46))
371 * 2. * TMath::Cos(2. * (phi - rp)));
372 if (TMath::Abs(id) == 211) weight *= 1.3; //pion flow
373 else if (TMath::Abs(id) == 2212) weight *= 1.0; //proton flow
379 //____________________________________________________________________
381 AliSPDMCTrackDensity::Print(Option_t* /*option*/) const
383 char ind[gROOT->GetDirLevel()+1];
384 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
385 ind[gROOT->GetDirLevel()] = '\0';
386 std::cout << ind << ClassName() << ": " << GetName() << '\n'
388 << ind << " Only primary tracks: " << fUseOnlyPrimary << '\n'
389 << ind << " R range: [" << fMinR << ',' << fMaxR
391 << ind << " Z range: [" << fMinZ << ',' << fMaxZ
392 << "]" << std::noboolalpha << std::endl;
396 //____________________________________________________________________