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 }