Simple example macro to make some analyzis for the ESD tracks - see comments (Not...
[u/mrichter/AliRoot.git] / TPC / AnalyzeESDtracks.C
1 //
2 // Process ESD tracks  - 
3 // Extract TPC tracks  - write them to tree
4 //
5 /*
6   .L AnalyzeESDtracks.C+
7   .L AliGenInfo.C+
8   AnalyzESDtracks(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 //AnalyzESDtracks(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 "TBranch.h"
62 #include "TTreeStream.h"
63 #include "TEventList.h"
64 #include "TCut.h"
65 #include "TFitter.h"
66 #include  "TMatrixD.h"
67
68 #include  "AliLog.h"
69 #include  "AliMagF.h"
70
71 #include  "AliESD.h"
72 #include  "AliESDfriend.h"
73 #include  "AliESDtrack.h"
74 #include  "AliTracker.h"
75 #include  "AliTPCseed.h"
76 #include  "AliTPCclusterMI.h"
77 #include  "TTreeStream.h"
78 #include  "TF1.h"
79 #include  "TGraph.h"
80 #include  "AliSignalProcesor.h"
81 #include  "TCanvas.h"
82
83
84 void AnalyzESDtracks(Int_t run){
85   //
86   // output redirect 
87   TTreeSRedirector  cstream("TPCtracks.root");
88   //
89   // dummy magnetic field
90   AliMagF mag("aaa","aaa",1,1,10);
91   AliTracker::SetFieldMap(&mag,kTRUE);
92   TFile f("AliESDs.root");
93   TTree * tree =(TTree*)f.Get("esdTree"); 
94   AliESD * esd =0;
95   tree->SetBranchAddress("ESD",&esd);
96   AliESDfriend *evf=0;
97   tree->AddFriend("esdFriendTree","AliESDfriends.root");
98   tree->SetBranchAddress("ESDfriend",&evf);
99
100   Int_t nevents = tree->GetEntries();
101   TClonesArray *clusters = new TClonesArray("AliTPCclusterMI",160);
102   for (Int_t irow=0; irow<160; irow++){
103     new ((*clusters)[irow])  AliTPCclusterMI;   // iitial dummy clusters
104   }
105   
106   for (Int_t ievent=0; ievent<nevents; ievent++){
107     tree->GetEntry(ievent);
108     if (!esd) continue;
109     if (!evf) continue;
110     esd->SetESDfriend(evf); //Attach the friend to the ESD      
111     for (Int_t itrack =0; itrack<esd->GetNumberOfTracks(); itrack++){
112       // Int_t itrack = 0;      
113       if (esd->GetTrack(itrack)->GetFriendTrack()==0) continue;
114       AliESDtrack * etrack = esd->GetTrack(itrack);
115       AliESDfriendTrack * ftrack = (AliESDfriendTrack *)esd->GetTrack(itrack)->GetFriendTrack();
116       AliTPCseed * seed =  (AliTPCseed*)(ftrack->GetCalibObject(0));
117       if (!seed) continue;  
118       for (Int_t irow=0; irow<160; irow++){
119         if (seed->GetClusterFast(irow)){
120           AliTPCclusterMI * cl = new ((*clusters)[irow])  AliTPCclusterMI(*(seed->GetClusterFast(irow)));
121           cl->SetLabel(itrack,0);
122         }
123         else{
124           AliTPCclusterMI * cl = (AliTPCclusterMI*)clusters->At(irow);
125           cl->SetX(0); cl->SetY(0); cl->SetZ(0); cl->SetQ(0); cl->SetLabel(-1,0);
126         }
127       }
128       Float_t dEdx = seed->GetdEdx();
129       Float_t dEdxI = seed->CookdEdx(0.05,0.6,0,77);
130       Float_t dEdxO = seed->CookdEdx(0.05,0.6,78,155);
131       Int_t ncl = seed->GetNumberOfClusters();
132       cstream<<"Tracks"<<
133         "Run="<<run<<
134         "Event="<<ievent<<
135         "dEdx="<<dEdx<<
136         "dEdxI="<<dEdxI<<
137         "dEdxO="<<dEdxO<<
138         "Track.="<<seed<<
139         "Etrack.="<<etrack<<
140         "Cl.="<<clusters<<
141         "\n";
142     }  
143   }
144 }
145
146
147 void FitSignals(TTree * treeB, TCut cut="Max-Median>150&&RMS06<2&&abs(Median-Mean09)<0.5"){
148   AliSignalProcesor proc;
149   TF1 * f1 = proc.GetAsymGauss();
150   TTreeSRedirector cstream("FitSignal.root");
151   TFile *f = cstream.GetFile();
152
153   char lname[100];
154   sprintf(lname,"Fit%s", cut.GetTitle());
155   TEventList *list = new TEventList(lname,lname);
156   sprintf(lname,">>Fit%s", cut.GetTitle());
157   treeB->Draw(lname,cut);
158   treeB->SetEventList(list);
159   for (Int_t ievent=0; ievent<list->GetN(); ievent++){
160     char ename[100];
161     sprintf(ename,"Fit%d", ievent);
162     Double_t nsample = treeB->Draw("Graph.fY-Mean09:Graph.fX","","",1,ievent);
163     Double_t * signal  = treeB->GetV1();
164     Double_t * time  = treeB->GetV2();
165     Double_t maxpos =0;
166     Double_t max = 0;
167     for (Int_t ipos = 0; ipos<nsample; ipos++){
168       if (signal[ipos]>max){
169         max    = signal[ipos];
170         maxpos = ipos;
171       }
172     }
173
174     Int_t first = TMath::Max(maxpos-10,0.);
175     Int_t last  = TMath::Min(maxpos+60, nsample);
176     //
177     f->cd();
178     TH1F his(ename,ename,last-first,first,last);
179     for (Int_t ipos=0; ipos<last-first; ipos++){
180       his.SetBinContent(ipos+1,signal[ipos+first]);
181     }
182     treeB->Draw("Sector:Row:Pad","","",1,ievent);
183     Double_t sector = treeB->GetV1()[0];
184     Double_t row    = treeB->GetV2()[0];
185     Double_t pad    = treeB->GetV3()[0];
186     //    TGraph  graph(last-first,&time[first],&signal[first]);
187     f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
188     //    TH1F * his = (TH1F*)graph.GetHistogram();
189     his.Fit(f1,"q");
190     his.Write(ename);
191     gPad->Clear();
192     his.Draw();
193     gPad->Update();
194     Double_t params[6];
195     for (Int_t ipar=0; ipar<6; ipar++) params[ipar] = f1->GetParameters()[ipar];
196     Double_t chi2 = TFitter::GetFitter()->Chisquare(6,params);
197     TMatrixD cov(6,6);
198     cov.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
199     //
200     // tail cancellation
201     //
202     Double_t x0[1000];
203     Double_t x1[1000];
204     Double_t x2[1000];
205     for (Int_t ipos=0; ipos<last-first; ipos++){
206       x0[ipos] = signal[ipos+first];
207     }
208     proc.TailCancelationALTRO1(x0,x1,0.85*0.339,0.09,last-first);
209     proc.TailCancelationALTRO1(x1,x2,0.85,0.789,last-first);
210     //
211     sprintf(ename,"Cancel1_%d", ievent);
212     TH1F his1(ename,ename,last-first,first,last);
213     for (Int_t ipos=0; ipos<last-first; ipos++){
214       his1.SetBinContent(ipos+1,x1[ipos]);
215     }
216     his1.Write(ename);
217     sprintf(ename,"Cancel2_%d", ievent);
218     TH1F his2(ename,ename,last-first,first,last);
219     for (Int_t ipos=0; ipos<last-first; ipos++){
220       his2.SetBinContent(ipos+1,x1[ipos]);
221     }
222     f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
223     his2.Fit(f1,"q");
224     his2.Write(ename);
225     Double_t params2[6];
226     for (Int_t ipar=0; ipar<6; ipar++) params2[ipar] = f1->GetParameters()[ipar];
227     Double_t chi22 = TFitter::GetFitter()->Chisquare(6,params2);    
228     TMatrixD cov2(6,6);
229     cov2.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
230
231     TGraph gr0(last-first, &time[first],x0);
232     TGraph gr1(last-first, &time[first],x1);
233     TGraph gr2(last-first, &time[first],x2);
234     //
235     cstream<<"Fit"<<
236       "Sector="<<sector<<
237       "Row="<<row<<
238       "Pad="<<pad<<
239       "First="<<first<<
240       "Max="<<max<<
241       "MaxPos="<<maxpos<<
242       "chi2="<<chi2<<
243       "chi22="<<chi22<<
244       "Cov="<<&cov<<
245       "Cov2="<<&cov2<<
246       "gr0.="<<&gr0<<
247       "gr1.="<<&gr1<<
248       "gr2.="<<&gr2<<
249       "p0="<<params[0]<<
250       "p1="<<params[1]<<
251       "p2="<<params[2]<<
252       "p3="<<params[3]<<
253       "p4="<<params[4]<<
254       "p5="<<params[5]<<
255       "p02="<<params2[0]<<
256       "p12="<<params2[1]<<
257       "p22="<<params2[2]<<
258       "p32="<<params2[3]<<
259       "p42="<<params2[4]<<
260       "p52="<<params2[5]<<
261       "\n";
262     //    delete his;
263   }
264
265 }