FASTSIM/AliQuenchingWeights.cxx

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