]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AnalyzeESDtracks.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Data base generation class for configuration entries (Haavard)
[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+
a2411279 8 AnalyzeESDtracks(567); // process tracks
eaad8c01 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
a2411279 26//AnalyzeESDtracks(567);
eaad8c01 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"
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"
663296cb 68#include "TRobustEstimator.h"
b78c817f 69#include "TTimeStamp.h"
eaad8c01 70
71#include "AliLog.h"
72#include "AliMagF.h"
73
74#include "AliESD.h"
75#include "AliESDfriend.h"
76#include "AliESDtrack.h"
77#include "AliTracker.h"
78#include "AliTPCseed.h"
79#include "AliTPCclusterMI.h"
663296cb 80#include "AliTPCParamSR.h"
81#include "AliTPCROC.h"
82
83
eaad8c01 84#include "TTreeStream.h"
85#include "TF1.h"
86#include "TGraph.h"
87#include "AliSignalProcesor.h"
88#include "TCanvas.h"
89
90
663296cb 91void FitSignals(TTree * treeB, TCut cut="Max-Median>150&&RMS06<2&&abs(Median-Mean09)<0.5", Int_t maxS=100);
ba6ac751 92void LaserCalib(TTreeSRedirector & cstream, TTree * chain, Float_t tmin, Float_t tmax, Float_t fraction=0.7);
de31eda2 93
a82f6396 94void AnalyzeESDtracks(Int_t run){
eaad8c01 95 //
96 // output redirect
a2411279 97 TTreeSRedirector * pcstream = new TTreeSRedirector("TPCtracks.root");
98 TTreeSRedirector &cstream = *pcstream;
eaad8c01 99 //
100 // dummy magnetic field
101 AliMagF mag("aaa","aaa",1,1,10);
102 AliTracker::SetFieldMap(&mag,kTRUE);
103 TFile f("AliESDs.root");
104 TTree * tree =(TTree*)f.Get("esdTree");
105 AliESD * esd =0;
106 tree->SetBranchAddress("ESD",&esd);
107 AliESDfriend *evf=0;
108 tree->AddFriend("esdFriendTree","AliESDfriends.root");
109 tree->SetBranchAddress("ESDfriend",&evf);
110
111 Int_t nevents = tree->GetEntries();
112 TClonesArray *clusters = new TClonesArray("AliTPCclusterMI",160);
113 for (Int_t irow=0; irow<160; irow++){
114 new ((*clusters)[irow]) AliTPCclusterMI; // iitial dummy clusters
115 }
116
117 for (Int_t ievent=0; ievent<nevents; ievent++){
118 tree->GetEntry(ievent);
119 if (!esd) continue;
120 if (!evf) continue;
121 esd->SetESDfriend(evf); //Attach the friend to the ESD
122 for (Int_t itrack =0; itrack<esd->GetNumberOfTracks(); itrack++){
123 // Int_t itrack = 0;
124 if (esd->GetTrack(itrack)->GetFriendTrack()==0) continue;
125 AliESDtrack * etrack = esd->GetTrack(itrack);
126 AliESDfriendTrack * ftrack = (AliESDfriendTrack *)esd->GetTrack(itrack)->GetFriendTrack();
127 AliTPCseed * seed = (AliTPCseed*)(ftrack->GetCalibObject(0));
128 if (!seed) continue;
129 for (Int_t irow=0; irow<160; irow++){
130 if (seed->GetClusterFast(irow)){
131 AliTPCclusterMI * cl = new ((*clusters)[irow]) AliTPCclusterMI(*(seed->GetClusterFast(irow)));
132 cl->SetLabel(itrack,0);
133 }
134 else{
135 AliTPCclusterMI * cl = (AliTPCclusterMI*)clusters->At(irow);
136 cl->SetX(0); cl->SetY(0); cl->SetZ(0); cl->SetQ(0); cl->SetLabel(-1,0);
137 }
138 }
139 Float_t dEdx = seed->GetdEdx();
140 Float_t dEdxI = seed->CookdEdx(0.05,0.6,0,77);
141 Float_t dEdxO = seed->CookdEdx(0.05,0.6,78,155);
142 Int_t ncl = seed->GetNumberOfClusters();
143 cstream<<"Tracks"<<
144 "Run="<<run<<
a82f6396 145 "Ncl="<<ncl<<
eaad8c01 146 "Event="<<ievent<<
147 "dEdx="<<dEdx<<
148 "dEdxI="<<dEdxI<<
149 "dEdxO="<<dEdxO<<
150 "Track.="<<seed<<
151 "Etrack.="<<etrack<<
152 "Cl.="<<clusters<<
153 "\n";
154 }
155 }
a2411279 156 delete pcstream;
a82f6396 157 //
158 // Fit signal part
159 //
160 TFile fs("TPCsignal.root");
161 TTree *treeB =(TTree*)fs.Get("SignalB");
663296cb 162 //
163 // Fit central electrode part
164 //
165 TTreeSRedirector * pcestream = new TTreeSRedirector("TimeRoot.root");
166 TTree * treece = (TTree*)fs.Get("Signalce");
167 if (tree) {
b78c817f 168 LaserCalib(*pcestream, treece, 800,1000, 0.7);
663296cb 169 delete pcestream;
170 }
ba6ac751 171 FitSignals(treeB,"Max-Median>150&&RMS06<1.0&&RMS09<1.5&&abs(Median-Mean09)<0.2&&abs(Mean06-Mean09)<0.2",1000);
663296cb 172 //
eaad8c01 173}
174
175
663296cb 176void FitSignals(TTree * treeB, TCut cut, Int_t max){
eaad8c01 177 AliSignalProcesor proc;
178 TF1 * f1 = proc.GetAsymGauss();
179 TTreeSRedirector cstream("FitSignal.root");
180 TFile *f = cstream.GetFile();
181
182 char lname[100];
183 sprintf(lname,"Fit%s", cut.GetTitle());
184 TEventList *list = new TEventList(lname,lname);
185 sprintf(lname,">>Fit%s", cut.GetTitle());
186 treeB->Draw(lname,cut);
187 treeB->SetEventList(list);
663296cb 188 Int_t nFits=0;
eaad8c01 189 for (Int_t ievent=0; ievent<list->GetN(); ievent++){
663296cb 190 if (nFits>max) break;
191 if (nFits%50==0) printf("%d\n",nFits);
eaad8c01 192 char ename[100];
193 sprintf(ename,"Fit%d", ievent);
194 Double_t nsample = treeB->Draw("Graph.fY-Mean09:Graph.fX","","",1,ievent);
195 Double_t * signal = treeB->GetV1();
196 Double_t * time = treeB->GetV2();
197 Double_t maxpos =0;
198 Double_t max = 0;
199 for (Int_t ipos = 0; ipos<nsample; ipos++){
200 if (signal[ipos]>max){
201 max = signal[ipos];
202 maxpos = ipos;
203 }
204 }
205
206 Int_t first = TMath::Max(maxpos-10,0.);
207 Int_t last = TMath::Min(maxpos+60, nsample);
208 //
209 f->cd();
210 TH1F his(ename,ename,last-first,first,last);
211 for (Int_t ipos=0; ipos<last-first; ipos++){
212 his.SetBinContent(ipos+1,signal[ipos+first]);
213 }
214 treeB->Draw("Sector:Row:Pad","","",1,ievent);
215 Double_t sector = treeB->GetV1()[0];
216 Double_t row = treeB->GetV2()[0];
217 Double_t pad = treeB->GetV3()[0];
218 // TGraph graph(last-first,&time[first],&signal[first]);
219 f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
220 // TH1F * his = (TH1F*)graph.GetHistogram();
221 his.Fit(f1,"q");
222 his.Write(ename);
223 gPad->Clear();
224 his.Draw();
225 gPad->Update();
226 Double_t params[6];
227 for (Int_t ipar=0; ipar<6; ipar++) params[ipar] = f1->GetParameters()[ipar];
228 Double_t chi2 = TFitter::GetFitter()->Chisquare(6,params);
229 TMatrixD cov(6,6);
230 cov.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
231 //
232 // tail cancellation
233 //
234 Double_t x0[1000];
235 Double_t x1[1000];
236 Double_t x2[1000];
237 for (Int_t ipos=0; ipos<last-first; ipos++){
238 x0[ipos] = signal[ipos+first];
239 }
240 proc.TailCancelationALTRO1(x0,x1,0.85*0.339,0.09,last-first);
241 proc.TailCancelationALTRO1(x1,x2,0.85,0.789,last-first);
242 //
243 sprintf(ename,"Cancel1_%d", ievent);
244 TH1F his1(ename,ename,last-first,first,last);
245 for (Int_t ipos=0; ipos<last-first; ipos++){
246 his1.SetBinContent(ipos+1,x1[ipos]);
247 }
248 his1.Write(ename);
249 sprintf(ename,"Cancel2_%d", ievent);
250 TH1F his2(ename,ename,last-first,first,last);
251 for (Int_t ipos=0; ipos<last-first; ipos++){
252 his2.SetBinContent(ipos+1,x1[ipos]);
253 }
254 f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
255 his2.Fit(f1,"q");
256 his2.Write(ename);
257 Double_t params2[6];
258 for (Int_t ipar=0; ipar<6; ipar++) params2[ipar] = f1->GetParameters()[ipar];
259 Double_t chi22 = TFitter::GetFitter()->Chisquare(6,params2);
260 TMatrixD cov2(6,6);
261 cov2.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
262
263 TGraph gr0(last-first, &time[first],x0);
264 TGraph gr1(last-first, &time[first],x1);
265 TGraph gr2(last-first, &time[first],x2);
266 //
267 cstream<<"Fit"<<
268 "Sector="<<sector<<
269 "Row="<<row<<
270 "Pad="<<pad<<
271 "First="<<first<<
272 "Max="<<max<<
273 "MaxPos="<<maxpos<<
274 "chi2="<<chi2<<
275 "chi22="<<chi22<<
276 "Cov="<<&cov<<
277 "Cov2="<<&cov2<<
278 "gr0.="<<&gr0<<
279 "gr1.="<<&gr1<<
280 "gr2.="<<&gr2<<
281 "p0="<<params[0]<<
282 "p1="<<params[1]<<
283 "p2="<<params[2]<<
284 "p3="<<params[3]<<
285 "p4="<<params[4]<<
286 "p5="<<params[5]<<
287 "p02="<<params2[0]<<
288 "p12="<<params2[1]<<
289 "p22="<<params2[2]<<
290 "p32="<<params2[3]<<
291 "p42="<<params2[4]<<
292 "p52="<<params2[5]<<
293 "\n";
294 // delete his;
663296cb 295 nFits++;
eaad8c01 296 }
297
298}
a82f6396 299
300
663296cb 301void LaserCalib(TTreeSRedirector & cstream, TTree * chain, Float_t tmin, Float_t tmax, Float_t fraction){
302 //
303 //
304 //
ba6ac751 305 const Double_t kMaxDelta=10;
663296cb 306 AliTPCParamSR param;
307 param.Update();
308 TFile fce("TPCsignal.root");
309 TTree * treece =(TTree*)fce.Get("Signalce");
310 if (chain) treece=chain;
311 //
312 TBranch * brsector = treece->GetBranch("Sector");
313 TBranch * brpad = treece->GetBranch("Pad");
314 TBranch * brrow = treece->GetBranch("Row");
b78c817f 315 TBranch * brTimeStamp = treece->GetBranch("TimeStamp");
663296cb 316 //
317 TBranch * brtime = treece->GetBranch("Time");
318 TBranch * brrms = treece->GetBranch("RMS06");
319 TBranch * brmax = treece->GetBranch("Max");
320 TBranch * brsum = treece->GetBranch("Qsum");
321
322 Int_t sector, pad, row=0;
323 Double_t time=0, rms=0, qMax=0, qSum=0;
b78c817f 324 UInt_t timeStamp=0;
663296cb 325 brsector->SetAddress(&sector);
326 brrow->SetAddress(&row);
327 brpad->SetAddress(&pad);
b78c817f 328 brTimeStamp->SetAddress(&timeStamp);
663296cb 329
330 brtime->SetAddress(&time);
331 brrms->SetAddress(&rms);
332 brmax->SetAddress(&qMax);
333 brsum->SetAddress(&qSum);
334
335
336 brsector->GetEntry(0);
337 //
338 Int_t firstSector = sector;
339 Int_t lastSector = sector;
340 Int_t fentry = 0;
341 Int_t lentry = 0;
342 //
343 // find time offset for differnt events
344 //
345 Int_t count = 0;
346 Double_t padTimes[500000];
347 TRobustEstimator restim;
348 Double_t meanS[72], sigmaS[72];
349 Int_t firstS[72], lastS[72];
350 Double_t sectorsID[72];
351 for (Int_t isector=0;isector<72; isector++){
352 firstS[isector]=-1;
353 lastS[isector] =-1;
354 };
ba6ac751 355 TH1F hisT("hisT","hisT",100,tmin,tmax);
356 treece->Draw("Time>>hisT","");
357 Float_t cbin = hisT.GetBinCenter(hisT.GetMaximumBin());
663296cb 358
359 for (Int_t ientry=0; ientry<treece->GetEntriesFast(); ientry++){
360 treece->GetEvent(ientry);
361 //
362 if (sector!=lastSector && sector==firstSector){
363 //if (sector!=lastSector){
364 lentry = ientry;
b78c817f 365 TTimeStamp stamp(timeStamp);
366 stamp.Print();
663296cb 367 printf("\nEvent\t%d\tFirst\t%d\tLast\t%d\t%d\n",count, fentry, lentry, lentry-fentry);
368 //
369 //
370 Int_t ngood=0;
371 for (Int_t ientry2=fentry; ientry2<lentry; ientry2++){
372 // brtime->GetEvent(ientry2);
373 // brsector->GetEvent(ientry2);
374 treece->GetEvent(ientry2);
ba6ac751 375 if (time>tmin&&time<tmax && TMath::Abs(time-cbin)<kMaxDelta){
663296cb 376 padTimes[ngood]=time;
377 ngood++;
378 if (firstS[sector]<0) firstS[sector]= ngood;
379 if (firstS[sector]>=0) lastS[sector] = ngood;
380 }
381 }
382 //
383 //
384 Double_t mean,sigma;
ba6ac751 385 restim.EvaluateUni(ngood,padTimes,mean, sigma,int(float(ngood)*fraction));
663296cb 386 printf("Event\t%d\t%f\t%f\n",count, mean, sigma);
387 for (Int_t isector=0; isector<72; isector++){
388 sectorsID[isector]=sector;
389 if (firstS[isector]>=0 &&lastS[isector]>=0 && lastS[isector]>firstS[isector] ){
390 Int_t ngoodS = lastS[isector]-firstS[isector];
391 restim.EvaluateUni(ngoodS, &padTimes[firstS[isector]],meanS[isector],
ba6ac751 392 sigmaS[isector],int(float(ngoodS)*fraction));
663296cb 393 }
394 }
395 TGraph graphM(72,sectorsID,meanS);
396 TGraph graphS(72,sectorsID,sigmaS);
397 cstream<<"TimeS"<<
ba6ac751 398 "CBin="<<cbin<<
663296cb 399 "Event="<<count<<
ba6ac751 400 "GM="<<&graphM<<
401 "GS="<<&graphS<<
663296cb 402 "\n";
403
404
405 for (Int_t ientry2=fentry; ientry2<lentry-1; ientry2++){
406 treece->GetEvent(ientry2);
407 Double_t x = param.GetPadRowRadii(sector,row);
408 Int_t maxpad = AliTPCROC::Instance()->GetNPads(sector,row);
409 Double_t y = (pad - 2.5 - 0.5*maxpad)*param.GetPadPitchWidth(sector);
410 Double_t alpha = TMath::DegToRad()*(10.+20.*(sector%18));
411 Double_t gx = x*TMath::Cos(alpha)-y*TMath::Sin(alpha);
412 Double_t gy = y*TMath::Cos(alpha)+x*TMath::Sin(alpha);
413
414 Int_t npadS = lastS[sector]-firstS[sector];
415 cstream<<"Time"<<
416 "Event="<<count<<
b78c817f 417 "TimeStamp="<<timeStamp<<
ba6ac751 418 "CBin="<<cbin<<
663296cb 419 "x="<<x<<
420 "y="<<y<<
421 "gx="<<gx<<
422 "gy="<<gy<<
423 "Sector="<<sector<<
424 "Row="<<row<<
425 "Pad="<<pad<<
426 "Time="<<time<<
427 "RMS="<<rms<<
428 "Time0="<<mean<<
429 "Sigma0="<<sigma<<
430 "TimeS0="<<meanS[sector]<<
431 "SigmaS0="<<sigmaS[sector]<<
432 "npad0="<<ngood<<
433 "npadS="<<npadS<<
434 "Max="<<qMax<<
435 "Sum="<<qSum<<
436 "\n";
437 }
438 treece->GetEvent(ientry);
439 fentry = ientry;
440 count++;
441 for (Int_t isector=0;isector<72; isector++){
442 firstS[isector]=-1;
443 lastS[isector] =-1;
444 }
445 }
446 lastSector=sector;
447 }
448}
449
450
451
a82f6396 452
453TChain *MakeChainCL(Int_t first, Int_t last){
454 TChain *chaincl = new TChain("TreeR","TreeR");
455 //
456 char fname[100];
457 for (Int_t i=first;i<last; i++){
458 if (i>0) sprintf(fname,"TPC.RecPoints%d.root/Event%d/TreeR",i,i);
459 if (i==0) sprintf(fname,"TPC.RecPoints.root/Event%d/TreeR",i);
460 chaincl->Add(fname);
461 }
462 return chaincl;
463}
464
465TTree* GetTree(Int_t ievent){
466 char fname[100];
467 char tname[100];
468 if (ievent>0) sprintf(fname,"TPC.RecPoints%d.root",ievent);
469 if (ievent==0) sprintf(fname,"TPC.RecPoints.root");
470 sprintf(tname,"Event%d/TreeR",ievent);
471 TFile * f = new TFile(fname);
472 TTree * tree = (TTree*)f->Get(tname);
473 return tree;
474
475}