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