// 24 luglio 2007 //---------------------------- // anna.parravicini@cnao.it //---------------------------- /* routine per plottare distribuzione particelle nello spazio fasi ed ellisse di Twiss >>>> Difference with previous versions: *the reference ellipse is plotted together with the RMS ellipse *geometrical emittance plot added *TProfile instead of TH2F, in order to have the average z-entry rather than the sum of z-entries. *x-bins [mm] are centred on the slit median axis width. namely, the bin is drawn from xo-(slit step)/2 to xo+(slit step)/2 *main parameters written on the graphs *wrt v4, aggiungo griglia ai pad e banner per output *wrt v5, bunner spostato a sinistra e scritte molto piu' cose, plottati dati prima del centraggio *wrt v6-7, aggiungo isocurva relative a Ithreshold sul grafico emitt geom. scrivo nel bunner la percentuale di corrente di fascio contenuta nelle ellisse rms e reference *wrt v8, modificato per lavorare con emittGeom_v1 >>>> To be done: * aggiungere logo CNAO * si potrebbero cercare i valori estremi presenti negli array con corrente nn nulla e settare i range dei grafici in modo automatico * decidere come settare il valore dell'errore sull'emittanza geometrica //----------------------------------------------------------*/ #include "TGraph.h" #include "TStyle.h" #include "TMath.h" #include "TCanvas.h" #include "TTree.h" #include "TProfile2D.h" //user's guide, page 44. See prova_TH2F_TProfile2D.C per esempio. int Main(){ // === VARIABLES DECLARATION =================================================== Int_t Nentries; //events number per branch Float_t pos[100000];//[mm] - array da riempire con Branches Float_t div[100000];//[mrad] - array da riempire con Branches Float_t charge[100000];//[arbitrary units] - array da riempire con Branches Float_t pos_t=0;//variabili d'appoggio per leggere Tree Float_t div_t=0; Float_t charge_t=0; Double_t Xmax;//extreme values for x and xp arrays Double_t Xmin; Double_t XPmax; Double_t XPmin; Double_t Imax; Double_t Imin; //RMS ellipse parameters Double_t alpha; Double_t beta; Double_t gamma; Double_t emi; Double_t limit; Int_t n=44; Double_t ell[88]; Double_t x[88]; Double_t step; //particles barycenter Double_t x_center;//[mm] Double_t xp_center;//[mrad] Int_t center; //0=data not centred; 1=data centred //bias and noise thresholds Int_t ThrB; // % of max entry current Int_t ThrZ; // % of max entry current //reference ellipse parameters Double_t alpha_ref; Double_t beta_ref; Double_t gamma_ref; Double_t emi_ref; Double_t limit_ref; Double_t ell_ref[88]; Double_t x_ref[88]; Double_t step_ref; //inside beam percentage Float_t rms_IN; //beam charge percentage inside RMS ellipse Float_t ref_IN;//beam charge percentage inside reference ellipse //geometrical ellipse threshold Double_t Ithreshold; Bool_t overloop; //in case it doesn't reach the required accuracy Float_t P; //we want compute 90% emittance Float_t Peff; //effective percentage of the emittance computation Double_t emiP;//[mm*mrad] geometrical emittance output value Double_t emiRes;//[mm*mrad] geometrical emittance error Double_t maxDiscr;//[mm*mrad] maximum discrepancy among required (P) and computed (Peff) emittance //scanning parameters Int_t sl_Nsteps;//number of slit steps Float_t sl_step;//slit step amplitude Int_t ws_points;// ws samples per profile, at constant velocity //measurement name string name; string outname; string bunner; //general settings gStyle->SetOptStat(kFALSE);//not to plot the statistic box gStyle->SetPalette(1); //gStyle->SetOptDate(31);//rimane nascosta sotto ai pad, non si vede // === CREO UN CANVAS in cui fare i plot e lo divido in 3 ======================= TCanvas *c=new TCanvas ("c", "Phase Space Plots",1200, 800); c->Divide(2,2); // ==== CARICO PUNTI SP FASI DA FILE (Il file deve stare in C:\root) ============= // --- from file to the tree------------------------------- // per fare histo 2D da file con 3 colonne, di cui la terza e' il "peso" // da dare alla cella TTree *t=new TTree ("t", "albero per caricare file di 3 colonne"); t->ReadFile("meas_distribution.txt", "x:xp:I"); Nentries=(Int_t)t->GetEntries(); //numero di eventi per ramo Xmax=t->GetMaximum("x"); Xmin=t->GetMinimum("x"); XPmax=t->GetMaximum ("xp"); XPmin=t->GetMinimum ("xp"); Imax=t->GetMaximum ("I"); Imin=t->GetMinimum ("I"); // --- from tree to the arrays------------------------------ //se mi bastano histo e non voglio grafico, //non e' necessario leggere l'albero e riempire gli array //before reading branches, we must give them an address t->SetBranchAddress("x",&pos_t); t->SetBranchAddress("xp",&div_t); t->SetBranchAddress("I",&charge_t); //fill the array from the branches for (Int_t j=0;jGetEntry(j); pos[j]=pos_t; div[j]=div_t; charge[j]=charge_t; } // === CARICO PARAMETRI DI TWISS e preparo i punti per disegnare l'ellisse ========== ifstream infile("twisspara.txt"); Double_t param[24]; // open the twisspara.txt file to understand the list of parameters written there for (Int_t i=0;i<22;i++) infile>>param[i]; infile>>name; name=name.substr(13); infile>>center; //RMS ellipse parameters alpha=param[0]; beta=param[1]; gamma=(1+alpha*alpha)/beta; emi=param[2]; limit=sqrt(emi*beta); //reference ellipse parameters alpha_ref=param[3]; beta_ref=param[4]; gamma_ref=(1+alpha_ref*alpha_ref)/beta_ref; emi_ref=param[5]; limit_ref=sqrt(emi_ref*beta_ref); //particles barycenter x_center=param[6]; xp_center=param[7]; //bias and noise thresolds Double_t thr; thr=param[18]; ThrB= (int) thr*100; thr=param[19]; ThrZ= (int) thr*100; //beam percentage inside ellipses rms_IN=param[20];//beam current percentage inside rms ellipse ref_IN=param[21];//insidea reference ellipse //scanning settings sl_Nsteps=param[9]; sl_step=param[10]; ws_points=param[11]; //geometrical ellipse Ithreshold=param[8]; overloop=param[12]; if (overloop) cout<<"Geom emitt is far from user's required accuracy."<Fill(pos[j],div[j],charge[j]); // h3d->GetYaxis()->SetRange(ws_points/3,ws_points*2/3); //range di bin da plottare h3d->GetXaxis()->SetTitle("#font[42]{x[mm]}"); h3d->GetYaxis()->SetTitle("#font[42]{xp[mrad]}"); h3d->Draw("lego2z fb bb"); //fb=remove front surface; bb=remove back surface of the cartesian axis // === NOT CENTRED DATA ================================================== c->cd(2); gPad->SetPad(0.34,0.5,0.64,1); gPad->SetGrid(); //gpad e' il puntatore al pad attivo TProfile2D *hnocenter = new TProfile2D("hnocenter","#font[42]{Phase space distribution (before centring data)}", sl_Nsteps,Xmin-sl_step/2+x_center, Xmax+sl_step/2+x_center,ws_points, XPmin+xp_center, XPmax+0.001*XPmax+xp_center); for (int j=0;jFill(pos[j]+x_center,div[j]+xp_center,charge[j]); // hnocenter->GetYaxis()->SetRange(ws_points/3,ws_points*2/3); //range di bin da plottare hnocenter->GetXaxis()->SetTitle("#font[42]{x[mm]}"); hnocenter->GetYaxis()->SetTitle("#font[42]{xp[mrad]}"); gStyle->SetOptStat(kFALSE);//not to plot the statistic box //hnocenter->Draw("cont1z"); //questo opzione riporta le isolinee, senza etichetta di valore hnocenter->Draw("colz"); /*** SCATTER PLOT (not very useful at now) // SCATTER PLOT 3D (non usa i colori, quindi non si capisce molto) TH3F *h3 = new TH3F("h3","#font[42]{Phase space scatter plot}", sl_Nsteps,Xmin, Xmax,ws_points, XPmin, XPmax, 100, Imin, Imax); for (int j=0;jFill(pos[j],div[j],charge[j]); h3->GetYaxis()->SetRange(ws_points/3,ws_points*2/3); //range di bin da plottare h3->GetXaxis()->SetTitle("#font[42]{x[mm]}"); h3->GetYaxis()->SetTitle("#font[42]{xp[mrad]}"); h3->GetZaxis()->SetTitle("#font[42]{I[a.u.]}"); h3->SetMarkerStyle(20);//pag 111 user's guide h3->SetMarkerSize(0.5);//pag 111 user's guide h3->Draw(""); */ // === PHASE SPACE TOP-VIEW WITH RMS AND REFERENCE ELLIPSES ================================================== c->cd(3); gPad->SetPad(0.02,0,0.32,0.5); gPad->SetGrid(); //gpad e' il puntatore al pad attivo //gPad->SetFrameFillColor(33); // histo 2D, with weigth I for each point. // If more points have the same x and xp, their I-values are sum up; it // follows that a I value larger than Imax is displayed: it's due to the binning! TProfile2D *h = new TProfile2D("h","#font[42]{Statistical emittance and ref ellipse}", sl_Nsteps,Xmin-sl_step/2, Xmax+sl_step/2,ws_points, XPmin, XPmax+0.001*XPmax); for (int j=0;jFill(pos[j],div[j],charge[j]); // h->GetYaxis()->SetRange(ws_points/3,ws_points*2/3); //range di bin da plottare h->GetXaxis()->SetTitle("#font[42]{x[mm]}"); h->GetYaxis()->SetTitle("#font[42]{xp[mrad]}"); gStyle->SetOptStat(kFALSE);//not to plot the statistic box h->Draw("colz"); // add the legend leg = new TLegend(0.47,0.11,0.9,0.25); leg->AddEntry(gEll,"RMS ellipse","l");//"l" use the line that markes the gEll graph leg->AddEntry(gEllRef,"Reference ellipse","l"); leg->Draw(); //Add statistical emitt value on the plot TPaveText *ptrms=new TPaveText (0.15,0.81,0.9,0.91,"NDC"); //NDC: normalized coordinates (x1,y1,x2,y2) char text[100]; sprintf (text,"#font[42]{#splitline{alpha RMS= %.3g - beta RMS= %.3g m}{Emitt RMS= %.3g #pi#upointmm#upointmrad}}", alpha, beta,emi); TText *trms=ptrms->AddText(text); trms->SetTextColor(1); trms->SetTextSize(0.04); trms->SetTextAlign(12);//=H and V centred gEll->Draw("C");//non riportiamo gli assi ("A") perche' gia' presenti. gEllRef->Draw("C"); ptrms->Draw(); // === GEOMETRICAL EMITTANCE PLOT ============================================================================ c->cd(4); gPad->SetPad(0.34,0,0.64,0.5); gPad->SetGrid(); // param[8]=Ithreshold, namely, all the bins with current >Ithreshold // contribute to the geometrical emittance, at the required percentage. // I would like a top-view of the 3D phase-space plot, with a different fill color // for bins over and under threshold // Only pints with I >Ithreshold are taken into account char hgeom_name[55]; sprintf(hgeom_name,"#font[42]{Geometrical emittance surface @ %.0f%%}",P); TProfile2D *hgeom = new TProfile2D("hgeom",hgeom_name,sl_Nsteps,Xmin-sl_step/2, Xmax+sl_step/2,ws_points, XPmin, XPmax+0.001*XPmax); for (int j=0;jIthreshold) hgeom->Fill(pos[j],div[j],charge[j]); } //Add geom emitt value on the plot TPaveText *ptgeom=new TPaveText (0.25,0.81,0.9,0.91,"NDC"); //NDC: normalized coordinates (x1,y1,x2,y2) char text[100]; sprintf (text,"#font[42]{#splitline{Geom Emitt=(%.2f#pm%.2f)#pi#upointmm#upointmrad}{(effective percentage: %.2f%%)}}",emiP,emiRes,Peff); TText *tgeom=ptgeom->AddText(text); tgeom->SetTextColor(1); tgeom->SetTextSize(0.04); tgeom->SetTextAlign(12);//=H and V centred //Add warning in case of bad geom emitt computation TPaveText *ptwarn=new TPaveText (0.15,0.15,0.7,0.3,"NDC"); //NDC: normalized coordinates (x1,y1,x2,y2) ptwarn->SetBorderSize(1); char text[100]; sprintf (text,"#font[102]{#splitline{Geom emitt accuracy far}{from user request}}"); TText *twarn=ptwarn->AddText(text); twarn->SetTextColor(2); twarn->SetTextSize(0.04); twarn->SetTextAlign(12);//=H and V centred //Add contour relative to Ithreshold //define contours value Double_t contours[1]; contours[0]=Ithreshold; //firstly, we must plot the histo to extract the contour hgeom->SetContour(1, contours);//se funziona, mettere qui direttamente Ithreshold hgeom->Draw("cont z list"); //questo disegna le isolinee, non le numera pero' c->Update();//e' necessario plottare il grafico per estrarre le line di livello //get contours TObjArray *cont_number = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours"); TList *contLevel = NULL; TGraph * curv = NULL; Int_t nGraphs =0; Int_t TotalConts =0; if (cont_number==NULL){ printf ("no contours extracted \n"); TotalConts=0; return; } else { TotalConts=cont_number->GetSize(); //TotalConts e' il numero di "contorni/isolinee" che trova nell'ultimo histo disegnato } printf ("Total contours =%d \n", TotalConts); for (int j=0; jAt(j); //contLevel e' la lista di grafici (linee disgiunte) che formano il j-th contorno (isolinea) printf("Contour %d has %d Graphs\n", j, contLevel->GetSize()); // contLevel->GetSize() dice il numero di grafici (linee disgiunte) che formano il j-th contorno } nGraphs=0; //TCanvas *cprova =new TCanvas ("c1", "Contour List", 610,0,600,600); TProfile2D *hgeom2 = new TProfile2D("hgeom",hgeom_name,sl_Nsteps,Xmin-sl_step/2, Xmax+sl_step/2,ws_points, XPmin, XPmax+0.001*XPmax); for (int j=0;jIthreshold) hgeom2->Fill(pos[j],div[j],charge[j]); } hgeom2->Draw("colz"); //disegna il grafico su cui saranno plottate le isolinee // hgeom2->GetYaxis()->SetRange(ws_points/3,ws_points*2/3); //range di bin da plottare hgeom2->GetXaxis()->SetTitle("#font[42]{x[mm]}"); hgeom2->GetYaxis()->SetTitle("#font[42]{xp[mrad]}"); //write isocurves labels TLatex la; la.SetTextSize(0.035); la.SetTextFont(42); char val[20]; Double_t xo,yo,zo; for (int j=0; jAt(j); //contLevel e' la lista dei grafici/curve dell'isolinea jth zo=contours[j]; printf ("Isocurve passing by the contour Z=%f\n", zo); cout<<"has "<GetSize()<<" graphs"<First(); //scrive in "curv" il primo grafico (linea disgiunta) che forma il contorno for (int i=0;iGetSize();i++){//cicla fino al numero di grafici che formano l'isolinea curv->GetPoint(0,xo,yo); nGraphs++; printf ("\t Graph: %d -- %d elements \n", nGraphs,curv->GetN()); //cosa fa GetN? Conta le entrate che stanno sulla linea "curv"? curv->SetLineWidth(2.5); curv->Draw("C"); sprintf(val, "%.3g", zo); //la.DrawLatex(xo,yo,val); // Here above line plots isocurve value on the plot itself. Most of the time // it creates confusion with no help...so I keep the line commented at now curv=(TGraph*)contLevel->After(curv);// prende il grafico successivo della TList }//end i-loop }//end j-loop // add the legend leg2 = new TLegend(0.40,0.11,0.9,0.20); leg2->SetTextSize(0.035); leg2->AddEntry(curv,"Current threshold contour","l");//"l" use the line that markes the gEll graph leg2->Draw(); c->Update(); ptgeom->Draw(); if (overloop) ptwarn->Draw(); // === OUTPUT BANNERS ================================================== // Add to the canvas the TPave where writing the measurement information c->cd(); //let's number the pads TPaveText *pad1=new TPaveText (0,0.5,0.019,1,"NDC"); //NDC: normalized coordinates (x1,y1,x2,y2) pad1->SetBorderSize(1); pad1->SetFillColor(29); TText *tpad1=pad1->AddText("PLOT 1"); tpad1->SetTextColor(1); tpad1->SetTextSize(0.025);//pixel size= textSize*canvasWidth tpad1->SetTextFont(102); tpad1->SetTextAlign(22); tpad1->SetTextAngle(90); pad1->Draw(); TPaveText *pad2=new TPaveText (0.32,0.5,0.339,1,"NDC"); pad2->SetBorderSize(1); pad2->SetFillColor(29); TText *tpad2=pad2->AddText("PLOT 2"); tpad2->SetTextColor(1); tpad2->SetTextSize(0.025);//pixel size= textSize*canvasWidth tpad2->SetTextFont(102); tpad2->SetTextAlign(22); tpad2->SetTextAngle(90); pad2->Draw(); TPaveText *pad3=new TPaveText (0,0,0.019,0.5,"NDC"); pad3->SetBorderSize(1); pad3->SetFillColor(29); TText *tpad3=pad3->AddText("PLOT 3"); tpad3->SetTextColor(1); tpad3->SetTextSize(0.025);//pixel size= textSize*canvasWidth tpad3->SetTextFont(102); tpad3->SetTextAlign(22); tpad3->SetTextAngle(90); pad3->Draw(); TPaveText *pad4=new TPaveText (0.32,0,0.339,0.5,"NDC"); pad4->SetBorderSize(1); pad4->SetFillColor(29); TText *tpad4=pad4->AddText("PLOT 4"); tpad4->SetTextColor(1); tpad4->SetTextSize(0.025);//pixel size= textSize*canvasWidth tpad4->SetTextFont(102); tpad4->SetTextAlign(22); tpad4->SetTextAngle(90); pad4->Draw(); //.................. //let's write the output information TPaveText *ptbunner=new TPaveText (0.64,0,1,1,"NDC"); ptbunner->SetBorderSize(1); ptbunner->SetFillColor(29); TText *t1=ptbunner->AddText("EMITTANCE MEASUREMENT PROGRAM"); t1->SetTextColor(1); t1->SetTextSize(0.03);//pixel size= textSize*canvasWidth t1->SetTextFont(102); t1->SetTextAlign(22); TText *t2=ptbunner->AddText(""); t2->SetTextColor(1); t2->SetTextSize(0.015); t2->SetTextFont(102); t2->SetTextAlign(23); bunner="#splitline{Input file name:}{"; bunner.append(name); bunner.append ("}"); TText *t3=ptbunner->AddText(bunner.c_str()); t3->SetTextColor(1); t3->SetTextSize(0.02); t3->SetTextFont(102); t3->SetTextAlign(13); char elab[200]; if (center) sprintf (elab,"#splitline{Data elaboration parameters:}{Bias thr %d%%, Noise thr %d%%, centred}", ThrB, ThrZ); else sprintf (elab,"#splitline{Data elaboration parameters:}{Bias thr %d%%, Noise thr %d%%, not-centred}", ThrB, ThrZ); TText *t4=ptbunner->AddText(elab); t4->SetTextColor(1); t4->SetTextSize(0.02);//pixel size= textSize*canvasWidth t4->SetTextFont(102); t4->SetTextAlign(13); TText *t10=ptbunner->AddText("#diamond PLOT 1"); t10->SetTextColor(1); t10->SetTextSize(0.025); t10->SetTextFont(102); t10->SetTextAlign(13); char description101[200]; sprintf (description101,"#splitline{3D plot of elaborated data used to compute}{RMS and Geom emittance.}"); TText *t101=ptbunner->AddText(description101); t101->SetTextColor(1); t101->SetTextSize(0.02); t101->SetTextFont(102); t101->SetTextAlign(13); TText *t20=ptbunner->AddText("#diamond PLOT 2"); t20->SetTextColor(1); t20->SetTextSize(0.025); t20->SetTextFont(102); t20->SetTextAlign(13); char description201[200]; sprintf (description201,"#splitline{Top-view of not centred data.}{Barycenter: x=%.3f mm, xp=%.3f mrad}",x_center, xp_center); TText *t201=ptbunner->AddText(description201); t201->SetTextColor(1); t201->SetTextSize(0.02); t201->SetTextFont(102); t201->SetTextAlign(13); TText *t30=ptbunner->AddText("#diamond PLOT 3"); t30->SetTextColor(1); t30->SetTextSize(0.025); t30->SetTextFont(102); t30->SetTextAlign(13); char description301[200]; sprintf (description301,"#splitline{Top-view of elaborated data. RMS and}{reference ellipses plotted on it.}"); TText *t301=ptbunner->AddText(description301); t301->SetTextColor(1); t301->SetTextSize(0.02); t301->SetTextFont(102); t301->SetTextAlign(13); char description302[200]; sprintf (description302,"#splitline{Reference ellipse: beam charge inside %d%%,}{alpha=%.3g, beta=%.3g m, emitt=%.3g #pi#upointmm#upointmrad}", ref_IN,alpha_ref, beta_ref, emi_ref); TText *t302=ptbunner->AddText(description302); t302->SetTextColor(1); t302->SetTextSize(0.02); t302->SetTextFont(102); t302->SetTextAlign(13); char description303[200]; sprintf (description303,"#splitline{RMS ellipse: beam charge inside %d%%, }{alpha=%.3g, beta=%.3g m, emitt=%.3g #pi#upointmm#upointmrad}",rms_IN, alpha, beta, emi); TText *t303=ptbunner->AddText(description303); t303->SetTextColor(1); t303->SetTextSize(0.02); t303->SetTextFont(102); t303->SetTextAlign(13); TText *t40=ptbunner->AddText("#diamond PLOT 4"); t40->SetTextColor(1); t40->SetTextSize(0.025); t40->SetTextFont(102); t40->SetTextAlign(13); char description401[200]; sprintf (description401,"#splitline{Top-view of elaborated data contributing}{to the Geom emitt computation.}"); TText *t401=ptbunner->AddText(description401); t401->SetTextColor(1); t401->SetTextSize(0.02); t401->SetTextFont(102); t401->SetTextAlign(13); char description402[200]; sprintf (description402,"#splitline{Required=%.1f%%, Effective=%.1f%%}{Current threshold=%.3g uA}", P,Peff,Ithreshold); TText *t402=ptbunner->AddText(description402); t402->SetTextColor(1); t402->SetTextSize(0.02); t402->SetTextFont(102); t402->SetTextAlign(13); char description403[200]; sprintf (description403,"#splitline{Maximum discrepancy between %.1f%% and %.1f%%}{geom emitt is %.2f#pi#upointmm#upointmrad}", P,Peff,maxDiscr); TText *t403=ptbunner->AddText(description403); t403->SetTextColor(1); t403->SetTextSize(0.02); t403->SetTextFont(102); t403->SetTextAlign(13); char description404[200]; sprintf (description404,"Geom Emitt at %.1f%%=(%.2f#pm%.2f)#pi#upointmm#upointmrad", Peff,emiP,emiRes); TText *t404=ptbunner->AddText(description404); t404->SetTextColor(1); t404->SetTextSize(0.02); t404->SetTextFont(102); t404->SetTextAlign(13); ptbunner->Draw(); // === SAVE THE OUTPUT =========================================================================== c->cd(); outname="C:\\measOutput\\"; outname.append(name); outname.append("_plot.gif"); c->SaveAs(outname.c_str()); return 0; }