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