[u/mrichter/AliRoot.git] / TPC / AnalyzeESDtracks.C

//
// Process ESD tracks  - 
// Extract TPC tracks  - write them to tree
//
/*
  .L AnalyzeESDtracks.C+
  .L AliGenInfo.C+
  AnalyzESDtracks(567);   // process tracks
  // Tracks are written to the file "TPCtracks.root"
  // Now yo can analyze it
  TFile fesd("AliESDs.root");
  TTree * treeE = (TTree*)fesd.Get("esdTree");
  TFile f("TPCtracks.root")
  TTree * tree =(TTree*)f.Get("Tracks");
  AliComparisonDraw comp;
  comp->fTree = tree;

*/

// 

/*
.L AnalyzeESDtracks.C+
.L AliGenInfo.C+

//AnalyzESDtracks(567);


TFile fesd("AliESDs.root");
TTree * treeE = (TTree*)fesd.Get("esdTree");
TFile f("TPCtracks.root")
TTree * tree =(TTree*)f.Get("Tracks");
AliComparisonDraw comp;
comp->fTree = tree;

TFile fs("TPCsignal.root");
TTree *treeB =(TTree*)fs.Get("SignalB");
TTree *treeN =(TTree*)fs.Get("SignalN");
TTree *treeS =(TTree*)fs.Get("Signal");
TTree *treef =(TTree*)fs.Get("Fit");


FitSignals(treeB,"Max-Median>100&&RMS06<2.5")
TFile ffit("FitSignal.root");
TTree * treeF = (TTree*)ffit->Get("Fit");


TChain chaincl("TreeR","TreeR")
chaincl.Add("TPC.RecPoints.root/Event0/TreeR")
chaincl.Add("TPC.RecPoints1.root/Event1/TreeR")
chaincl.Add("TPC.RecPoints2.root/Event2/TreeR")
chaincl.Add("TPC.RecPoints3.root/Event3/TreeR")

*/


#include "TObject.h"
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TBranch.h"
#include "TTreeStream.h"
#include "TEventList.h"
#include "TCut.h"
#include "TFitter.h"
#include  "TMatrixD.h"

#include  "AliLog.h"
#include  "AliMagF.h"

#include  "AliESD.h"
#include  "AliESDfriend.h"
#include  "AliESDtrack.h"
#include  "AliTracker.h"
#include  "AliTPCseed.h"
#include  "AliTPCclusterMI.h"
#include  "TTreeStream.h"
#include  "TF1.h"
#include  "TGraph.h"
#include  "AliSignalProcesor.h"
#include  "TCanvas.h"


void AnalyzeESDtracks(Int_t run){
  //
  // output redirect 
  TTreeSRedirector  cstream("TPCtracks.root");
  //
  // dummy magnetic field
  AliMagF mag("aaa","aaa",1,1,10);
  AliTracker::SetFieldMap(&mag,kTRUE);
  TFile f("AliESDs.root");
  TTree * tree =(TTree*)f.Get("esdTree"); 
  AliESD * esd =0;
  tree->SetBranchAddress("ESD",&esd);
  AliESDfriend *evf=0;
  tree->AddFriend("esdFriendTree","AliESDfriends.root");
  tree->SetBranchAddress("ESDfriend",&evf);

  Int_t nevents = tree->GetEntries();
  TClonesArray *clusters = new TClonesArray("AliTPCclusterMI",160);
  for (Int_t irow=0; irow<160; irow++){
    new ((*clusters)[irow])  AliTPCclusterMI;   // iitial dummy clusters
  }
  
  for (Int_t ievent=0; ievent<nevents; ievent++){
    tree->GetEntry(ievent);
    if (!esd) continue;
    if (!evf) continue;
    esd->SetESDfriend(evf); //Attach the friend to the ESD      
    for (Int_t itrack =0; itrack<esd->GetNumberOfTracks(); itrack++){
      // Int_t itrack = 0;      
      if (esd->GetTrack(itrack)->GetFriendTrack()==0) continue;
      AliESDtrack * etrack = esd->GetTrack(itrack);
      AliESDfriendTrack * ftrack = (AliESDfriendTrack *)esd->GetTrack(itrack)->GetFriendTrack();
      AliTPCseed * seed =  (AliTPCseed*)(ftrack->GetCalibObject(0));
      if (!seed) continue;  
      for (Int_t irow=0; irow<160; irow++){
	if (seed->GetClusterFast(irow)){
	  AliTPCclusterMI * cl = new ((*clusters)[irow])  AliTPCclusterMI(*(seed->GetClusterFast(irow)));
	  cl->SetLabel(itrack,0);
	}
	else{
	  AliTPCclusterMI * cl = (AliTPCclusterMI*)clusters->At(irow);
	  cl->SetX(0); cl->SetY(0); cl->SetZ(0); cl->SetQ(0); cl->SetLabel(-1,0);
	}
      }
      Float_t dEdx = seed->GetdEdx();
      Float_t dEdxI = seed->CookdEdx(0.05,0.6,0,77);
      Float_t dEdxO = seed->CookdEdx(0.05,0.6,78,155);
      Int_t ncl = seed->GetNumberOfClusters();
      cstream<<"Tracks"<<
	"Run="<<run<<
	"Ncl="<<ncl<<
	"Event="<<ievent<<
	"dEdx="<<dEdx<<
	"dEdxI="<<dEdxI<<
	"dEdxO="<<dEdxO<<
	"Track.="<<seed<<
	"Etrack.="<<etrack<<
	"Cl.="<<clusters<<
	"\n";
    }  
  }
  //
  // Fit signal part
  //
  TFile fs("TPCsignal.root");
  TTree *treeB =(TTree*)fs.Get("SignalB");
  FitSignals(treeB,"Max-Median>150&&RMS06<2.5");
}


void FitSignals(TTree * treeB, TCut cut="Max-Median>150&&RMS06<2&&abs(Median-Mean09)<0.5"){
  AliSignalProcesor proc;
  TF1 * f1 = proc.GetAsymGauss();
  TTreeSRedirector cstream("FitSignal.root");
  TFile *f = cstream.GetFile();

  char lname[100];
  sprintf(lname,"Fit%s", cut.GetTitle());
  TEventList *list = new TEventList(lname,lname);
  sprintf(lname,">>Fit%s", cut.GetTitle());
  treeB->Draw(lname,cut);
  treeB->SetEventList(list);
  for (Int_t ievent=0; ievent<list->GetN(); ievent++){
    char ename[100];
    sprintf(ename,"Fit%d", ievent);
    Double_t nsample = treeB->Draw("Graph.fY-Mean09:Graph.fX","","",1,ievent);
    Double_t * signal  = treeB->GetV1();
    Double_t * time  = treeB->GetV2();
    Double_t maxpos =0;
    Double_t max = 0;
    for (Int_t ipos = 0; ipos<nsample; ipos++){
      if (signal[ipos]>max){
	max    = signal[ipos];
	maxpos = ipos;
      }
    }

    Int_t first = TMath::Max(maxpos-10,0.);
    Int_t last  = TMath::Min(maxpos+60, nsample);
    //
    f->cd();
    TH1F his(ename,ename,last-first,first,last);
    for (Int_t ipos=0; ipos<last-first; ipos++){
      his.SetBinContent(ipos+1,signal[ipos+first]);
    }
    treeB->Draw("Sector:Row:Pad","","",1,ievent);
    Double_t sector = treeB->GetV1()[0];
    Double_t row    = treeB->GetV2()[0];
    Double_t pad    = treeB->GetV3()[0];
    //    TGraph  graph(last-first,&time[first],&signal[first]);
    f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
    //    TH1F * his = (TH1F*)graph.GetHistogram();
    his.Fit(f1,"q");
    his.Write(ename);
    gPad->Clear();
    his.Draw();
    gPad->Update();
    Double_t params[6];
    for (Int_t ipar=0; ipar<6; ipar++) params[ipar] = f1->GetParameters()[ipar];
    Double_t chi2 = TFitter::GetFitter()->Chisquare(6,params);
    TMatrixD cov(6,6);
    cov.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
    //
    // tail cancellation
    //
    Double_t x0[1000];
    Double_t x1[1000];
    Double_t x2[1000];
    for (Int_t ipos=0; ipos<last-first; ipos++){
      x0[ipos] = signal[ipos+first];
    }
    proc.TailCancelationALTRO1(x0,x1,0.85*0.339,0.09,last-first);
    proc.TailCancelationALTRO1(x1,x2,0.85,0.789,last-first);
    //
    sprintf(ename,"Cancel1_%d", ievent);
    TH1F his1(ename,ename,last-first,first,last);
    for (Int_t ipos=0; ipos<last-first; ipos++){
      his1.SetBinContent(ipos+1,x1[ipos]);
    }
    his1.Write(ename);
    sprintf(ename,"Cancel2_%d", ievent);
    TH1F his2(ename,ename,last-first,first,last);
    for (Int_t ipos=0; ipos<last-first; ipos++){
      his2.SetBinContent(ipos+1,x1[ipos]);
    }
    f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
    his2.Fit(f1,"q");
    his2.Write(ename);
    Double_t params2[6];
    for (Int_t ipar=0; ipar<6; ipar++) params2[ipar] = f1->GetParameters()[ipar];
    Double_t chi22 = TFitter::GetFitter()->Chisquare(6,params2);    
    TMatrixD cov2(6,6);
    cov2.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());

    TGraph gr0(last-first, &time[first],x0);
    TGraph gr1(last-first, &time[first],x1);
    TGraph gr2(last-first, &time[first],x2);
    //
    cstream<<"Fit"<<
      "Sector="<<sector<<
      "Row="<<row<<
      "Pad="<<pad<<
      "First="<<first<<
      "Max="<<max<<
      "MaxPos="<<maxpos<<
      "chi2="<<chi2<<
      "chi22="<<chi22<<
      "Cov="<<&cov<<
      "Cov2="<<&cov2<<
      "gr0.="<<&gr0<<
      "gr1.="<<&gr1<<
      "gr2.="<<&gr2<<
      "p0="<<params[0]<<
      "p1="<<params[1]<<
      "p2="<<params[2]<<
      "p3="<<params[3]<<
      "p4="<<params[4]<<
      "p5="<<params[5]<<
      "p02="<<params2[0]<<
      "p12="<<params2[1]<<
      "p22="<<params2[2]<<
      "p32="<<params2[3]<<
      "p42="<<params2[4]<<
      "p52="<<params2[5]<<
      "\n";
    //    delete his;
  }

}


TChain *MakeChainCL(Int_t first, Int_t last){
  TChain *chaincl = new TChain("TreeR","TreeR");
  //
  char fname[100];
  for (Int_t i=first;i<last; i++){
    if (i>0) sprintf(fname,"TPC.RecPoints%d.root/Event%d/TreeR",i,i);
    if (i==0) sprintf(fname,"TPC.RecPoints.root/Event%d/TreeR",i);
    chaincl->Add(fname);
  }
  return chaincl;
}

TTree* GetTree(Int_t ievent){
  char fname[100];
  char tname[100];
  if (ievent>0) sprintf(fname,"TPC.RecPoints%d.root",ievent);
  if (ievent==0) sprintf(fname,"TPC.RecPoints.root");
  sprintf(tname,"Event%d/TreeR",ievent);
  TFile * f  = new TFile(fname);
  TTree * tree = (TTree*)f->Get(tname);
  return tree;

}
Commit	Line	Data
eaad8c01	1	//
	2	// Process ESD tracks -
	3	// Extract TPC tracks - write them to tree
	4	//
	5	/*
	6	.L AnalyzeESDtracks.C+
	7	.L AliGenInfo.C+
	8	AnalyzESDtracks(567); // process tracks
	9	// Tracks are written to the file "TPCtracks.root"
	10	// Now yo can analyze it
	11	TFile fesd("AliESDs.root");
	12	TTree * treeE = (TTree*)fesd.Get("esdTree");
	13	TFile f("TPCtracks.root")
	14	TTree * tree =(TTree*)f.Get("Tracks");
	15	AliComparisonDraw comp;
	16	comp->fTree = tree;
	17
	18	*/
	19
	20	//
	21
	22	/*
	23	.L AnalyzeESDtracks.C+
	24	.L AliGenInfo.C+
	25
	26	//AnalyzESDtracks(567);
	27
	28
	29	TFile fesd("AliESDs.root");
	30	TTree * treeE = (TTree*)fesd.Get("esdTree");
	31	TFile f("TPCtracks.root")
	32	TTree * tree =(TTree*)f.Get("Tracks");
	33	AliComparisonDraw comp;
	34	comp->fTree = tree;
	35
	36	TFile fs("TPCsignal.root");
	37	TTree treeB =(TTree)fs.Get("SignalB");
	38	TTree treeN =(TTree)fs.Get("SignalN");
	39	TTree treeS =(TTree)fs.Get("Signal");
	40	TTree treef =(TTree)fs.Get("Fit");
	41
	42
	43	FitSignals(treeB,"Max-Median>100&&RMS06<2.5")
	44	TFile ffit("FitSignal.root");
	45	TTree * treeF = (TTree*)ffit->Get("Fit");
	46
	47
	48	TChain chaincl("TreeR","TreeR")
	49	chaincl.Add("TPC.RecPoints.root/Event0/TreeR")
	50	chaincl.Add("TPC.RecPoints1.root/Event1/TreeR")
	51	chaincl.Add("TPC.RecPoints2.root/Event2/TreeR")
	52	chaincl.Add("TPC.RecPoints3.root/Event3/TreeR")
	53
	54	*/
	55
	56
	57
	58	#include "TObject.h"
	59	#include "TFile.h"
	60	#include "TTree.h"
a82f6396	61	#include "TChain.h"
eaad8c01	62	#include "TBranch.h"
	63	#include "TTreeStream.h"
	64	#include "TEventList.h"
	65	#include "TCut.h"
	66	#include "TFitter.h"
	67	#include "TMatrixD.h"
	68
	69	#include "AliLog.h"
	70	#include "AliMagF.h"
	71
	72	#include "AliESD.h"
	73	#include "AliESDfriend.h"
	74	#include "AliESDtrack.h"
	75	#include "AliTracker.h"
	76	#include "AliTPCseed.h"
	77	#include "AliTPCclusterMI.h"
	78	#include "TTreeStream.h"
	79	#include "TF1.h"
	80	#include "TGraph.h"
	81	#include "AliSignalProcesor.h"
	82	#include "TCanvas.h"
	83
	84
a82f6396	85	void AnalyzeESDtracks(Int_t run){
eaad8c01	86	//
	87	// output redirect
	88	TTreeSRedirector cstream("TPCtracks.root");
	89	//
	90	// dummy magnetic field
	91	AliMagF mag("aaa","aaa",1,1,10);
	92	AliTracker::SetFieldMap(&mag,kTRUE);
	93	TFile f("AliESDs.root");
	94	TTree * tree =(TTree*)f.Get("esdTree");
	95	AliESD * esd =0;
	96	tree->SetBranchAddress("ESD",&esd);
	97	AliESDfriend *evf=0;
	98	tree->AddFriend("esdFriendTree","AliESDfriends.root");
	99	tree->SetBranchAddress("ESDfriend",&evf);
	100
	101	Int_t nevents = tree->GetEntries();
	102	TClonesArray *clusters = new TClonesArray("AliTPCclusterMI",160);
	103	for (Int_t irow=0; irow<160; irow++){
	104	new ((*clusters)[irow]) AliTPCclusterMI; // iitial dummy clusters
	105	}
	106
	107	for (Int_t ievent=0; ievent<nevents; ievent++){
	108	tree->GetEntry(ievent);
	109	if (!esd) continue;
	110	if (!evf) continue;
	111	esd->SetESDfriend(evf); //Attach the friend to the ESD
	112	for (Int_t itrack =0; itrack<esd->GetNumberOfTracks(); itrack++){
	113	// Int_t itrack = 0;
	114	if (esd->GetTrack(itrack)->GetFriendTrack()==0) continue;
	115	AliESDtrack * etrack = esd->GetTrack(itrack);
	116	AliESDfriendTrack * ftrack = (AliESDfriendTrack *)esd->GetTrack(itrack)->GetFriendTrack();
	117	AliTPCseed * seed = (AliTPCseed*)(ftrack->GetCalibObject(0));
	118	if (!seed) continue;
	119	for (Int_t irow=0; irow<160; irow++){
	120	if (seed->GetClusterFast(irow)){
	121	AliTPCclusterMI * cl = new ((clusters)[irow]) AliTPCclusterMI((seed->GetClusterFast(irow)));
	122	cl->SetLabel(itrack,0);
	123	}
	124	else{
	125	AliTPCclusterMI * cl = (AliTPCclusterMI*)clusters->At(irow);
	126	cl->SetX(0); cl->SetY(0); cl->SetZ(0); cl->SetQ(0); cl->SetLabel(-1,0);
	127	}
	128	}
	129	Float_t dEdx = seed->GetdEdx();
	130	Float_t dEdxI = seed->CookdEdx(0.05,0.6,0,77);
	131	Float_t dEdxO = seed->CookdEdx(0.05,0.6,78,155);
	132	Int_t ncl = seed->GetNumberOfClusters();
	133	cstream<<"Tracks"<<
	134	"Run="<<run<<
a82f6396	135	"Ncl="<<ncl<<
eaad8c01	136	"Event="<<ievent<<
	137	"dEdx="<<dEdx<<
	138	"dEdxI="<<dEdxI<<
	139	"dEdxO="<<dEdxO<<
	140	"Track.="<<seed<<
	141	"Etrack.="<<etrack<<
	142	"Cl.="<<clusters<<
	143	"\n";
	144	}
	145	}
a82f6396	146	//
	147	// Fit signal part
	148	//
	149	TFile fs("TPCsignal.root");
	150	TTree treeB =(TTree)fs.Get("SignalB");
	151	FitSignals(treeB,"Max-Median>150&&RMS06<2.5");
eaad8c01	152	}
	153
	154
	155	void FitSignals(TTree * treeB, TCut cut="Max-Median>150&&RMS06<2&&abs(Median-Mean09)<0.5"){
	156	AliSignalProcesor proc;
	157	TF1 * f1 = proc.GetAsymGauss();
	158	TTreeSRedirector cstream("FitSignal.root");
	159	TFile *f = cstream.GetFile();
	160
	161	char lname[100];
	162	sprintf(lname,"Fit%s", cut.GetTitle());
	163	TEventList *list = new TEventList(lname,lname);
	164	sprintf(lname,">>Fit%s", cut.GetTitle());
	165	treeB->Draw(lname,cut);
	166	treeB->SetEventList(list);
	167	for (Int_t ievent=0; ievent<list->GetN(); ievent++){
	168	char ename[100];
	169	sprintf(ename,"Fit%d", ievent);
	170	Double_t nsample = treeB->Draw("Graph.fY-Mean09:Graph.fX","","",1,ievent);
	171	Double_t * signal = treeB->GetV1();
	172	Double_t * time = treeB->GetV2();
	173	Double_t maxpos =0;
	174	Double_t max = 0;
	175	for (Int_t ipos = 0; ipos<nsample; ipos++){
	176	if (signal[ipos]>max){
	177	max = signal[ipos];
	178	maxpos = ipos;
	179	}
	180	}
	181
	182	Int_t first = TMath::Max(maxpos-10,0.);
	183	Int_t last = TMath::Min(maxpos+60, nsample);
	184	//
	185	f->cd();
	186	TH1F his(ename,ename,last-first,first,last);
	187	for (Int_t ipos=0; ipos<last-first; ipos++){
	188	his.SetBinContent(ipos+1,signal[ipos+first]);
	189	}
	190	treeB->Draw("Sector:Row:Pad","","",1,ievent);
	191	Double_t sector = treeB->GetV1()[0];
	192	Double_t row = treeB->GetV2()[0];
	193	Double_t pad = treeB->GetV3()[0];
	194	// TGraph graph(last-first,&time[first],&signal[first]);
	195	f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
	196	// TH1F * his = (TH1F*)graph.GetHistogram();
	197	his.Fit(f1,"q");
	198	his.Write(ename);
	199	gPad->Clear();
	200	his.Draw();
	201	gPad->Update();
	202	Double_t params[6];
	203	for (Int_t ipar=0; ipar<6; ipar++) params[ipar] = f1->GetParameters()[ipar];
	204	Double_t chi2 = TFitter::GetFitter()->Chisquare(6,params);
	205	TMatrixD cov(6,6);
	206	cov.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
	207	//
	208	// tail cancellation
	209	//
	210	Double_t x0[1000];
	211	Double_t x1[1000];
	212	Double_t x2[1000];
	213	for (Int_t ipos=0; ipos<last-first; ipos++){
	214	x0[ipos] = signal[ipos+first];
	215	}
216	proc.TailCancelationALTRO1(x0,x1,0.85*0.339,0.09,last-first);
217	proc.TailCancelationALTRO1(x1,x2,0.85,0.789,last-first);
218	//
219	sprintf(ename,"Cancel1_%d", ievent);
220	TH1F his1(ename,ename,last-first,first,last);
221	for (Int_t ipos=0; ipos<last-first; ipos++){
222	his1.SetBinContent(ipos+1,x1[ipos]);
223	}
224	his1.Write(ename);
225	sprintf(ename,"Cancel2_%d", ievent);
226	TH1F his2(ename,ename,last-first,first,last);
227	for (Int_t ipos=0; ipos<last-first; ipos++){
228	his2.SetBinContent(ipos+1,x1[ipos]);
229	}
230	f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
231	his2.Fit(f1,"q");
232	his2.Write(ename);
233	Double_t params2[6];
234	for (Int_t ipar=0; ipar<6; ipar++) params2[ipar] = f1->GetParameters()[ipar];
235	Double_t chi22 = TFitter::GetFitter()->Chisquare(6,params2);
236	TMatrixD cov2(6,6);
237	cov2.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
238
239	TGraph gr0(last-first, &time[first],x0);
240	TGraph gr1(last-first, &time[first],x1);
241	TGraph gr2(last-first, &time[first],x2);
242	//
243	cstream<<"Fit"<<
244	"Sector="<<sector<<
245	"Row="<<row<<
246	"Pad="<<pad<<
247	"First="<<first<<
248	"Max="<<max<<
249	"MaxPos="<<maxpos<<
250	"chi2="<<chi2<<
251	"chi22="<<chi22<<
252	"Cov="<<&cov<<
253	"Cov2="<<&cov2<<
254	"gr0.="<<&gr0<<
255	"gr1.="<<&gr1<<
256	"gr2.="<<&gr2<<
257	"p0="<<params[0]<<
258	"p1="<<params[1]<<
259	"p2="<<params[2]<<
260	"p3="<<params[3]<<
261	"p4="<<params[4]<<
262	"p5="<<params[5]<<
263	"p02="<<params2[0]<<
264	"p12="<<params2[1]<<
265	"p22="<<params2[2]<<
266	"p32="<<params2[3]<<
267	"p42="<<params2[4]<<
268	"p52="<<params2[5]<<
269	"\n";
270	// delete his;
271	}
272
273	}
a82f6396	274
	275
	276
	277	TChain *MakeChainCL(Int_t first, Int_t last){
	278	TChain *chaincl = new TChain("TreeR","TreeR");
	279	//
	280	char fname[100];
	281	for (Int_t i=first;i<last; i++){
	282	if (i>0) sprintf(fname,"TPC.RecPoints%d.root/Event%d/TreeR",i,i);
	283	if (i==0) sprintf(fname,"TPC.RecPoints.root/Event%d/TreeR",i);
	284	chaincl->Add(fname);
	285	}
	286	return chaincl;
	287	}
	288
	289	TTree* GetTree(Int_t ievent){
	290	char fname[100];
	291	char tname[100];
	292	if (ievent>0) sprintf(fname,"TPC.RecPoints%d.root",ievent);
	293	if (ievent==0) sprintf(fname,"TPC.RecPoints.root");
	294	sprintf(tname,"Event%d/TreeR",ievent);
	295	TFile * f = new TFile(fname);
	296	TTree * tree = (TTree*)f->Get(tname);
	297	return tree;
	298
	299	}