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