FASTSIM/AliQuenchingWeights.cxx

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