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