FASTSIM/AliQuenchingWeights.cxx

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