TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGPP / TPC / macros / MakeReportTPC.C
1 /*
2    //gSystem->AddIncludePath("-I$ALICE_ROOT/PWGPP/TPC");
3    //.L $ALICE_ROOT/PWGPP/TPC/macros/MakeReportTPC.C+
4    //
5    //MakeReportTPC();
6
7    gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
8    gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
9    AliXRDPROOFtoolkit tool;
10    TChain * chain = tool.MakeChain("summaryTPCQA.txt","tpcQA",0,200000);
11
12    TTree *tree = chain->CopyTree("1");
13 */
14
15
16 #if !defined(__CINT__) || defined(__MAKECINT__)
17 #include "TSystem.h"
18 #include "TMath.h"
19 #include "TVectorD.h"
20 #include "TFile.h"
21 #include "TGrid.h"
22 #include "TF1.h"
23 #include "TH1.h"
24 #include "TH2.h"
25 #include "TProfile.h"
26 #include "THnSparse.h"
27 #include "TTreeStream.h"
28 #include "AliPerformanceTPC.h"
29 #include "AliPerformanceDEdx.h"
30 #endif
31
32 TObject *pTPCObject=0;
33 TObject *pTPCObjectGain=0;
34 TTreeSRedirector  *pcstream=0;
35 void Init();
36 void ReadObjects(const char * path=0);
37 void MakeReportTPC(Int_t run, const char *path=0 );
38 void AnalyzeNCL();
39 void AnalyzeDrift();
40 void AnalyzeGain();
41
42 void Init(){
43   //
44   //
45   //
46   gSystem->Load("libANALYSIS.so");
47   gSystem->Load("libANALYSISalice.so");
48   gSystem->Load("libTender.so");
49   gSystem->Load("libCORRFW.so");
50   gSystem->Load("libPWG0base.so");
51   gSystem->Load("libPWG0dep.so");
52   gSystem->Load("libPWG0selectors.so");
53   gSystem->Load("libPWGPP.so");
54   gSystem->Load("libTPCcalib"); 
55 //   gSystem->Setenv("alien_CLOSE_SE","ALICE::GSI::SE");
56 //   TGrid::Connect("alien://",0,0,"t"); 
57 //   gSystem->Setenv("alien_CLOSE_SE","ALICE::GSI::SE");    
58 }
59
60
61 void ReadObjects(const char * path){
62   //
63   //
64   //
65   TFile *f =0;
66   if (path) f=TFile::Open(Form("%s/TPC.PerformanceQAQA.root",path));
67   if (!path) f=TFile::Open("TPC.Performance.root");
68   if (!f){
69     printf("File %s not available\n", path);
70     return;
71   }
72   TList *list = (TList*)f->Get("TPC");
73   if (!list){
74     printf("QA %s not available\n", path);
75     return;
76   }
77   pTPCObject= (AliPerformanceTPC*)list->FindObject("AliPerformanceTPC");
78   pTPCObjectGain = (AliPerformanceDEdx*)list->FindObject("AliPerformanceDEdxTPCInner");
79 }
80
81
82 void MakeReportTPC(Int_t run, const char *path ){
83   //
84   // make a tpcQA report
85   //  typical variables exported to the tree
86   //
87   Init();
88   ReadObjects(path);
89   
90   pcstream= new TTreeSRedirector("TPC.PerformanceSummary.root");
91   (*pcstream)<<"tpcQA"<<"run="<<run;
92   AnalyzeNCL();
93   AnalyzeDrift();
94   AnalyzeGain();
95   (*pcstream)<<"tpcQA"<<"\n";
96   delete pcstream;
97 }
98
99
100 void AnalyzeNCL(){
101   //
102   // NCL statistic
103   //    
104   // variables:
105   static Double_t meanTPCnclF=0;
106   static Double_t rmsTPCnclF=0;
107   static Double_t slopeATPCnclF=0;
108   static Double_t slopeCTPCnclF=0;
109   static Double_t slopeATPCnclFErr=0;
110   static Double_t slopeCTPCnclFErr=0;
111   static Double_t meanTPCncl=0;
112   static Double_t rmsTPCncl=0;
113   static Double_t slopeATPCncl=0;
114   static Double_t slopeCTPCncl=0;
115   static Double_t slopeATPCnclErr=0;
116   static Double_t slopeCTPCnclErr=0;
117   AliPerformanceTPC * pTPC=  (AliPerformanceTPC *)pTPCObject;
118
119   TH1* his1D=0;
120   TProfile* hprof=0;
121   static TF1 *fpol1 = new TF1("fpol1","pol1");
122   //
123   // all clusters
124   // eta cut - +-1
125   // pt cut  - +-0.250 GeV
126   pTPC->GetTPCTrackHisto()->GetAxis(5)->SetRangeUser(-1.,1.);
127   pTPC->GetTPCTrackHisto()->GetAxis(7)->SetRangeUser(0.25,10);
128   his1D = pTPC->GetTPCTrackHisto()->Projection(0);
129   meanTPCncl= his1D->GetMean();
130   rmsTPCncl= his1D->GetRMS();
131   delete his1D;
132   hprof = pTPC->GetTPCTrackHisto()->Projection(0,5)->ProfileX();
133   hprof->Fit(fpol1,"QNR","QNR",0,0.8);
134   slopeATPCncl= fpol1->GetParameter(1);
135   slopeATPCnclErr= fpol1->GetParError(1);
136   hprof->Fit(fpol1,"QNR","QNR",-0.8,0.0);
137   slopeCTPCncl= fpol1->GetParameter(1);
138   slopeCTPCnclErr= fpol1->GetParameter(1);
139   delete hprof;
140   //
141   // findable clusters
142   //
143   pTPC->GetTPCTrackHisto()->GetAxis(2)->SetRangeUser(0.4,1.1);
144   his1D = pTPC->GetTPCTrackHisto()->Projection(2);
145   meanTPCnclF= his1D->GetMean();
146   rmsTPCnclF= his1D->GetRMS();
147   delete his1D;
148   his1D = pTPC->GetTPCTrackHisto()->Projection(2,5)->ProfileX();
149   his1D->Fit(fpol1,"QNR","QNR",0,0.8);
150   slopeATPCnclF= fpol1->GetParameter(1);
151   slopeATPCnclFErr= fpol1->GetParError(1);
152   his1D->Fit(fpol1,"QNR","QNR",-0.8,0.0);
153   slopeCTPCnclF= fpol1->GetParameter(1);
154   slopeCTPCnclFErr= fpol1->GetParameter(1);
155   delete his1D;
156   printf("Cluster QA report\n");
157   printf("meanTPCnclF=\t%f\n",meanTPCnclF);
158   printf("rmsTPCnclF=\t%f\n",rmsTPCnclF);
159   printf("slopeATPCnclF=\t%f\n",slopeATPCnclF);
160   printf("slopeCTPCnclF=\t%f\n",slopeCTPCnclF);
161   printf("meanTPCncl=\t%f\n",meanTPCncl);
162   printf("rmsTPCncl=\t%f\n",rmsTPCncl);
163   printf("slopeATPCncl=\t%f\n",slopeATPCncl);
164   printf("slopeCTPCncl=\t%f\n",slopeCTPCncl);
165   //
166   // dump results to the tree
167   //
168   (*pcstream)<<"tpcQA"<<
169     "meanTPCnclF="<<meanTPCnclF <<   
170     "rmsTPCnclF="<<rmsTPCnclF <<
171     "slopeATPCnclF="<< slopeATPCnclF<<
172     "slopeCTPCnclF="<< slopeCTPCnclF<<
173     "slopeATPCnclFErr="<< slopeATPCnclFErr<<
174     "slopeCTPCnclFErr="<< slopeCTPCnclFErr<<
175     "meanTPCncl="<<meanTPCncl <<
176     "rmsTPCncl="<< rmsTPCncl<<
177     "slopeATPCncl="<< slopeATPCncl<<
178     "slopeCTPCncl="<< slopeCTPCncl<<
179     "slopeATPCnclErr="<< slopeATPCnclErr<<
180     "slopeCTPCnclErr="<< slopeCTPCnclErr;
181
182     
183 }
184
185 void AnalyzeDrift(){
186   //
187   // Analyze drift velocity imperfections
188   //
189   // variables:
190   static Double_t offsetdZA=0;
191   static Double_t slopedZA=0;
192   static Double_t offsetdZC=0;
193   static Double_t slopedZC=0;
194   static Double_t offsetdZAErr=0;
195   static Double_t slopedZAErr=0;
196   static Double_t offsetdZCErr=0;
197   static Double_t slopedZCErr=0;
198   static Double_t offsetdZAchi2=0;
199   static Double_t slopedZAchi2=0;
200   static Double_t offsetdZCchi2=0;
201   static Double_t slopedZCchi2=0;
202   AliPerformanceTPC * pTPC=  (AliPerformanceTPC *)pTPCObject;
203
204   TH1* his1D=0;
205   TH2* his2D=0;
206   static TF1 *fpol1 = new TF1("fpol1","pol1");
207   TObjArray arrayFit;
208   his2D = pTPC->GetTPCTrackHisto()->Projection(4,5);
209   his2D->FitSlicesY(0,0,-1,10,"QNR",&arrayFit);
210   delete his2D;
211   his1D = (TH1*) arrayFit.At(1);
212   his1D->Fit(fpol1,"QNRROB=0.8","QNR",-0.8,-0.1);
213   offsetdZC=fpol1->GetParameter(0);
214   slopedZC=fpol1->GetParameter(1);
215   //offsetdZCchi2=fpol1->GetChisquare();
216   slopedZCchi2=fpol1->GetChisquare();
217   //
218   his1D->Fit(fpol1,"QNRROB=0.8","QNR",0.1,0.8);
219   offsetdZA=fpol1->GetParameter(0);
220   slopedZA=fpol1->GetParameter(1);
221   offsetdZAErr=fpol1->GetParError(0);
222   slopedZAErr=fpol1->GetParError(1);
223   offsetdZAchi2=fpol1->GetChisquare();
224   slopedZAchi2=fpol1->GetChisquare();
225   //
226   printf("Drift velocity QA report\n");
227   printf("offsetdZA\t%f\n",offsetdZA);
228   printf("slopedZA\t%f\n",slopedZA);
229   printf("offsetdZC\t%f\n",offsetdZC);
230   printf("slopedZC\t%f\n",slopedZC);
231   //
232   // dump drift QA values
233   //
234   (*pcstream)<<"tpcQA"<<
235     "offsetdZA="<< offsetdZA<<
236     "slopedZA="<< slopedZA<<
237     "offsetdZC="<< offsetdZC<<
238     "slopedZC="<<slopedZC<<
239     //
240     "offsetdZAErr="<< offsetdZAErr<<
241     "slopedZAErr="<< slopedZAErr<<
242     "offsetdZCErr="<< offsetdZCErr<<
243     "slopedZCErr="<<slopedZCErr<<
244     //
245     "offsetdZAchi2="<< offsetdZAchi2<<
246     "slopedZAchi2="<< slopedZAchi2<<
247     "offsetdZCchi2="<< offsetdZCchi2<<
248     "slopedZCchi2="<<slopedZCchi2;
249   
250 }
251
252
253 void AnalyzeGain(){
254   //
255   // gain statistic
256   //
257   AliPerformanceDEdx * pTPCgain =  (AliPerformanceDEdx *)pTPCObjectGain;
258   static TVectorD MIPvsSector(36);
259   static TVectorD sector(36);
260   static Float_t MIPmean = 0;
261   static Float_t MIPresolution = 0;
262   static Float_t attachSlopeC = 0;
263   static Float_t attachSlopeA = 0;
264
265   MIPvsSector.Zero();
266   //
267   // select MIP particles
268   //
269   if (pTPCgain){
270     pTPCgain->GetDeDxHisto()->GetAxis(7)->SetRangeUser(0.4,0.55);
271     pTPCgain->GetDeDxHisto()->GetAxis(0)->SetRangeUser(35,60);
272     pTPCgain->GetDeDxHisto()->GetAxis(6)->SetRangeUser(80,160);
273     pTPCgain->GetDeDxHisto()->GetAxis(5)->SetRangeUser(-1,1);
274     //
275     // MIP position and resolution
276     //
277     TF1 gausFit("gausFit","gaus");
278     TH1 * hisProj1D = pTPCgain->GetDeDxHisto()->Projection(0);
279     hisProj1D->Fit(&gausFit,"QN","QN");
280     delete hisProj1D;
281     MIPmean = gausFit.GetParameter(1);
282     MIPresolution = 0;
283     if (MIPmean!=0) MIPresolution = gausFit.GetParameter(2)/MIPmean;
284     //
285     // MIP position vs. dip angle (attachment)
286     //
287     pTPCgain->GetDeDxHisto()->GetAxis(5)->SetRangeUser(-3,0);
288     TH2* his2D = pTPCgain->GetDeDxHisto()->Projection(0,5);
289     TF1 * fpol = new TF1("fpol","pol1");
290     TObjArray arrayFit;
291     his2D->FitSlicesY(0,0,-1,10,"QN",&arrayFit);
292     delete his2D;
293     TH1 * his1D = (TH1*) arrayFit.At(1);
294     his1D->Fit(fpol,"QNROB=0.8","QNR",-1,0);
295     attachSlopeC = fpol->GetParameter(1);
296     //
297     pTPCgain->GetDeDxHisto()->GetAxis(5)->SetRangeUser(0,3);
298     TH2* his2DA = pTPCgain->GetDeDxHisto()->Projection(0,5);
299     TF1 * fpolA = new TF1("fpolA","pol1");
300     TObjArray arrayFitA;
301     his2DA->FitSlicesY(0,0,-1,10,"QN",&arrayFit);
302     delete his2DA;
303     TH1 * his1DA = (TH1*) arrayFit.At(1);
304     his1DA->Fit(fpolA,"QNROB=0.8","QN",0,1);
305     attachSlopeA = fpolA->GetParameter(1);
306     //
307     // MIP position vs. sector
308     //
309     pTPCgain->GetDeDxHisto()->GetAxis(5)->SetRangeUser(-3,0); // C side
310     for(Int_t i = 0; i < 18; i++) { // loop over sectors; correct mapping to be checked!
311       TH1* his1D=0;
312       Float_t phiLow = -TMath::Pi() + i*(20./360.)*(2*TMath::Pi());
313       Float_t phiUp  = -TMath::Pi() + (i+1)*(20./360.)*(2*TMath::Pi());
314       pTPCgain->GetDeDxHisto()->GetAxis(1)->SetRangeUser(phiLow,phiUp);
315       his1D = pTPCgain->GetDeDxHisto()->Projection(0);
316       TF1 gausFunc("gausFunc","gaus");
317       his1D->Fit(&gausFunc, "QN");
318       MIPvsSector(i) = gausFunc.GetParameter(1);
319       sector(i)=i;
320       delete his1D;
321     }
322     //
323     pTPCgain->GetDeDxHisto()->GetAxis(5)->SetRangeUser(0,3); // A side
324     for(Int_t i = 0; i < 18; i++) { // loop over sectors; correct mapping to be checked!
325       TH1* his1D=0;
326       Float_t phiLow = -TMath::Pi() + i*(20./360.)*(2*TMath::Pi());
327       Float_t phiUp  = -TMath::Pi() + (i+1)*(20./360.)*(2*TMath::Pi());
328       pTPCgain->GetDeDxHisto()->GetAxis(1)->SetRangeUser(phiLow,phiUp);
329       his1D = pTPCgain->GetDeDxHisto()->Projection(0);
330       TF1 gausFunc("gausFunc","gaus");
331       his1D->Fit(&gausFunc, "QN");
332       MIPvsSector(i+18) = gausFunc.GetParameter(1);
333       sector(i+18)=i+18;
334       delete his1D;
335     }
336   }
337   //
338   printf("Gain QA report\n");
339   printf("MIP mean\t%f\n",MIPmean);
340   printf("MIP resolution\t%f\n",MIPresolution);
341   printf("MIPslopeA\t%f\n",attachSlopeA);
342   printf("MIPslopeC\t%f\n",attachSlopeC);
343   //
344   (*pcstream)<<"tpcQA"<<
345     "MIPattachSlopeC="<<attachSlopeC<<
346     "MIPattachSlopeA="<<attachSlopeA<<
347     "MIPresolution="<<MIPresolution<<
348     "MIPvsSector.="<<&MIPvsSector<<
349     "sector.="<<&sector<<
350     "MIPmean="<<MIPmean;
351
352 }