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