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