FASTSIM/AliQuenchingWeights.cxx

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