Remove Warnings.
[u/mrichter/AliRoot.git] / FASTSIM / AliQuenchingWeights.cxx
CommitLineData
b6d061b7 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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/* $Id$ */
17
18//----------------------------------------------------------------------------
19// Implementation of the class to calculate the parton energy loss
20// Based on the "BDMPS" quenching weights by C.A.Salgado and U.A.Wiedemann
21//
22// References:
23// C.A.Salgado and U.A.Wiedemann, Phys.Rev.D68 (2003) 014008 [hep-ph/0302184]
24// A.Dainese, Eur.Phys.J.C, in press, [nucl-ex/0312005]
25//
7a76a12e 26//
b6d061b7 27// Origin: C. Loizides constantinos.loizides@cern.ch
28// A. Dainese andrea.dainese@pd.infn.it
7a76a12e 29//
b6d061b7 30//----------------------------------------------------------------------------
31
32#include <Riostream.h>
33#include <TH1F.h>
34#include <TH2F.h>
35#include <TCanvas.h>
36#include <TGraph.h>
37#include <TROOT.h>
38#include <TSystem.h>
39#include <TLegend.h>
40#include "AliQuenchingWeights.h"
41
42ClassImp(AliQuenchingWeights)
43
44// conversion from fm to GeV^-1: 1 fm = fmGeV GeV^-1
7a76a12e 45const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197;
b6d061b7 46
47// maximum value of R
7a76a12e 48const Double_t AliQuenchingWeights::fgkRMax = 1.e6;
b6d061b7 49
50// counter for histogram labels
7a76a12e 51Int_t AliQuenchingWeights::fgCounter = 0;
b6d061b7 52
53AliQuenchingWeights::AliQuenchingWeights()
54 : TObject()
55{
7a76a12e 56 //default constructor
57
b6d061b7 58 fTablesLoaded=kFALSE;
59 fMultSoft=kTRUE;
60 fHistos=0;
61 SetMu();
62 SetQTransport();
63 SetECMethod();
64 SetLengthMax();
65 fLengthMaxOld=0;
7a76a12e 66 fInstanceNumber=fgCounter++;
b6d061b7 67}
68
69AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a)
70 : TObject()
71{
7a76a12e 72 //copy constructor
73
b6d061b7 74 fTablesLoaded=kFALSE;
75 fHistos=0;
76 fLengthMaxOld=0;
77 fMultSoft=a.GetMultSoft();;
78 fMu=a.GetMu();
79 fQTransport=a.GetQTransport();
80 fECMethod=(kECMethod)a.GetECMethod();
81 fLengthMax=a.GetLengthMax();
7a76a12e 82 fInstanceNumber=fgCounter++;
b6d061b7 83 //Missing in the class is the pathname
84 //to the tables, can be added if needed
85}
86
87AliQuenchingWeights::~AliQuenchingWeights()
88{
89 Reset();
90}
91
92void AliQuenchingWeights::Reset()
93{
7a76a12e 94 //reset tables if there were used
95
b6d061b7 96 if(!fHistos) return;
97 for(Int_t l=0;l<fLengthMaxOld;l++){
98 delete fHistos[0][l];
99 delete fHistos[1][l];
100 }
101 delete[] fHistos;
102 fHistos=0;
103 fLengthMaxOld=0;
104}
105
106void AliQuenchingWeights::SetECMethod(kECMethod type)
107{
7a76a12e 108 //set energy constraint method
109
b6d061b7 110 fECMethod=type;
111 if(fECMethod==kDefault)
112 Info("SetECMethod","Energy Constraint Method set to DEFAULT:\nIf (sampled energy loss > parton energy) then sampled energy loss = parton energy.");
113 else
114 Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy.");
115}
116
117Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
118{
119 // read in tables for multiple scattering approximation
120 // path to continuum and to discrete part
121
122 fTablesLoaded = kFALSE;
123 fMultSoft=kTRUE;
124
125 Char_t fname[1024];
126 sprintf(fname,"%s",gSystem->ExpandPathName(contall));
bb545331 127 //PH ifstream fincont(fname);
128 fstream fincont(fname,ios::in);
129#if defined(__HP_aCC) || defined(__DECCXX)
130 if(!fincont.rdbuf()->is_open()) return -1;
131#else
b6d061b7 132 if(!fincont.is_open()) return -1;
bb545331 133#endif
b6d061b7 134
135 Int_t nn=0; //quarks
136 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
137 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
7a76a12e 138 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>fcaq[13][nn]>>
139 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>fcaq[18][nn]>>
140 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>fcaq[23][nn]>>
141 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>fcaq[28][nn]>>
142 fcaq[29][nn]>>fcaq[30][nn]>>fcaq[31][nn]>>fcaq[32][nn]>>fcaq[33][nn])
b6d061b7 143 {
144 nn++;
145 if(nn==261) break;
146 }
147
148 nn=0; //gluons
149 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
150 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
7a76a12e 151 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>fcag[13][nn]>>
152 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>fcag[18][nn]>>
153 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>fcag[23][nn]>>
154 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>fcag[28][nn]>>
155 fcag[29][nn]>>fcag[30][nn]>>fcag[31][nn]>>fcag[32][nn]>>fcag[33][nn])
156 {
b6d061b7 157 nn++;
158 if(nn==261) break;
159 }
160 fincont.close();
161
162 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
bb545331 163 //PH ifstream findisc(fname);
164 fstream findisc(fname,ios::in);
165#if defined(__HP_aCC) || defined(__DECCXX)
166 if(!findisc.rdbuf()->is_open()) return -1;
167#else
b6d061b7 168 if(!findisc.is_open()) return -1;
bb545331 169#endif
b6d061b7 170
171 nn=0; //quarks
172 while(findisc>>frrr[nn]>>fdaq[nn]) {
173 nn++;
7a76a12e 174 if(nn==34) break;
b6d061b7 175 }
176 nn=0; //gluons
177 while(findisc>>frrrg[nn]>>fdag[nn]) {
178 nn++;
7a76a12e 179 if(nn==34) break;
b6d061b7 180 }
181 findisc.close();
182 fTablesLoaded = kTRUE;
183 return 0;
184}
185
186/*
187C***************************************************************************
188C Quenching Weights for Multiple Soft Scattering
189C February 10, 2003
190C
191C Refs:
192C
193C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
194C
195C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
196C
197C
198C This package contains quenching weights for gluon radiation in the
199C multiple soft scattering approximation.
200C
201C swqmult returns the quenching weight for a quark (ipart=1) or
202C a gluon (ipart=2) traversing a medium with transport coeficient q and
203C length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where
204C wc=0.5*q*L^2 and w is the energy radiated. The output values are
205C the continuous and discrete (prefactor of the delta function) parts
206C of the quenching weights.
207C
208C In order to use this routine, the files cont.all and disc.all need to be
209C in the working directory.
210C
211C An initialization of the tables is needed by doing call initmult before
212C using swqmult.
213C
214C Please, send us any comment:
215C
7a76a12e 216C urs.wiedemann@cern.ch
217C carlos.salgado@cern.ch
218C
219C
b6d061b7 220C-------------------------------------------------------------------
221
222 SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete)
223*
7a76a12e 224 REAL*8 xx(400), daq(34), caq(34,261), rrr(34)
b6d061b7 225 COMMON /dataqua/ xx, daq, caq, rrr
226*
7a76a12e 227 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
b6d061b7 228 COMMON /dataglu/ xxg, dag, cag, rrrg
229
230 REAL*8 rrrr,xxxx, continuous, discrete
231 REAL*8 rrin, xxin
232 INTEGER nrlow, nrhigh, nxlow, nxhigh
233 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
234 REAL*8 xfraclow, xfrachigh
235 REAL*8 clow, chigh
236*
7a76a12e 237
238 continuous=0.d0
239 discrete=0.d0
240
b6d061b7 241 rrin = rrrr
242 xxin = xxxx
243*
7a76a12e 244 do 666, nr=1,34
b6d061b7 245 if (rrin.lt.rrr(nr)) then
246 rrhigh = rrr(nr)
247 else
248 rrhigh = rrr(nr-1)
249 rrlow = rrr(nr)
250 nrlow = nr
251 nrhigh = nr-1
252 goto 665
253 endif
254 666 enddo
255 665 continue
256*
257 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
258 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
7a76a12e 259 if (rrin.gt.10000d0) then
260 rfraclow = dlog(rrhigh/rrin)/dlog(rrhigh/rrlow)
261 rfrachigh = dlog(rrin/rrlow)/dlog(rrhigh/rrlow)
262 endif
263*
264 if (ipart.eq.1.and.rrin.ge.rrr(1)) then
265 nrlow=1
266 nrhigh=1
267 rfraclow=1
268 rfrachigh=0
269 endif
270
271 if (ipart.ne.1.and.rrin.ge.rrrg(1)) then
272 nrlow=1
273 nrhigh=1
274 rfraclow=1
275 rfrachigh=0
276 endif
277
278 if (xxxx.ge.xx(261)) go to 245
279
280 nxlow = int(xxin/0.01) + 1
281 nxhigh = nxlow + 1
282 xfraclow = (xx(nxhigh)-xxin)/0.01
283 xfrachigh = (xxin - xx(nxlow))/0.01
b6d061b7 284*
285 if(ipart.eq.1) then
286 clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh)
287 chigh = xfraclow*caq(nrhigh,nxlow)+xfrachigh*caq(nrhigh,nxhigh)
288 else
289 clow = xfraclow*cag(nrlow,nxlow)+xfrachigh*cag(nrlow,nxhigh)
290 chigh = xfraclow*cag(nrhigh,nxlow)+xfrachigh*cag(nrhigh,nxhigh)
291 endif
292
293 continuous = rfraclow*clow + rfrachigh*chigh
294
7a76a12e 295245 continue
296
b6d061b7 297 if(ipart.eq.1) then
298 discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh)
299 else
300 discrete = rfraclow*dag(nrlow) + rfrachigh*dag(nrhigh)
301 endif
302*
303 END
304
305 subroutine initmult
7a76a12e 306 REAL*8 xxq(400), daq(34), caq(34,261), rrr(34)
b6d061b7 307 COMMON /dataqua/ xxq, daq, caq, rrr
308*
7a76a12e 309 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
b6d061b7 310 COMMON /dataglu/ xxg, dag, cag, rrrg
311*
7a76a12e 312 OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90)
b6d061b7 313 do 110 nn=1,261
314 read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn),
315 + caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn),
316 + caq(9,nn), caq(10,nn), caq(11,nn), caq(12,nn),
317 + caq(13,nn),
318 + caq(14,nn), caq(15,nn), caq(16,nn), caq(17,nn),
319 + caq(18,nn),
320 + caq(19,nn), caq(20,nn), caq(21,nn), caq(22,nn),
321 + caq(23,nn),
322 + caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn),
323 + caq(28,nn),
7a76a12e 324 + caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn),
325 + caq(33,nn), caq(34,nn)
b6d061b7 326 110 continue
327 do 111 nn=1,261
328 read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn),
329 + cag(4,nn), cag(5,nn), cag(6,nn), cag(7,nn), cag(8,nn),
330 + cag(9,nn), cag(10,nn), cag(11,nn), cag(12,nn),
331 + cag(13,nn),
332 + cag(14,nn), cag(15,nn), cag(16,nn), cag(17,nn),
333 + cag(18,nn),
334 + cag(19,nn), cag(20,nn), cag(21,nn), cag(22,nn),
335 + cag(23,nn),
336 + cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn),
337 + cag(28,nn),
7a76a12e 338 + cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn),
339 + cag(33,nn), cag(34,nn)
b6d061b7 340 111 continue
341 close(20)
342*
7a76a12e 343 OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91)
344 do 112 nn=1,34
b6d061b7 345 read (21,*) rrr(nn), daq(nn)
346 112 continue
7a76a12e 347 do 113 nn=1,34
b6d061b7 348 read (21,*) rrrg(nn), dag(nn)
349 113 continue
350 close(21)
351*
352 goto 888
353 90 PRINT*, 'input - output error'
354 91 PRINT*, 'input - output error #2'
355 888 continue
356
357 end
358
7a76a12e 359
b6d061b7 360=======================================================================
361
362 Adapted to ROOT macro by A. Dainese - 13/07/2003
7a76a12e 363 Ported to class by C. Loizides - 12/02/2004
364 New version for extended R values added - 06/03/2004
b6d061b7 365*/
366
367Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
368 Double_t &continuous,Double_t &discrete) const
369{
370 // Calculate Multiple Scattering approx.
371 // weights for given parton type,
372 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
373
7a76a12e 374 //set result to zero
375 continuous=0.;
376 discrete=0.;
377
b6d061b7 378 // read-in data before first call
379 if(!fTablesLoaded){
380 Error("CalcMult","Tables are not loaded.");
381 return -1;
382 }
383 if(!fMultSoft){
384 Error("CalcMult","Tables are not loaded for Multiple Scattering.");
385 return -1;
386 }
387
388 Double_t rrin = rrrr;
389 Double_t xxin = xxxx;
390
7a76a12e 391 if(xxin>fxx[260]) return -1;
b6d061b7 392 Int_t nxlow = (Int_t)(xxin/0.01) + 1;
393 Int_t nxhigh = nxlow + 1;
394 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.01;
395 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01;
396
397 //why this?
7a76a12e 398 if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD
b6d061b7 399 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
400
401 Int_t nrlow=0,nrhigh=0;
402 Double_t rrhigh=0,rrlow=0;
7a76a12e 403 for(Int_t nr=1; nr<=34; nr++) {
b6d061b7 404 if(rrin<frrr[nr-1]) {
405 rrhigh = frrr[nr-1];
406 } else {
407 rrhigh = frrr[nr-1-1];
408 rrlow = frrr[nr-1];
409 nrlow = nr;
410 nrhigh = nr-1;
411 break;
412 }
413 }
414
415 rrin = rrrr; // AD
416
417 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
418 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
419
7a76a12e 420 if(rrin>1.e4){
421 rfraclow = TMath::Log2(rrhigh/rrin)/TMath::Log2(rrhigh/rrlow);
422 rfrachigh = TMath::Log2(rrin/rrlow)/TMath::Log2(rrhigh/rrlow);
423 }
424 if((ipart==1) && (rrin>=frrr[0]))
425 {
426 nrlow=1;
427 nrhigh=1;
428 rfraclow=1.;
429 rfrachigh=0.;
430 }
431 if((ipart==2) && (rrin>=frrrg[0]))
432 {
433 nrlow=1;
434 nrhigh=1;
435 rfraclow=1.;
436 rfrachigh=0.;
437 }
438
b6d061b7 439 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
440
441 Double_t clow=0,chigh=0;
442 if(ipart==1) {
443 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
444 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
445 } else {
446 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
447 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
448 }
449
450 continuous = rfraclow*clow + rfrachigh*chigh;
451 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
452 // rfraclow,clow,rfrachigh,chigh,continuous);
453
454 if(ipart==1) {
455 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
456 } else {
457 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
458 }
459
460 return 0;
461}
462
463Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *discall)
464{
465 // read in tables for Single Hard Approx.
466 // path to continuum and to discrete part
467
468 fTablesLoaded = kFALSE;
469 fMultSoft=kFALSE;
470
471 Char_t fname[1024];
472 sprintf(fname,"%s",gSystem->ExpandPathName(contall));
bb545331 473 //PH ifstream fincont(fname);
474 fstream fincont(fname,ios::in);
475#if defined(__HP_aCC) || defined(__DECCXX)
476 if(!fincont.rdbuf()->is_open()) return -1;
477#else
b6d061b7 478 if(!fincont.is_open()) return -1;
bb545331 479#endif
b6d061b7 480
481 Int_t nn=0; //quarks
482 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
483 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
484 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>
485 fcaq[13][nn]>>
486 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>
487 fcaq[18][nn]>>
488 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>
489 fcaq[23][nn]>>
490 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>
491 fcaq[28][nn]>>
492 fcaq[29][nn])
493 {
494 nn++;
495 if(nn==261) break;
496 }
497
498 nn=0; //gluons
499 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
500 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
501 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>
502 fcag[13][nn]>>
503 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>
504 fcag[18][nn]>>
505 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>
506 fcag[23][nn]>>
507 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>
508 fcag[28][nn]>>
509 fcag[29][nn]) {
510 nn++;
511 if(nn==261) break;
512 }
513 fincont.close();
514
515 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
bb545331 516 //PH ifstream findisc(fname);
517 fstream findisc(fname,ios::in);
518#if defined(__HP_aCC) || defined(__DECCXX)
519 if(!findisc.rdbuf()->is_open()) return -1;
520#else
b6d061b7 521 if(!findisc.is_open()) return -1;
bb545331 522#endif
b6d061b7 523
524 nn=0; //quarks
525 while(findisc>>frrr[nn]>>fdaq[nn]) {
526 nn++;
527 if(nn==30) break;
528 }
529 nn=0; //gluons
530 while(findisc>>frrrg[nn]>>fdag[nn]) {
531 nn++;
532 if(nn==30) break;
533 }
534 findisc.close();
535
536 fTablesLoaded = kTRUE;
537 return 0;
538}
539
540/*
541C***************************************************************************
542C Quenching Weights for Single Hard Scattering
543C February 20, 2003
544C
545C Refs:
546C
547C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
548C
549C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
550C
551C
552C This package contains quenching weights for gluon radiation in the
553C single hard scattering approximation.
554C
555C swqlin returns the quenching weight for a quark (ipart=1) or
556C a gluon (ipart=2) traversing a medium with Debye screening mass mu and
557C length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where
558C wc=0.5*mu^2*L and w is the energy radiated. The output values are
559C the continuous and discrete (prefactor of the delta function) parts
560C of the quenching weights.
561C
562C In order to use this routine, the files contlin.all and disclin.all
563C need to be in the working directory.
564C
565C An initialization of the tables is needed by doing call initlin before
566C using swqlin.
567C
568C Please, send us any comment:
569C
570C urs.wiedemann@cern.ch
571C carlos.salgado@cern.ch
572C
573C
574C-------------------------------------------------------------------
575
576
577 SUBROUTINE swqlin(ipart,rrrr,xxxx,continuous,discrete)
578*
579 REAL*8 xx(400), dalq(30), calq(30,261), rrr(30)
580 COMMON /datalinqua/ xx, dalq, calq, rrr
581*
582 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
583 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
584
585 REAL*8 rrrr,xxxx, continuous, discrete
586 REAL*8 rrin, xxin
587 INTEGER nrlow, nrhigh, nxlow, nxhigh
588 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
589 REAL*8 xfraclow, xfrachigh
590 REAL*8 clow, chigh
591*
592 rrin = rrrr
593 xxin = xxxx
594*
595 nxlow = int(xxin/0.038) + 1
596 nxhigh = nxlow + 1
597 xfraclow = (xx(nxhigh)-xxin)/0.038
598 xfrachigh = (xxin - xx(nxlow))/0.038
599*
600 do 666, nr=1,30
601 if (rrin.lt.rrr(nr)) then
602 rrhigh = rrr(nr)
603 else
604 rrhigh = rrr(nr-1)
605 rrlow = rrr(nr)
606 nrlow = nr
607 nrhigh = nr-1
608 goto 665
609 endif
610 666 enddo
611 665 continue
612*
613 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
614 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
615*
616 if(ipart.eq.1) then
617 clow = xfraclow*calq(nrlow,nxlow)+xfrachigh*calq(nrlow,nxhigh)
618 chigh = xfraclow*calq(nrhigh,nxlow)+xfrachigh*calq(nrhigh,nxhigh)
619 else
620 clow = xfraclow*calg(nrlow,nxlow)+xfrachigh*calg(nrlow,nxhigh)
621 chigh = xfraclow*calg(nrhigh,nxlow)+xfrachigh*calg(nrhigh,nxhigh)
622 endif
623
624 continuous = rfraclow*clow + rfrachigh*chigh
625
626 if(ipart.eq.1) then
627 discrete = rfraclow*dalq(nrlow) + rfrachigh*dalq(nrhigh)
628 else
629 discrete = rfraclow*dalg(nrlow) + rfrachigh*dalg(nrhigh)
630 endif
631*
632 END
633
634 subroutine initlin
635 REAL*8 xxlq(400), dalq(30), calq(30,261), rrr(30)
636 COMMON /datalinqua/ xxlq, dalq, calq, rrr
637*
638 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
639 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
640*
641 OPEN(UNIT=20,FILE='contlin.all',STATUS='OLD',ERR=90)
642 do 110 nn=1,261
643 read (20,*) xxlq(nn), calq(1,nn), calq(2,nn), calq(3,nn),
644 + calq(4,nn), calq(5,nn), calq(6,nn), calq(7,nn), calq(8,nn),
645 + calq(9,nn), calq(10,nn), calq(11,nn), calq(12,nn),
646 + calq(13,nn),
647 + calq(14,nn), calq(15,nn), calq(16,nn), calq(17,nn),
648 + calq(18,nn),
649 + calq(19,nn), calq(20,nn), calq(21,nn), calq(22,nn),
650 + calq(23,nn),
651 + calq(24,nn), calq(25,nn), calq(26,nn), calq(27,nn),
652 + calq(28,nn),
653 + calq(29,nn), calq(30,nn)
654 110 continue
655 do 111 nn=1,261
656 read (20,*) xxlg(nn), calg(1,nn), calg(2,nn), calg(3,nn),
657 + calg(4,nn), calg(5,nn), calg(6,nn), calg(7,nn), calg(8,nn),
658 + calg(9,nn), calg(10,nn), calg(11,nn), calg(12,nn),
659 + calg(13,nn),
660 + calg(14,nn), calg(15,nn), calg(16,nn), calg(17,nn),
661 + calg(18,nn),
662 + calg(19,nn), calg(20,nn), calg(21,nn), calg(22,nn),
663 + calg(23,nn),
664 + calg(24,nn), calg(25,nn), calg(26,nn), calg(27,nn),
665 + calg(28,nn),
666 + calg(29,nn), calg(30,nn)
667 111 continue
668 close(20)
669*
670 OPEN(UNIT=21,FILE='disclin.all',STATUS='OLD',ERR=91)
671 do 112 nn=1,30
672 read (21,*) rrr(nn), dalq(nn)
673 112 continue
674 do 113 nn=1,30
675 read (21,*) rrrlg(nn), dalg(nn)
676 113 continue
677 close(21)
678*
679 goto 888
680 90 PRINT*, 'input - output error'
681 91 PRINT*, 'input - output error #2'
682 888 continue
683
684 end
685
686=======================================================================
687
688 Ported to class by C. Loizides - 17/02/2004
689
690*/
691
692Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx,
693 Double_t &continuous,Double_t &discrete) const
694{
695 // calculate Single Hard approx.
696 // weights for given parton type,
697 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
698
699 // read-in data before first call
700 if(!fTablesLoaded){
701 Error("CalcMult","Tables are not loaded.");
702 return -1;
703 }
704 if(!fMultSoft){
705 Error("CalcMult","Tables are not loaded for Single Hard Scattering.");
706 return -1;
707 }
708
709 Double_t rrin = rrrr;
710 Double_t xxin = xxxx;
711
712 Int_t nxlow = (Int_t)(xxin/0.038) + 1;
713 Int_t nxhigh = nxlow + 1;
714 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.038;
715 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.038;
716
717 //why this?
718 if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD
719 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
720
721 Int_t nrlow=0,nrhigh=0;
722 Double_t rrhigh=0,rrlow=0;
723 for(Int_t nr=1; nr<=30; nr++) {
724 if(rrin<frrr[nr-1]) {
725 rrhigh = frrr[nr-1];
726 } else {
727 rrhigh = frrr[nr-1-1];
728 rrlow = frrr[nr-1];
729 nrlow = nr;
730 nrhigh = nr-1;
731 break;
732 }
733 }
734
735 rrin = rrrr; // AD
736
737 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
738 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
739
740 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
741
742 Double_t clow=0,chigh=0;
743 if(ipart==1) {
744 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
745 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
746 } else {
747 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
748 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
749 }
750
751 continuous = rfraclow*clow + rfrachigh*chigh;
752 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
753 // rfraclow,clow,rfrachigh,chigh,continuous);
754
755 if(ipart==1) {
756 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
757 } else {
758 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
759 }
760
761 return 0;
762}
763
764Int_t AliQuenchingWeights::CalcMult(Int_t ipart,
765 Double_t w,Double_t qtransp,Double_t length,
766 Double_t &continuous,Double_t &discrete) const
767{
7a76a12e 768 // Calculate Multiple Scattering approx.
769 // weights for given parton type,
770 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
771
b6d061b7 772 Double_t wc=CalcWC(qtransp,length);
773 Double_t rrrr=CalcR(wc,length);
774 Double_t xxxx=w/wc;
775 return CalcMult(ipart,rrrr,xxxx,continuous,discrete);
776}
777
778Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
779 Double_t w,Double_t mu,Double_t length,
780 Double_t &continuous,Double_t &discrete) const
781{
7a76a12e 782 // calculate Single Hard approx.
783 // weights for given parton type,
784 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
785
b6d061b7 786 Double_t wcbar=CalcWCbar(mu,length);
787 Double_t rrrr=CalcR(wcbar,length);
788 Double_t xxxx=w/wcbar;
789 return CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
790}
791
792Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const
793{
7a76a12e 794 //calculate R value and
795 //check if it is less then maximum
796
797 Double_t R = wc*l*fgkConvFmToInvGeV;
798 if(R>fgkRMax) {
799 Warning("CalcR","Value of R = %.2f; should be less than %.2f",R,fgkRMax);
b6d061b7 800 return -R;
801 }
802 return R;
803}
804
805Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const
806{
807 // return DeltaE for MS or SH scattering
808 // for given parton type, length and energy
809 // Dependant on ECM (energy constraint method)
810 // e is used to determine where to set bins to zero.
811
812 if(!fHistos){
813 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
814 return -1000.;
815 }
816 if((ipart<1) || (ipart>2)) {
817 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
818 return -1000;
819 }
820
821 Int_t l=(Int_t)length;
822 if((length-(Double_t)l)>0.5) l++;
823 if(l<=0) return 0.;
824 if(l>fLengthMax) l=fLengthMax;
825
826 if(fECMethod==kReweight){
827 TH1F *dummy=new TH1F(*fHistos[ipart-1][l-1]);
828 dummy->SetName("dummy");
829 for(Int_t bin=dummy->FindBin(e)+1;bin<=1100;bin++){
830 dummy->SetBinContent(bin,0.);
831 }
832 dummy->Scale(1./dummy->Integral());
833 Double_t ret=dummy->GetRandom();
834 delete dummy;
835 return ret;
836 //****** !! ALTERNATIVE WAY OF DOING IT !! ***
837 //Double_t ret = 2.*e;
838 //while(ret>e) ret=fHistos[ipart-1][l-1]->GetRandom();
839 //return ret;
840 //********************************************
841 } else { //kDefault
842 Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
843 if(ret>e) return e;
844 return ret;
845 }
846}
847
848Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const
849{
850 //return quenched parton energy
851 //for given parton type, length and energy
852
853 Double_t loss=GetELossRandom(ipart,length,e);
854 return e-loss;
855}
856
e99e3ed5 857Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const
b6d061b7 858{
7a76a12e 859 // return DeltaE for MS or SH scattering
860 // for given parton type, length distribution and energy
861 // Dependant on ECM (energy constraint method)
862 // e is used to determine where to set bins to zero.
863
b6d061b7 864 if(!hell){
865 Warning("GetELossRandom","Pointer to length distribution is NULL.");
866 return 0.;
867 }
868 Double_t ell=hell->GetRandom();
869 return GetELossRandom(ipart,ell,e);
870}
871
872Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e) const
873{
874 //return quenched parton energy
875 //for given parton type, length distribution and energy
876
877 Double_t loss=GetELossRandom(ipart,hell,e);
878 return e-loss;
879}
880
881Int_t AliQuenchingWeights::SampleEnergyLoss()
882{
883 // Has to be called to fill the histograms
884 //
885 // For stored values fQTransport loop over
886 // particle type and length = 1 to fMaxLength (fm)
887 // to fill energy loss histos
888 //
889 // Take histogram of continuous weights
890 // Take discrete_weight
891 // If discrete_weight > 1, put all channels to 0, except channel 1
892 // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral
893
894 // read-in data before first call
895 if(!fTablesLoaded){
896 Error("CalcMult","Tables are not loaded.");
897 return -1;
898 }
899
900 if(fMultSoft) {
901 Int_t lmax=CalcLengthMax(fQTransport);
902 if(fLengthMax>lmax){
7a76a12e 903 Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
b6d061b7 904 fLengthMax=lmax;
905 }
906 } else {
907 Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested.");
908 }
909
910 Reset();
911 fHistos=new TH1F**[2];
912 fHistos[0]=new TH1F*[fLengthMax];
913 fHistos[1]=new TH1F*[fLengthMax];
914 fLengthMaxOld=fLengthMax; //remember old value in case
915 //user wants to reset
916
917 Int_t medvalue=0;
918 Char_t meddesc[100];
919 if(fMultSoft) {
920 medvalue=(Int_t)(fQTransport*1000.);
921 sprintf(meddesc,"MS");
922 } else {
923 medvalue=(Int_t)(fMu*1000.);
924 sprintf(meddesc,"SH");
925 }
926
927 for(Int_t ipart=1;ipart<=2;ipart++){
928 for(Int_t l=1;l<=fLengthMax;l++){
929
930 Char_t hname[100];
931 sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
932 Double_t wc = CalcWC(l);
933 fHistos[ipart-1][l-1] = new TH1F(hname,hname,1100,0.,1.1*wc);
934 fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
935 fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
936 fHistos[ipart-1][l-1]->SetLineColor(4);
937
938 Double_t rrrr = CalcR(wc,l);
939 Double_t discrete=0.;
940 // loop on histogram channels
941 for(Int_t bin=1; bin<=1100; bin++) {
942 Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
943 Double_t continuous;
944 CalcMult(ipart,rrrr,xxxx,continuous,discrete);
945 fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
946 }
947 // add discrete part to distribution
948 if(discrete>=1.)
949 for(Int_t bin=2;bin<=1100;bin++)
950 fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
951 else {
952 Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,1100);
953 fHistos[ipart-1][l-1]->Fill(0.,val);
954 }
955 Double_t hint=fHistos[ipart-1][l-1]->Integral(1,1100);
956 fHistos[ipart-1][l-1]->Scale(1./hint);
957 }
958 }
959 return 0;
960}
961
962const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Int_t l) const
963{
7a76a12e 964 //return quenching histograms
965 //for ipart and length
966
b6d061b7 967 if(!fHistos){
968 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
969 return 0;
970 }
971 if((ipart<1) || (ipart>2)) {
972 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
973 return 0;
974 }
975
976 if(l<=0) return 0;
977 if(l>fLengthMax) l=fLengthMax;
978
979 return fHistos[ipart-1][l-1];
980}
981
982TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const
983{
984 // ipart = 1 for quark, 2 for gluon
985 // medval a) qtransp = transport coefficient (GeV^2/fm)
986 // b) mu = Debye mass (GeV)
987 // length = path length in medium (fm)
988 // Get from SW tables:
989 // - discrete weight (probability to have NO energy loss)
990 // - continuous weight, as a function of dE
991 // compute up to max dE = 1.1*wc
992 // (as an histogram named hContQW_<ipart>_<medval*1000>_<l>
993 // and range [0,1.1*wc] (1100 channels)
994
995 Double_t wc = 0;
996 Char_t meddesc[100];
997 if(fMultSoft) {
998 wc=CalcWC(medval,length);
999 sprintf(meddesc,"MS");
1000 } else {
1001 wc=CalcWCbar(medval,length);
1002 sprintf(meddesc,"SH");
1003 }
1004
1005 Char_t hname[100];
1006 sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart,
1007 (Int_t)(medval*1000.),(Int_t)length);
1008
1009 TH1F *hist = new TH1F("hist",hname,1100,0.,1.1*wc);
1010 hist->SetXTitle("#Delta E [GeV]");
1011 hist->SetYTitle("p(#Delta E)");
1012 hist->SetLineColor(4);
1013
1014 Double_t rrrr = CalcR(wc,length);
1015 // loop on histogram channels
1016 for(Int_t bin=1; bin<=1100; bin++) {
1017 Double_t xxxx = hist->GetBinCenter(bin)/wc;
1018 Double_t continuous,discrete;
1019 Int_t ret=0;
1020 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1021 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1022 if(ret!=0){
1023 delete hist;
1024 return 0;
1025 };
1026 hist->SetBinContent(bin,continuous);
1027 }
1028 return hist;
1029}
1030
1031TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const
1032{
1033 // ipart = 1 for quark, 2 for gluon
1034 // medval a) qtransp = transport coefficient (GeV^2/fm)
1035 // b) mu = Debye mass (GeV)
1036 // length = path length in medium (fm)
1037 // Get from SW tables:
1038 // - discrete weight (probability to have NO energy loss)
1039 // - continuous weight, as a function of dE
1040 // compute up to max dE = 1.1*wc
1041 // (as an histogram named hContQW_<ipart>_<medval*1000>_<l>
1042 // and range [0,1.1] (1100 channels)
1043
1044 Double_t wc = 0;
1045 Char_t meddesc[100];
1046 if(fMultSoft) {
1047 wc=CalcWC(medval,length);
1048 sprintf(meddesc,"MS");
1049 } else {
1050 wc=CalcWCbar(medval,length);
1051 sprintf(meddesc,"SH");
1052 }
1053
1054 Char_t hname[100];
1055 sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart,
1056 (Int_t)(medval*1000.),(Int_t)length);
1057
1058 TH1F *histx = new TH1F("histx",hname,1100,0.,1.1);
1059 histx->SetXTitle("x = #Delta E/#omega_{c}");
1060 if(fMultSoft)
1061 histx->SetYTitle("p(#Delta E/#omega_{c})");
1062 else
1063 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1064 histx->SetLineColor(4);
1065
1066 Double_t rrrr = CalcR(wc,length);
1067 // loop on histogram channels
1068 for(Int_t bin=1; bin<=1100; bin++) {
1069 Double_t xxxx = histx->GetBinCenter(bin);
1070 Double_t continuous,discrete;
1071 Int_t ret=0;
1072 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1073 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1074 if(ret!=0){
1075 delete histx;
1076 return 0;
1077 };
1078 histx->SetBinContent(bin,continuous);
1079 }
1080 return histx;
1081}
1082
e99e3ed5 1083TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
b6d061b7 1084{
7a76a12e 1085 //compute energy loss histogram for
1086 //parton type, medium value, length and energy
1087
b6d061b7 1088 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1089 if(fMultSoft){
1090 dummy->SetQTransport(medval);
1091 dummy->InitMult();
1092 } else {
1093 dummy->SetMu(medval);
1094 dummy->InitSingleHard();
1095 }
1096 dummy->SampleEnergyLoss();
1097
1098 Char_t name[100];
1099 Char_t hname[100];
1100 if(ipart==1){
1101 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1102 sprintf(hname,"hLossQuarks");
1103 } else {
1104 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1105 sprintf(hname,"hLossGluons");
1106 }
1107
1108 TH1F *h = new TH1F(hname,name,250,0,250);
1109 for(Int_t i=0;i<100000;i++){
1110 //if(i % 1000 == 0) cout << "." << flush;
1111 Double_t loss=dummy->GetELossRandom(ipart,l,e);
1112 h->Fill(loss);
1113 }
1114
1115 h->SetStats(kTRUE);
1116
1117 delete dummy;
1118 return h;
1119}
1120
e99e3ed5 1121TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
b6d061b7 1122{
7a76a12e 1123 //compute energy loss histogram for
1124 //parton type, medium value,
1125 //length distribution and energy
1126
b6d061b7 1127 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1128 if(fMultSoft){
1129 dummy->SetQTransport(medval);
1130 dummy->InitMult();
1131 } else {
1132 dummy->SetMu(medval);
1133 dummy->InitSingleHard();
1134 }
1135 dummy->SampleEnergyLoss();
1136
1137 Char_t name[100];
1138 Char_t hname[100];
1139 if(ipart==1){
1140 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1141 sprintf(hname,"hLossQuarks");
1142 } else {
1143 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1144 sprintf(hname,"hLossGluons");
1145 }
1146
1147 TH1F *h = new TH1F(hname,name,250,0,250);
1148 for(Int_t i=0;i<100000;i++){
1149 //if(i % 1000 == 0) cout << "." << flush;
1150 Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
1151 h->Fill(loss);
1152 }
1153
1154 h->SetStats(kTRUE);
1155
1156 delete dummy;
1157 return h;
1158}
1159
1160void AliQuenchingWeights::PlotDiscreteWeights(Int_t len) const
1161{
7a76a12e 1162 //plot discrete weights for given length
1163
b6d061b7 1164 TCanvas *c;
1165 if(fMultSoft)
1166 c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
1167 else
1168 c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
1169 c->cd();
1170
1171 TH2F *hframe = new TH2F("hdisc","",2,0,5.1,2,0,1);
1172 hframe->SetStats(0);
1173 if(fMultSoft)
1174 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1175 else
1176 hframe->SetXTitle("#mu [GeV]");
1177 hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
1178 hframe->Draw();
1179
1180 TGraph *gq=new TGraph(20);
1181 Int_t i=0;
1182 if(fMultSoft) {
1183 for(Double_t q=0.05;q<=5.05;q+=0.25){
1184 Double_t disc,cont;
1185 CalcMult(1,1.0, q, len, cont, disc);
1186 gq->SetPoint(i,q,disc);i++;
1187 }
1188 } else {
1189 for(Double_t m=0.05;m<=5.05;m+=0.25){
1190 Double_t disc,cont;
1191 CalcSingleHard(1,1.0, m, len, cont, disc);
1192 gq->SetPoint(i,m,disc);i++;
1193 }
1194 }
1195 gq->SetMarkerStyle(20);
1196 gq->Draw("pl");
1197
1198 TGraph *gg=new TGraph(20);
1199 i=0;
1200 if(fMultSoft){
1201 for(Double_t q=0.05;q<=5.05;q+=0.25){
1202 Double_t disc,cont;
1203 CalcMult(2,1.0, q, 5., cont, disc);
1204 gg->SetPoint(i,q,disc);i++;
1205 }
1206 } else {
1207 for(Double_t m=0.05;m<=5.05;m+=0.25){
1208 Double_t disc,cont;
1209 CalcSingleHard(2,1.0, m, 5., cont, disc);
1210 gg->SetPoint(i,m,disc);i++;
1211 }
1212 }
1213 gg->SetMarkerStyle(24);
1214 gg->Draw("pl");
1215
1216 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1217 l1a->SetFillStyle(0);
1218 l1a->SetBorderSize(0);
1219 Char_t label[100];
1220 sprintf(label,"L = %d fm",len);
1221 l1a->AddEntry(gq,label,"");
1222 l1a->AddEntry(gq,"quark","pl");
1223 l1a->AddEntry(gg,"gluon","pl");
1224 l1a->Draw();
1225
1226 c->Update();
1227}
1228
1229void AliQuenchingWeights::PlotContWeights(Int_t itype,Int_t ell) const
1230{
7a76a12e 1231 //plot continous weights for
1232 //given parton type and length
1233
b6d061b7 1234 Float_t medvals[3];
1235 Char_t title[1024];
1236 Char_t name[1024];
1237 if(fMultSoft) {
1238 if(itype==1)
1239 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1240 else if(itype==2)
1241 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1242 else return;
1243 medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1244 sprintf(name,"ccont-ms-%d",itype);
1245 } else {
1246 if(itype==1)
1247 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1248 else if(itype==2)
1249 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1250 else return;
1251 medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1252 sprintf(name,"ccont-ms-%d",itype);
1253 }
1254
1255 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1256 c->cd();
1257 TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1258 h1->SetName("h1");
1259 h1->SetTitle(title);
1260 h1->SetStats(0);
1261 h1->SetLineColor(1);
1262 h1->DrawCopy();
1263 TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1264 h2->SetName("h2");
1265 h2->SetLineColor(2);
1266 h2->DrawCopy("SAME");
1267 TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1268 h3->SetName("h3");
1269 h3->SetLineColor(3);
1270 h3->DrawCopy("SAME");
1271
1272 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1273 l1a->SetFillStyle(0);
1274 l1a->SetBorderSize(0);
1275 Char_t label[100];
1276 sprintf(label,"L = %d fm",ell);
1277 l1a->AddEntry(h1,label,"");
1278 if(fMultSoft) {
1279 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1280 l1a->AddEntry(h1,label,"pl");
1281 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1282 l1a->AddEntry(h2,label,"pl");
1283 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1284 l1a->AddEntry(h3,label,"pl");
1285 } else {
1286 sprintf(label,"#mu = %.1f GeV",medvals[0]);
1287 l1a->AddEntry(h1,label,"pl");
1288 sprintf(label,"#mu = %.1f GeV",medvals[1]);
1289 l1a->AddEntry(h2,label,"pl");
1290 sprintf(label,"#mu = %.1f GeV",medvals[2]);
1291 l1a->AddEntry(h3,label,"pl");
1292 }
1293 l1a->Draw();
1294
1295 c->Update();
1296}
1297
1298void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const
1299{
7a76a12e 1300 //plot continous weights for
1301 //given parton type and medium value
1302
b6d061b7 1303 Char_t title[1024];
1304 Char_t name[1024];
1305 if(fMultSoft) {
1306 if(itype==1)
1307 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1308 else if(itype==2)
1309 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1310 else return;
1311 sprintf(name,"ccont2-ms-%d",itype);
1312 } else {
1313 if(itype==1)
1314 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1315 else if(itype==2)
1316 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1317 else return;
1318 sprintf(name,"ccont2-sh-%d",itype);
1319 }
1320 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1321 c->cd();
1322 TH1F *h1=ComputeQWHisto(itype,medval,8);
1323 h1->SetName("h1");
1324 h1->SetTitle(title);
1325 h1->SetStats(0);
1326 h1->SetLineColor(1);
1327 h1->DrawCopy();
1328 TH1F *h2=ComputeQWHisto(itype,medval,5);
1329 h2->SetName("h2");
1330 h2->SetLineColor(2);
1331 h2->DrawCopy("SAME");
1332 TH1F *h3=ComputeQWHisto(itype,medval,2);
1333 h3->SetName("h3");
1334 h3->SetLineColor(3);
1335 h3->DrawCopy("SAME");
1336
1337 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1338 l1a->SetFillStyle(0);
1339 l1a->SetBorderSize(0);
1340 Char_t label[100];
1341 if(fMultSoft)
1342 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval);
1343 else
1344 sprintf(label,"#mu = %.1f GeV",medval);
1345
1346 l1a->AddEntry(h1,label,"");
1347 l1a->AddEntry(h1,"L = 8 fm","pl");
1348 l1a->AddEntry(h2,"L = 5 fm","pl");
1349 l1a->AddEntry(h3,"L = 2 fm","pl");
1350 l1a->Draw();
1351
1352 c->Update();
1353}
1354
e99e3ed5 1355void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e) const
b6d061b7 1356{
7a76a12e 1357 //plot average energy loss for given length
1358 //and parton energy
1359
b6d061b7 1360 if(!fTablesLoaded){
1361 Error("CalcMult","Tables are not loaded.");
1362 return;
1363 }
1364
1365 Char_t title[1024];
1366 Char_t name[1024];
1367 if(fMultSoft){
1368 sprintf(title,"Average Loss for Multiple Scattering");
1369 sprintf(name,"cavgelossms");
1370 } else {
1371 sprintf(title,"Average Loss for Single Hard Scattering");
1372 sprintf(name,"cavgelosssh");
1373 }
1374
1375 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1376 c->cd();
1377 TH2F *hframe = new TH2F("avgloss",title,2,0,5.1,2,0,100);
1378 hframe->SetStats(0);
1379 if(fMultSoft)
1380 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1381 else
1382 hframe->SetXTitle("#mu [GeV]");
1383 hframe->SetYTitle("<E_{loss}> [GeV]");
1384 hframe->Draw();
1385
1386 TGraph *gq=new TGraph(20);
1387 Int_t i=0;
1388 for(Double_t v=0.05;v<=5.05;v+=0.25){
1389 TH1F *dummy=ComputeELossHisto(1,v,len,e);
1390 Double_t avgloss=dummy->GetMean();
1391 gq->SetPoint(i,v,avgloss);i++;
1392 delete dummy;
1393 }
1394 gq->SetMarkerStyle(20);
1395 gq->Draw("pl");
1396
1397 TGraph *gg=new TGraph(20);
1398 i=0;
1399 for(Double_t v=0.05;v<=5.05;v+=0.25){
1400 TH1F *dummy=ComputeELossHisto(2,v,len,e);
1401 Double_t avgloss=dummy->GetMean();
1402 gg->SetPoint(i,v,avgloss);i++;
1403 delete dummy;
1404 }
1405 gg->SetMarkerStyle(24);
1406 gg->Draw("pl");
1407
1408 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1409 l1a->SetFillStyle(0);
1410 l1a->SetBorderSize(0);
1411 Char_t label[100];
1412 sprintf(label,"L = %d fm",len);
1413 l1a->AddEntry(gq,label,"");
1414 l1a->AddEntry(gq,"quark","pl");
1415 l1a->AddEntry(gg,"gluon","pl");
1416 l1a->Draw();
1417
1418 c->Update();
1419}
1420
e99e3ed5 1421void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
b6d061b7 1422{
7a76a12e 1423 //plot average energy loss for given
1424 //length distribution and parton energy
1425
b6d061b7 1426 if(!fTablesLoaded){
1427 Error("CalcMult","Tables are not loaded.");
1428 return;
1429 }
1430
1431 Char_t title[1024];
1432 Char_t name[1024];
1433 if(fMultSoft){
1434 sprintf(title,"Average Loss for Multiple Scattering");
1435 sprintf(name,"cavgelossms2");
1436 } else {
1437 sprintf(title,"Average Loss for Single Hard Scattering");
1438 sprintf(name,"cavgelosssh2");
1439 }
1440
1441 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1442 c->cd();
1443 TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1444 hframe->SetStats(0);
1445 if(fMultSoft)
1446 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1447 else
1448 hframe->SetXTitle("#mu [GeV]");
1449 hframe->SetYTitle("<E_{loss}> [GeV]");
1450 hframe->Draw();
1451
1452 TGraph *gq=new TGraph(20);
1453 Int_t i=0;
1454 for(Double_t v=0.05;v<=5.05;v+=0.25){
1455 TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1456 Double_t avgloss=dummy->GetMean();
1457 gq->SetPoint(i,v,avgloss);i++;
1458 delete dummy;
1459 }
1460 gq->SetMarkerStyle(20);
1461 gq->Draw("pl");
1462
1463 TGraph *gg=new TGraph(20);
1464 i=0;
1465 for(Double_t v=0.05;v<=5.05;v+=0.25){
1466 TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1467 Double_t avgloss=dummy->GetMean();
1468 gg->SetPoint(i,v,avgloss);i++;
1469 delete dummy;
1470 }
1471 gg->SetMarkerStyle(24);
1472 gg->Draw("pl");
1473
1474 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1475 l1a->SetFillStyle(0);
1476 l1a->SetBorderSize(0);
1477 Char_t label[100];
1478 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1479 l1a->AddEntry(gq,label,"");
1480 l1a->AddEntry(gq,"quark","pl");
1481 l1a->AddEntry(gg,"gluon","pl");
1482 l1a->Draw();
1483
1484 c->Update();
1485}
1486
1487void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Int_t len) const
1488{
7a76a12e 1489 //plot relative energy loss for given
1490 //length and parton energy versus pt
1491
b6d061b7 1492 if(!fTablesLoaded){
1493 Error("CalcMult","Tables are not loaded.");
1494 return;
1495 }
1496
1497 Char_t title[1024];
1498 Char_t name[1024];
1499 if(fMultSoft){
1500 sprintf(title,"Relative Loss for Multiple Scattering");
1501 sprintf(name,"cavgelossvsptms");
1502 } else {
1503 sprintf(title,"Relative Loss for Single Hard Scattering");
1504 sprintf(name,"cavgelossvsptsh");
1505 }
1506
1507 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1508 c->cd();
1509 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
1510 hframe->SetStats(0);
1511 hframe->SetXTitle("p_{T} [GeV]");
1512 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
1513 hframe->Draw();
1514
1515 TGraph *gq=new TGraph(40);
1516 Int_t i=0;
1517 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1518 TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
1519 Double_t avgloss=dummy->GetMean();
1520 gq->SetPoint(i,pt,avgloss/pt);i++;
1521 delete dummy;
1522 }
1523 gq->SetMarkerStyle(20);
1524 gq->Draw("pl");
1525
1526 TGraph *gg=new TGraph(40);
1527 i=0;
1528 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1529 TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
1530 Double_t avgloss=dummy->GetMean();
1531 gg->SetPoint(i,pt,avgloss/pt);i++;
1532 delete dummy;
1533 }
1534 gg->SetMarkerStyle(24);
1535 gg->Draw("pl");
1536
1537 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1538 l1a->SetFillStyle(0);
1539 l1a->SetBorderSize(0);
1540 Char_t label[100];
1541 sprintf(label,"L = %d fm",len);
1542 l1a->AddEntry(gq,label,"");
1543 l1a->AddEntry(gq,"quark","pl");
1544 l1a->AddEntry(gg,"gluon","pl");
1545 l1a->Draw();
1546
1547 c->Update();
1548}
1549
1550void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
1551{
7a76a12e 1552 //plot relative energy loss for given
1553 //length distribution and parton energy versus pt
1554
b6d061b7 1555 if(!fTablesLoaded){
1556 Error("CalcMult","Tables are not loaded.");
1557 return;
1558 }
1559
1560 Char_t title[1024];
1561 Char_t name[1024];
1562 if(fMultSoft){
1563 sprintf(title,"Relative Loss for Multiple Scattering");
1564 sprintf(name,"cavgelossvsptms2");
1565 } else {
1566 sprintf(title,"Relative Loss for Single Hard Scattering");
1567 sprintf(name,"cavgelossvsptsh2");
1568 }
1569 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1570 c->cd();
1571 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
1572 hframe->SetStats(0);
1573 hframe->SetXTitle("p_{T} [GeV]");
1574 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
1575 hframe->Draw();
1576
1577 TGraph *gq=new TGraph(40);
1578 Int_t i=0;
1579 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1580 TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
1581 Double_t avgloss=dummy->GetMean();
1582 gq->SetPoint(i,pt,avgloss/pt);i++;
1583 delete dummy;
1584 }
1585 gq->SetMarkerStyle(20);
1586 gq->Draw("pl");
1587
1588 TGraph *gg=new TGraph(40);
1589 i=0;
1590 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1591 TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
1592 Double_t avgloss=dummy->GetMean();
1593 gg->SetPoint(i,pt,avgloss/pt);i++;
1594 delete dummy;
1595 }
1596 gg->SetMarkerStyle(24);
1597 gg->Draw("pl");
1598
1599 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1600 l1a->SetFillStyle(0);
1601 l1a->SetBorderSize(0);
1602 Char_t label[100];
1603 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1604 l1a->AddEntry(gq,label,"");
1605 l1a->AddEntry(gq,"quark","pl");
1606 l1a->AddEntry(gg,"gluon","pl");
1607 l1a->Draw();
1608
1609 c->Update();
1610}
1611