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