StRoot: StFtpcCalibMaker/macros/vertex.C Source File

00001 // $Id: vertex.C,v 1.3 2007/01/19 08:23:07 jcs Exp $
00002 //
00003 // $Log: vertex.C,v $
00004 // Revision 1.3  2007/01/19 08:23:07  jcs
00005 // replace true, false with kTRUE, kFALSE
00006 //
00007 // Revision 1.2  2006/03/15 15:14:06  jcs
00008 // add lines for listing CVS update info
00009 //
00010 
00011 // vertex-extrapolation (d.h. plot - program for new Tree-Structure !!!)
00012 
00013 struct HIT 
00014 {
00015   Float_t x,y,z;
00016   Float_t rad,phi;
00017   Float_t raderror,phierror;
00018 };
00019 
00020 struct CLUSTER
00021 {
00022   Float_t timepos,padpos,timesigma,padsigma;
00023   Float_t peakheight, charge;
00024   Int_t timebin,pad;
00025   Int_t padlength,timelength;
00026   Int_t row,sec;
00027   Int_t flag;
00028   Int_t numpeaks;
00029 };
00030 
00031 struct EVENT
00032 {
00033   Float_t run;
00034   Int_t nevent;
00035 };
00036 
00037 struct TEVENT
00038 {
00039   Float_t run;
00040   Int_t nevent;
00041 };
00042 
00043 
00044 struct TCLUSTER
00045 {
00046   Int_t row,sec, padlength, timelength;
00047   Float_t peakheight, charge;
00048   Int_t ntracks;
00049   Float_t padpos, timepos;
00050   Float_t padpossigma, timepossigma; 
00051 };
00052 
00053 struct THIT
00054 {
00055   Float_t x,y,z;
00056   Float_t ex,ey,ez;
00057   Float_t globResX,globResY,globResPhi,globResR;
00058   Float_t primResX,primResY,primResPhi,primResR;
00059 };
00060 
00061 struct TREVENT
00062 {
00063   Float_t run;
00064   Int_t nevent;
00065 };
00066 
00067 struct TRACK
00068 {
00069   Float_t px,py,pz;
00070   Float_t eta,p,pt;
00071   Int_t npoints;
00072   Int_t charge;
00073   Int_t type;
00074   Int_t sec;
00075 };
00076 
00077 struct MVERTEX
00078 {
00079   Float_t x,y,z;
00080 };
00081 
00082 CLUSTER cluster;
00083 HIT hit;
00084 TCLUSTER tcluster;
00085 THIT thit;
00086 TRACK track;
00087 EVENT event;
00088 TEVENT tevent;
00089 TREVENT trevent;
00090 MVERTEX mvertex;
00091 
00092 
00093 void vertex(TString eingabe,int evt)
00094 {
00095 
00096   TBranch *bhit, *bevent, *bcluster;
00097   TBranch *bthit, *btcluster, *btevent;
00098   TBranch *btrevent, *btrack,*btrvertex;
00099   TDirectory *histdir, *vertexdir;
00100 
00101   TTree *dtree;
00102   TTree *dttree;
00103   TTree *dtrtree;
00104 
00105   TStyle *plain  = new TStyle("Plain","Plain Style (no colors/fill areas)");
00106   plain->SetTitleOffset(1.25);
00107   plain->SetCanvasBorderMode(0);
00108   plain->SetPadBorderMode(0);
00109   plain->SetPadColor(10);
00110   plain->SetCanvasColor(10);
00111   plain->SetTitleColor(0);
00112   plain->SetStatColor(0);
00113   plain->SetPalette(1);
00114   plain->SetOptStat(11);
00115   plain->SetOptFit();
00116 
00117   TString histoname;
00118 
00119   cout<<"Vertex Analysis started..."<<endl;
00120   cout<<endl;
00121 
00122   TFile *f=new TFile(eingabe+".root");
00123 
00124   eingabe +="-vertex";
00125   TPostScript *fps=new TPostScript(eingabe+".ps",112);
00126 
00127   // get trees in file !!!
00128 
00129   /*
00130   dtree=(TTree*) f->Get("cl");
00131   bhit=dtree->GetBranch("hit");
00132   bhit->SetAddress(&hit);
00133   bcluster=dtree->GetBranch("cluster");
00134   bcluster->SetAddress(&cluster);
00135   bevent=dtree->GetBranch("event");
00136   bevent->SetAddress(&event);
00137   */
00138 
00139   dtrtree=(TTree*) f->Get("tr");
00140   btrevent=dtrtree->GetBranch("event");
00141   btrevent->SetAddress(&trevent);
00142   btrack=dtrtree->GetBranch("track");
00143   btrack->SetAddress(&track);
00144   btrvertex=dtrtree->GetBranch("vertex");
00145   btrvertex->SetAddress(&mvertex);
00146 
00147   dttree=(TTree*) f->Get("clot");
00148   btcluster=dttree->GetBranch("cluster");
00149   btcluster->SetAddress(&tcluster);
00150   bthit=dttree->GetBranch("hit");
00151   bthit->SetAddress(&thit);
00152   btevent=dttree->GetBranch("event");
00153   btevent->SetAddress(&tevent);
00154 
00155   TCanvas *c1 = new TCanvas("vertex","vertex",100,10,800,600);
00156   TH2F *hr = new TH2F("hr","",100,-275,275,100,0,30);
00157   hr->SetTitle("Vertex Extrapolation (linear)");
00158   hr->GetYaxis()->SetTitleOffset(1.1);
00159   hr->GetYaxis()->SetTitle("rad");
00160   hr->GetXaxis()->SetTitle("z");
00161   hr->DrawCopy(); 
00162 
00163   TH1F *hvertex=new TH1F("hvertex","z-position of used vertex",201,-100,100);
00164   TH1F *vertexeast=new TH1F("vertexeast","z-position of extrapolated vertex east",101,-50,50);
00165   TH1F *vertexwest=new TH1F("vertexwest","z-position of extrapolated vertex west",101,-50,50);
00166   TH1F *vertexall=new TH1F("vertexall","z-position of extrapolated vertex",101,-50,50);
00167 
00168   Int_t maxentries1 = (Int_t)btcluster->GetEntries();
00169   cout<<"Process Cluster-on-Track-Tree with "<<maxentries1<<" clusters..."<<endl; 
00170 
00171   cout<<endl;
00172   cout<<"Investigate Event "<<evt<<"..."<<endl;
00173   Float_t mean_vertex=0;
00174   Float_t mean_vertex_west=0;
00175   Float_t mean_vertex_east=0;
00176 
00177   Float_t radius[11],z[11];
00178 
00179   Bool_t EndOfEvent=kFALSE;
00180   Int_t i=0;
00181   Int_t ntracksold=0;
00182   Int_t ntrack=0;
00183   Int_t ntrackeast=0;
00184   Int_t ntrackwest=0;
00185 
00186   TF1 *linfit=new TF1("linfit","[0]*x+[1]",-275,275);
00187   linfit->SetParameters(1,10);
00188   linfit->SetLineWidth(0.5);
00189 
00190   for (int k=0;k<=maxentries1;k++)
00191     {
00192       btcluster->GetEntry(k);
00193       bthit->GetEntry(k);
00194       btevent->GetEntry(k);
00195       
00196       if (EndOfEvent) break;
00197       if (tevent.nevent > evt) 
00198        {
00199          if (i==0) break;
00200          tevent.nevent=evt;
00201          EndOfEvent=kTRUE;
00202        }
00203       if (tevent.nevent==evt)
00204         {
00205           if (tcluster.ntracks != ntracksold)
00206             {
00207 
00208               TGraph *ver=new TGraph(i,z,radius);
00209               ver->SetMarkerStyle(22);
00210               ver->SetMarkerSize(0.75);
00211               ver->Draw("P");
00212 
00213               if (z[0]<0)
00214                 {
00215                   ntrackeast++;
00216                   ver->SetMarkerColor(2);
00217                   linfit->SetLineColor(2);
00218                 }
00219               else
00220                 {
00221                   ntrackwest++;
00222                   ver->SetMarkerColor(3);
00223                   linfit->SetLineColor(3);
00224                 }
00225 
00226               ver->Fit(linfit,"QR+"); 
00227 
00228               cout<<"Event "<<tevent.nevent<<" Track "<<ntrack<<" Vertex = "<<-linfit->GetParameter(1)/linfit->GetParameter(0)<<endl;
00229               ntrack++;
00230               mean_vertex=(mean_vertex+(-linfit->GetParameter(1)/linfit->GetParameter(0)));
00231               vertexall->Fill(-linfit->GetParameter(1)/linfit->GetParameter(0));
00232               if (z[0]<0)
00233                 {
00234                   mean_vertex_east=mean_vertex_east+-linfit->GetParameter(1)/linfit->GetParameter(0);
00235                   vertexeast->Fill(-linfit->GetParameter(1)/linfit->GetParameter(0));
00236                 }
00237               else
00238                 {
00239                   mean_vertex_west=mean_vertex_west+-linfit->GetParameter(1)/linfit->GetParameter(0);   
00240                   vertexwest->Fill(-linfit->GetParameter(1)/linfit->GetParameter(0));
00241                 }
00242               i=0;
00243 
00244             }
00245           radius[i]=sqrt(thit.x*thit.x+thit.y*thit.y);
00246           z[i]=thit.z;
00247           cout<<"DEBUG: ntrack = "<<tcluster.ntracks<<" cluster "<<i<<" x,y,z = "<<thit.x<<" "<<thit.y<<" "<<thit.z<<" radius = "<<radius[i]<<endl;
00248           i++;
00249             
00250           ntracksold=tcluster.ntracks;
00251         }
00252     }
00253   cout<<"---------------------------"<<endl;
00254   cout<<"Mean_vertex = "<<mean_vertex/ntrack<<endl;
00255   cout<<"Mean_vertex_east = "<<mean_vertex_east/ntrackeast<<endl;
00256   cout<<"Mean_vertex_west = "<<mean_vertex_west/ntrackwest<<endl;
00257   cout<<"---------------------------"<<endl;
00258 
00259   TCanvas *c2 = new TCanvas("vertex2","vertex2",150,60,850,650);
00260 
00261   vertexall->DrawCopy();
00262   vertexwest->SetLineColor(3);vertexwest->DrawCopy("same");
00263   vertexeast->SetLineColor(2);vertexeast->DrawCopy("same");
00264   TLegend *leg = new TLegend(0.7,0.6,0.95,0.8);
00265   leg->SetFillColor(10);
00266   leg->SetTextSize(0.03);
00267   leg->AddEntry(vertexall,"vertex East&West","L");
00268   leg->AddEntry(vertexeast,"vertex East","L");
00269   leg->AddEntry(vertexwest,"vertex West","L");
00270   leg->Draw();
00271   //cout<<vertexwest->GetEntries()+ vertexeast->GetEntries()<<endl;
00272   Int_t maxentries2 = (Int_t)btrack->GetEntries();
00273   cout<<endl;
00274   cout<<"Process Track-Tree with "<<maxentries2<<" Tracks..."<<endl;  
00275   
00276   Int_t oldevent=0;
00277   for (int k=0;k<=maxentries2;k++)
00278   //for (int k=0;k<=1;k++)
00279     {
00280       btrack->GetEntry(k);
00281       btrvertex->GetEntry(k);
00282       btrevent->GetEntry(k);
00283       if (trevent.nevent!=oldevent)
00284         {
00285           hvertex->Fill(mvertex.z);
00286           if (trevent.nevent==evt)
00287             {
00288               cout<<"---------------------------"<<endl;
00289               cout<<"Used Vertex (Event "<<evt<<")  x = "<<mvertex.x<<" y = "<<mvertex.y<<" z = "<<mvertex.z<<endl;
00290               cout<<"---------------------------"<<endl;
00291             }
00292         }
00293       oldevent=trevent.nevent;
00294     }
00295  
00296   cout<<endl;
00297   //cout<<"Analyse beendet !"<<endl;
00298   cout<<"Histograms in ROOT-file : "<<eingabe+".root"<<endl;
00299 
00300   //ana->Write();
00301   //ana->Close();
00302   //f->Close();
00303 
00304   //return 0;
00305   //c1->Delete();
00306 }
00307