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