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 //            Origin:  C. Loizides   constantinos.loizides@cern.ch
  27 //                     A. Dainese    andrea.dainese@pd.infn.it
  28 //----------------------------------------------------------------------------
  29
  30 #include <Riostream.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TCanvas.h>
  34 #include <TGraph.h>
  35 #include <TROOT.h>
  36 #include <TSystem.h>
  37 #include <TLegend.h>
  38 #include "AliQuenchingWeights.h"
  39
  40 ClassImp(AliQuenchingWeights)
  41
  42 // conversion from fm to GeV^-1: 1 fm = fmGeV GeV^-1
  43 const Double_t AliQuenchingWeights::gkConvFmToInvGeV = 1./0.197;
  44
  45 // maximum value of R
  46 const Double_t AliQuenchingWeights::gkRMax = 40000.;
  47
  48 // counter for histogram labels
  49 Int_t AliQuenchingWeights::gCounter = 0;
  50
  51 AliQuenchingWeights::AliQuenchingWeights()
  52                    : TObject()
  53 {
  54   fTablesLoaded=kFALSE;
  55   fMultSoft=kTRUE;
  56   fHistos=0;
  57   SetMu();
  58   SetQTransport();
  59   SetECMethod();
  60   SetLengthMax();
  61   fLengthMaxOld=0;
  62   fInstanceNumber=gCounter++;
  63 }
  64
  65 AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a)
  66                    : TObject()
  67 {
  68   fTablesLoaded=kFALSE;
  69   fHistos=0;
  70   fLengthMaxOld=0;
  71   fMultSoft=a.GetMultSoft();;
  72   fMu=a.GetMu();
  73   fQTransport=a.GetQTransport();
  74   fECMethod=(kECMethod)a.GetECMethod();
  75   fLengthMax=a.GetLengthMax();
  76   fInstanceNumber=gCounter++;
  77   //Missing in the class is the pathname
  78   //to the tables, can be added if needed
  79 }
  80
  81 AliQuenchingWeights::~AliQuenchingWeights()
  82 {
  83   Reset();
  84 }
  85
  86 void AliQuenchingWeights::Reset()
  87 {
  88   if(!fHistos) return;
  89   for(Int_t l=0;l<fLengthMaxOld;l++){
  90     delete fHistos[0][l];
  91     delete fHistos[1][l];
  92   }
  93   delete[] fHistos;
  94   fHistos=0;
  95   fLengthMaxOld=0;
  96 }
  97
  98 void AliQuenchingWeights::SetECMethod(kECMethod type)
  99 {
 100   fECMethod=type;
 101   if(fECMethod==kDefault)
 102     Info("SetECMethod","Energy Constraint Method set to DEFAULT:\nIf (sampled energy loss > parton energy) then sampled energy loss = parton energy.");
 103   else
 104     Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy.");
 105 }
 106
 107 Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
 108 {
 109   // read in tables for multiple scattering approximation
 110   // path to continuum and to discrete part
 111
 112   fTablesLoaded = kFALSE;
 113   fMultSoft=kTRUE;
 114
 115   Char_t fname[1024];
 116   sprintf(fname,"%s",gSystem->ExpandPathName(contall));
 117   //PH  ifstream fincont(fname);
 118   fstream fincont(fname,ios::in);
 119 #if defined(__HP_aCC) || defined(__DECCXX)
 120   if(!fincont.rdbuf()->is_open()) return -1;
 121 #else
 122   if(!fincont.is_open()) return -1;
 123 #endif
 124
 125   Int_t nn=0; //quarks
 126   while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
 127         fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
 128         fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>
 129         fcaq[13][nn]>>
 130         fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>
 131         fcaq[18][nn]>>
 132         fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>
 133         fcaq[23][nn]>>
 134         fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>
 135         fcaq[28][nn]>>
 136         fcaq[29][nn])
 137     {
 138       nn++;
 139       if(nn==261) break;
 140     }
 141
 142   nn=0;       //gluons
 143   while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
 144         fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
 145         fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>
 146         fcag[13][nn]>>
 147         fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>
 148         fcag[18][nn]>>
 149         fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>
 150         fcag[23][nn]>>
 151         fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>
 152         fcag[28][nn]>>
 153         fcag[29][nn]) {
 154     nn++;
 155     if(nn==261) break;
 156   }
 157   fincont.close();
 158
 159   sprintf(fname,"%s",gSystem->ExpandPathName(discall));
 160   //PH  ifstream findisc(fname);
 161   fstream findisc(fname,ios::in);
 162 #if defined(__HP_aCC) || defined(__DECCXX)
 163   if(!findisc.rdbuf()->is_open()) return -1;
 164 #else
 165   if(!findisc.is_open()) return -1;
 166 #endif
 167
 168   nn=0; //quarks
 169   while(findisc>>frrr[nn]>>fdaq[nn]) {
 170     nn++;
 171     if(nn==30) break;
 172   }
 173   nn=0; //gluons
 174   while(findisc>>frrrg[nn]>>fdag[nn]) {
 175     nn++;
 176     if(nn==30) break;
 177   }
 178   findisc.close();
 179   fTablesLoaded = kTRUE;
 180   return 0;
 181 }
 182
 183 /*
 184 C***************************************************************************
 185 C       Quenching Weights for Multiple Soft Scattering
 186 C               February 10, 2003
 187 C
 188 C       Refs:
 189 C
 190 C  Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
 191 C
 192 C  Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
 193 C
 194 C
 195 C   This package contains quenching weights for gluon radiation in the
 196 C   multiple soft scattering approximation.
 197 C
 198 C   swqmult returns the quenching weight for a quark (ipart=1) or
 199 C   a gluon (ipart=2) traversing a medium with transport coeficient q and
 200 C   length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where
 201 C   wc=0.5*q*L^2 and w is the energy radiated. The output values are
 202 C   the continuous and discrete (prefactor of the delta function) parts
 203 C   of the quenching weights.
 204 C
 205 C   In order to use this routine, the files cont.all and disc.all need to be
 206 C   in the working directory.
 207 C
 208 C   An initialization of the tables is needed by doing call initmult before
 209 C   using swqmult.
 210 C
 211 C   Please, send us any comment:
 212 C
 213 C       urs.wiedemann@cern.ch
 214 C       carlos.salgado@cern.ch
 215 C
 216 C
 217 C-------------------------------------------------------------------
 218
 219       SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete)
 220 *
 221       REAL*8           xx(400), daq(30), caq(30,261), rrr(30)
 222       COMMON /dataqua/    xx, daq, caq, rrr
 223 *
 224       REAL*8           xxg(400), dag(30), cag(30,261), rrrg(30)
 225       COMMON /dataglu/    xxg, dag, cag, rrrg
 226
 227       REAL*8           rrrr,xxxx, continuous, discrete
 228       REAL*8           rrin, xxin
 229       INTEGER          nrlow, nrhigh, nxlow, nxhigh
 230       REAL*8           rrhigh, rrlow, rfraclow, rfrachigh
 231       REAL*8           xfraclow, xfrachigh
 232       REAL*8           clow, chigh
 233 *
 234       rrin = rrrr
 235       xxin = xxxx
 236 *
 237       nxlow = int(xxin/0.005) + 1
 238       nxhigh = nxlow + 1
 239       xfraclow = (xx(nxhigh)-xxin)/0.005
 240       xfrachigh = (xxin - xx(nxlow))/0.005
 241 *
 242       do 666, nr=1,30
 243          if (rrin.lt.rrr(nr)) then
 244             rrhigh = rrr(nr)
 245          else
 246             rrhigh = rrr(nr-1)
 247             rrlow = rrr(nr)
 248             nrlow = nr
 249             nrhigh = nr-1
 250             goto 665
 251          endif
 252  666     enddo
 253  665     continue
 254 *
 255       rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
 256       rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
 257 *
 258       if(ipart.eq.1) then
 259       clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh)
 260       chigh = xfraclow*caq(nrhigh,nxlow)+xfrachigh*caq(nrhigh,nxhigh)
 261       else
 262       clow = xfraclow*cag(nrlow,nxlow)+xfrachigh*cag(nrlow,nxhigh)
 263       chigh = xfraclow*cag(nrhigh,nxlow)+xfrachigh*cag(nrhigh,nxhigh)
 264       endif
 265
 266       continuous = rfraclow*clow + rfrachigh*chigh
 267
 268       if(ipart.eq.1) then
 269       discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh)
 270       else
 271       discrete = rfraclow*dag(nrlow) + rfrachigh*dag(nrhigh)
 272       endif
 273 *
 274       END
 275
 276       subroutine initmult
 277       REAL*8           xxq(400), daq(30), caq(30,261), rrr(30)
 278       COMMON /dataqua/    xxq, daq, caq, rrr
 279 *
 280       REAL*8           xxg(400), dag(30), cag(30,261), rrrg(30)
 281       COMMON /dataglu/    xxg, dag, cag, rrrg
 282 *
 283       OPEN(UNIT=20,FILE='cont.all',STATUS='OLD',ERR=90)
 284       do 110 nn=1,261
 285       read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn),
 286      +     caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn),
 287      +     caq(9,nn), caq(10,nn), caq(11,nn), caq(12,nn),
 288      +     caq(13,nn),
 289      +     caq(14,nn), caq(15,nn), caq(16,nn), caq(17,nn),
 290      +     caq(18,nn),
 291      +     caq(19,nn), caq(20,nn), caq(21,nn), caq(22,nn),
 292      +     caq(23,nn),
 293      +     caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn),
 294      +     caq(28,nn),
 295      +     caq(29,nn), caq(30,nn)
 296  110     continue
 297       do 111 nn=1,261
 298       read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn),
 299      +     cag(4,nn), cag(5,nn), cag(6,nn), cag(7,nn), cag(8,nn),
 300      +     cag(9,nn), cag(10,nn), cag(11,nn), cag(12,nn),
 301      +     cag(13,nn),
 302      +     cag(14,nn), cag(15,nn), cag(16,nn), cag(17,nn),
 303      +     cag(18,nn),
 304      +     cag(19,nn), cag(20,nn), cag(21,nn), cag(22,nn),
 305      +     cag(23,nn),
 306      +     cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn),
 307      +     cag(28,nn),
 308      +     cag(29,nn), cag(30,nn)
 309  111     continue
 310       close(20)
 311 *
 312       OPEN(UNIT=21,FILE='disc.all',STATUS='OLD',ERR=91)
 313       do 112 nn=1,30
 314       read (21,*) rrr(nn), daq(nn)
 315  112     continue
 316       do 113 nn=1,30
 317       read (21,*) rrrg(nn), dag(nn)
 318  113     continue
 319       close(21)
 320 *
 321       goto 888
 322  90   PRINT*, 'input - output error'
 323  91   PRINT*, 'input - output error #2'
 324  888  continue
 325
 326       end
 327
 328 =======================================================================
 329
 330    Adapted to ROOT macro by A. Dainese - 13/07/2003
 331    ported to class by C. Loizides - 12/02/2004
 332 */
 333
 334 Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
 335                               Double_t &continuous,Double_t &discrete) const
 336 {
 337   // Calculate Multiple Scattering approx.
 338   // weights for given parton type,
 339   // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
 340
 341   // read-in data before first call
 342   if(!fTablesLoaded){
 343     Error("CalcMult","Tables are not loaded.");
 344     return -1;
 345   }
 346   if(!fMultSoft){
 347     Error("CalcMult","Tables are not loaded for Multiple Scattering.");
 348     return -1;
 349   }
 350
 351   Double_t rrin = rrrr;
 352   Double_t xxin = xxxx;
 353
 354   Int_t nxlow     = (Int_t)(xxin/0.01) + 1;
 355   Int_t nxhigh    = nxlow + 1;
 356   Double_t xfraclow  = (fxx[nxhigh-1]-xxin)/0.01;
 357   Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01;
 358
 359   //why this?
 360   if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD
 361   if(rrin>=frrr[0])  rrin = 0.95*frrr[0];  // AD
 362
 363   Int_t nrlow=0,nrhigh=0;
 364   Double_t rrhigh=0,rrlow=0;
 365   for(Int_t nr=1; nr<=30; nr++) {
 366     if(rrin<frrr[nr-1]) {
 367       rrhigh = frrr[nr-1];
 368     } else {
 369       rrhigh = frrr[nr-1-1];
 370       rrlow  = frrr[nr-1];
 371       nrlow  = nr;
 372       nrhigh = nr-1;
 373       break;
 374     }
 375   }
 376
 377   rrin = rrrr; // AD
 378
 379   Double_t rfraclow  = (rrhigh-rrin)/(rrhigh-rrlow);
 380   Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
 381
 382   //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
 383
 384   Double_t clow=0,chigh=0;
 385   if(ipart==1) {
 386     clow  = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
 387     chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
 388   } else {
 389     clow  = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
 390     chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
 391   }
 392
 393   continuous = rfraclow*clow + rfrachigh*chigh;
 394   //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
 395   //     rfraclow,clow,rfrachigh,chigh,continuous);
 396
 397   if(ipart==1) {
 398     discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
 399   } else {
 400     discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
 401   }
 402
 403   return 0;
 404 }
 405
 406 Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *discall)
 407 {
 408   // read in tables for Single Hard Approx.
 409   // path to continuum and to discrete part
 410
 411   fTablesLoaded = kFALSE;
 412   fMultSoft=kFALSE;
 413
 414   Char_t fname[1024];
 415   sprintf(fname,"%s",gSystem->ExpandPathName(contall));
 416   //PH  ifstream fincont(fname);
 417   fstream fincont(fname,ios::in);
 418 #if defined(__HP_aCC) || defined(__DECCXX)
 419   if(!fincont.rdbuf()->is_open()) return -1;
 420 #else
 421   if(!fincont.is_open()) return -1;
 422 #endif
 423
 424   Int_t nn=0; //quarks
 425   while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
 426         fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
 427         fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>
 428         fcaq[13][nn]>>
 429         fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>
 430         fcaq[18][nn]>>
 431         fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>
 432         fcaq[23][nn]>>
 433         fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>
 434         fcaq[28][nn]>>
 435         fcaq[29][nn])
 436     {
 437       nn++;
 438       if(nn==261) break;
 439     }
 440
 441   nn=0;       //gluons
 442   while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
 443         fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
 444         fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>
 445         fcag[13][nn]>>
 446         fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>
 447         fcag[18][nn]>>
 448         fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>
 449         fcag[23][nn]>>
 450         fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>
 451         fcag[28][nn]>>
 452         fcag[29][nn]) {
 453     nn++;
 454     if(nn==261) break;
 455   }
 456   fincont.close();
 457
 458   sprintf(fname,"%s",gSystem->ExpandPathName(discall));
 459   //PH  ifstream findisc(fname);
 460   fstream findisc(fname,ios::in);
 461 #if defined(__HP_aCC) || defined(__DECCXX)
 462   if(!findisc.rdbuf()->is_open()) return -1;
 463 #else
 464   if(!findisc.is_open()) return -1;
 465 #endif
 466
 467   nn=0; //quarks
 468   while(findisc>>frrr[nn]>>fdaq[nn]) {
 469     nn++;
 470     if(nn==30) break;
 471   }
 472   nn=0; //gluons
 473   while(findisc>>frrrg[nn]>>fdag[nn]) {
 474     nn++;
 475     if(nn==30) break;
 476   }
 477   findisc.close();
 478
 479   fTablesLoaded = kTRUE;
 480   return 0;
 481 }
 482
 483 /*
 484 C***************************************************************************
 485 C       Quenching Weights for Single Hard Scattering
 486 C               February 20, 2003
 487 C
 488 C       Refs:
 489 C
 490 C  Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
 491 C
 492 C  Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
 493 C
 494 C
 495 C   This package contains quenching weights for gluon radiation in the
 496 C   single hard scattering approximation.
 497 C
 498 C   swqlin returns the quenching weight for a quark (ipart=1) or
 499 C   a gluon (ipart=2) traversing a medium with Debye screening mass mu and
 500 C   length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where
 501 C   wc=0.5*mu^2*L and w is the energy radiated. The output values are
 502 C   the continuous and discrete (prefactor of the delta function) parts
 503 C   of the quenching weights.
 504 C
 505 C   In order to use this routine, the files contlin.all and disclin.all
 506 C   need to be in the working directory.
 507 C
 508 C   An initialization of the tables is needed by doing call initlin before
 509 C   using swqlin.
 510 C
 511 C   Please, send us any comment:
 512 C
 513 C       urs.wiedemann@cern.ch
 514 C       carlos.salgado@cern.ch
 515 C
 516 C
 517 C-------------------------------------------------------------------
 518
 519
 520       SUBROUTINE swqlin(ipart,rrrr,xxxx,continuous,discrete)
 521 *
 522       REAL*8           xx(400), dalq(30), calq(30,261), rrr(30)
 523       COMMON /datalinqua/    xx, dalq, calq, rrr
 524 *
 525       REAL*8           xxlg(400), dalg(30), calg(30,261), rrrlg(30)
 526       COMMON /datalinglu/    xxlg, dalg, calg, rrrlg
 527
 528       REAL*8           rrrr,xxxx, continuous, discrete
 529       REAL*8           rrin, xxin
 530       INTEGER          nrlow, nrhigh, nxlow, nxhigh
 531       REAL*8           rrhigh, rrlow, rfraclow, rfrachigh
 532       REAL*8           xfraclow, xfrachigh
 533       REAL*8           clow, chigh
 534 *
 535       rrin = rrrr
 536       xxin = xxxx
 537 *
 538       nxlow = int(xxin/0.038) + 1
 539       nxhigh = nxlow + 1
 540       xfraclow = (xx(nxhigh)-xxin)/0.038
 541       xfrachigh = (xxin - xx(nxlow))/0.038
 542 *
 543       do 666, nr=1,30
 544          if (rrin.lt.rrr(nr)) then
 545             rrhigh = rrr(nr)
 546          else
 547             rrhigh = rrr(nr-1)
 548             rrlow = rrr(nr)
 549             nrlow = nr
 550             nrhigh = nr-1
 551             goto 665
 552          endif
 553  666     enddo
 554  665     continue
 555 *
 556       rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
 557       rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
 558 *
 559       if(ipart.eq.1) then
 560       clow = xfraclow*calq(nrlow,nxlow)+xfrachigh*calq(nrlow,nxhigh)
 561       chigh = xfraclow*calq(nrhigh,nxlow)+xfrachigh*calq(nrhigh,nxhigh)
 562       else
 563       clow = xfraclow*calg(nrlow,nxlow)+xfrachigh*calg(nrlow,nxhigh)
 564       chigh = xfraclow*calg(nrhigh,nxlow)+xfrachigh*calg(nrhigh,nxhigh)
 565       endif
 566
 567       continuous = rfraclow*clow + rfrachigh*chigh
 568
 569       if(ipart.eq.1) then
 570       discrete = rfraclow*dalq(nrlow) + rfrachigh*dalq(nrhigh)
 571       else
 572       discrete = rfraclow*dalg(nrlow) + rfrachigh*dalg(nrhigh)
 573       endif
 574 *
 575       END
 576
 577       subroutine initlin
 578       REAL*8           xxlq(400), dalq(30), calq(30,261), rrr(30)
 579       COMMON /datalinqua/    xxlq, dalq, calq, rrr
 580 *
 581       REAL*8           xxlg(400), dalg(30), calg(30,261), rrrlg(30)
 582       COMMON /datalinglu/    xxlg, dalg, calg, rrrlg
 583 *
 584       OPEN(UNIT=20,FILE='contlin.all',STATUS='OLD',ERR=90)
 585       do 110 nn=1,261
 586       read (20,*) xxlq(nn), calq(1,nn), calq(2,nn), calq(3,nn),
 587      +     calq(4,nn), calq(5,nn), calq(6,nn), calq(7,nn), calq(8,nn),
 588      +     calq(9,nn), calq(10,nn), calq(11,nn), calq(12,nn),
 589      +     calq(13,nn),
 590      +     calq(14,nn), calq(15,nn), calq(16,nn), calq(17,nn),
 591      +     calq(18,nn),
 592      +     calq(19,nn), calq(20,nn), calq(21,nn), calq(22,nn),
 593      +     calq(23,nn),
 594      +     calq(24,nn), calq(25,nn), calq(26,nn), calq(27,nn),
 595      +     calq(28,nn),
 596      +     calq(29,nn), calq(30,nn)
 597  110     continue
 598       do 111 nn=1,261
 599       read (20,*) xxlg(nn), calg(1,nn), calg(2,nn), calg(3,nn),
 600      +     calg(4,nn), calg(5,nn), calg(6,nn), calg(7,nn), calg(8,nn),
 601      +     calg(9,nn), calg(10,nn), calg(11,nn), calg(12,nn),
 602      +     calg(13,nn),
 603      +     calg(14,nn), calg(15,nn), calg(16,nn), calg(17,nn),
 604      +     calg(18,nn),
 605      +     calg(19,nn), calg(20,nn), calg(21,nn), calg(22,nn),
 606      +     calg(23,nn),
 607      +     calg(24,nn), calg(25,nn), calg(26,nn), calg(27,nn),
 608      +     calg(28,nn),
 609      +     calg(29,nn), calg(30,nn)
 610  111     continue
 611       close(20)
 612 *
 613       OPEN(UNIT=21,FILE='disclin.all',STATUS='OLD',ERR=91)
 614       do 112 nn=1,30
 615       read (21,*) rrr(nn), dalq(nn)
 616  112     continue
 617       do 113 nn=1,30
 618       read (21,*) rrrlg(nn), dalg(nn)
 619  113     continue
 620       close(21)
 621 *
 622       goto 888
 623  90   PRINT*, 'input - output error'
 624  91   PRINT*, 'input - output error #2'
 625  888  continue
 626
 627       end
 628
 629 =======================================================================
 630
 631    Ported to class by C. Loizides - 17/02/2004
 632
 633 */
 634
 635 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx,
 636                                     Double_t &continuous,Double_t &discrete) const
 637 {
 638   // calculate Single Hard approx.
 639   // weights for given parton type,
 640   // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
 641
 642   // read-in data before first call
 643   if(!fTablesLoaded){
 644     Error("CalcMult","Tables are not loaded.");
 645     return -1;
 646   }
 647   if(!fMultSoft){
 648     Error("CalcMult","Tables are not loaded for Single Hard Scattering.");
 649     return -1;
 650   }
 651
 652   Double_t rrin = rrrr;
 653   Double_t xxin = xxxx;
 654
 655   Int_t nxlow     = (Int_t)(xxin/0.038) + 1;
 656   Int_t nxhigh    = nxlow + 1;
 657   Double_t xfraclow  = (fxx[nxhigh-1]-xxin)/0.038;
 658   Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.038;
 659
 660   //why this?
 661   if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD
 662   if(rrin>=frrr[0])  rrin = 0.95*frrr[0];  // AD
 663
 664   Int_t nrlow=0,nrhigh=0;
 665   Double_t rrhigh=0,rrlow=0;
 666   for(Int_t nr=1; nr<=30; nr++) {
 667     if(rrin<frrr[nr-1]) {
 668       rrhigh = frrr[nr-1];
 669     } else {
 670       rrhigh = frrr[nr-1-1];
 671       rrlow  = frrr[nr-1];
 672       nrlow  = nr;
 673       nrhigh = nr-1;
 674       break;
 675     }
 676   }
 677
 678   rrin = rrrr; // AD
 679
 680   Double_t rfraclow  = (rrhigh-rrin)/(rrhigh-rrlow);
 681   Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
 682
 683   //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
 684
 685   Double_t clow=0,chigh=0;
 686   if(ipart==1) {
 687     clow  = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
 688     chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
 689   } else {
 690     clow  = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
 691     chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
 692   }
 693
 694   continuous = rfraclow*clow + rfrachigh*chigh;
 695   //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
 696   //     rfraclow,clow,rfrachigh,chigh,continuous);
 697
 698   if(ipart==1) {
 699     discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
 700   } else {
 701     discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
 702   }
 703
 704   return 0;
 705 }
 706
 707 Int_t AliQuenchingWeights::CalcMult(Int_t ipart,
 708                               Double_t w,Double_t qtransp,Double_t length,
 709                               Double_t &continuous,Double_t &discrete) const
 710 {
 711   Double_t wc=CalcWC(qtransp,length);
 712   Double_t rrrr=CalcR(wc,length);
 713   Double_t xxxx=w/wc;
 714   return CalcMult(ipart,rrrr,xxxx,continuous,discrete);
 715 }
 716
 717 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
 718                                     Double_t w,Double_t mu,Double_t length,
 719                                     Double_t &continuous,Double_t &discrete) const
 720 {
 721   Double_t wcbar=CalcWCbar(mu,length);
 722   Double_t rrrr=CalcR(wcbar,length);
 723   Double_t xxxx=w/wcbar;
 724   return CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
 725 }
 726
 727 Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const
 728 {
 729   Double_t R = wc*l*gkConvFmToInvGeV;
 730   if(R>gkRMax) {
 731     Warning("CalcR","Value of R = %.2f; should be less than %.2f",R,gkRMax);
 732     return -R;
 733   }
 734   return R;
 735 }
 736
 737 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const
 738 {
 739   // return DeltaE for MS or SH scattering
 740   // for given parton type, length and energy
 741   // Dependant on ECM (energy constraint method)
 742   // e is used to determine where to set bins to zero.
 743
 744   if(!fHistos){
 745     Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
 746     return -1000.;
 747   }
 748   if((ipart<1) || (ipart>2)) {
 749     Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
 750     return -1000;
 751   }
 752
 753   Int_t l=(Int_t)length;
 754   if((length-(Double_t)l)>0.5) l++;
 755   if(l<=0) return 0.;
 756   if(l>fLengthMax) l=fLengthMax;
 757
 758   if(fECMethod==kReweight){
 759     TH1F *dummy=new TH1F(*fHistos[ipart-1][l-1]);
 760     dummy->SetName("dummy");
 761     for(Int_t bin=dummy->FindBin(e)+1;bin<=1100;bin++){
 762       dummy->SetBinContent(bin,0.);
 763     }
 764     dummy->Scale(1./dummy->Integral());
 765     Double_t ret=dummy->GetRandom();
 766     delete dummy;
 767     return ret;
 768     //****** !! ALTERNATIVE WAY OF DOING IT !! ***
 769     //Double_t ret = 2.*e;
 770     //while(ret>e) ret=fHistos[ipart-1][l-1]->GetRandom();
 771     //return ret;
 772     //********************************************
 773   } else { //kDefault
 774     Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
 775     if(ret>e) return e;
 776     return ret;
 777   }
 778 }
 779
 780 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const
 781 {
 782   //return quenched parton energy
 783   //for given parton type, length and energy
 784
 785   Double_t loss=GetELossRandom(ipart,length,e);
 786   return e-loss;
 787 }
 788
 789 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const
 790 {
 791   if(!hell){
 792     Warning("GetELossRandom","Pointer to length distribution is NULL.");
 793     return 0.;
 794   }
 795   Double_t ell=hell->GetRandom();
 796   return GetELossRandom(ipart,ell,e);
 797 }
 798
 799 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e)  const
 800 {
 801   //return quenched parton energy
 802   //for given parton type, length distribution and energy
 803
 804   Double_t loss=GetELossRandom(ipart,hell,e);
 805   return e-loss;
 806 }
 807
 808 Int_t AliQuenchingWeights::SampleEnergyLoss()
 809 {
 810   // Has to be called to fill the histograms
 811   //
 812   // For stored values fQTransport loop over
 813   // particle type and length = 1 to fMaxLength (fm)
 814   // to fill energy loss histos
 815   //
 816   //    Take histogram of continuous weights
 817   //    Take discrete_weight
 818   //    If discrete_weight > 1, put all channels to 0, except channel 1
 819   //    Fill channel 1 with discrete_weight/(1-discrete_weight)*integral
 820
 821   // read-in data before first call
 822   if(!fTablesLoaded){
 823     Error("CalcMult","Tables are not loaded.");
 824     return -1;
 825   }
 826
 827   if(fMultSoft) {
 828     Int_t lmax=CalcLengthMax(fQTransport);
 829     if(fLengthMax>lmax){
 830       Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,gkRMax);
 831       fLengthMax=lmax;
 832     }
 833   } else {
 834       Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested.");
 835   }
 836
 837   Reset();
 838   fHistos=new TH1F**[2];
 839   fHistos[0]=new TH1F*[fLengthMax];
 840   fHistos[1]=new TH1F*[fLengthMax];
 841   fLengthMaxOld=fLengthMax; //remember old value in case
 842                             //user wants to reset
 843
 844   Int_t medvalue=0;
 845   Char_t meddesc[100];
 846   if(fMultSoft) {
 847     medvalue=(Int_t)(fQTransport*1000.);
 848     sprintf(meddesc,"MS");
 849   } else {
 850     medvalue=(Int_t)(fMu*1000.);
 851     sprintf(meddesc,"SH");
 852   }
 853
 854   for(Int_t ipart=1;ipart<=2;ipart++){
 855     for(Int_t l=1;l<=fLengthMax;l++){
 856
 857       Char_t hname[100];
 858       sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
 859       Double_t wc = CalcWC(l);
 860       fHistos[ipart-1][l-1] = new TH1F(hname,hname,1100,0.,1.1*wc);
 861       fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
 862       fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
 863       fHistos[ipart-1][l-1]->SetLineColor(4);
 864
 865       Double_t rrrr = CalcR(wc,l);
 866       Double_t discrete=0.;
 867       // loop on histogram channels
 868       for(Int_t bin=1; bin<=1100; bin++) {
 869         Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
 870         Double_t continuous;
 871         CalcMult(ipart,rrrr,xxxx,continuous,discrete);
 872         fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
 873       }
 874       // add discrete part to distribution
 875       if(discrete>=1.)
 876         for(Int_t bin=2;bin<=1100;bin++)
 877           fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
 878       else {
 879         Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,1100);
 880         fHistos[ipart-1][l-1]->Fill(0.,val);
 881       }
 882       Double_t hint=fHistos[ipart-1][l-1]->Integral(1,1100);
 883       fHistos[ipart-1][l-1]->Scale(1./hint);
 884     }
 885   }
 886   return 0;
 887 }
 888
 889 const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Int_t l) const
 890 {
 891   if(!fHistos){
 892     Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
 893     return 0;
 894   }
 895   if((ipart<1) || (ipart>2)) {
 896     Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
 897     return 0;
 898   }
 899
 900   if(l<=0) return 0;
 901   if(l>fLengthMax) l=fLengthMax;
 902
 903   return fHistos[ipart-1][l-1];
 904 }
 905
 906 TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const
 907 {
 908   // ipart = 1 for quark, 2 for gluon
 909   // medval a) qtransp = transport coefficient (GeV^2/fm)
 910   //        b) mu      = Debye mass (GeV)
 911   // length = path length in medium (fm)
 912   // Get from SW tables:
 913   // - discrete weight (probability to have NO energy loss)
 914   // - continuous weight, as a function of dE
 915   // compute up to max dE = 1.1*wc
 916   //   (as an histogram named hContQW_<ipart>_<medval*1000>_<l>
 917   //   and range [0,1.1*wc] (1100 channels)
 918
 919   Double_t wc = 0;
 920   Char_t meddesc[100];
 921   if(fMultSoft) {
 922     wc=CalcWC(medval,length);
 923     sprintf(meddesc,"MS");
 924   } else {
 925     wc=CalcWCbar(medval,length);
 926     sprintf(meddesc,"SH");
 927   }
 928
 929   Char_t hname[100];
 930   sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart,
 931                 (Int_t)(medval*1000.),(Int_t)length);
 932
 933   TH1F *hist = new TH1F("hist",hname,1100,0.,1.1*wc);
 934   hist->SetXTitle("#Delta E [GeV]");
 935   hist->SetYTitle("p(#Delta E)");
 936   hist->SetLineColor(4);
 937
 938   Double_t rrrr = CalcR(wc,length);
 939   // loop on histogram channels
 940   for(Int_t bin=1; bin<=1100; bin++) {
 941     Double_t xxxx = hist->GetBinCenter(bin)/wc;
 942     Double_t continuous,discrete;
 943     Int_t ret=0;
 944     if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
 945     else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
 946     if(ret!=0){
 947       delete hist;
 948       return 0;
 949     };
 950     hist->SetBinContent(bin,continuous);
 951   }
 952   return hist;
 953 }
 954
 955 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const
 956 {
 957   // ipart = 1 for quark, 2 for gluon
 958   // medval a) qtransp = transport coefficient (GeV^2/fm)
 959   //        b) mu      = Debye mass (GeV)
 960   // length = path length in medium (fm)
 961   // Get from SW tables:
 962   // - discrete weight (probability to have NO energy loss)
 963   // - continuous weight, as a function of dE
 964   // compute up to max dE = 1.1*wc
 965   //   (as an histogram named hContQW_<ipart>_<medval*1000>_<l>
 966   //   and range [0,1.1] (1100 channels)
 967
 968   Double_t wc = 0;
 969   Char_t meddesc[100];
 970   if(fMultSoft) {
 971     wc=CalcWC(medval,length);
 972     sprintf(meddesc,"MS");
 973   } else {
 974     wc=CalcWCbar(medval,length);
 975     sprintf(meddesc,"SH");
 976   }
 977
 978   Char_t hname[100];
 979   sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart,
 980                 (Int_t)(medval*1000.),(Int_t)length);
 981
 982   TH1F *histx = new TH1F("histx",hname,1100,0.,1.1);
 983   histx->SetXTitle("x = #Delta E/#omega_{c}");
 984   if(fMultSoft)
 985     histx->SetYTitle("p(#Delta E/#omega_{c})");
 986   else
 987     histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
 988   histx->SetLineColor(4);
 989
 990   Double_t rrrr = CalcR(wc,length);
 991   // loop on histogram channels
 992   for(Int_t bin=1; bin<=1100; bin++) {
 993     Double_t xxxx = histx->GetBinCenter(bin);
 994     Double_t continuous,discrete;
 995     Int_t ret=0;
 996     if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
 997     else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
 998     if(ret!=0){
 999       delete histx;
1000       return 0;
1001     };
1002     histx->SetBinContent(bin,continuous);
1003   }
1004   return histx;
1005 }
1006
1007 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
1008 {
1009   AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1010   if(fMultSoft){
1011     dummy->SetQTransport(medval);
1012     dummy->InitMult();
1013   } else {
1014     dummy->SetMu(medval);
1015     dummy->InitSingleHard();
1016   }
1017   dummy->SampleEnergyLoss();
1018
1019   Char_t name[100];
1020   Char_t hname[100];
1021   if(ipart==1){
1022     sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1023     sprintf(hname,"hLossQuarks");
1024   } else {
1025     sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1026     sprintf(hname,"hLossGluons");
1027   }
1028
1029   TH1F *h = new TH1F(hname,name,250,0,250);
1030   for(Int_t i=0;i<100000;i++){
1031     //if(i % 1000 == 0) cout << "." << flush;
1032     Double_t loss=dummy->GetELossRandom(ipart,l,e);
1033     h->Fill(loss);
1034   }
1035
1036   h->SetStats(kTRUE);
1037
1038   delete dummy;
1039   return h;
1040 }
1041
1042 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
1043 {
1044   AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1045   if(fMultSoft){
1046     dummy->SetQTransport(medval);
1047     dummy->InitMult();
1048   } else {
1049     dummy->SetMu(medval);
1050     dummy->InitSingleHard();
1051   }
1052   dummy->SampleEnergyLoss();
1053
1054   Char_t name[100];
1055   Char_t hname[100];
1056   if(ipart==1){
1057     sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1058     sprintf(hname,"hLossQuarks");
1059   } else {
1060     sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1061     sprintf(hname,"hLossGluons");
1062   }
1063
1064   TH1F *h = new TH1F(hname,name,250,0,250);
1065   for(Int_t i=0;i<100000;i++){
1066     //if(i % 1000 == 0) cout << "." << flush;
1067     Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
1068     h->Fill(loss);
1069   }
1070
1071   h->SetStats(kTRUE);
1072
1073   delete dummy;
1074   return h;
1075 }
1076
1077 void AliQuenchingWeights::PlotDiscreteWeights(Int_t len) const
1078 {
1079   TCanvas *c;
1080   if(fMultSoft)
1081     c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
1082   else
1083     c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
1084   c->cd();
1085
1086   TH2F *hframe = new TH2F("hdisc","",2,0,5.1,2,0,1);
1087   hframe->SetStats(0);
1088   if(fMultSoft)
1089     hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1090   else
1091     hframe->SetXTitle("#mu [GeV]");
1092   hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
1093   hframe->Draw();
1094
1095   TGraph *gq=new TGraph(20);
1096   Int_t i=0;
1097   if(fMultSoft) {
1098     for(Double_t q=0.05;q<=5.05;q+=0.25){
1099       Double_t disc,cont;
1100       CalcMult(1,1.0, q, len, cont, disc);
1101       gq->SetPoint(i,q,disc);i++;
1102     }
1103   } else {
1104     for(Double_t m=0.05;m<=5.05;m+=0.25){
1105       Double_t disc,cont;
1106       CalcSingleHard(1,1.0, m, len, cont, disc);
1107       gq->SetPoint(i,m,disc);i++;
1108     }
1109   }
1110   gq->SetMarkerStyle(20);
1111   gq->Draw("pl");
1112
1113   TGraph *gg=new TGraph(20);
1114   i=0;
1115   if(fMultSoft){
1116     for(Double_t q=0.05;q<=5.05;q+=0.25){
1117       Double_t disc,cont;
1118       CalcMult(2,1.0, q, 5., cont, disc);
1119       gg->SetPoint(i,q,disc);i++;
1120     }
1121   } else {
1122     for(Double_t m=0.05;m<=5.05;m+=0.25){
1123       Double_t disc,cont;
1124       CalcSingleHard(2,1.0, m, 5., cont, disc);
1125       gg->SetPoint(i,m,disc);i++;
1126     }
1127   }
1128   gg->SetMarkerStyle(24);
1129   gg->Draw("pl");
1130
1131   TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1132   l1a->SetFillStyle(0);
1133   l1a->SetBorderSize(0);
1134   Char_t label[100];
1135   sprintf(label,"L = %d fm",len);
1136   l1a->AddEntry(gq,label,"");
1137   l1a->AddEntry(gq,"quark","pl");
1138   l1a->AddEntry(gg,"gluon","pl");
1139   l1a->Draw();
1140
1141   c->Update();
1142 }
1143
1144 void AliQuenchingWeights::PlotContWeights(Int_t itype,Int_t ell) const
1145 {
1146   Float_t medvals[3];
1147   Char_t title[1024];
1148   Char_t name[1024];
1149   if(fMultSoft) {
1150     if(itype==1)
1151       sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1152     else if(itype==2)
1153       sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1154     else return;
1155     medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1156     sprintf(name,"ccont-ms-%d",itype);
1157   } else {
1158     if(itype==1)
1159       sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1160     else if(itype==2)
1161       sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1162     else return;
1163     medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1164     sprintf(name,"ccont-ms-%d",itype);
1165   }
1166
1167   TCanvas *c = new TCanvas(name,title,0,0,500,400);
1168   c->cd();
1169   TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1170   h1->SetName("h1");
1171   h1->SetTitle(title);
1172   h1->SetStats(0);
1173   h1->SetLineColor(1);
1174   h1->DrawCopy();
1175   TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1176   h2->SetName("h2");
1177   h2->SetLineColor(2);
1178   h2->DrawCopy("SAME");
1179   TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1180   h3->SetName("h3");
1181   h3->SetLineColor(3);
1182   h3->DrawCopy("SAME");
1183
1184   TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1185   l1a->SetFillStyle(0);
1186   l1a->SetBorderSize(0);
1187   Char_t label[100];
1188   sprintf(label,"L = %d fm",ell);
1189   l1a->AddEntry(h1,label,"");
1190   if(fMultSoft) {
1191     sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1192     l1a->AddEntry(h1,label,"pl");
1193     sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1194     l1a->AddEntry(h2,label,"pl");
1195     sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1196     l1a->AddEntry(h3,label,"pl");
1197   } else {
1198     sprintf(label,"#mu = %.1f GeV",medvals[0]);
1199     l1a->AddEntry(h1,label,"pl");
1200     sprintf(label,"#mu = %.1f GeV",medvals[1]);
1201     l1a->AddEntry(h2,label,"pl");
1202     sprintf(label,"#mu = %.1f GeV",medvals[2]);
1203     l1a->AddEntry(h3,label,"pl");
1204   }
1205   l1a->Draw();
1206
1207   c->Update();
1208 }
1209
1210 void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const
1211 {
1212   Char_t title[1024];
1213   Char_t name[1024];
1214   if(fMultSoft) {
1215     if(itype==1)
1216       sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1217     else if(itype==2)
1218       sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1219     else return;
1220     sprintf(name,"ccont2-ms-%d",itype);
1221   } else {
1222     if(itype==1)
1223       sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1224     else if(itype==2)
1225       sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1226     else return;
1227     sprintf(name,"ccont2-sh-%d",itype);
1228   }
1229   TCanvas *c = new TCanvas(name,title,0,0,500,400);
1230   c->cd();
1231   TH1F *h1=ComputeQWHisto(itype,medval,8);
1232   h1->SetName("h1");
1233   h1->SetTitle(title);
1234   h1->SetStats(0);
1235   h1->SetLineColor(1);
1236   h1->DrawCopy();
1237   TH1F *h2=ComputeQWHisto(itype,medval,5);
1238   h2->SetName("h2");
1239   h2->SetLineColor(2);
1240   h2->DrawCopy("SAME");
1241   TH1F *h3=ComputeQWHisto(itype,medval,2);
1242   h3->SetName("h3");
1243   h3->SetLineColor(3);
1244   h3->DrawCopy("SAME");
1245
1246   TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1247   l1a->SetFillStyle(0);
1248   l1a->SetBorderSize(0);
1249   Char_t label[100];
1250   if(fMultSoft)
1251     sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval);
1252   else
1253     sprintf(label,"#mu = %.1f GeV",medval);
1254
1255   l1a->AddEntry(h1,label,"");
1256   l1a->AddEntry(h1,"L = 8 fm","pl");
1257   l1a->AddEntry(h2,"L = 5 fm","pl");
1258   l1a->AddEntry(h3,"L = 2 fm","pl");
1259   l1a->Draw();
1260
1261   c->Update();
1262 }
1263
1264 void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e)  const
1265 {
1266   if(!fTablesLoaded){
1267     Error("CalcMult","Tables are not loaded.");
1268     return;
1269   }
1270
1271   Char_t title[1024];
1272   Char_t name[1024];
1273   if(fMultSoft){
1274     sprintf(title,"Average Loss for Multiple Scattering");
1275     sprintf(name,"cavgelossms");
1276   } else {
1277     sprintf(title,"Average Loss for Single Hard Scattering");
1278     sprintf(name,"cavgelosssh");
1279   }
1280
1281   TCanvas *c = new TCanvas(name,title,0,0,500,400);
1282   c->cd();
1283   TH2F *hframe = new TH2F("avgloss",title,2,0,5.1,2,0,100);
1284   hframe->SetStats(0);
1285   if(fMultSoft)
1286     hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1287   else
1288     hframe->SetXTitle("#mu [GeV]");
1289   hframe->SetYTitle("<E_{loss}> [GeV]");
1290   hframe->Draw();
1291
1292   TGraph *gq=new TGraph(20);
1293   Int_t i=0;
1294   for(Double_t v=0.05;v<=5.05;v+=0.25){
1295     TH1F *dummy=ComputeELossHisto(1,v,len,e);
1296     Double_t avgloss=dummy->GetMean();
1297     gq->SetPoint(i,v,avgloss);i++;
1298     delete dummy;
1299   }
1300   gq->SetMarkerStyle(20);
1301   gq->Draw("pl");
1302
1303   TGraph *gg=new TGraph(20);
1304   i=0;
1305   for(Double_t v=0.05;v<=5.05;v+=0.25){
1306     TH1F *dummy=ComputeELossHisto(2,v,len,e);
1307     Double_t avgloss=dummy->GetMean();
1308     gg->SetPoint(i,v,avgloss);i++;
1309     delete dummy;
1310   }
1311   gg->SetMarkerStyle(24);
1312   gg->Draw("pl");
1313
1314   TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1315   l1a->SetFillStyle(0);
1316   l1a->SetBorderSize(0);
1317   Char_t label[100];
1318   sprintf(label,"L = %d fm",len);
1319   l1a->AddEntry(gq,label,"");
1320   l1a->AddEntry(gq,"quark","pl");
1321   l1a->AddEntry(gg,"gluon","pl");
1322   l1a->Draw();
1323
1324   c->Update();
1325 }
1326
1327 void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
1328 {
1329   if(!fTablesLoaded){
1330     Error("CalcMult","Tables are not loaded.");
1331     return;
1332   }
1333
1334   Char_t title[1024];
1335   Char_t name[1024];
1336   if(fMultSoft){
1337     sprintf(title,"Average Loss for Multiple Scattering");
1338     sprintf(name,"cavgelossms2");
1339   } else {
1340     sprintf(title,"Average Loss for Single Hard Scattering");
1341     sprintf(name,"cavgelosssh2");
1342   }
1343
1344   TCanvas *c = new TCanvas(name,title,0,0,500,400);
1345   c->cd();
1346   TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1347   hframe->SetStats(0);
1348   if(fMultSoft)
1349     hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1350   else
1351     hframe->SetXTitle("#mu [GeV]");
1352   hframe->SetYTitle("<E_{loss}> [GeV]");
1353   hframe->Draw();
1354
1355   TGraph *gq=new TGraph(20);
1356   Int_t i=0;
1357   for(Double_t v=0.05;v<=5.05;v+=0.25){
1358     TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1359     Double_t avgloss=dummy->GetMean();
1360     gq->SetPoint(i,v,avgloss);i++;
1361     delete dummy;
1362   }
1363   gq->SetMarkerStyle(20);
1364   gq->Draw("pl");
1365
1366   TGraph *gg=new TGraph(20);
1367   i=0;
1368   for(Double_t v=0.05;v<=5.05;v+=0.25){
1369     TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1370     Double_t avgloss=dummy->GetMean();
1371     gg->SetPoint(i,v,avgloss);i++;
1372     delete dummy;
1373   }
1374   gg->SetMarkerStyle(24);
1375   gg->Draw("pl");
1376
1377   TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1378   l1a->SetFillStyle(0);
1379   l1a->SetBorderSize(0);
1380   Char_t label[100];
1381   sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1382   l1a->AddEntry(gq,label,"");
1383   l1a->AddEntry(gq,"quark","pl");
1384   l1a->AddEntry(gg,"gluon","pl");
1385   l1a->Draw();
1386
1387   c->Update();
1388 }
1389
1390 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Int_t len) const
1391 {
1392   if(!fTablesLoaded){
1393     Error("CalcMult","Tables are not loaded.");
1394     return;
1395   }
1396
1397   Char_t title[1024];
1398   Char_t name[1024];
1399   if(fMultSoft){
1400     sprintf(title,"Relative Loss for Multiple Scattering");
1401     sprintf(name,"cavgelossvsptms");
1402   } else {
1403     sprintf(title,"Relative Loss for Single Hard Scattering");
1404     sprintf(name,"cavgelossvsptsh");
1405   }
1406
1407   TCanvas *c = new TCanvas(name,title,0,0,500,400);
1408   c->cd();
1409   TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
1410   hframe->SetStats(0);
1411   hframe->SetXTitle("p_{T} [GeV]");
1412   hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
1413   hframe->Draw();
1414
1415   TGraph *gq=new TGraph(40);
1416   Int_t i=0;
1417   for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1418     TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
1419     Double_t avgloss=dummy->GetMean();
1420     gq->SetPoint(i,pt,avgloss/pt);i++;
1421     delete dummy;
1422   }
1423   gq->SetMarkerStyle(20);
1424   gq->Draw("pl");
1425
1426   TGraph *gg=new TGraph(40);
1427   i=0;
1428   for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1429     TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
1430     Double_t avgloss=dummy->GetMean();
1431     gg->SetPoint(i,pt,avgloss/pt);i++;
1432     delete dummy;
1433   }
1434   gg->SetMarkerStyle(24);
1435   gg->Draw("pl");
1436
1437   TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1438   l1a->SetFillStyle(0);
1439   l1a->SetBorderSize(0);
1440   Char_t label[100];
1441   sprintf(label,"L = %d fm",len);
1442   l1a->AddEntry(gq,label,"");
1443   l1a->AddEntry(gq,"quark","pl");
1444   l1a->AddEntry(gg,"gluon","pl");
1445   l1a->Draw();
1446
1447   c->Update();
1448 }
1449
1450 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
1451 {
1452   if(!fTablesLoaded){
1453     Error("CalcMult","Tables are not loaded.");
1454     return;
1455   }
1456
1457   Char_t title[1024];
1458   Char_t name[1024];
1459   if(fMultSoft){
1460     sprintf(title,"Relative Loss for Multiple Scattering");
1461     sprintf(name,"cavgelossvsptms2");
1462   } else {
1463     sprintf(title,"Relative Loss for Single Hard Scattering");
1464     sprintf(name,"cavgelossvsptsh2");
1465   }
1466   TCanvas *c = new TCanvas(name,title,0,0,500,400);
1467   c->cd();
1468   TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
1469   hframe->SetStats(0);
1470   hframe->SetXTitle("p_{T} [GeV]");
1471   hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
1472   hframe->Draw();
1473
1474   TGraph *gq=new TGraph(40);
1475   Int_t i=0;
1476   for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1477     TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
1478     Double_t avgloss=dummy->GetMean();
1479     gq->SetPoint(i,pt,avgloss/pt);i++;
1480     delete dummy;
1481   }
1482   gq->SetMarkerStyle(20);
1483   gq->Draw("pl");
1484
1485   TGraph *gg=new TGraph(40);
1486   i=0;
1487   for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1488     TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
1489     Double_t avgloss=dummy->GetMean();
1490     gg->SetPoint(i,pt,avgloss/pt);i++;
1491     delete dummy;
1492   }
1493   gg->SetMarkerStyle(24);
1494   gg->Draw("pl");
1495
1496   TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1497   l1a->SetFillStyle(0);
1498   l1a->SetBorderSize(0);
1499   Char_t label[100];
1500   sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1501   l1a->AddEntry(gq,label,"");
1502   l1a->AddEntry(gq,"quark","pl");
1503   l1a->AddEntry(gg,"gluon","pl");
1504   l1a->Draw();
1505
1506   c->Update();
1507 }
1508