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