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