]>
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+ | |
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 | ||
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" | |
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 | 91 | void FitSignals(TTree * treeB, TCut cut="Max-Median>150&&RMS06<2&&abs(Median-Mean09)<0.5", Int_t maxS=100); |
ba6ac751 | 92 | void LaserCalib(TTreeSRedirector & cstream, TTree * chain, Float_t tmin, Float_t tmax, Float_t fraction=0.7); |
de31eda2 | 93 | |
a82f6396 | 94 | void 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 | 176 | void 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 | 301 | void 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(§or); |
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 | |
453 | TChain *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 | ||
465 | TTree* 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 | } |