[u/mrichter/AliRoot.git] / TOF / AliTOFanalyzeMatchingV2Phenix.C

#ifndef __CINT__
#include "TFile.h"
#include "TTask.h"
#include "TTree.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TFolder.h"
#include "TBenchmark.h"       
#include "AliTOFTrackV2.h"
#endif


Int_t AliTOFanalyzeMatchingV2Phenix() 
{
  // Read TPC and TRD reconstructed tracks and matched 
  // with TOF digits
  // origin: F. Pierella | pierella@bo.infn.it
  // Analysis done a la TOF in PHENIX

  gROOT->Reset();
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  gStyle->SetPalette(1);
  gStyle->SetCanvasColor(33);
  gStyle->SetFrameFillColor(18);
  

  gBenchmark->Start("AliTOFanalyzeMatchingV2");


  if(gAlice){
    delete gAlice;
    gAlice=0;
  }
  else{
    // Dynamically link some shared libs
    if(gClassTable->GetID("AliRun") < 0) {
      gROOT->LoadMacro("$ALICE_ROOT/macros/loadlibs.C");
      loadlibs();
    } // end if
  }

  // ==================>histos
  TH2F *hpvsm     = new TH2F("hpvsm"," ",170,-0.5,1.2,400,-4.,4.);
  hpvsm->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
  hpvsm->GetYaxis()->SetTitle("charge*momentum [GeV/c]");

  TH2F *hdEdXvsp  = new TH2F("hdEdXvsp"," ",400,-4.,4.,1000,0.,1000.);
  hdEdXvsp->GetYaxis()->SetTitle("dE/dX [au]");
  hdEdXvsp->GetXaxis()->SetTitle("charge*momentum [GeV/c]");

  TH2F *hdEdXvsabsp  = new TH2F("hdEdXvsabsp"," ",400,0.,4.,100,0.,1000.);
  hdEdXvsabsp->GetYaxis()->SetTitle("dE/dX [au]");
  hdEdXvsabsp->GetXaxis()->SetTitle("momentum [GeV/c]");

  TH2F *hinvpvstof= new TH2F("hinvpvstof"," ",300,10.,40.,400,-4.,4.);
  hinvpvstof->GetYaxis()->SetTitle("charge/momentum [(GeV/c)^{-1}]");
  hinvpvstof->GetXaxis()->SetTitle("Time of flight [ns]");

  TH2F *htofvsp   = new TH2F("htofvsp"," ",400,0.,4.,300,10.,40.);
  htofvsp->GetYaxis()->SetTitle("Time of flight [ns]");
  htofvsp->GetXaxis()->SetTitle("momentum [GeV/c]");

  TH2F *hsqmvsp   = new TH2F("hsqmvsp"," ",400,-4.,4.,130,-0.1,1.2);
  hsqmvsp->GetYaxis()->SetTitle("mass^{2} (GeV/c^{2})^{2}");
  hsqmvsp->GetXaxis()->SetTitle("charge*momentum [GeV/c]");


  TH2F *hinvbgvsp   = new TH2F("hinvbgvsp"," ",400,-4.,4.,50,0.,5);
  hinvbgvsp->GetXaxis()->SetTitle("charge*momentum [GeV/c]");
  hinvbgvsp->GetYaxis()->SetTitle("1/(#beta #gamma)");

  TH2F *hpvstofde     = new TH2F("hpvstofde"," ",350,-5,30.,400,-4.,4.);
  hpvstofde->GetXaxis()->SetTitle("Time-of-flight difference from electron [ns]");
  hpvstofde->GetYaxis()->SetTitle("charge*momentum [GeV/c]");


  TH2F *hinvpvstofde     = new TH2F("hinvpvstofde"," ",350,-5,30.,400,-4.,4.);
  hinvpvstofde->GetXaxis()->SetTitle("Time-of-flight difference from electron [ns]");
  hinvpvstofde->GetYaxis()->SetTitle("charge/momentum [(GeV/c)^{-1}]");
  // =======================================> done


  // chain in the case of more events
  TChain ch("T");
  ch.Add("tofTracks*.root"); // use wildcards

  TClonesArray *arr = new TClonesArray("AliTOFTrackV2"); 
  T->GetBranch("tracks")->SetAutoDelete(kFALSE);               
  T->SetBranchAddress("tracks",&arr);                  
  Int_t nentries = (Int_t)(T->GetEntries());           
  cout << "number of events " << nentries << endl;
  for (Int_t ev=0;ev<nentries;ev++) {
    arr->Clear();
    T->GetEntry(ev);
    Int_t ntracks = arr->GetEntriesFast();   
    cout << ntracks << endl;
    for (Int_t i=0;i<ntracks;i++) {
      AliTOFTrackV2 *toftrack = (AliTOFTrackV2*)arr->At(i);
      if(toftrack->GetTrackLabel()<0) continue; // reject fake

      Float_t momTPC=toftrack->GetPTPC(); // signed momentum
      if(toftrack->GetPdgCode()<0) momTPC=-momTPC;
      Int_t pdgCode=toftrack->GetPdgCode();
      //if(TMath::Abs(pdgCode)!=211) continue;
      hdEdXvsp->Fill(momTPC,toftrack->GetdEdX());
      hdEdXvsabsp->Fill(TMath::Abs(momTPC),toftrack->GetdEdX());

      Int_t matchStatus=toftrack->GetMatchingStatus();
      Float_t tof=toftrack->GetTof();
      if(matchStatus==3 || matchStatus==4){ // track has a tof
	htofvsp->Fill(TMath::Abs(momTPC),tof);
	hinvpvstof->Fill(tof,1./momTPC);
      } // if(matchStatus==3 || matchStatus==4)

      Float_t length=toftrack->GetLength();
      // starting mass calculation/ only possible with track length
      if(length>0. && (matchStatus==3 || matchStatus==4)) // to be skipped when track length will be in the Kalman 
	{
	  Float_t squareMass=momTPC*momTPC*((29.9792*tof/length)*(29.9792*tof/length)-1);
	  hsqmvsp->Fill(momTPC,TMath::Abs(squareMass));
	  Float_t dummy=squareMass;
	  if(dummy<0) dummy=-dummy;
	  Float_t mass=TMath::Sqrt(dummy);
	  if(squareMass<0) mass=-mass;
	  hpvsm->Fill(mass,momTPC);
	  Float_t invbg=TMath::Sqrt(dummy)/TMath::Abs(momTPC);
	  hinvbgvsp->Fill(momTPC,invbg);
	  Float_t betaEl=TMath::Abs(momTPC)/TMath::Sqrt(momTPC*momTPC+0.000510999*0.000510999);
	  Float_t tofEl=length/(betaEl*29.9792); // [ns]
	  hpvstofde->Fill((tof-tofEl),momTPC);
	  hinvpvstofde->Fill((tof-tofEl),1./momTPC);
	} // if(toftrack->GetLength()>0.)

    } // end loop on tracks
  } // end loop on 'events'
  
  gStyle->SetPalette(1);
  gStyle->SetCanvasColor(10);
  gStyle->SetFrameFillColor(10);

  //TPaveLabel pl;

  TCanvas* c1 = new TCanvas("c1", "TOF a la PHENIX (I)",210, 210, 740, 690);
  c1->SetFillColor(10);
  c1->cd();
  gPad->SetGrid();
  c1->SetHighLightColor(2);
  c1->SetFillColor(10);
  c1->SetBorderSize(2);

  c1->Divide(2,2);
  c1->cd(1);
  gPad->SetGridx();
  gPad->SetGridy();
  hpvsm->Draw("contz");
  
  //Float_t x1=0.67, y1=0.875, x2=0.85, y2=0.95;
  /*hpvsm->Draw("contz");*/ 
  //pl.DrawPaveLabel(x1,y1,x2,y2,"CONTZ","brNDC");
  
  
  c1->cd(2);
  gPad->SetGridx();
  gPad->SetGridy();
  hdEdXvsp->Draw();
  c1->cd(3);
  gPad->SetGridx();
  gPad->SetGridy();
  hinvpvstof->Draw();
  c1->cd(4);
  gPad->SetGridx();
  gPad->SetGridy();
  htofvsp->Draw();

  TCanvas* c2 = new TCanvas("c2", "TOF a la PHENIX (II)", 210, 210, 740, 690);
  c2->SetFillColor(10);
  c2->cd();

  //gStyle->SetOptStat(0);
  c2->SetHighLightColor(2);
  c2->SetFillColor(10);
  c2->SetBorderSize(2);
  c2->SetGridy();
  c2->Divide(2,2);


  c2->cd(1);
  gPad->SetGridx();
  gPad->SetGridy();
  hsqmvsp->Draw("contz");

  c2->cd(2);
  gPad->SetGridx();
  gPad->SetGridy();
  hinvbgvsp->Draw();


  c2->cd(3);
  gPad->SetGridx();
  gPad->SetGridy();
  gPad->SetLogx();
  hdEdXvsabsp->Draw();


  c2->cd(4);
  gPad->SetGridx();
  gPad->SetGridy();
  hpvstofde->Draw();


  TCanvas* c2b = new TCanvas("c2b", " ", 210, 210, 740, 690);
  c2b->SetFillColor(10);
  c2b->cd();
  //gStyle->SetOptStat(0);
  c2b->SetHighLightColor(2);
  c2b->SetFillColor(10);
  c2b->SetBorderSize(2);
  c2b->SetGridy();
  c2b->Divide(2,2);
  c2b->cd(1);
  gPad->SetGridx();
  gPad->SetGridy();
  hinvpvstofde->SetFillColor(33);
  hinvpvstofde->SetFillStyle(3001);
  hinvpvstofde->Draw();
  c2b->cd(2);
  gPad->SetGridx();
  gPad->SetGridy();
  hpvstofde->Draw();


  TCanvas* c3 = new TCanvas("c3", "TOF a la PHENIX (II)", 210, 210, 740, 690);
  c3->SetFillColor(10);
  c3->cd();

  //gStyle->SetOptStat(0);
  c3->SetHighLightColor(2);
  c3->SetFillColor(10);
  c3->SetBorderSize(2);
  c3->SetGridy();
  c3->Divide(2,2);
  c3->cd(1); 
  gPad->SetGridx();
  gPad->SetGridy();
  hpvsm->Draw("contz");

  c3->cd(2);
  gPad->SetGridx();
  gPad->SetGridy();
  hsqmvsp->Draw("contz");

  c3->cd(3);
  gPad->SetGridx();
  gPad->SetGridy();
  gPad->SetLogx();
  hdEdXvsabsp->Draw();

}
Commit	Line	Data
	1	#ifndef __CINT__
	2	#include "TFile.h"
	3	#include "TTask.h"
	4	#include "TTree.h"
	5	#include "TSystem.h"
	6	#include "TROOT.h"
	7	#include "TFolder.h"
	8	#include "TBenchmark.h"
	9	#include "AliTOFTrackV2.h"
	10	#endif
	11
	12
	13	Int_t AliTOFanalyzeMatchingV2Phenix()
	14	{
	15	// Read TPC and TRD reconstructed tracks and matched
	16	// with TOF digits
	17	// origin: F. Pierella \| pierella@bo.infn.it
	18	// Analysis done a la TOF in PHENIX
	19
	20	gROOT->Reset();
	21	gStyle->SetOptStat(0);
	22	gStyle->SetOptTitle(0);
	23	gStyle->SetPalette(1);
	24	gStyle->SetCanvasColor(33);
	25	gStyle->SetFrameFillColor(18);
	26
	27
	28	gBenchmark->Start("AliTOFanalyzeMatchingV2");
	29
	30
	31	if(gAlice){
	32	delete gAlice;
	33	gAlice=0;
	34	}
	35	else{
	36	// Dynamically link some shared libs
	37	if(gClassTable->GetID("AliRun") < 0) {
	38	gROOT->LoadMacro("$ALICE_ROOT/macros/loadlibs.C");
	39	loadlibs();
	40	} // end if
	41	}
	42
	43	// ==================>histos
	44	TH2F *hpvsm = new TH2F("hpvsm"," ",170,-0.5,1.2,400,-4.,4.);
	45	hpvsm->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
	46	hpvsm->GetYaxis()->SetTitle("charge*momentum [GeV/c]");
	47
	48	TH2F *hdEdXvsp = new TH2F("hdEdXvsp"," ",400,-4.,4.,1000,0.,1000.);
	49	hdEdXvsp->GetYaxis()->SetTitle("dE/dX [au]");
	50	hdEdXvsp->GetXaxis()->SetTitle("charge*momentum [GeV/c]");
	51
	52	TH2F *hdEdXvsabsp = new TH2F("hdEdXvsabsp"," ",400,0.,4.,100,0.,1000.);
	53	hdEdXvsabsp->GetYaxis()->SetTitle("dE/dX [au]");
	54	hdEdXvsabsp->GetXaxis()->SetTitle("momentum [GeV/c]");
	55
	56	TH2F *hinvpvstof= new TH2F("hinvpvstof"," ",300,10.,40.,400,-4.,4.);
	57	hinvpvstof->GetYaxis()->SetTitle("charge/momentum [(GeV/c)^{-1}]");
	58	hinvpvstof->GetXaxis()->SetTitle("Time of flight [ns]");
	59
	60	TH2F *htofvsp = new TH2F("htofvsp"," ",400,0.,4.,300,10.,40.);
	61	htofvsp->GetYaxis()->SetTitle("Time of flight [ns]");
	62	htofvsp->GetXaxis()->SetTitle("momentum [GeV/c]");
	63
	64	TH2F *hsqmvsp = new TH2F("hsqmvsp"," ",400,-4.,4.,130,-0.1,1.2);
	65	hsqmvsp->GetYaxis()->SetTitle("mass^{2} (GeV/c^{2})^{2}");
	66	hsqmvsp->GetXaxis()->SetTitle("charge*momentum [GeV/c]");
	67
	68
	69	TH2F *hinvbgvsp = new TH2F("hinvbgvsp"," ",400,-4.,4.,50,0.,5);
	70	hinvbgvsp->GetXaxis()->SetTitle("charge*momentum [GeV/c]");
	71	hinvbgvsp->GetYaxis()->SetTitle("1/(#beta #gamma)");
	72
	73	TH2F *hpvstofde = new TH2F("hpvstofde"," ",350,-5,30.,400,-4.,4.);
	74	hpvstofde->GetXaxis()->SetTitle("Time-of-flight difference from electron [ns]");
	75	hpvstofde->GetYaxis()->SetTitle("charge*momentum [GeV/c]");
	76
	77
	78	TH2F *hinvpvstofde = new TH2F("hinvpvstofde"," ",350,-5,30.,400,-4.,4.);
	79	hinvpvstofde->GetXaxis()->SetTitle("Time-of-flight difference from electron [ns]");
	80	hinvpvstofde->GetYaxis()->SetTitle("charge/momentum [(GeV/c)^{-1}]");
	81	// =======================================> done
	82
	83
	84
	85	// chain in the case of more events
	86	TChain ch("T");
	87	ch.Add("tofTracks*.root"); // use wildcards
	88
	89	TClonesArray *arr = new TClonesArray("AliTOFTrackV2");
	90	T->GetBranch("tracks")->SetAutoDelete(kFALSE);
	91	T->SetBranchAddress("tracks",&arr);
	92	Int_t nentries = (Int_t)(T->GetEntries());
	93	cout << "number of events " << nentries << endl;
	94	for (Int_t ev=0;ev<nentries;ev++) {
	95	arr->Clear();
	96	T->GetEntry(ev);
	97	Int_t ntracks = arr->GetEntriesFast();
	98	cout << ntracks << endl;
	99	for (Int_t i=0;i<ntracks;i++) {
	100	AliTOFTrackV2 toftrack = (AliTOFTrackV2)arr->At(i);
	101	if(toftrack->GetTrackLabel()<0) continue; // reject fake
	102
	103	Float_t momTPC=toftrack->GetPTPC(); // signed momentum
	104	if(toftrack->GetPdgCode()<0) momTPC=-momTPC;
	105	Int_t pdgCode=toftrack->GetPdgCode();
	106	//if(TMath::Abs(pdgCode)!=211) continue;
	107	hdEdXvsp->Fill(momTPC,toftrack->GetdEdX());
	108	hdEdXvsabsp->Fill(TMath::Abs(momTPC),toftrack->GetdEdX());
	109
	110	Int_t matchStatus=toftrack->GetMatchingStatus();
	111	Float_t tof=toftrack->GetTof();
	112	if(matchStatus==3 \|\| matchStatus==4){ // track has a tof
	113	htofvsp->Fill(TMath::Abs(momTPC),tof);
	114	hinvpvstof->Fill(tof,1./momTPC);
	115	} // if(matchStatus==3 \|\| matchStatus==4)
	116
	117	Float_t length=toftrack->GetLength();
	118	// starting mass calculation/ only possible with track length
	119	if(length>0. && (matchStatus==3 \|\| matchStatus==4)) // to be skipped when track length will be in the Kalman
	120	{
	121	Float_t squareMass=momTPCmomTPC((29.9792tof/length)(29.9792*tof/length)-1);
	122	hsqmvsp->Fill(momTPC,TMath::Abs(squareMass));
	123	Float_t dummy=squareMass;
	124	if(dummy<0) dummy=-dummy;
	125	Float_t mass=TMath::Sqrt(dummy);
	126	if(squareMass<0) mass=-mass;
	127	hpvsm->Fill(mass,momTPC);
	128	Float_t invbg=TMath::Sqrt(dummy)/TMath::Abs(momTPC);
	129	hinvbgvsp->Fill(momTPC,invbg);
	130	Float_t betaEl=TMath::Abs(momTPC)/TMath::Sqrt(momTPCmomTPC+0.0005109990.000510999);
	131	Float_t tofEl=length/(betaEl*29.9792); // [ns]
	132	hpvstofde->Fill((tof-tofEl),momTPC);
	133	hinvpvstofde->Fill((tof-tofEl),1./momTPC);
	134	} // if(toftrack->GetLength()>0.)
	135
	136	} // end loop on tracks
	137	} // end loop on 'events'
	138
	139	gStyle->SetPalette(1);
	140	gStyle->SetCanvasColor(10);
	141	gStyle->SetFrameFillColor(10);
	142
	143	//TPaveLabel pl;
	144
	145	TCanvas* c1 = new TCanvas("c1", "TOF a la PHENIX (I)",210, 210, 740, 690);
	146	c1->SetFillColor(10);
	147	c1->cd();
	148	gPad->SetGrid();
	149	c1->SetHighLightColor(2);
	150	c1->SetFillColor(10);
	151	c1->SetBorderSize(2);
	152
	153	c1->Divide(2,2);
	154	c1->cd(1);
	155	gPad->SetGridx();
	156	gPad->SetGridy();
	157	hpvsm->Draw("contz");
	158
	159	//Float_t x1=0.67, y1=0.875, x2=0.85, y2=0.95;
	160	/hpvsm->Draw("contz");/
	161	//pl.DrawPaveLabel(x1,y1,x2,y2,"CONTZ","brNDC");
	162
	163
	164
	165	c1->cd(2);
	166	gPad->SetGridx();
	167	gPad->SetGridy();
	168	hdEdXvsp->Draw();
	169	c1->cd(3);
	170	gPad->SetGridx();
	171	gPad->SetGridy();
	172	hinvpvstof->Draw();
	173	c1->cd(4);
	174	gPad->SetGridx();
	175	gPad->SetGridy();
	176	htofvsp->Draw();
	177
	178	TCanvas* c2 = new TCanvas("c2", "TOF a la PHENIX (II)", 210, 210, 740, 690);
	179	c2->SetFillColor(10);
	180	c2->cd();
	181
	182	//gStyle->SetOptStat(0);
	183	c2->SetHighLightColor(2);
	184	c2->SetFillColor(10);
	185	c2->SetBorderSize(2);
	186	c2->SetGridy();
	187	c2->Divide(2,2);
	188
	189
	190	c2->cd(1);
	191	gPad->SetGridx();
	192	gPad->SetGridy();
	193	hsqmvsp->Draw("contz");
	194
	195	c2->cd(2);
	196	gPad->SetGridx();
	197	gPad->SetGridy();
	198	hinvbgvsp->Draw();
	199
	200
	201	c2->cd(3);
	202	gPad->SetGridx();
	203	gPad->SetGridy();
	204	gPad->SetLogx();
	205	hdEdXvsabsp->Draw();
	206
	207
	208	c2->cd(4);
	209	gPad->SetGridx();
	210	gPad->SetGridy();
	211	hpvstofde->Draw();
	212
	213
	214	TCanvas* c2b = new TCanvas("c2b", " ", 210, 210, 740, 690);
	215	c2b->SetFillColor(10);
	216	c2b->cd();
	217	//gStyle->SetOptStat(0);
	218	c2b->SetHighLightColor(2);
	219	c2b->SetFillColor(10);
	220	c2b->SetBorderSize(2);
	221	c2b->SetGridy();
	222	c2b->Divide(2,2);
	223	c2b->cd(1);
	224	gPad->SetGridx();
	225	gPad->SetGridy();
	226	hinvpvstofde->SetFillColor(33);
	227	hinvpvstofde->SetFillStyle(3001);
	228	hinvpvstofde->Draw();
	229	c2b->cd(2);
	230	gPad->SetGridx();
	231	gPad->SetGridy();
	232	hpvstofde->Draw();
	233
	234
	235
	236	TCanvas* c3 = new TCanvas("c3", "TOF a la PHENIX (II)", 210, 210, 740, 690);
	237	c3->SetFillColor(10);
	238	c3->cd();
	239
	240	//gStyle->SetOptStat(0);
	241	c3->SetHighLightColor(2);
	242	c3->SetFillColor(10);
	243	c3->SetBorderSize(2);
	244	c3->SetGridy();
	245	c3->Divide(2,2);
	246	c3->cd(1);
	247	gPad->SetGridx();
	248	gPad->SetGridy();
	249	hpvsm->Draw("contz");
	250
	251	c3->cd(2);
	252	gPad->SetGridx();
	253	gPad->SetGridy();
	254	hsqmvsp->Draw("contz");
	255
	256	c3->cd(3);
	257	gPad->SetGridx();
	258	gPad->SetGridy();
	259	gPad->SetLogx();
	260	hdEdXvsabsp->Draw();
	261
	262	}