FASTSIM/AliQuenchingWeights.cxx

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