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