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