FASTSIM/AliQuenchingWeights.cxx

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