1 /// \file AnalyzeESDtracks.C
3 /// \brief Process ESD tracks - extract TPC tracks - write them to tree
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;
23 .L AnalyzeESDtracks.C+
26 //AnalyzeESDtracks(567);
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;
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");
43 FitSignals(treeB,"Max-Median>100&&RMS06<2.5")
44 TFile ffit("FitSignal.root");
45 TTree * treeF = (TTree*)ffit->Get("Fit");
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")
63 #include "TTreeStream.h"
64 #include "TEventList.h"
68 #include "TRobustEstimator.h"
69 #include "TTimeStamp.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"
84 #include "TTreeStream.h"
87 #include "AliSignalProcesor.h"
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);
94 void AnalyzeESDtracks(Int_t run){
96 TTreeSRedirector * pcstream = new TTreeSRedirector("TPCtracks.root");
97 TTreeSRedirector &cstream = *pcstream;
99 TFile f("AliESDs.root");
100 TTree * tree =(TTree*)f.Get("esdTree");
102 tree->SetBranchAddress("ESD",&esd);
104 tree->AddFriend("esdFriendTree","AliESDfriends.root");
105 tree->SetBranchAddress("ESDfriend",&evf);
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
113 for (Int_t ievent=0; ievent<nevents; ievent++){
114 tree->GetEntry(ievent);
117 esd->SetESDfriend(evf); //Attach the friend to the ESD
118 for (Int_t itrack =0; itrack<esd->GetNumberOfTracks(); itrack++){
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));
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);
131 AliTPCclusterMI * cl = (AliTPCclusterMI*)clusters->At(irow);
132 cl->SetX(0); cl->SetY(0); cl->SetZ(0); cl->SetQ(0); cl->SetLabel(-1,0);
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();
156 TFile fs("TPCsignal.root");
157 TTree *treeB =(TTree*)fs.Get("SignalB");
159 // Fit central electrode part
161 TTreeSRedirector * pcestream = new TTreeSRedirector("TimeRoot.root");
162 TTree * treece = (TTree*)fs.Get("Signalce");
164 LaserCalib(*pcestream, treece, 800,1000, 0.7);
167 FitSignals(treeB,"Max-Median>150&&RMS06<1.0&&RMS09<1.5&&abs(Median-Mean09)<0.2&&abs(Mean06-Mean09)<0.2",1000);
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();
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);
185 for (Int_t ievent=0; ievent<list->GetN(); ievent++){
186 if (nFits>max) break;
187 if (nFits%50==0) printf("%d\n",nFits);
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();
195 for (Int_t ipos = 0; ipos<nsample; ipos++){
196 if (signal[ipos]>max){
202 Int_t first = TMath::Max(maxpos-10,0.);
203 Int_t last = TMath::Min(maxpos+60, nsample);
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]);
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();
223 for (Int_t ipar=0; ipar<6; ipar++) params[ipar] = f1->GetParameters()[ipar];
224 Double_t chi2 = TFitter::GetFitter()->Chisquare(6,params);
226 cov.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
233 for (Int_t ipos=0; ipos<last-first; ipos++){
234 x0[ipos] = signal[ipos+first];
236 proc.TailCancelationALTRO1(x0,x1,0.85*0.339,0.09,last-first);
237 proc.TailCancelationALTRO1(x1,x2,0.85,0.789,last-first);
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]);
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]);
250 f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
254 for (Int_t ipar=0; ipar<6; ipar++) params2[ipar] = f1->GetParameters()[ipar];
255 Double_t chi22 = TFitter::GetFitter()->Chisquare(6,params2);
257 cov2.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
259 TGraph gr0(last-first, &time[first],x0);
260 TGraph gr1(last-first, &time[first],x1);
261 TGraph gr2(last-first, &time[first],x2);
297 void LaserCalib(TTreeSRedirector & cstream, TTree * chain, Float_t tmin, Float_t tmax, Float_t fraction){
300 const Double_t kMaxDelta=10;
303 TFile fce("TPCsignal.root");
304 TTree * treece =(TTree*)fce.Get("Signalce");
305 if (chain) treece=chain;
307 TBranch * brsector = treece->GetBranch("Sector");
308 TBranch * brpad = treece->GetBranch("Pad");
309 TBranch * brrow = treece->GetBranch("Row");
310 TBranch * brTimeStamp = treece->GetBranch("TimeStamp");
312 TBranch * brtime = treece->GetBranch("Time");
313 TBranch * brrms = treece->GetBranch("RMS06");
314 TBranch * brmax = treece->GetBranch("Max");
315 TBranch * brsum = treece->GetBranch("Qsum");
317 Int_t sector, pad, row=0;
318 Double_t time=0, rms=0, qMax=0, qSum=0;
320 brsector->SetAddress(§or);
321 brrow->SetAddress(&row);
322 brpad->SetAddress(&pad);
323 brTimeStamp->SetAddress(&timeStamp);
325 brtime->SetAddress(&time);
326 brrms->SetAddress(&rms);
327 brmax->SetAddress(&qMax);
328 brsum->SetAddress(&qSum);
331 brsector->GetEntry(0);
333 Int_t firstSector = sector;
334 Int_t lastSector = sector;
338 // find time offset for differnt events
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++){
350 TH1F hisT("hisT","hisT",100,tmin,tmax);
351 treece->Draw("Time>>hisT","");
352 Float_t cbin = hisT.GetBinCenter(hisT.GetMaximumBin());
354 for (Int_t ientry=0; ientry<treece->GetEntriesFast(); ientry++){
355 treece->GetEvent(ientry);
357 if (sector!=lastSector && sector==firstSector){
358 //if (sector!=lastSector){
360 TTimeStamp stamp(timeStamp);
362 printf("\nEvent\t%d\tFirst\t%d\tLast\t%d\t%d\n",count, fentry, lentry, lentry-fentry);
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;
373 if (firstS[sector]<0) firstS[sector]= ngood;
374 if (firstS[sector]>=0) lastS[sector] = ngood;
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));
390 TGraph graphM(72,sectorsID,meanS);
391 TGraph graphS(72,sectorsID,sigmaS);
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);
409 Int_t npadS = lastS[sector]-firstS[sector];
412 "TimeStamp="<<timeStamp<<
425 "TimeS0="<<meanS[sector]<<
426 "SigmaS0="<<sigmaS[sector]<<
433 treece->GetEvent(ientry);
436 for (Int_t isector=0;isector<72; isector++){
448 TChain *MakeChainCL(Int_t first, Int_t last){
449 TChain *chaincl = new TChain("TreeR","TreeR");
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);
460 TTree* GetTree(Int_t ievent){
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);