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, ztot, ntot;
   8  void spectator(Float_t, Int_t*, Int_t*);
   9
  10  TH2F *hsp = new TH2F("hsp","Spectator protons vs b",100,0.,20.,100,0.,150.);
  11  hsp -> SetXTitle("b (fm)");
  12  hsp -> SetYTitle("Num. of spectator protons");
  13
  14  TH2F *hsn = new TH2F("hsn","Spectator neutrons vs b",100,0.,20.,100,0.,150.);
  15  hsn -> SetXTitle("b (fm)");
  16  hsn -> SetYTitle("Num. of spectator neutrons");
  17
  18  TH2F *hFragp = new TH2F("hFragp","Free spectator protons vs b",100,0.,20.,100,0.,150.);
  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",100,0.,20.,100,0.,150.);
  23  hFragn -> SetXTitle("b (fm)");
  24  hFragn -> SetYTitle("Num. of free spectator neutrons");
  25
  26  TF1 *funb = new  TF1("funb","x",0,15);
  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    Int_t NFrag = gallio->GetFragmentNum();
  42 //
  43 // Attach neutrons
  44    ztot=0;
  45    ntot=0;
  46    gallio->AttachNeutrons();
  47    nspectpfree = nspectp-ztot-2*nalpha;
  48    nspectnfree = nspectn-ntot-2*nalpha;
  49    hFragp -> Fill(b,nspectpfree);
  50    hFragn -> Fill(b,nspectnfree);
  51 //
  52 // Print
  53    if(debug == 1){
  54      printf("\n b = %f\n",b);
  55      printf("\n #spect p = %d, #spect n = %d\n",nspectp,nspectn);
  56      printf("\n #spect p free = %d, #spect n free = %d\n",nspectpfree,nspectnfree);
  57      printf("\n #fragments = %f\n",NFrag);
  58      for(Int_t i=0; i<NFrag; i++){
  59         printf("\n      ZZ[%d] = %d, NN[%d] = %d\n",i,zz[i],i,nn[i]);
  60      }
  61      printf("\n NAlpha = %d, Ztot = %d, Ntot = %d\n\n",nalpha,ztot,ntot);
  62    }
  63    delete gallio;
  64   } //Event loop
  65
  66    TCanvas *c1 = new TCanvas("c1","Fragmentation",0,10,580,700);
  67    c1->cd();
  68    c1->SetFillColor(29);
  69    c1->Divide(2,2);
  70    c1->cd(1);
  71    hsp -> Draw("box");
  72    c1->cd(2);
  73    hsn -> Draw("box");
  74    c1->cd(3);
  75    hFragp -> Draw("box");
  76    c1->cd(4);
  77    hFragn -> Draw("box");
  78
  79 }
  80
  81
  82 void spectator(Float_t b, Int_t* NSpectp, Int_t* NSpectn)
  83 {
  84    Float_t SppvsB[6] = {3.633,-1.518,1.360,-.4899e-1,-.2398e-2,.1066e-3};
  85    Float_t SpnvsB[6] = {5.639,-1.685,1.803,-.3129e-1,-.6618e-2,.2352e-3};
  86    Float_t Sigmap[4] = {.5668,-.2200e-1,.3657e-3,-.2201e-5};
  87    Float_t Sigman[4] = {.4185,-.9798e-2,.1052e-3,-.4238e-6};
  88
  89    Float_t rnsp = SppvsB[0]+SppvsB[1]*b+SppvsB[2]*(b*b)+SppvsB[3]*(b*b*b)+
  90                   SppvsB[4]*(b*b*b*b)+SppvsB[5]*(b*b*b*b*b);
  91    Float_t rnsn = SpnvsB[0]+SpnvsB[1]*b+SpnvsB[2]*(b*b)+SpnvsB[3]*(b*b*b)+
  92                   SpnvsB[4]*(b*b*b*b)+SpnvsB[5]*(b*b*b*b*b);
  93    Float_t snsp = Sigmap[0]+Sigmap[1]*rnsp+Sigmap[2]*(rnsp*rnsp)+Sigmap[3]*
  94                   (rnsp*rnsp*rnsp);
  95    Float_t snsn = Sigman[0]+Sigman[1]*rnsn+Sigman[2]*(rnsn*rnsn)+Sigman[3]*
  96                   (rnsn*rnsn*rnsn);
  97
  98    snsp = snsp*rnsp;
  99    snsn = snsn*rnsn;
 100
 101    Float_t xgaup = gRandom->Gaus(0.0,1.0);
 102    snsp = snsp*xgaup;
 103    Float_t xgaun = gRandom->Gaus(0.0,1.0);
 104    snsn = snsn*xgaun;
 105    rnsp=rnsp+snsp;
 106    rnsn=rnsn+snsn;
 107
 108    *NSpectp = Int_t(rnsp);
 109    *NSpectn = Int_t(rnsn);
 110
 111 }