Updated macro
authorcoppedis <coppedis@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Sep 2010 07:43:55 +0000 (07:43 +0000)
committercoppedis <coppedis@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Sep 2010 07:43:55 +0000 (07:43 +0000)
ZDC/Fragment.C

index 9661e41..9963f14 100644 (file)
@@ -4,26 +4,26 @@ void Fragment(Int_t nev=0, Int_t debug=0)
  
  Float_t b;
  Int_t nspectp, nspectn, nspectpfree, nspectnfree; 
- Int_t zz[100], nn[100], nalpha, ztot, ntot;
+ Int_t zz[100], nn[100], nalpha, NFrag, ztot, ntot, ndeu;
  void spectator(Float_t, Int_t*, Int_t*);
  
- TH2F *hsp = new TH2F("hsp","Spectator protons vs b",100,0.,20.,100,0.,150.);
+ TH2F *hsp = new TH2F("hsp","Spectator protons vs b",50,0.,20.,90,0.,90.);
  hsp -> SetXTitle("b (fm)");
  hsp -> SetYTitle("Num. of spectator protons");
  
- TH2F *hsn = new TH2F("hsn","Spectator neutrons vs b",100,0.,20.,100,0.,150.);
+ TH2F *hsn = new TH2F("hsn","Spectator neutrons vs b",50,0.,20.,130,0.,130.);
  hsn -> SetXTitle("b (fm)");
  hsn -> SetYTitle("Num. of spectator neutrons");
  
- TH2F *hFragp = new TH2F("hFragp","Free spectator protons vs b",100,0.,20.,100,0.,150.);
+ TH2F *hFragp = new TH2F("hFragp","Free spectator protons vs b",50,0.,20.,60,0.,60.);
  hFragp -> SetXTitle("b (fm)");
  hFragp -> SetYTitle("Num. of free spectator protons");
  
- TH2F *hFragn = new TH2F("hFragn","Free spectator neutrons vs b",100,0.,20.,100,0.,150.);
+ TH2F *hFragn = new TH2F("hFragn","Free spectator neutrons vs b",50,0.,20.,80,0.,80.);
  hFragn -> SetXTitle("b (fm)");
  hFragn -> SetYTitle("Num. of free spectator neutrons");
 
- TF1 *funb = new  TF1("funb","x",0,15);
+ TF1 *funb = new  TF1("funb","x",0,20);
  for(Int_t ievent=0; ievent<nev; ievent++){  
    if(ievent%1000==0){printf("Processing event %d\n",ievent);}
    b= Float_t(funb->GetRandom()); 
@@ -38,43 +38,86 @@ void Fragment(Int_t nev=0, Int_t debug=0)
 //
 // Generation of fragments
    gallio->GenerateIMF();
-   Int_t NFrag = gallio->GetFragmentNum();
+   nalpha = gallio->GetNalpha();
+   NFrag = gallio->GetFragmentNum();
 //
-// Attach neutrons
-   ztot=0;
-   ntot=0;
+   // Attach neutrons
+   ztot = gallio->GetZtot();
+   ntot = gallio->GetNtot();
    gallio->AttachNeutrons();
-   nspectpfree = nspectp-ztot-2*nalpha;
-   nspectnfree = nspectn-ntot-2*nalpha;
-   hFragp -> Fill(b,nspectpfree);
-   hFragn -> Fill(b,nspectnfree);
+   //
+   nspectpfree = (nspectp-ztot-2*nalpha);
+   nspectnfree = (nspectn-ntot-2*nalpha);
+   
+   // Removing deuterons
+   ndeu = (Int_t) (nspectnfree*gallio->DeuteronNumber());
+   nspectpfree -= ndeu;
+   nspectnfree -= ndeu;
+   //
+   hFragp -> Fill(b, nspectpfree);
+   hFragn -> Fill(b, nspectnfree);
 //
 // Print
    if(debug == 1){
-     printf("\n        b = %f\n",b);
-     printf("\n        #spect p = %d, #spect n = %d\n",nspectp,nspectn); 
-     printf("\n        #spect p free = %d, #spect n free = %d\n",nspectpfree,nspectnfree); 
-     printf("\n        #fragments = %f\n",NFrag);
-     for(Int_t i=0; i<NFrag; i++){
+     printf("\n        b = %f",b);
+     printf("  #spect p = %d, #spect n = %d\n",nspectp,nspectn); 
+     printf("  #spect p free = %d, #spect n free = %d\n",nspectpfree,nspectnfree); 
+     printf("  #fragments = %f ",NFrag);
+     /*for(Int_t i=0; i<NFrag; i++){
        printf("\n      ZZ[%d] = %d, NN[%d] = %d\n",i,zz[i],i,nn[i]);
-     }
-     printf("\n        NAlpha = %d, Ztot = %d, Ntot = %d\n\n",nalpha,ztot,ntot);
+     }*/
+     printf("  NAlpha = %d, Ztot = %d, Ntot = %d\n\n",nalpha,ztot,ntot);
    }
    delete gallio;
   } //Event loop
+  
+  TProfile *profsp = hsp->ProfileX("profsp",-1,-1,"s");
+  TProfile *profsn = hsn->ProfileX("profsn",-1,-1,"s");
+  TProfile *profFragp = hFragp->ProfileX("profFragp",-1,-1,"s");
+  TProfile *profFragn = hFragn->ProfileX("profFragn",-1,-1,"s");
 
-   TCanvas *c1 = new TCanvas("c1","Fragmentation",0,10,580,700);
-   c1->cd();
-   c1->SetFillColor(29);
-   c1->Divide(2,2);
-   c1->cd(1);
-   hsp -> Draw("box");   
-   c1->cd(2);
-   hsn -> Draw("box");    
-   c1->cd(3);
-   hFragp -> Draw("box");   
-   c1->cd(4);
-   hFragn -> Draw("box");    
+   
+  //***********************************************************
+  // #### ROOT initialization
+  gROOT->Reset();
+  gStyle->SetCanvasColor(10);
+  gStyle->SetFrameFillColor(10);
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(0);
+  //
+  TCanvas *c1 = new TCanvas("c1","Fragmentation I",0,0,700,700);
+  c1->cd();
+  c1->SetFillColor(kAzure+7);
+  c1->Divide(2,2);
+  c1->cd(1);
+  hsp -> Draw("colz");   
+  c1->cd(2);
+  hsn -> Draw("colz");    
+  c1->cd(3);
+  hFragp -> Draw("colz");   
+  c1->cd(4);
+  hFragn -> Draw("colz"); 
+  //
+  c1->Print("Fragment1.gif");   
+  //
+  TCanvas *c2 = new TCanvas("c2","Fragmentation II",300,10,700,700);
+  c2->cd();
+  c2->SetFillColor(kAzure+10);
+  c2->Divide(2,2);
+  c2->cd(1);
+  profsp ->SetLineColor(kRed);  profsp ->SetLineWidth(2);
+  profsp -> Draw("colz");   
+  c2->cd(2);
+  profsn ->SetLineColor(kRed);  profsn ->SetLineWidth(2);
+  profsn -> Draw("colz");    
+  c2->cd(3);
+  profFragp -> Draw("colz");   
+  profFragp ->SetLineColor(kAzure);  profFragp ->SetLineWidth(2);
+  c2->cd(4);
+  profFragn -> Draw("colz");    
+  profFragn ->SetLineColor(kAzure);  profFragn ->SetLineWidth(2);
+  //
+  c1->Print("Fragment2.gif");   
 
 }