CMake: Changing DA detector name to DA0 to match DAQ settings
[u/mrichter/AliRoot.git] / ZDC / Fragment.C
1 void Fragment(Int_t nev=0, Int_t debug=0)
2 {
3  gROOT->Reset();
4  
5  Float_t b;
6  Int_t nspectp, nspectn, nspectpfree, nspectnfree; 
7  Int_t zz[100], nn[100], nalpha, NFrag, ztot, ntot, ndeu;
8  void spectator(Float_t, Int_t*, Int_t*);
9  
10  TH2F *hsp = new TH2F("hsp","Spectator protons vs b",50,0.,20.,90,0.,90.);
11  hsp -> SetXTitle("b (fm)");
12  hsp -> SetYTitle("Num. of spectator protons");
13  
14  TH2F *hsn = new TH2F("hsn","Spectator neutrons vs b",50,0.,20.,130,0.,130.);
15  hsn -> SetXTitle("b (fm)");
16  hsn -> SetYTitle("Num. of spectator neutrons");
17  
18  TH2F *hFragp = new TH2F("hFragp","Free spectator protons vs b",50,0.,20.,60,0.,60.);
19  hFragp -> SetXTitle("b (fm)");
20  hFragp -> SetYTitle("Num. of free spectator protons");
21  
22  TH2F *hFragn = new TH2F("hFragn","Free spectator neutrons vs b",50,0.,20.,80,0.,80.);
23  hFragn -> SetXTitle("b (fm)");
24  hFragn -> SetYTitle("Num. of free spectator neutrons");
25
26  TF1 *funb = new  TF1("funb","x",0,20);
27  for(Int_t ievent=0; ievent<nev; ievent++){  
28    if(ievent%1000==0){printf("Processing event %d\n",ievent);}
29    b= Float_t(funb->GetRandom()); 
30    spectator(b, &nspectp, &nspectn);
31    hsp -> Fill(b,nspectp);
32    hsn -> Fill(b,nspectn);
33    AliZDCFragment *gallio = new AliZDCFragment(b);
34    for(Int_t j=0; j<=99; j++){
35       zz[j] =0;
36       nn[j] =0;
37    }
38 //
39 // Generation of fragments
40    gallio->GenerateIMF();
41    nalpha = gallio->GetNalpha();
42    NFrag = gallio->GetFragmentNum();
43 //
44    // Attach neutrons
45    ztot = gallio->GetZtot();
46    ntot = gallio->GetNtot();
47    gallio->AttachNeutrons();
48    //
49    nspectpfree = (nspectp-ztot-2*nalpha);
50    nspectnfree = (nspectn-ntot-2*nalpha);
51    
52    // Removing deuterons
53    ndeu = (Int_t) (nspectnfree*gallio->DeuteronNumber());
54    nspectpfree -= ndeu;
55    nspectnfree -= ndeu;
56    //
57    hFragp -> Fill(b, nspectpfree);
58    hFragn -> Fill(b, nspectnfree);
59 //
60 // Print
61    if(debug == 1){
62      printf("\n b = %f",b);
63      printf("   #spect p = %d, #spect n = %d\n",nspectp,nspectn); 
64      printf("   #spect p free = %d, #spect n free = %d\n",nspectpfree,nspectnfree); 
65      printf("   #fragments = %f ",NFrag);
66      /*for(Int_t i=0; i<NFrag; i++){
67         printf("\n      ZZ[%d] = %d, NN[%d] = %d\n",i,zz[i],i,nn[i]);
68      }*/
69      printf("   NAlpha = %d, Ztot = %d, Ntot = %d\n\n",nalpha,ztot,ntot);
70    }
71    delete gallio;
72   } //Event loop
73   
74   TProfile *profsp = hsp->ProfileX("profsp",-1,-1,"s");
75   TProfile *profsn = hsn->ProfileX("profsn",-1,-1,"s");
76   TProfile *profFragp = hFragp->ProfileX("profFragp",-1,-1,"s");
77   TProfile *profFragn = hFragn->ProfileX("profFragn",-1,-1,"s");
78
79    
80   //***********************************************************
81   // #### ROOT initialization
82   gROOT->Reset();
83   gStyle->SetCanvasColor(10);
84   gStyle->SetFrameFillColor(10);
85   gStyle->SetPalette(1);
86   gStyle->SetOptStat(0);
87   //
88   TCanvas *c1 = new TCanvas("c1","Fragmentation I",0,0,700,700);
89   c1->cd();
90   c1->SetFillColor(kAzure+7);
91   c1->Divide(2,2);
92   c1->cd(1);
93   hsp -> Draw("colz");   
94   c1->cd(2);
95   hsn -> Draw("colz");    
96   c1->cd(3);
97   hFragp -> Draw("colz");   
98   c1->cd(4);
99   hFragn -> Draw("colz"); 
100   //
101   c1->Print("Fragment1.gif");   
102   //
103   TCanvas *c2 = new TCanvas("c2","Fragmentation II",300,10,700,700);
104   c2->cd();
105   c2->SetFillColor(kAzure+10);
106   c2->Divide(2,2);
107   c2->cd(1);
108   profsp ->SetLineColor(kRed);  profsp ->SetLineWidth(2);
109   profsp -> Draw("colz");   
110   c2->cd(2);
111   profsn ->SetLineColor(kRed);  profsn ->SetLineWidth(2);
112   profsn -> Draw("colz");    
113   c2->cd(3);
114   profFragp -> Draw("colz");   
115   profFragp ->SetLineColor(kAzure);  profFragp ->SetLineWidth(2);
116   c2->cd(4);
117   profFragn -> Draw("colz");    
118   profFragn ->SetLineColor(kAzure);  profFragn ->SetLineWidth(2);
119   //
120   c1->Print("Fragment2.gif");   
121
122 }
123
124
125 void spectator(Float_t b, Int_t* NSpectp, Int_t* NSpectn)
126 {
127    Float_t SppvsB[6] = {3.633,-1.518,1.360,-.4899e-1,-.2398e-2,.1066e-3};
128    Float_t SpnvsB[6] = {5.639,-1.685,1.803,-.3129e-1,-.6618e-2,.2352e-3};
129    Float_t Sigmap[4] = {.5668,-.2200e-1,.3657e-3,-.2201e-5};
130    Float_t Sigman[4] = {.4185,-.9798e-2,.1052e-3,-.4238e-6};
131    
132    Float_t rnsp = SppvsB[0]+SppvsB[1]*b+SppvsB[2]*(b*b)+SppvsB[3]*(b*b*b)+
133                   SppvsB[4]*(b*b*b*b)+SppvsB[5]*(b*b*b*b*b);
134    Float_t rnsn = SpnvsB[0]+SpnvsB[1]*b+SpnvsB[2]*(b*b)+SpnvsB[3]*(b*b*b)+
135                   SpnvsB[4]*(b*b*b*b)+SpnvsB[5]*(b*b*b*b*b);
136    Float_t snsp = Sigmap[0]+Sigmap[1]*rnsp+Sigmap[2]*(rnsp*rnsp)+Sigmap[3]*
137                   (rnsp*rnsp*rnsp); 
138    Float_t snsn = Sigman[0]+Sigman[1]*rnsn+Sigman[2]*(rnsn*rnsn)+Sigman[3]*
139                   (rnsn*rnsn*rnsn); 
140    
141    snsp = snsp*rnsp;
142    snsn = snsn*rnsn;
143    
144    Float_t xgaup = gRandom->Gaus(0.0,1.0);
145    snsp = snsp*xgaup;
146    Float_t xgaun = gRandom->Gaus(0.0,1.0);
147    snsn = snsn*xgaun;
148    rnsp=rnsp+snsp;
149    rnsn=rnsn+snsn;
150    
151    *NSpectp = Int_t(rnsp);
152    *NSpectn = Int_t(rnsn);
153
154 }