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