FASTSIM/AliQuenchingWeights.cxx

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