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