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