First commit.

diff --git a/LHC/multi.C b/LHC/multi.C
new file mode 100644 (file)
index 0000000..f2f3796
--- /dev/null
+++ b/LHC/multi.C
@@ -0,0 +1,187 @@
+Double_t fs1(Double_t *x, Double_t* dum);
+Double_t fs2(Double_t *x, Double_t* dum);
+Double_t fs3(Double_t *x, Double_t* dum);
+
+static Int_t n;        // array size
+static Double_t* as1;  // dsigma/db
+static Double_t* as2;  // <N>(b)
+static Double_t* ab;   // b
+// |  PbPb           |    geom.   |      hard
+// |  0   -  5   fm |   0 - 10 % |   0    -  42.8 %
+// |  5   -  8.6 fm |  10 - 30 % |  42.8  -  83.5 %
+// |  8.6 - 11.2 fm |  30 - 50 % |  83.5  -  96.8 %
+// | 11.2 - 13.2 fm |  50 - 70 % |  96.8  -  99.6 %
+// | 13.2 - 15.0 fm |  70 - 90 % |  99.6  -  99.9 %
+// | 15.0 -  oo  fm |  90 -100 % |  99.9  - 100.0 %
+ 
+void multi()
+{
+    const Float_t drift = 90;
+    
+//  open file and get graph objects
+    TFile*  file = new TFile("graph.root", "read");
+    TGraph* gs1 =  gAlice = (TGraph*)(file->Get("Graph;1"));
+    TGraph* gs2 =  gAlice = (TGraph*)(file->Get("Graph;2"));
+// initialize arrays
+    n   = gs1->GetN();
+    ab  = gs1->GetX();
+    as1 = gs1->GetY();
+    as2 = gs2->GetY();
+//
+// function
+//
+// dsigma_mb/db 
+    TF1* tf1 = new TF1("tf1", fs1, 0., ab[n-1], 0);
+// dsigma_hard/db
+    TF1* tf2 = new TF1("tf2", fs3, 0., ab[n-1], 0);
+
+    Float_t db = 20./100.;
+    for (Int_t i = 0; i< 100; i++) 
+    {
+       Float_t bb = Float_t(i)*db;
+       Float_t s  = tf1->Integral(0,bb);
+       Float_t sh = tf2->Integral(0,bb);
+       printf("\n %d %f %f %f %f" , i, bb, s, s/7.769*100., sh/4.327*100.);
+    }
+    
+    
+// histogram for random distribution
+    TH1F* hs1 = new TH1F("hs1", "dSigma/db", n, 0,ab[n-1]);
+    hs1->Eval(tf1);
+//    hs1->Draw();
+    TH1F* hmult  = new TH1F("hmult",  "Multiplicity", 100, 0, 10000);
+    TH1F* hmulte = new TH1F("hmulte", "Extra Multiplicity", 100, 0, 10000);
+//     
+    Float_t t1, t2;
+    Double_t b, mult, dm, multe;
+    Int_t i;
+    Double_t *par;
+    Double_t scale = 8000./as2[0];
+    Double_t dT    = 125./2.;
+    
+    for (Int_t i=0; i < 1000000; i++)
+    {
+       // impact parameters
+       b = hs1->GetRandom();
+       // multiplicity
+        mult = scale*(Float_t) fs2(&b,par);
+       mult += gRandom->Gaus(0, TMath::Sqrt(mult));
+       
+       hmult->Fill(Float_t(mult));
+       t1 = 0;
+       t2 = 0;
+       multe = 0;
+       Double_t arg;
+/*     
+       while (1) {
+           arg=0.;
+           while(arg==0.) arg = gRandom->Rndm();
+           t1 -= dT*TMath::Log(arg); // (musec)   
+           if (t1 > drift) break;
+           b = hs1->GetRandom();
+           mult = scale*(Float_t) fs2(&b,par);
+           mult += gRandom->Gaus(0, TMath::Sqrt(mult));
+//         multe+=mult*(drift-t1)/drift;
+           multe+=mult;
+       }
+
+       while (1) {
+           arg=0.;
+           while(arg==0.) arg = gRandom->Rndm();
+           t2 -= dT*TMath::Log(arg); // (musec)   
+           if (t2 > drift) break;
+           b = hs1->GetRandom();
+           mult = scale*(Float_t) fs2(&b,par);
+           mult += gRandom->Gaus(0, TMath::Sqrt(mult));
+//         multe+=mult*(drift-t2)/drift;
+           multe+=mult;
+       }
+*/
+       hmulte->Fill(Float_t(multe));
+    }
+    Double_t scale = 1./hmult->Integral();
+    hmult->Scale(scale/100.);
+    printf("\n scale %f \n ", scale);
+    
+    
+    TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
+//    c1->Divide(1,2);
+//    c1->cd(1);
+    hmult->Draw();
+//    c1->cd(2);
+//   hmulte->Draw();
+    
+    
+    TCanvas *c2 = new TCanvas("c2","Canvas 2",400,10,600,700);
+    c2->Divide(1,2);
+    c2->cd(1);
+    tf1->Draw();
+    tf1->GetHistogram()->SetXTitle("b (fm)");
+    tf1->GetHistogram()->SetYTitle("d#sigma /db (barn/fm)");    
+    c2->cd(2);
+    tf2->Draw();
+    tf2->GetHistogram()->SetXTitle("b (fm)");
+    tf2->GetHistogram()->SetYTitle("d#sigma /db (a.u.)");    
+    c2->Update();
+    out = new TFile("mult.root", "recreate");
+    hmult->Write();
+    out->Close();
+    
+}
+
+
+Double_t fs1(Double_t* x, Double_t* dum)
+{
+    Double_t xx = x[0];
+    Int_t i;
+    Double_t y;
+    
+    for (i=0; i<n; i++) {
+       
+       if (xx <ab[i]) {
+           y = as1[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as1[i]-as1[i-1]);
+           break;
+       }
+    }
+    return y;
+}
+
+
+Double_t fs2(Double_t *x, Double_t* dum)
+{
+    Double_t xx = x[0];
+    Int_t i;
+    for (i=0; i<n; i++) {
+       if (xx <ab[i]) {
+           Double_t y = as2[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as2[i]-as2[i-1]);
+           break;
+       }
+    }
+    return y;
+}
+
+Double_t fs3(Double_t *x, Double_t* dum)
+{
+    Double_t xx = x[0];
+    Int_t i;
+    Double_t y1, y2;
+    
+    for (i=0; i<n; i++) {
+       if (xx <ab[i]) {
+           y1 = as1[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as1[i]-as1[i-1]);
+           break;
+       }
+    }
+
+    for (i=0; i<n; i++) {
+       if (xx <ab[i]) {
+           Double_t y2 = as2[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as2[i]-as2[i-1]);
+           break;
+       }
+    }
+
+
+    return y1*y2;
+}
+
+