LHC/multi.C

   1 Double_t fs1(Double_t *x, Double_t* dum);
   2 Double_t fs2(Double_t *x, Double_t* dum);
   3 Double_t fs3(Double_t *x, Double_t* dum);
   4
   5 static Int_t n;        // array size
   6 static Double_t* as1;  // dsigma/db
   7 static Double_t* as2;  // <N>(b)
   8 static Double_t* ab;   // b
   9 // |  PbPb           |    geom.   |      hard
  10 // |  0   -  5   fm |   0 - 10 % |   0    -  42.8 %
  11 // |  5   -  8.6 fm |  10 - 30 % |  42.8  -  83.5 %
  12 // |  8.6 - 11.2 fm |  30 - 50 % |  83.5  -  96.8 %
  13 // | 11.2 - 13.2 fm |  50 - 70 % |  96.8  -  99.6 %
  14 // | 13.2 - 15.0 fm |  70 - 90 % |  99.6  -  99.9 %
  15 // | 15.0 -  oo  fm |  90 -100 % |  99.9  - 100.0 %
  16
  17 void multi()
  18 {
  19     const Float_t drift = 90;
  20
  21 //  open file and get graph objects
  22     TFile*  file = new TFile("DsigmaDb.root", "read");
  23     TGraph* gs1 =  gAlice = (TGraph*)(file->Get("Graph;1"));
  24     TGraph* gs2 =  gAlice = (TGraph*)(file->Get("Graph;2"));
  25 // initialize arrays
  26     n   = gs1->GetN();
  27     ab  = gs1->GetX();
  28     as1 = gs1->GetY();
  29     as2 = gs2->GetY();
  30 //
  31 // function
  32 //
  33 // dsigma_mb/db
  34     TF1* tf1 = new TF1("tf1", fs1, 0., ab[n-1], 0);
  35 // dsigma_hard/db
  36     TF1* tf2 = new TF1("tf2", fs3, 0., ab[n-1], 0);
  37
  38     Float_t db = 20./100.;
  39     for (Int_t i = 0; i< 100; i++)
  40     {
  41         Float_t bb = Float_t(i)*db;
  42         Float_t s  = tf1->Integral(0,bb);
  43         Float_t sh = tf2->Integral(0,bb);
  44         printf("\n %d %f %f %f %f" , i, bb, s, s/7.769*100., sh/4.327*100.);
  45     }
  46
  47
  48 // histogram for random distribution
  49     TH1F* hs1 = new TH1F("hs1", "dSigma/db", n, 0,ab[n-1]);
  50     hs1->Eval(tf1);
  51 //    hs1->Draw();
  52     TH1F* hmult  = new TH1F("hmult",  "Multiplicity", 100, 0, 10000);
  53     TH1F* hmulte = new TH1F("hmulte", "Extra Multiplicity", 100, 0, 10000);
  54 //
  55     Float_t t1, t2;
  56     Double_t b, mult, dm, multe;
  57     Int_t i;
  58     Double_t *par;
  59     Double_t scale = 8000./as2[0];
  60     Double_t dT    = 125./2.;
  61
  62     for (Int_t i=0; i < 1000000; i++)
  63     {
  64         // impact parameters
  65         b = hs1->GetRandom();
  66         // multiplicity
  67         mult = scale*(Float_t) fs2(&b,par);
  68         mult += gRandom->Gaus(0, TMath::Sqrt(mult));
  69
  70         hmult->Fill(Float_t(mult));
  71         t1 = 0;
  72         t2 = 0;
  73         multe = 0;
  74         Double_t arg;
  75 /*
  76         while (1) {
  77             arg=0.;
  78             while(arg==0.) arg = gRandom->Rndm();
  79             t1 -= dT*TMath::Log(arg); // (musec)
  80             if (t1 > drift) break;
  81             b = hs1->GetRandom();
  82             mult = scale*(Float_t) fs2(&b,par);
  83             mult += gRandom->Gaus(0, TMath::Sqrt(mult));
  84 //          multe+=mult*(drift-t1)/drift;
  85             multe+=mult;
  86         }
  87
  88         while (1) {
  89             arg=0.;
  90             while(arg==0.) arg = gRandom->Rndm();
  91             t2 -= dT*TMath::Log(arg); // (musec)
  92             if (t2 > drift) break;
  93             b = hs1->GetRandom();
  94             mult = scale*(Float_t) fs2(&b,par);
  95             mult += gRandom->Gaus(0, TMath::Sqrt(mult));
  96 //          multe+=mult*(drift-t2)/drift;
  97             multe+=mult;
  98         }
  99 */
 100         hmulte->Fill(Float_t(multe));
 101     }
 102     Double_t scale = 1./hmult->Integral();
 103     hmult->Scale(scale/100.);
 104     printf("\n scale %f \n ", scale);
 105
 106
 107     TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
 108 //    c1->Divide(1,2);
 109 //    c1->cd(1);
 110     hmult->Draw();
 111 //    c1->cd(2);
 112 //   hmulte->Draw();
 113
 114
 115     TCanvas *c2 = new TCanvas("c2","Canvas 2",400,10,600,700);
 116     c2->Divide(1,2);
 117     c2->cd(1);
 118     tf1->Draw();
 119     tf1->GetHistogram()->SetXTitle("b (fm)");
 120     tf1->GetHistogram()->SetYTitle("d#sigma /db (barn/fm)");
 121     c2->cd(2);
 122     tf2->Draw();
 123     tf2->GetHistogram()->SetXTitle("b (fm)");
 124     tf2->GetHistogram()->SetYTitle("d#sigma /db (a.u.)");
 125     c2->Update();
 126     out = new TFile("mult.root", "recreate");
 127     hmult->Write();
 128     out->Close();
 129
 130 }
 131
 132
 133 Double_t fs1(Double_t* x, Double_t* dum)
 134 {
 135     Double_t xx = x[0];
 136     Int_t i;
 137     Double_t y;
 138
 139     for (i=0; i<n; i++) {
 140
 141         if (xx <ab[i]) {
 142             y = as1[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as1[i]-as1[i-1]);
 143             break;
 144         }
 145     }
 146     return y;
 147 }
 148
 149
 150 Double_t fs2(Double_t *x, Double_t* dum)
 151 {
 152     Double_t xx = x[0];
 153     Int_t i;
 154     for (i=0; i<n; i++) {
 155         if (xx <ab[i]) {
 156             Double_t y = as2[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as2[i]-as2[i-1]);
 157             break;
 158         }
 159     }
 160     return y;
 161 }
 162
 163 Double_t fs3(Double_t *x, Double_t* dum)
 164 {
 165     Double_t xx = x[0];
 166     Int_t i;
 167     Double_t y1, y2;
 168
 169     for (i=0; i<n; i++) {
 170         if (xx <ab[i]) {
 171             y1 = as1[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as1[i]-as1[i-1]);
 172             break;
 173         }
 174     }
 175
 176     for (i=0; i<n; i++) {
 177         if (xx <ab[i]) {
 178             Double_t y2 = as2[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as2[i]-as2[i-1]);
 179             break;
 180         }
 181     }
 182
 183
 184     return y1*y2;
 185 }
 186
 187