]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/Upgrade/macros/spaceChargeFluctuation.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding function to fit q(z) normalization parameters of distortion maps
[u/mrichter/AliRoot.git] / TPC / Upgrade / macros / spaceChargeFluctuation.C
CommitLineData
c0172f82 1/*
7212f8c3 2 Macro to study the space charge fluctuations.
3 3 functions using the ToyMC + analytical fomula to describe given MC results
4
5 function to histogram space charge using the raw data ana anlyzing them
6 To use given function - CPU conusming therefore batch farms used
7 See $ALICE_ROOT/TPC/Upgrade/macros/spaceChargeFluctuation.sh macro to see example to run the code
8
9
10
11 .x $HOME/NimStyle.C
ae4cea45 12 .x $HOME/rootlogon.C
c0172f82 13 .L $ALICE_ROOT/TPC/Upgrade/macros/spaceChargeFluctuation.C+
7212f8c3 14
2ff913b6 15
c0172f82 16
17 */
18#include "TMath.h"
19#include "TRandom.h"
20#include "TTreeStream.h"
21#include "TVectorD.h"
22#include "TCanvas.h"
23#include "TStopwatch.h"
24#include "AliTPCParam.h"
25#include "AliTPCcalibDB.h"
26#include "AliTPCAltroMapping.h"
27#include "AliAltroRawStream.h"
28#include "AliSysInfo.h"
29#include "AliTPCRawStreamV3.h"
30#include "AliCDBManager.h"
31#include "TGeoGlobalMagField.h"
32#include "AliMagF.h"
33#include "AliRawReaderRoot.h"
34#include "AliRawReader.h"
35#include "TH3.h"
36#include "TH2.h"
76dbc35e 37#include "AliTPCCalPad.h"
38#include "AliTPCCalROC.h"
39#include "TChain.h"
40#include "AliXRDPROOFtoolkit.h"
41#include "TLegend.h"
42#include "TCut.h"
0eff8c53 43#include "TGraphErrors.h"
44#include "TStatToolkit.h"
c0172f82 45
ae4cea45 46#include "AliDCSSensor.h"
47#include "AliCDBEntry.h"
48#include "AliDCSSensorArray.h"
7212f8c3 49#include "TStyle.h"
50#include "AliTPCSpaceCharge3D.h"
4b955de3 51#include "AliExternalTrackParam.h"
633a08cb 52#include "AliTrackerBase.h"
53#include "TDatabasePDG.h"
511fd208 54#include "TROOT.h"
4b955de3 55//
56// constants
57//
58Double_t omegaTau=0.325;
7212f8c3 59//
60// Function declaration
61//
62// TOY MC
63void spaceChargeFluctuationToyMC(Int_t nframes, Double_t interactionRate);
64void spaceChargeFluctuationToyDraw();
65void spaceChargeFluctuationToyDrawSummary();
c0172f82 66
7212f8c3 67//
68// RAW data analysis
69//
76dbc35e 70TH1 * GenerateMapRawIons(Int_t useGain,const char *fileName="raw.root", const char *outputName="histo.root", Int_t maxEvents=25);
0eff8c53 71void DoMerge();
7212f8c3 72void AnalyzeMaps1D(); // make nice plots
73void MakeFluctuationStudy3D(Int_t nhistos, Int_t nevents, Int_t niter);
74TH3D * NormalizeHistoQ(TH3D * hisInput, Bool_t normEpsilon);
4b955de3 75//
7212f8c3 76TH3D * PermutationHistoZ(TH3D * hisInput, Double_t deltaZ);
77TH3D * PermutationHistoPhi(TH3D * hisInput, Double_t deltaPhi);
4b955de3 78TH3D * PermutationHistoLocalPhi(TH3D * hisInput, Int_t deltaPhi);
79void MakeSpaceChargeFluctuationScan(Double_t scale, Int_t nfiles, Int_t sign);
7212f8c3 80void DrawFluctuationdeltaZ(Int_t stat=0, Double_t norm=10000);
81void DrawFluctuationSector(Int_t stat=0, Double_t norm=10000);
5a998ace 82void MakeLocalDistortionPlotsGlobalFitPolDrift(Float_t xmin, Float_t xmax);
83void MakeLocalDistortionPlots(Int_t npoints=20000, Int_t npointsZ=10000);
84
c0172f82 85
4b955de3 86void spaceChargeFluctuation(Int_t mode=0, Float_t arg0=0, Float_t arg1=0, Float_t arg2=0){
c0172f82 87 //
7212f8c3 88 // function called from the shell script
c0172f82 89 //
7212f8c3 90 gRandom->SetSeed(0);
76dbc35e 91 if (mode==0) GenerateMapRawIons(arg0);
0eff8c53 92 if (mode==1) DoMerge();
7212f8c3 93 if (mode==2) spaceChargeFluctuationToyMC(arg0,arg1);
94 if (mode==3) MakeFluctuationStudy3D(10000, arg0, arg1);
4b955de3 95 if (mode==4) MakeSpaceChargeFluctuationScan(arg0,arg1,arg2); // param: scale, nfiles, sign Bz
7212f8c3 96 if (mode==5) {
97 DrawFluctuationdeltaZ(arg0,arg1);
98 DrawFluctuationSector(arg0,arg1);
99 }
5a998ace 100 if (mode==6) {
101 MakeLocalDistortionPlotsGlobalFitPolDrift(arg0,arg1);
102 }
103 if (mode==7) {
104 MakeLocalDistortionPlots(arg0,arg1);
105
106 }
c0172f82 107}
06ab15a8 108
109
110Double_t RndmdNchdY(Double_t s){
111 //
112 // dNch/deta - last 2 points inventeted (to find it somewhere ?)
113 //
114 // http://arxiv.org/pdf/1012.1657v2.pdf - table 1. ALICE PbPb
115 // Scaled according s^0.15
116 // http://arxiv.org/pdf/1210.3615v2.pdf
117 // This we can cite.
118 // Usage example::
119 /*
120 TH1F his550("his550","his550",1000,0,3000)
121 for (Int_t i=0; i<300000; i++) his550->Fill(RndmdNchdY(5.5));
122 his550->Draw();
123 TF1 f1("f1","[0]*x^(-(0.00001+abs([1])))",1,2000)
124 f1.SetParameters(1,-1)
125 his550->Fit("f1","","",10,3000);
126 TH1F his276("his276","his276",1000,0,3000)
2ff913b6 127 for (Int_t i=0; i<300000; i++) his276->Fill(RndmdNchdY(2.76));
06ab15a8 128 his276->Draw();
129
130 */
131 static TSpline3 * spline276=0;
2ff913b6 132 const Double_t sref=2.76; // reference s
133
06ab15a8 134 if (!spline276){
2ff913b6 135 // Refence multiplicities for 2.76 TeV
06ab15a8 136 // multplicity from archive except of the last point was set to 0
2ff913b6 137 //
138 const Double_t mult[20]={1601, 1294, 966, 649, 426, 261, 149, 76, 35, 0.001};
139 const Double_t cent[20]={2.5, 7.5, 15, 25, 35, 45, 55, 65, 75, 100.};
06ab15a8 140 TGraphErrors * gr = new TGraphErrors(10,cent,mult);
141 spline276 = new TSpline3("spline276",gr);
142 }
2ff913b6 143 Double_t norm = TMath::Power((s*s)/(sref*sref),0.15);
06ab15a8 144 spline276->Eval(gRandom->Rndm()*100.);
145 return spline276->Eval(gRandom->Rndm()*100.)*norm;
146}
147
148
149
150
151
4b955de3 152void pileUpToyMC(Int_t nframes){
153 //
154 //
155 //
156 /*
157 Int)t nframes=1000;
158 */
159 TTreeSRedirector *pcstream = new TTreeSRedirector("pileup.root","recreate");
4b955de3 160 Double_t central = 2350;
161 Double_t pmean=5;
162 TVectorD vectorT(nframes);
163 //
164 for (Int_t irate=1; irate<10; irate++){
165 printf("rate\t%d\n",irate);
166 for (Int_t iframe=0; iframe<nframes; iframe++){
167 if (iframe%100000==0)printf("iframe=%d\n",iframe);
168 Int_t ntracksAll=0;
169 Int_t nevents=gRandom->Poisson(irate);
170 Int_t ntracks=0; // to be taken from the MB primary distribution
171 Bool_t hasCentral=0;
172 for (Int_t ievent=0; ievent<nevents; ievent++){
06ab15a8 173 ntracks=RndmdNchdY(5.5);
4b955de3 174 ntracksAll+=ntracks;
175 if (ntracks>central) hasCentral = kTRUE;
176 }
177 (*pcstream)<<"pileupFrame"<<
178 "rate="<<irate<<
179 "nevents="<<nevents<<
180 "ntracks="<<ntracks<<
181 "ntracksAll="<<ntracksAll<<
182 "hasCentral"<<hasCentral<<
183 "\n";
184 vectorT[iframe]=ntracksAll;
185 }
186 Double_t mean = TMath::Mean(nframes, vectorT.GetMatrixArray());
187 Double_t rms = TMath::RMS(nframes, vectorT.GetMatrixArray());
188 Double_t median = TMath::Median(nframes, vectorT.GetMatrixArray());
189 Double_t ord90 = TMath::KOrdStat(nframes,vectorT.GetMatrixArray() , Int_t(nframes*0.90));
190 Double_t ord95 = TMath::KOrdStat(nframes,vectorT.GetMatrixArray() , Int_t(nframes*0.95));
191 Double_t ord99 = TMath::KOrdStat(nframes,vectorT.GetMatrixArray() , Int_t(nframes*0.99));
192 Double_t ord999 = TMath::KOrdStat(nframes,vectorT.GetMatrixArray() , Int_t(nframes*0.999));
193 Double_t ord9999 = TMath::KOrdStat(nframes,vectorT.GetMatrixArray() , Int_t(nframes*0.9999));
194 (*pcstream)<<"pileup"<<
195 "rate="<<irate<<
196 "mean="<<mean<<
197 "rms="<<rms<<
198 "median="<<median<<
199 "ord90="<<ord90<<
200 "ord95="<<ord95<<
201 "ord99="<<ord99<<
202 "ord999="<<ord999<<
203 "ord9999="<<ord9999<<
204 "\n";
205 }
206 delete pcstream;
207 // Draw
208 pcstream = new TTreeSRedirector("pileup.root","update");
209 TTree * treeStat = (TTree*)(pcstream->GetFile()->Get("pileup"));
210 TTree * treeFrame = (TTree*)(pcstream->GetFile()->Get("pileupFrame"));
211 Int_t mentries = treeStat->Draw("ord999","1","goff");
212 Double_t maximum = TMath::MaxElement(mentries, treeStat->GetV1());
213 const char * names[6]={"mean","median","ord90","ord95","ord99","ord999"};
214 const char * titles[6]={"Mean","Median","90 %","95 %","99 %","99.9 %"};
215 const Int_t mcolors[6]={1,2,3,4,6,7};
216 //
217 //
218 TF1 * f1 = new TF1("f1","[0]*x+[1]*sqrt(x)");
219 Double_t par0=0;
220 //
221 TCanvas * canvasMult = new TCanvas("canvasCumul","canvasCumul");
222 canvasMult->SetLeftMargin(0.13);
223 TLegend * legend= new TLegend(0.14,0.6,0.45,0.89, "Effective dN_{ch}/d#eta");
224 TGraphErrors *graphs[6]={0};
225 for (Int_t igr=0; igr<6; igr++){
226 graphs[igr] = TStatToolkit::MakeGraphErrors(treeStat,Form("%s:rate",names[igr]),"1",21+(igr%5),mcolors[igr],0);
227 graphs[igr]->SetMinimum(0);
228 graphs[igr]->GetYaxis()->SetTitleOffset(1.3);
229 graphs[igr]->SetMaximum(maximum*1.1);
230 graphs[igr]->GetXaxis()->SetTitle("<N_{ev}>");
231 graphs[igr]->GetYaxis()->SetTitle("dN_{ch}/d#eta");
232 TF1 * f2 = new TF1("f2","[0]*x+[1]*sqrt(x)");
233 f2->SetLineColor(mcolors[igr]);
234 f2->SetLineWidth(0.5);
235 if (igr>0) f2->FixParameter(0,par0);
236 graphs[igr]->Fit(f2,"","");
237 if (igr==0) par0=f2->GetParameter(0);
238 if (igr==0) graphs[igr]->Draw("ap");
239 graphs[igr]->Draw("p");
240 legend->AddEntry(graphs[igr], titles[igr],"p");
241 }
242 legend->SetBorderSize(0);
243 legend->Draw();
244
245 canvasMult->SaveAs("effectiveMult.pdf");
246 canvasMult->SaveAs("effectiveMult.png");
247 gStyle->SetOptStat(0);
248 TH2F * hisMult = new TH2F("ntracksNevent","ntracksnevents",9,1,10,100,0,2*maximum);
249 {
250 treeFrame->Draw("ntracksAll:rate>>ntracksNevent","","colz");
251 hisMult->GetXaxis()->SetTitle("<N_{ev}>");
252 hisMult->GetYaxis()->SetTitle("dN_{ch}/d#eta");
253 hisMult->GetYaxis()->SetTitleOffset(1.3);
254 hisMult->Draw("colz");
255 }
256 canvasMult->SaveAs("effectiveMultColz.pdf");
257 canvasMult->SaveAs("effectiveMultColz.png");
258 //
259 //
260 //
261 TH2F * hisMult5 = new TH2F("ntracksNevent5","ntracksnEvents5",9,1,10,100,0,maximum);
262 {
263 treeFrame->Draw("ntracksAll:nevents>>ntracksNevent5","abs(rate-5)<0.5","colz");
264 hisMult5->GetXaxis()->SetTitle("N_{ev}");
265 hisMult5->GetYaxis()->SetTitle("dN_{ch}/d#eta");
266 hisMult5->GetYaxis()->SetTitleOffset(1.3);
267 hisMult5->Draw("colz");
268 }
269 canvasMult->SaveAs("effectiveMultF5.pdf");
270 canvasMult->SaveAs("effectiveMultF5.png");
2ff913b6 271
272 {
273 gStyle->SetOptFit(1);
274 gStyle->SetOptStat(1);
275 gStyle->SetOptTitle(1);
276 TCanvas * canvasMultH = new TCanvas("canvasCumulH","canvasCumulH",700,700);
277 canvasMultH->Divide(1,2);
278 canvasMultH->cd(1);
279 TH1F his550("his550","his550",1000,0,3000);
280 TH1F his276("his276","his276",1000,0,3000);
281 for (Int_t i=0; i<300000; i++) his550.Fill(RndmdNchdY(5.5));
282 for (Int_t i=0; i<300000; i++) his276.Fill(RndmdNchdY(2.76));
283 TF1 f1("f1","[0]*x^(-(0.00001+abs([1])))",1,2000);
284 f1.SetParameters(1,-1);
285 his550.GetXaxis()->SetTitle("dN_{ch}/d#eta");
286 his276.GetXaxis()->SetTitle("dN_{ch}/d#eta");
287 his550.Fit("f1","","",10,3000);
288 his276.Fit("f1","","",10,3000);
289 canvasMultH->cd(1)->SetLogx(1);
290 canvasMultH->cd(1)->SetLogy(1);
291 his550.Draw();
292 canvasMultH->cd(2)->SetLogx(1);
293 canvasMultH->cd(2)->SetLogy(1);
294 his276.Draw("");
295 canvasMultH->SaveAs("dNchdEta.pdf");
296 }
4b955de3 297 delete pcstream;
298}
c0172f82 299
7212f8c3 300void spaceChargeFluctuationToyMC(Int_t nframes, Double_t interactionRate){
c0172f82 301 //
7212f8c3 302 // Toy MC to generate space charge fluctuation, to estimate the fluctuation of the integral space charge in part of the
303 // TPC
304 // Parameters:
305 // nframes - number of frames to simulate
306 // 1. Make a toy simulation part for given setup
307 // 2. Make a summary plots for given setups - see function spaceChargeFluctuationToyMCDraw()
c0172f82 308 //
309 TTreeSRedirector *pcstream = new TTreeSRedirector("spaceChargeFluctuation.root","recreate");
4e715689 310 Double_t driftTime=0.1;
4e715689 311 Double_t eventMean=interactionRate*driftTime;
312 Double_t trackMean=500;
313 Double_t FPOT=1.0, EEND=3000;
314 Double_t EEXPO=0.8567;
315 const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
c0172f82 316
317 for (Int_t iframe=0; iframe<nframes; iframe++){
4e715689 318 printf("iframe=%d\n",iframe);
c0172f82 319 Int_t nevents=gRandom->Poisson(interactionRate*driftTime);
320 Int_t ntracksAll=0;
4e715689 321 TVectorD vecTracksPhi180(180);
322 TVectorD vecTracksPhi36(36);
323 TVectorD vecEPhi180(180);
324 TVectorD vecEPhi36(36);
325 Double_t dESum=0;
c0172f82 326 for (Int_t ievent=0; ievent<nevents; ievent++){
4e715689 327 Int_t ntracks=gRandom->Exp(trackMean); // to be taken from the MB primary distribution
328 Float_t RAN = gRandom->Rndm();
329 ntracks=TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/2.;
c0172f82 330 ntracksAll+=ntracks;
4e715689 331 for (Int_t itrack=0; itrack<ntracks; itrack++){
332 Double_t phi = gRandom->Rndm();
333 vecTracksPhi180(Int_t(phi*180))+=1;
334 vecTracksPhi36(Int_t(phi*36)) +=1;
335 // simplified MC to get track length including loopers
336 Double_t theta= gRandom->Rndm();
337 Double_t pt = gRandom->Exp(0.5)+0.05;
338 Double_t crv = TMath::Abs(5*kB2C/pt); //GetC(b); // bz*kB2C/pt;
339 Double_t deltaPhi=0;
340 if (TMath::Abs(2*crv*(245-85)/2.) <1.) deltaPhi=TMath::ASin(crv*(245-85)/2.);
341 else
342 deltaPhi=TMath::Pi();
343 Double_t dE=deltaPhi/crv;
344 Double_t xloop=1;
345 if (1./crv<250) {
346 xloop = TMath::Min(1./(TMath::Abs(theta)+0.0001),10.);
347 if (xloop<1) xloop=1;
348 }
349 dESum+=xloop*dE;
350 if (itrack==0) (*pcstream)<<"track"<<
351 "pt="<<pt<<
352 "crv="<<crv<<
353 "theta="<<theta<<
354 "dE="<<dE<<
355 "xloop="<<xloop<<
356 "\n";
357
358 vecEPhi180(Int_t(phi*180)) +=dE*xloop;
359 vecEPhi36(Int_t(phi*36)) +=dE*xloop;
c0172f82 360 }
4e715689 361 (*pcstream)<<"event"<<
362 "ntracks="<<ntracks<<
363 "nevents="<<nevents<<
364 "\n";
c0172f82 365 }
366 (*pcstream)<<"ntracks"<<
7212f8c3 367 "rate="<<interactionRate<< // interaction rate
368 "eventMean="<<eventMean<< // mean number of events per frame
369 "trackMean="<<trackMean<< // assumed mean of the tracks per event
370 //
4e715689 371 "nevents="<<nevents<< // number of events withing time frame
372 "ntracksAll="<<ntracksAll<< // number of tracks within time frame
373 "dESum="<<dESum<< // sum of the energy loss
374 "vecTracksPhi36.="<<&vecTracksPhi36<< // number of tracks in phi bin (36 bins) within time frame
375 "vecTracksPhi180.="<<&vecTracksPhi180<< // number of tracks in phi bin (180 bins) within time frame
376 "vecEPhi36.="<<&vecEPhi36<< // number of tracks in phi bin (36 bins) within time frame
377 "vecEPhi180.="<<&vecEPhi180<< // number of tracks in phi bin (180 bins) within time frame
c0172f82 378 "\n";
379 }
380 delete pcstream;
7212f8c3 381 spaceChargeFluctuationToyDraw();
382}
383
384
385void spaceChargeFluctuationToyDraw(){
4e715689 386 //
387 // Toy MC to simulate the space charge integral fluctuation
7212f8c3 388 // Draw function for given setup
389 // for MC generation part see : void spaceChargeFluctuationToyMC
390 TTreeSRedirector *pcstream = new TTreeSRedirector("spaceChargeFluctuation.root","update");
391 TFile * f = pcstream->GetFile();
4e715689 392 TTree * treeStat = (TTree*)f->Get("ntracks");
393 TTree * treedE = (TTree*)f->Get("track");
394 TTree * treeEv = (TTree*)f->Get("event");
7212f8c3 395
4e715689 396 Int_t nentries=treedE->Draw("dE*xloop","1","",1000000);
7212f8c3 397
4e715689 398 Double_t meandE=TMath::Mean(nentries,treedE->GetV1());
399 Double_t rmsdE=TMath::RMS(nentries,treedE->GetV1());
400 treeStat->SetAlias("meandE",Form("(%f+0)",meandE));
401 treeStat->SetAlias("rmsdE",Form("(%f+0)",rmsdE));
402 nentries=treeEv->Draw("ntracks","1","",1000000);
403 Double_t meanT=TMath::Mean(nentries,treeEv->GetV1());
404 Double_t rmsT=TMath::RMS(nentries,treeEv->GetV1());
405 treeStat->SetAlias("tracksMean",Form("(%f+0)",meanT));
406 treeStat->SetAlias("tracksRMS",Form("(%f+0)",rmsT));
7212f8c3 407 nentries = treeStat->Draw("eventMean","","");
408 Double_t meanEvents =TMath::Mean(nentries,treeStat->GetV1());
409 treeStat->SetMarkerStyle(21);
410 treeStat->SetMarkerSize(0.4);
4e715689 411 //
7212f8c3 412 const Int_t kColors[6]={1,2,3,4,6,7};
413 const Int_t kStyle[6]={20,21,24,25,24,25};
414 const char * htitles[6]={"Events","Tracks","Tracks #phi region (1/180)","Q #phi region (1/180)", "Tracks #phi region (1/36)","Q #phi region (1/36)"};
415 const char * hnames[6]={"Events","Tracks","TracksPhi180","QPhi180", "TracksPhi36","QPhi36"};
416
417 TH1* hisFluc[6]={0};
418 TH1* hisPull[6]={0};
419 TVectorD *vecFitFluc[6]={0};
420 TVectorD *vecFitFlucPull[6]={0};
421 //
422 // histograms
423 //
424 treeStat->Draw("nevents/eventMean>>hisEv(100,0.85,1.15)","");
425 hisFluc[0]=(TH1*)treeStat->GetHistogram()->Clone();
426 treeStat->Draw("ntracksAll/(eventMean*tracksMean)>>hisTrackAll(100,0.85,1.1)","","same");
427 hisFluc[1]=(TH1*)treeStat->GetHistogram()->Clone();
428 treeStat->Draw("vecTracksPhi180.fElements/(eventMean*tracksMean/180)>>hisTrackSector(100,0.85,1.1)","1/180","same");
429 hisFluc[2]=(TH1*)treeStat->GetHistogram()->Clone();
430 treeStat->Draw("vecEPhi180.fElements/(eventMean*tracksMean*meandE/180)>>hisdESector(100,0.85,1.1)","1/180","same");
431 hisFluc[3]=(TH1*)treeStat->GetHistogram()->Clone();
432 treeStat->Draw("vecTracksPhi36.fElements/(eventMean*tracksMean/36)>>hisTrackSector36(100,0.85,1.1)","1/36","same");
433 hisFluc[4]=(TH1*)treeStat->GetHistogram()->Clone();
434 treeStat->Draw("vecEPhi36.fElements/(eventMean*tracksMean*meandE/36)>>hisdESector36(100,0.85,1.1)","1/36","same");
435 hisFluc[5]=(TH1*)treeStat->GetHistogram()->Clone();
436 //
437 // pulls
4e715689 438 //
7212f8c3 439 treeStat->Draw("((nevents/eventMean)-1)/sqrt(1/eventMean)>>pullEvent(100,-6,6)","","err"); //tracks All pull
440 hisPull[0]=(TH1*)treeStat->GetHistogram()->Clone();
441 treeStat->Draw("(ntracksAll/(eventMean*tracksMean)-1)/sqrt(1/eventMean+(tracksRMS/tracksMean)**2/eventMean)>>pullTrackAll(100,-6,6)","","err"); //tracks All pull
442 hisPull[1]=(TH1*)treeStat->GetHistogram()->Clone();
443 treeStat->Draw("(vecTracksPhi180.fElements/(eventMean*tracksMean/180.)-1)/sqrt(1/eventMean+(tracksRMS/tracksMean)**2/eventMean+180/(tracksMean*eventMean))>>pullTrack180(100,-6,6)","1/180","errsame"); //tracks spread
444 hisPull[2]=(TH1*)treeStat->GetHistogram()->Clone();
445 treeStat->Draw("(vecEPhi180.fElements/(eventMean*tracksMean*meandE/180)-1)/sqrt(1/eventMean+(tracksRMS/tracksMean)**2/eventMean+180/(tracksMean*eventMean)+180*(rmsdE/meandE)**2/(eventMean*tracksMean))>>hisPulldE180(100,-6,6)","1/180","errsame"); //dE spread
446 hisPull[3]=(TH1*)treeStat->GetHistogram()->Clone();
447 treeStat->Draw("(vecTracksPhi36.fElements/(eventMean*tracksMean/36.)-1)/sqrt(1/eventMean+(tracksRMS/tracksMean)**2/eventMean+36/(tracksMean*eventMean))>>pullTrack36(100,-6,6)","1/36","errsame"); //tracks spread
448 hisPull[4]=(TH1*)treeStat->GetHistogram()->Clone();
449 treeStat->Draw("(vecEPhi36.fElements/(eventMean*tracksMean*meandE/36)-1)/sqrt(1/eventMean+(tracksRMS/tracksMean)**2/eventMean+36/(tracksMean*eventMean)+36*(rmsdE/meandE)**2/(eventMean*tracksMean))>>hisPulldE36(100,-6,6)","1/36","errsame"); //dE spread
450 hisPull[5]=(TH1*)treeStat->GetHistogram()->Clone();
4e715689 451 //
7212f8c3 452
453 for (Int_t ihis=0; ihis<6; ihis++) {
454 vecFitFluc[ihis] = new TVectorD(3);
455 vecFitFlucPull[ihis] = new TVectorD(3);
456 TF1 * fg = new TF1(Form("fg%d",ihis),"gaus");
457 fg->SetLineWidth(0.5);
458 fg->SetLineColor(kColors[ihis]);
459 hisFluc[ihis]->Fit(fg,"","");
460 fg->GetParameters( vecFitFluc[ihis]->GetMatrixArray());
461 hisPull[ihis]->Fit(fg,"","");
462 fg->GetParameters( vecFitFlucPull[ihis]->GetMatrixArray());
463 hisFluc[ihis]->SetName(Form("Fluctuation%s",hnames[ihis]));
464 hisFluc[ihis]->SetTitle(Form("Fluctuation%s",htitles[ihis]));
465 hisPull[ihis]->SetName(Form("Pull%s",hnames[ihis]));
466 hisPull[ihis]->SetTitle(Form("Pull%s",htitles[ihis]));
467 }
468
469 gStyle->SetOptStat(0);
470 TCanvas * canvasQFluc= new TCanvas("SpaceChargeFluc","SpaceChargeFluc",600,700);
471 canvasQFluc->Divide(1,2);
472 canvasQFluc->cd(1);
473 TLegend * legendFluc = new TLegend(0.11,0.55,0.45,0.89,"Relative fluctuation");
474 TLegend * legendPull = new TLegend(0.11,0.55,0.45,0.89,"Fluctuation pulls");
475 for (Int_t ihis=0; ihis<6; ihis++){
476 hisFluc[ihis]->SetMarkerStyle(kStyle[ihis]);
477 hisFluc[ihis]->SetMarkerColor(kColors[ihis]);
478 hisFluc[ihis]->SetMarkerSize(0.8);
479 if (ihis==0) hisFluc[ihis]->Draw("err");
480 hisFluc[ihis]->Draw("errsame");
481 legendFluc->AddEntry(hisFluc[ihis],htitles[ihis]);
482 }
483 legendFluc->Draw();
484
485 canvasQFluc->cd(2);
486 for (Int_t ihis=0; ihis<6; ihis++){
487 hisPull[ihis]->SetMarkerStyle(kStyle[ihis]);
488 hisPull[ihis]->SetMarkerColor(kColors[ihis]);
489 hisPull[ihis]->SetMarkerSize(0.8);
490 if (ihis==0) hisPull[ihis]->Draw("err");
491 hisPull[ihis]->Draw("errsame");
492 legendPull->AddEntry(hisPull[ihis],htitles[ihis]);
493 }
494 legendPull->Draw();
4e715689 495 //
7212f8c3 496 for (Int_t ihis=0; ihis<6; ihis++){
497 hisFluc[ihis]->Write();
498 hisPull[ihis]->Write();
499 }
500 (*pcstream)<<"summary"<< // summary information for given setup
501 "meanEvents="<<meanEvents<< // mean number of events in the frame
502 "meandE="<<meandE<< // mean "energy loss" of track
503 "rmsdE="<<rmsdE<< // rms
504 "meanT="<<meanT<< // mean number of tracks per MB event
505 "rmsT="<<rmsT<< // rms of onumber of tracks
506 // // fit of the relative fluctuation
507 "vflucE.="<<vecFitFluc[0]<< // in events
508 "vflucEP.="<<vecFitFlucPull[0]<< // in events pull
509 "vflucTr.="<<vecFitFluc[1]<< // in tracks
510 "vflucTrP.="<<vecFitFlucPull[1]<<
511 //
512 "vflucTr180.="<<vecFitFluc[2]<<
513 "vflucTr180P.="<<vecFitFlucPull[2]<<
514 "vflucE180.="<<vecFitFluc[3]<<
515 "vflucE180P.="<<vecFitFlucPull[3]<<
516 //
517 "vflucTr36.="<<vecFitFluc[4]<<
518 "vflucTr36P.="<<vecFitFlucPull[4]<<
519 "vflucE36.="<<vecFitFluc[5]<<
520 "vflucE36P.="<<vecFitFlucPull[5]<<
521 "\n";
522 canvasQFluc->SaveAs("CanvasFluctuation.pdf");
523 canvasQFluc->SaveAs("CanvasFluctuation.png");
524 delete pcstream;
c0172f82 525
4e715689 526}
7212f8c3 527
528void spaceChargeFluctuationToyDrawSummary(){
4e715689 529 //
7212f8c3 530 // make a summary information plots using several runs with differnt mean IR setting
531 // Input:
532 // space.list - list of root files produced by spaceChargeFluctuationToyDraw
533 // Output:
534 // canvas saved in current directory
4e715689 535 //
7212f8c3 536 TChain * chain = AliXRDPROOFtoolkit::MakeChain("space.list","summary",0,100);
537 chain->SetMarkerStyle(21);
538 const Int_t kColors[6]={1,2,3,4,6,7};
539 const Int_t kStyle[6]={20,21,24,25,24,25};
540 const char * htitles[6]={"Events","Tracks","Tracks #phi region (1/180)","Q #phi region (1/180)", "Tracks #phi region (1/36)","Q #phi region (1/36)"};
541 // const char * hnames[6]={"Events","Tracks","TracksPhi180","QPhi180", "TracksPhi36","QPhi36"};
542 //
543 Double_t meanT,rmsT=0;
544 Double_t meandE,rmsdE=0;
545 Int_t entries = chain->Draw("meanT:rmsT:meandE:rmsdE","1","goff");
546 meanT =TMath::Mean(entries, chain->GetV1());
547 rmsT =TMath::Mean(entries, chain->GetV2());
548 meandE =TMath::Mean(entries, chain->GetV3());
549 rmsdE =TMath::Mean(entries, chain->GetV4());
550 //
551 //
552 //
553 TGraphErrors * graphs[6]={0};
554 TF1 * functions[6]={0};
555
556 graphs[5]=TStatToolkit::MakeGraphErrors(chain,"vflucE36.fElements[2]:meanEvents:0.025*vflucE36.fElements[2]","1",kStyle[5],kColors[5],1);
557 graphs[4]=TStatToolkit::MakeGraphErrors(chain,"vflucTr36.fElements[2]:meanEvents:0.025*vflucTr36.fElements[2]","1",kStyle[4],kColors[4],1);
558 graphs[3]=TStatToolkit::MakeGraphErrors(chain,"vflucE180.fElements[2]:meanEvents:0.025*vflucE180.fElements[2]","1",kStyle[3],kColors[3],1);
559 graphs[2]=TStatToolkit::MakeGraphErrors(chain,"vflucTr180.fElements[2]:meanEvents:0.025*vflucTr180.fElements[2]","1",kStyle[2],kColors[2],1);
560 graphs[1]=TStatToolkit::MakeGraphErrors(chain,"vflucTr.fElements[2]:meanEvents:0.025*vflucTr.fElements[2]","1",kStyle[1],kColors[1],1);
561 graphs[0]=TStatToolkit::MakeGraphErrors(chain,"vflucE.fElements[2]:meanEvents:0.025*vflucE.fElements[2]","1",kStyle[0],kColors[0],1);
562
563 functions[5]=new TF1("fe","sqrt(1+[0]**2+[1]/[2]+[1]*[3]**2/[2])/sqrt(x)",2000,200000);
564 functions[5]->SetParameters(rmsT/meanT,36.,meanT,rmsdE/meandE);
565 functions[4]=new TF1("fe","sqrt(1+[0]**2+[1]/[2])/sqrt(x)",2000,200000);
566 functions[4]->SetParameters(rmsT/meanT,36.,meanT,0);
567 functions[3]=new TF1("fe","sqrt(1+[0]**2+[1]/[2]+[1]*[3]**2/[2])/sqrt(x)",2000,200000);
568 functions[3]->SetParameters(rmsT/meanT,180.,meanT,rmsdE/meandE);
569 functions[2]=new TF1("fe","sqrt(1+[0]**2+[1]/[2])/sqrt(x)",2000,200000);
570 functions[2]->SetParameters(rmsT/meanT,180.,meanT,0);
571 functions[1]=new TF1("fe","sqrt(1+[0]**2)/sqrt(x)",2000,200000);
572 functions[1]->SetParameters(rmsT/meanT,0);
573 functions[0]=new TF1("fe","sqrt(1)/sqrt(x)",2000,200000);
574
575
576 TCanvas *canvasF= new TCanvas("fluc","fluc",600,500);
577 // TLegend *legend = new TLegend(0.5,0.65,0.89,0.89,"Relative fluctuation #sigma=#sqrt{1+#frac{#sigma_{T}^{2}}{#mu_{T}^{2}}}");
578 TLegend *legendF = new TLegend(0.45,0.5,0.89,0.89,"Relative fluctuation of charge");
579 for (Int_t ihis=0; ihis<4; ihis++){
580 graphs[ihis]->SetMinimum(0.00);
581 graphs[ihis]->SetMaximum(0.05);
582 if (ihis==0) graphs[ihis]->Draw("ap");
583 graphs[ihis]->GetXaxis()->SetTitle("events");
584 graphs[ihis]->GetXaxis()->SetNdivisions(507);
585 graphs[ihis]->GetYaxis()->SetTitle("#frac{#sigma}{#mu}");
586 graphs[ihis]->Draw("p");
587 legendF->AddEntry(graphs[ihis],htitles[ihis],"p");
588 if (functions[ihis]){
589 functions[ihis]->SetLineColor(kColors[ihis]);
590 functions[ihis]->SetLineWidth(0.5);
591 functions[ihis]->Draw("same");
592 }
593 }
594 legendF->Draw();
595 canvasF->SaveAs("spaceChargeFlucScan.pdf");
596 canvasF->SaveAs("spaceChargeFlucScan.png");
597
598 TCanvas *canvasF36= new TCanvas("fluc36","fluc36",600,500);
599 // TLegend *legend = new TLegend(0.5,0.65,0.89,0.89,"Relative fluctuation #sigma=#sqrt{1+#frac{#sigma_{T}^{2}}{#mu_{T}^{2}}}");
600 TLegend *legendF36 = new TLegend(0.45,0.5,0.89,0.89,"Relative fluctuation of charge");
601 for (Int_t ihis=0; ihis<6; ihis++){
602 if (ihis==2 || ihis==3) continue;
603 graphs[ihis]->SetMinimum(0.00);
604 graphs[ihis]->SetMaximum(0.05);
605 if (ihis==0) graphs[ihis]->Draw("ap");
606 graphs[ihis]->GetXaxis()->SetTitle("events");
607 graphs[ihis]->GetXaxis()->SetNdivisions(507);
608 graphs[ihis]->GetYaxis()->SetTitle("#frac{#sigma}{#mu}");
609 graphs[ihis]->Draw("p");
610 legendF36->AddEntry(graphs[ihis],htitles[ihis],"p");
611 if (functions[ihis]){
612 functions[ihis]->SetLineColor(kColors[ihis]);
613 functions[ihis]->SetLineWidth(0.5);
614 functions[ihis]->Draw("same");
615 }
616 }
617 legendF36->Draw();
618 canvasF36->SaveAs("spaceChargeFlucScan36.pdf");
619 canvasF36->SaveAs("spaceChargeFlucScan36.png");
620
621
622}
623
624
625
626void FitMultiplicity(const char * fname="mult_dist_pbpb.root"){
627 //
628 // Fit multiplicity distribution using as a power law in limited range
629 // const char * fname="mult_dist_pbpb.root"
630 TFile *fmult=TFile::Open(fname);
4e715689 631 TF1 f1("f1","[0]*(x+abs([2]))**(-abs([1]))",1,3000);
632 TH1* his = (TH1*) fmult->Get("mult_dist_PbPb_normalizedbywidth");
633 f1.SetParameters(his->GetEntries(),1,1);
634 his->Fit(&f1,"","",2,3000);
635
636 Double_t FPOT=1.0, EEND=3000, EEXPO= TMath::Abs(f1.GetParameter(1));
637 EEXPO=0.8567;
638 const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
639 TH1F *hisr= new TH1F("aaa","aaa",4000,0,4000);
640 for (Int_t i=0; i<400000; i++){
641 Float_t RAN = gRandom->Rndm();
642 hisr->Fill(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO));
643 }
4e715689 644}
c0172f82 645
c0172f82 646
647
648
76dbc35e 649TH1 * GenerateMapRawIons(Int_t useGainMap, const char *fileName, const char *outputName, Int_t maxEvents){
c0172f82 650 //
651 // Generate 3D maps of the space charge for the rad data maps
652 // different threshold considered
653 // Paramaters:
0eff8c53 654 // useGainMap - switch usage of the gain map
c0172f82 655 // fileName - name of input raw file
656 // outputName - name of output file with the space charge histograms
657 // maxEvents - grouping of the events
658 //
659 //
7212f8c3 660 gRandom->SetSeed(0); //set initial seed to be independent for different jobs
661
c0172f82 662 TTreeSRedirector * pcstream = new TTreeSRedirector(outputName, "recreate");
663 const char * ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";
664 AliCDBManager * man = AliCDBManager::Instance();
665 man->SetDefaultStorage(ocdbpath);
666 man->SetRun(0);
667 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 2.76/2.));
668 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
669 AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
76dbc35e 670 AliTPCCalPad * gain = AliTPCcalibDB::Instance()->GetDedxGainFactor();
671 AliTPCCalPad * noise = AliTPCcalibDB::Instance()->GetPadNoise();
672
c0172f82 673 TStopwatch timer;
674 timer.Start();
675 // arrays of space charges - different elements corresponds to different threshold to accumulate charge
76dbc35e 676 TH1D * hisQ1D[3]={0};
76dbc35e 677 TH1D * hisQ1DROC[3]={0};
7212f8c3 678 TH2D * hisQ2DRPhi[3]={0};
679 TH2D * hisQ2DRZ[3]={0};
680 TH2D * hisQ2DRPhiROC[3]={0};
681 TH2D * hisQ2DRZROC[3]={0};
682 TH3D * hisQ3D[3]={0}; // 3D maps space charge from drift volume
683 TH3D * hisQ3DROC[3]={0}; // 3D maps space charge from ROC
684
5c70ec78 685 Int_t nbinsRow=param->GetNRowLow()+param->GetNRowUp();
686 Double_t *xbins = new Double_t[nbinsRow+1];
687 xbins[0]=param->GetPadRowRadiiLow(0)-1; //underflow bin
688 for (Int_t ibin=0; ibin<param->GetNRowLow();ibin++) xbins[1+ibin]=param->GetPadRowRadiiLow(ibin);
689 for (Int_t ibin=0; ibin<param->GetNRowUp();ibin++) xbins[1+ibin+param->GetNRowLow()]=param->GetPadRowRadiiUp(ibin);
c0172f82 690 //
691 for (Int_t ith=0; ith<3; ith++){
692 char chname[100];
7212f8c3 693 // 1D
694 snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
695 hisQ1D[ith] = new TH1D(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
696 snprintf(chname,100,"hisQ1DROC_Th%d",2*ith+2);
697 hisQ1DROC[ith] = new TH1D(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
698 // 3D
c0172f82 699 snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
7212f8c3 700 hisQ3D[ith] = new TH3D(chname, chname,360, 0,TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
701 snprintf(chname,100,"hisQ3DROC_Th%d",2*ith+2);
702 hisQ3DROC[ith] = new TH3D(chname, chname,360, 0,TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
703 // 2D
704 snprintf(chname,100,"hisQ2DRPhi_Th%d",2*ith+2);
705 hisQ2DRPhi[ith] = new TH2D(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
76dbc35e 706 snprintf(chname,100,"hisQ2DRZ_Th%d",2*ith+2);
76dbc35e 707 hisQ2DRZ[ith] = new TH2D(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36), 125,-250,250);
c0172f82 708 //
7212f8c3 709 snprintf(chname,100,"hisQ2DRPhiROC_Th%d",2*ith+2);
710 hisQ2DRPhiROC[ith] = new TH2D(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
76dbc35e 711 snprintf(chname,100,"hisQ2DRZROC_Th%d",2*ith+2);
712 hisQ2DRZROC[ith] = new TH2D(chname,chname,param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36), 125,-250,250);
7212f8c3 713 //
5c70ec78 714 hisQ1D[ith]->GetXaxis()->Set(nbinsRow,xbins);
5c70ec78 715 hisQ1DROC[ith]->GetXaxis()->Set(nbinsRow,xbins);
7212f8c3 716 hisQ3D[ith]->GetYaxis()->Set(nbinsRow,xbins);
5c70ec78 717 hisQ3DROC[ith]->GetYaxis()->Set(nbinsRow,xbins);
c0172f82 718 //
7212f8c3 719 hisQ2DRPhi[ith]->GetYaxis()->Set(nbinsRow,xbins);
720 hisQ2DRZ[ith]->GetXaxis()->Set(nbinsRow,xbins);
721 hisQ2DRPhiROC[ith]->GetYaxis()->Set(nbinsRow,xbins);
722 hisQ2DRZROC[ith]->GetXaxis()->Set(nbinsRow,xbins);
723 //
c0172f82 724 hisQ1D[ith]->SetDirectory(0);
7212f8c3 725 hisQ1DROC[ith]->SetDirectory(0);
726 hisQ3D[ith]->SetDirectory(0);
c0172f82 727 hisQ3DROC[ith]->SetDirectory(0);
7212f8c3 728 //
729 hisQ2DRPhi[ith]->SetDirectory(0);
730 hisQ2DRZ[ith]->SetDirectory(0);
76dbc35e 731 hisQ2DRZROC[ith]->SetDirectory(0);
7212f8c3 732 hisQ2DRPhiROC[ith]->SetDirectory(0);
c0172f82 733 }
734 //
735 //
736 AliRawReader *reader = new AliRawReaderRoot(fileName);
737 reader->Reset();
738 AliAltroRawStream* stream = new AliAltroRawStream(reader);
739 stream->SelectRawData("TPC");
740 Int_t evtnr=0;
741 Int_t chunkNr=0;
742 //
743
744 while (reader->NextEvent()) {
0eff8c53 745 Double_t shiftZ= gRandom->Rndm()*250.;
c0172f82 746 //
747 if(evtnr>=maxEvents) {
748 chunkNr++;
c0172f82 749 pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
750 pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
751 for (Int_t ith=0; ith<3; ith++){
752 hisQ1D[ith]->Write(Form("His1DDrift_%d",ith));
7212f8c3 753 hisQ2DRPhi[ith]->Write(Form("His2DRPhiDrift_%d",ith));
754 hisQ2DRZ[ith]->Write(Form("His2DRZDrift_%d",ith));
c0172f82 755 hisQ3D[ith]->Write(Form("His3DDrift_%d",ith));
756 hisQ1DROC[ith]->Write(Form("His1DROC_%d",ith));
7212f8c3 757 hisQ2DRPhiROC[ith]->Write(Form("His2DRPhiROC_%d",ith));
76dbc35e 758 hisQ2DRZROC[ith]->Write(Form("His2DRZROC_%d",ith));
c0172f82 759 hisQ3DROC[ith]->Write(Form("His3DROC_%d",ith));
5c70ec78 760 (*pcstream)<<"histo"<<
76dbc35e 761 "events="<<evtnr<<
762 "useGain="<<useGainMap<<
5c70ec78 763 Form("hist1D_%d.=",ith*2+2)<<hisQ1D[ith]<<
7212f8c3 764 Form("hist2DRPhi_%d.=",ith*2+2)<<hisQ2DRPhi[ith]<<
76dbc35e 765 Form("hist2DRZ_%d.=",ith*2+2)<<hisQ2DRZ[ith]<<
766 Form("hist3D_%d.=",ith*2+2)<<hisQ3D[ith]<<
5c70ec78 767 Form("hist1DROC_%d.=",ith*2+2)<<hisQ1DROC[ith]<<
7212f8c3 768 Form("hist2DRPhiROC_%d.=",ith*2+2)<<hisQ2DRPhiROC[ith]<<
76dbc35e 769 Form("hist2DRZROC_%d.=",ith*2+2)<<hisQ2DRZROC[ith]<<
770 Form("hist3DROC_%d.=",ith*2+2)<<hisQ3DROC[ith];
5c70ec78 771 }
772 (*pcstream)<<"histo"<<"\n";
773 for (Int_t ith=0; ith<3; ith++){
c0172f82 774 hisQ1D[ith]->Reset();
7212f8c3 775 hisQ2DRPhi[ith]->Reset();
76dbc35e 776 hisQ2DRZ[ith]->Reset();
c0172f82 777 hisQ3D[ith]->Reset();
778 hisQ1DROC[ith]->Reset();
7212f8c3 779 hisQ2DRPhiROC[ith]->Reset();
76dbc35e 780 hisQ2DRZROC[ith]->Reset();
c0172f82 781 hisQ3DROC[ith]->Reset();
782 }
76dbc35e 783 evtnr=0;
c0172f82 784 }
785 cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
786 evtnr++;
787 AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
788 AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
789 //
790 while (input.NextDDL()){
76dbc35e 791 Int_t sector = input.GetSector();
792 AliTPCCalROC * gainROC =gain->GetCalROC(sector);
793 AliTPCCalROC * noiseROC =noise->GetCalROC(sector);
c0172f82 794 while ( input.NextChannel() ) {
795 Int_t row = input.GetRow();
c0172f82 796 Int_t pad = input.GetPad();
797 Int_t nPads = param->GetNPads(sector,row);
798 Double_t localX = param->GetPadRowRadii(sector,row);
7212f8c3 799 Double_t localY = (pad-nPads/2)*param->GetPadPitchWidth(sector);
c0172f82 800 Double_t localPhi= TMath::ATan2(localY,localX);
801 Double_t phi = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
c0172f82 802 Double_t padLength=param->GetPadPitchLength(sector,row);
76dbc35e 803 Double_t gainPad = gainROC->GetValue(row,pad);
804 Double_t noisePad = noiseROC->GetValue(row,pad);
c0172f82 805 //
806 while ( input.NextBunch() ){
807 Int_t startTbin = (Int_t)input.GetStartTimeBin();
808 Int_t bunchlength = (Int_t)input.GetBunchLength();
809 const UShort_t *sig = input.GetSignals();
810 Int_t aboveTh[3]={0};
811 for (Int_t i=0; i<bunchlength; i++){
76dbc35e 812 if (sig[i]<4*noisePad) continue;
c0172f82 813 for (Int_t ith=0; ith<3; ith++){
814 if (sig[i]>(ith*2)+2) aboveTh[ith]++;
815 }
816 }
817 for (Int_t ith=0; ith<3; ith++){
818 if (aboveTh[ith%3]>1){
819 for (Int_t i=0; i<bunchlength; i++){
820 //
821 // normalization
822 //
76dbc35e 823 Double_t zIonDrift =(param->GetZLength()-startTbin*param->GetZWidth());
c0172f82 824 zIonDrift+=shiftZ;
76dbc35e 825 Double_t signal=sig[i];
826 if (useGainMap) signal/=gainPad;
7212f8c3 827 Double_t shiftPhi = ((sector%36)<18) ? 0: TMath::TwoPi();
c0172f82 828 if (TMath::Abs(zIonDrift)<param->GetZLength()){
829 if ((sector%36)>=18) zIonDrift*=-1; // c side has opposite sign
76dbc35e 830 if (sector%36<18) hisQ1D[ith]->Fill(localX, signal/padLength);
7212f8c3 831 hisQ2DRPhi[ith]->Fill(phi+shiftPhi,localX, signal/padLength);
76dbc35e 832 hisQ2DRZ[ith]->Fill(localX, zIonDrift, signal/padLength);
833 hisQ3D[ith]->Fill(phi,localX,zIonDrift,signal/padLength);
c0172f82 834 }
835 //
836 Double_t zIonROC = ((sector%36)<18)? shiftZ: -shiftZ; // z position of the "ion disc" - A side C side opposite sign
76dbc35e 837 if (sector%36<18) hisQ1DROC[ith]->Fill(localX, signal/padLength);
7212f8c3 838 hisQ2DRPhiROC[ith]->Fill(phi+shiftPhi,localX, signal/padLength);
76dbc35e 839 hisQ2DRZROC[ith]->Fill(localX, zIonROC, signal/padLength);
840 hisQ3DROC[ith]->Fill(phi,localX,zIonROC,signal/padLength);
c0172f82 841 }
842 }
843 }
844 }
845 }
846 }
847 }
848 timer.Print();
849 delete pcstream;
850 return 0;
851}
76dbc35e 852
853
7212f8c3 854void DoMerge(){
76dbc35e 855 //
7212f8c3 856 // Merge results to the tree
76dbc35e 857 //
7212f8c3 858 TFile * fhisto = new TFile("histo.root","recreate");
859 TTree * tree = 0;
860 TChain *chain = AliXRDPROOFtoolkit::MakeChainRandom("histo.list","histo",0,100,1);
861 chain->SetBranchStatus("hist3DROC_6*",kFALSE);
862 chain->SetBranchStatus("hist3DROC_4*",kFALSE);
863 tree = chain->CopyTree("1");
864 tree->Write("histo");
865 delete fhisto;
866}
867
868
869
870
871void AnalyzeMaps1D(){
76dbc35e 872 //
7212f8c3 873 // Analyze space charge maps stored as s hitograms in trees
76dbc35e 874 //
7212f8c3 875 TFile * fhisto = new TFile("histo.root");
876 TTree * tree = (TTree*)fhisto->Get("histo");
76dbc35e 877 //
878 TH1 *his1Th[3]={0,0,0};
879 TF1 *fq1DStep= new TF1("fq1DStep","([0]+[1]*(x>134))/x**min(abs([2]),3)",85,245);
880 fq1DStep->SetParameters(1,-0.5,1);
881 tree->Draw("hist1DROC_2.fArray:hist1D_2.fXaxis.fXbins.fArray>>his(40,85,245)","","prof");
882 tree->GetHistogram()->Fit(fq1DStep);
883 // normalize step between the IROC-OROC
884 tree->SetAlias("normQ",Form("(1+%f*(hist1D_2.fXaxis.fXbins.fArray>136))",fq1DStep->GetParameter(1)/fq1DStep->GetParameter(0)));
885 //
886 {
887 Int_t entries= tree->Draw("hist1DROC_2.fArray/(events*normQ)","1","goff");
888 Double_t median=TMath::Median(entries,tree->GetV1());
889 TCut cut10Median = Form("hist1DROC_2.fArray/(events*normQ)<%f",10*median);
890 //
891 tree->Draw("hist1DROC_2.fArray/(events*normQ):hist1D_2.fXaxis.fXbins.fArray>>his1Th0(40,86,245)",cut10Median+"","prof");
892 his1Th[0] = tree->GetHistogram();
893 tree->Draw("hist1DROC_4.fArray/(events*normQ):hist1D_2.fXaxis.fXbins.fArray>>his1Th1(40,86,245)",cut10Median+"","prof");
894 his1Th[1] = tree->GetHistogram();
895 tree->Draw("hist1DROC_6.fArray/(events*normQ):hist1D_2.fXaxis.fXbins.fArray>>his1Th2(40,86,245)",cut10Median+"","prof");
896 his1Th[2]=tree->GetHistogram();
897 }
898 //
899 TCanvas *canvasR = new TCanvas("canvasR","canvasR",600,500);
900 canvasR->cd();
901 for (Int_t i=0; i<3; i++){
902 his1Th[i]->SetMarkerStyle(21);
903 his1Th[i]->SetMarkerColor(i+2);
904 fq1DStep->SetLineColor(i+2);
905 his1Th[i]->Fit(fq1DStep,"","");
906 his1Th[i]->GetXaxis()->SetTitle("r (cm)");
907 his1Th[i]->GetYaxis()->SetTitle("#frac{N_{el}}{N_{ev}}(ADC/cm)");
908 }
909 TLegend * legend = new TLegend(0.11,0.11,0.7,0.39,"1D space Charge map (ROC part) (z,phi integrated)");
910 for (Int_t i=0; i<3; i++){
911 his1Th[i]->SetMinimum(0);fq1DStep->SetLineColor(i+2);
912 his1Th[i]->Fit(fq1DStep,"qnr","qnr");
913 if (i==0) his1Th[i]->Draw("");
914 his1Th[i]->Draw("same");
915 legend->AddEntry(his1Th[i],Form("Thr=%d Slope=%2.2f",2*i+2,fq1DStep->GetParameter(2)));
916 }
917 legend->Draw();
7212f8c3 918 canvasR->SaveAs("spaceCharge1d.png");
919 canvasR->SaveAs("spaceCharge1d.eps");
76dbc35e 920 //
921 //
922 //
923}
0eff8c53 924void MakeFluctuationStudy3D(Int_t nhistos, Int_t nevents, Int_t niter){
925 //
926 //
927 //
0eff8c53 928 //
7212f8c3 929 // Input:
930 // nhistos - maximal number of histograms to be used for sum
931 // nevents - number of events to make a fluctuation studies
932 // niter - number of itterations
933 // Algortihm:
0eff8c53 934 // 1. Make a summary integral 3D/2D/1D maps
935 // 2. Create several maps with niter events - Poisson flucturation in n
936 // 3. Store results 3D maps in the tree (and also as histogram) current and mean
7212f8c3 937 //
938
939 TFile * fhisto = TFile::Open("histo.root");
76dbc35e 940 TTree * tree = (TTree*)fhisto->Get("histo");
0eff8c53 941 tree->SetCacheSize(10000000000);
76dbc35e 942
7212f8c3 943 TTreeSRedirector * pcstream = new TTreeSRedirector("fluctuation.root", "update");
944
945
76dbc35e 946 TH1D * his1DROC=0, * his1DROCSum=0, * his1DROCN=0;
947 TH1D * his1DDrift=0, * his1DDriftSum=0, * his1DDriftN=0 ;
7212f8c3 948 TH2D * his2DRPhiROC=0, * his2DRPhiROCSum=0, * his2DRPhiROCN=0;
0eff8c53 949 TH2D * his2DRZROC=0, * his2DRZROCSum=0, * his2DRZROCN=0;
7212f8c3 950 TH2D * his2DRPhiDrift=0, * his2DRPhiDriftSum=0, * his2DRPhiDriftN=0;
0eff8c53 951 TH2D * his2DRZDrift=0, * his2DRZDriftSum=0, * his2DRZDriftN=0;
76dbc35e 952 TH3D * his3DROC=0, * his3DROCSum=0, * his3DROCN=0;
953 TH3D * his3DDrift=0, * his3DDriftSum=0, * his3DDriftN=0;
954 //
0eff8c53 955 if (nhistos<0 || nhistos> tree->GetEntries()) nhistos = tree->GetEntries();
956 Int_t eventsPerChunk=0;
76dbc35e 957 tree->SetBranchAddress("hist1D_2.",&his1DDrift);
958 tree->SetBranchAddress("hist1DROC_2.",&his1DROC);
7212f8c3 959 tree->SetBranchAddress("hist2DRPhi_2.",&his2DRPhiDrift);
0eff8c53 960 tree->SetBranchAddress("hist2DRZ_2.",&his2DRZDrift);
7212f8c3 961 tree->SetBranchAddress("hist2DRPhiROC_2.",&his2DRPhiROC);
76dbc35e 962 tree->SetBranchAddress("hist3D_2.",&his3DDrift);
963 tree->SetBranchAddress("hist3DROC_2.",&his3DROC);
0eff8c53 964 tree->SetBranchAddress("hist2DRZROC_2.",&his2DRZROC);
965 tree->SetBranchAddress("events",&eventsPerChunk);
966 //
967 // 1. Make a summary integral 3D/2D/1D maps
76dbc35e 968 //
0eff8c53 969 Int_t neventsAll=0;
76dbc35e 970 for (Int_t i=0; i<nhistos; i++){
971 tree->GetEntry(i);
972 if (i%25==0) printf("%d\n",i);
973 if (his1DROCSum==0) his1DROCSum=new TH1D(*his1DROC);
974 if (his1DDriftSum==0) his1DDriftSum=new TH1D(*his1DDrift);
7212f8c3 975 if (his2DRPhiROCSum==0) his2DRPhiROCSum=new TH2D(*his2DRPhiROC);
0eff8c53 976 if (his2DRZROCSum==0) his2DRZROCSum=new TH2D(*his2DRZROC);
7212f8c3 977 if (his2DRPhiDriftSum==0) his2DRPhiDriftSum=new TH2D(*his2DRPhiDrift);
0eff8c53 978 if (his2DRZDriftSum==0) his2DRZDriftSum=new TH2D(*his2DRZDrift);
76dbc35e 979 if (his3DROCSum==0) his3DROCSum=new TH3D(*his3DROC);
980 if (his3DDriftSum==0) his3DDriftSum=new TH3D(*his3DDrift);
981 his1DROCSum->Add(his1DROC);
982 his1DDriftSum->Add(his1DDrift);
7212f8c3 983 his2DRPhiROCSum->Add(his2DRPhiROC);
0eff8c53 984 his2DRZROCSum->Add(his2DRZROC);
7212f8c3 985 his2DRPhiDriftSum->Add(his2DRPhiDrift);
0eff8c53 986 his2DRZDriftSum->Add(his2DRZDrift);
76dbc35e 987 his3DROCSum->Add(his3DROC);
988 his3DDriftSum->Add(his3DDrift);
0eff8c53 989 neventsAll+=eventsPerChunk;
76dbc35e 990 }
991 //
0eff8c53 992 // 2. Create several maps with niter events - Poisson flucturation in n
993 //
994 for (Int_t iter=0; iter<niter; iter++){
995 printf("Itteration=\t%d\n",iter);
0eff8c53 996 Int_t nchunks=gRandom->Poisson(nevents)/eventsPerChunk; // chunks with n typically 25 events
997 for (Int_t i=0; i<nchunks; i++){
0eff8c53 998 tree->GetEntry(gRandom->Rndm()*nhistos);
999 if (i%10==0) printf("%d\t%d\n",iter, i);
1000 if (his1DROCN==0) his1DROCN=new TH1D(*his1DROC);
1001 if (his1DDriftN==0) his1DDriftN=new TH1D(*his1DDrift);
7212f8c3 1002 if (his2DRPhiROCN==0) his2DRPhiROCN=new TH2D(*his2DRPhiROC);
1003 if (his2DRPhiDriftN==0) his2DRPhiDriftN=new TH2D(*his2DRPhiDrift);
0eff8c53 1004 if (his2DRZROCN==0) his2DRZROCN=new TH2D(*his2DRZROC);
1005 if (his2DRZDriftN==0) his2DRZDriftN=new TH2D(*his2DRZDrift);
1006 if (his3DROCN==0) his3DROCN=new TH3D(*his3DROC);
1007 if (his3DDriftN==0) his3DDriftN=new TH3D(*his3DDrift);
1008 his1DROCN->Add(his1DROC);
1009 his1DDriftN->Add(his1DDrift);
7212f8c3 1010 his2DRPhiROCN->Add(his2DRPhiROC);
0eff8c53 1011 his2DRZDriftN->Add(his2DRZDrift);
1012 his2DRZROCN->Add(his2DRZROC);
7212f8c3 1013 his2DRPhiDriftN->Add(his2DRPhiDrift);
0eff8c53 1014 his3DROCN->Add(his3DROC);
1015 his3DDriftN->Add(his3DDrift);
1016 }
1017 //
1018 // 3. Store results 3D maps in the tree (and also as histogram) current and mea
1019 //
1020 Int_t eventsUsed= nchunks*eventsPerChunk;
1021 (*pcstream)<<"fluctuation"<<
1022 "neventsAll="<<neventsAll<< // total number of event to define mean
1023 "nmean="<<nevents<< // mean number of events used
1024 "eventsUsed="<<eventsUsed<< // number of chunks used for one fluct. study
1025 //
1026 // 1,2,3D histogram per group and total
1027 "his1DROCN.="<<his1DROCN<<
1028 "his1DROCSum.="<<his1DROCSum<<
1029 "his1DDriftN.="<<his1DDriftN<<
1030 "his1DDriftSum.="<<his1DDriftSum<<
7212f8c3 1031 "his2DRPhiROCN.="<<his2DRPhiROCN<<
1032 "his2DRPhiROCSum.="<<his2DRPhiROCSum<<
1033 "his2DRPhiDriftN.="<<his2DRPhiDriftN<<
1034 "his2DRPhiDriftSum.="<<his2DRPhiDriftSum<<
0eff8c53 1035 "his2DRZROCN.="<<his2DRZROCN<<
1036 "his2DRZROCSum.="<<his2DRZROCSum<<
1037 "his2DRZDriftN.="<<his2DRZDriftN<<
1038 "his2DRZDriftSum.="<<his2DRZDriftSum<<
1039 "his3DROCN.="<<his3DROCN<<
1040 "his3DROCSum.="<<his3DROCSum<<
1041 "his3DDriftN.="<<his3DDriftN<<
1042 "his3DDriftSum.="<<his3DDriftSum<<
1043 "\n";
1044 pcstream->GetFile()->mkdir(Form("Fluc%d",iter));
1045 pcstream->GetFile()->cd(Form("Fluc%d",iter));
4e715689 1046 //
7212f8c3 1047 his2DRPhiROCN->Write("his2DRPhiROCN");
0eff8c53 1048 his2DRZROCN->Write("his2DRZROCN");
4e715689 1049 //
7212f8c3 1050 his2DRPhiROCSum->Write("his2DRPhiROCSum");
0eff8c53 1051 his2DRZROCSum->Write("his2DRZROCSum");
4e715689 1052 //
7212f8c3 1053 his2DRPhiDriftN->Write("his2DRPhiDriftN");
0eff8c53 1054 his2DRZDriftN->Write("his2DRZDriftN");
4e715689 1055 //
7212f8c3 1056 his2DRPhiDriftSum->Write("his2DRPhiDriftSum");
0eff8c53 1057 his2DRZDriftSum->Write("his2DRZDriftSum");
1058 //
1059 his3DROCN->Write("his3DROCN");
1060 his3DROCSum->Write("his3DROCSum");
1061 his3DDriftN->Write("his3DDriftN");
1062 his3DDriftSum->Write("his3DDriftSum");
1063
1064 his1DROCN->Reset();
1065 his1DDriftN->Reset();
7212f8c3 1066 his2DRPhiROCN->Reset();
0eff8c53 1067 his2DRZDriftN->Reset();
1068 his2DRZROCN->Reset();
7212f8c3 1069 his2DRPhiDriftN->Reset();
0eff8c53 1070 his3DROCN->Reset();
1071 his3DDriftN->Reset();
76dbc35e 1072 }
0eff8c53 1073
1074 delete pcstream;
76dbc35e 1075}
1076
1077
1078void DrawDCARPhiTrendTime(){
1079 //
0eff8c53 1080 // Macros to draw the DCA correlation with the luminosity (estimated from the occupancy)
76dbc35e 1081 //
0eff8c53 1082 // A side and c side 0 differnt behaviour -
1083 // A side - space charge effect
1084 // C side - space charge effect+ FC charging:
1085 // Variables to query from the QA/calibration DB - tree:
1086 // QA.TPC.CPass1.dcar_posA_0 -dca rphi in cm - offset
1087 // Calib.TPC.occQA.Sum() - luminosity is estimated using the mean occupancy per run
1088 //
76dbc35e 1089 TFile *fdb = TFile::Open("outAll.root");
0eff8c53 1090 if (!fdb) fdb = TFile::Open("http://www-alice.gsi.de/TPC/CPassMonitor/outAll.root");
76dbc35e 1091 TTree * tree = (TTree*)fdb->Get("joinAll");
1092 tree->SetCacheSize(100000000);
1093 tree->SetMarkerStyle(25);
1094
1095 //QA.TPC.CPass1.dcar_posA_0 QA.TPC.CPass1.dcar_posA_0_Err QA.TPC.CPass1.meanMult Calib.TPC.occQA. DAQ.L3_magnetCurrent
1096
1097 TGraphErrors * grA = TStatToolkit::MakeGraphErrors(tree,"QA.TPC.CPass1.dcar_posA_0:Calib.TPC.occQA.Sum()*sign(DAQ.L3_magnetCurrent):2*QA.TPC.CPass1.dcar_posA_0_Err","run>190000&&QA.TPC.CPass1.status==1",25,2,0.5);
1098 TGraphErrors * grC = TStatToolkit::MakeGraphErrors(tree,"QA.TPC.CPass1.dcar_posC_0:Calib.TPC.occQA.Sum()*sign(DAQ.L3_magnetCurrent):2*QA.TPC.CPass1.dcar_posC_0_Err","run>190000&&QA.TPC.CPass1.status==1",25,4,0.5);
1099 Double_t mean,rms;
1100 TStatToolkit::EvaluateUni(grA->GetN(),grA->GetY(), mean,rms,grA->GetN()*0.8);
1101 grA->SetMinimum(mean-5*rms);
1102 grA->SetMaximum(mean+3*rms);
1103
1104
1105 grA->GetXaxis()->SetTitle("occ*sign(bz)");
1106 grA->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
1107 grA->Draw("ap");
1108 grC->Draw("p");
1109 TLegend* legend = new TLegend(0.11,0.11,0.5,0.3,"DCA_{rphi} as function of IR (2013)" );
1110 legend->AddEntry(grA,"A side","p");
1111 legend->AddEntry(grC,"C side","p");
1112 legend->Draw();
76dbc35e 1113}
1114
1115
1116
1117void DrawOpenGate(){
1118 //
0eff8c53 1119 // Make nice plot to demonstrate the space charge effect in run with the open gating grid
1120 // For the moment the inmput is harwired - the CPass0 calibration data used
1121 // Make nice drawing (with axis labels):
1122 // To fix (longer term)
1123 // the distortion map to be recalculated - using gaussian fit (currently we use mean)
1124 // the histogram should be extended
76dbc35e 1125 TFile f("/hera/alice/alien/alice/data/2013/LHC13g/000197470/cpass0/OCDB/root_archive.zip#meanITSVertex.root");
1126 TFile fref("/hera/alice/alien/alice/data/2013/LHC13g/000197584/cpass0/OCDB/root_archive.zip#meanITSVertex.root");
1127 //
1128 TTree * treeTOFdy=(TTree*)f.Get("TOFdy");
1129 TTree * treeTOFdyRef=(TTree*)fref.Get("TOFdy");
1130 treeTOFdy->AddFriend(treeTOFdyRef,"R");
1131 treeTOFdy->SetMarkerStyle(25);
1132 TTree * treeITSdy=(TTree*)f.Get("ITSdy");
1133 TTree * treeITSdyRef=(TTree*)fref.Get("ITSdy");
1134 treeITSdy->AddFriend(treeITSdyRef,"R");
1135 treeITSdy->SetMarkerStyle(25);
1136 TTree * treeVertexdy=(TTree*)f.Get("Vertexdy");
1137 TTree * treeVertexdyRef=(TTree*)fref.Get("Vertexdy");
1138 treeVertexdy->AddFriend(treeVertexdyRef,"R");
1139 treeVertexdy->SetMarkerStyle(25);
0eff8c53 1140
1141 // treeITSdy->Draw("mean-R.mean:sector:abs(theta)","entries>50&&abs(snp)<0.1&&theta<0","colz")
76dbc35e 1142
0eff8c53 1143 treeITSdy->Draw("mean-R.mean:sector:abs(theta)","entries>50&&abs(snp)<0.1&&theta>0","colz");
c634db45 1144}
1145
0eff8c53 1146
ae4cea45 1147void DrawCurrent(const char * ocdb="/cvmfs/alice.gsi.de/alice/data/2013/OCDB/TPC/Calib/HighVoltage", Int_t run0=100000, Int_t run1=110000){
c634db45 1148 //
1149 //
1150 /*
1151 const char * ocdb="/cvmfs/alice.gsi.de/alice/data/2013/OCDB/TPC/Calib/HighVoltage";
1152 Int_t run0=197460;
1153 Int_t run1=197480;
1154 */
1155 const Int_t knpoints=100000;
1156 TVectorD vecTime(knpoints);
1157 TVectorD vecI(knpoints);
1158 Int_t npoints=0;
1159 for (Int_t irun=run0; irun<run1; irun++){
1160 TFile * f = TFile::Open(Form("%s/Run%d_%d_v1_s0.root",ocdb,irun,irun));
1161 if (!f) continue;
1162 AliCDBEntry * entry = (AliCDBEntry *)f->Get("AliCDBEntry");
1163 if (!entry) continue;
1164 AliDCSSensorArray * array = (AliDCSSensorArray *)entry->GetObject();
1165 if (!array) continue;
1166 AliDCSSensor * sensor = array->GetSensor("TPC_VHV_D_I_MON");
1167 //sensor->Draw(Form("%d",irun));
1168 TGraph *graph = sensor->GetGraph();
1169 for (Int_t ipoint=0; ipoint<graph->GetN(); ipoint++){
1170 vecTime[npoints]=sensor->GetStartTime()+graph->GetX()[ipoint]*3600;
1171 vecI[npoints]=graph->GetY()[ipoint];
1172 npoints++;
1173 }
1174 }
1175 TGraph * graph = new TGraph(npoints, vecTime.GetMatrixArray(), vecI.GetMatrixArray());
1176 graph->Draw("alp");
1177
76dbc35e 1178
1179}
4e715689 1180
1181
4b955de3 1182void MakeSpaceChargeFluctuationScan(Double_t scale, Int_t nfilesMerge, Int_t sign){
7212f8c3 1183 //
7212f8c3 1184 //
1185 // Input:
633a08cb 1186 // scale - scaling of the space charge (defaul 1 corrsponds to the epsilon ~ 5)
7212f8c3 1187 // nfilesMerge - amount of chunks to merge
1188 // - =0 all chunks used
1189 // <0 subset form full statistic
633a08cb 1190 // >0 subset from the limited (1000 mean) statistic
1191 // Output"
1192 // For given SC setups the distortion on the space point and track level characterezed
1193 // SpaceChargeFluc%d_%d.root - space point distortion maps
1194 // SpaceChargeTrackFluc%d%d.root - tracks distortion caused by space point distortion
1195 //
1196
4b955de3 1197 // Make fluctuation scan:
1198 // 1.) Shift of z disk - to show which granularity in time needed
1199 // 2.) Shift in sector - to show influence of the gass gain and epsilon
1200 // 3.) Smearing in phi - to define phi granularity needed
633a08cb 1201 // 4.) Rebin z - commented out (not delete it for the moment)
1202 // 5.) Rebin phi - commented out
4b955de3 1203
1204
1205 //
1206 // Some constant definition
1207 //
1208 Int_t nitteration=100; // number of itteration in the lookup
1209 Int_t fullNorm =10000; // normalization fro the full statistic
5a998ace 1210 gROOT->ProcessLine(".x $ALICE_ROOT/TPC/Upgrade/macros/ConfigOCDB.C\(1\)");
4b955de3 1211 //
1212 // Init magnetic field and OCDB
1213 //
1214
1215 Double_t bsign= sign;
1216 if (bsign>1) bsign=-1;
2ff913b6 1217 const Int_t nTracks=2000;
4b955de3 1218 const char *ocdb="local://$ALICE_ROOT/OCDB/";
1219 AliCDBManager::Instance()->SetDefaultStorage(ocdb);
1220 AliCDBManager::Instance()->SetRun(0);
1221 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", bsign, bsign, AliMagF::k5kG));
1222 //
1223
1224
1225 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("SpaceChargeFluc%d_%d.root",nfilesMerge,sign),"recreate");
1226 TTreeSRedirector *pcstreamTrack = new TTreeSRedirector(Form("SpaceChargeTrackFluc%d_%d.root",nfilesMerge,sign),"recreate");
7212f8c3 1227 TH1D *his1DROCN=0, *his1DROCSum=0;
1228 TH2D *his2DRPhiROCN=0, *his2DRPhiROCSum=0, *his2DRZROCN=0, *his2DRZROCSum=0;
1229 TH3D *his3DROCN=0, *his3DROCSum=0;
1230 TH1D *his1DROCNC=0, *his1DROCSumC=0;
1231 TH2D *his2DRPhiROCNC=0, *his2DRPhiROCSumC=0, *his2DRZROCNC=0, *his2DRZROCSumC=0;
1232 TH3D *his3DROCNC=0, *his3DROCSumC=0;
1233 TH1 * histos[8]={his1DROCN, his1DROCSum, his2DRPhiROCN, his2DRPhiROCSum, his2DRZROCN, his2DRZROCSum, his3DROCN, his3DROCSum};
1234 Int_t neventsAll=0, neventsAllC=0;
1235 Int_t neventsChunk=0, neventsChunkC=0;
1236 const Double_t ePerADC = 500.;
1237 const Double_t fgke0 = 8.854187817e-12;
1238 //
1239 //
1240 //
1241 const char *inputFile="fluctuation.root";
1242 TObjArray * fileList = (gSystem->GetFromPipe("cat fluctuation.list")).Tokenize("\n");
1243 if (fileList->GetEntries()==0) fileList->AddLast(new TObjString(inputFile));
1244 Int_t nfiles = fileList->GetEntries();
1245 Int_t indexPer[1000];
1246 Double_t numbersPer[10000];
1247 for (Int_t i=0; i<nfiles; i++) numbersPer[i]=gRandom->Rndm();
1248 TMath::Sort(nfiles, numbersPer,indexPer);
1249
1250 for (Int_t ifile=0; ifile<nfiles; ifile++){
1251 if (nfilesMerge>0 && ifile>=nfilesMerge) continue; // merge only limited amount if specified by argument
1252 TFile *fhistos = TFile::Open(fileList->At(indexPer[ifile])->GetName());
1253 if (!fhistos) continue;
1254 TTree * treeHis = (TTree*)fhistos->Get("fluctuation");
1255 if (!treeHis) { printf("file %s does not exist or tree does not exist\n",fileList->At(ifile)->GetName()); continue;}
1256 Int_t nchunks=treeHis->GetEntries();
1257 Int_t chunk=nchunks*gRandom->Rndm();
1258 treeHis->SetBranchAddress("his1DROCN.",&his1DROCNC);
1259 treeHis->SetBranchAddress("his1DROCSum.",&his1DROCSumC);
1260 treeHis->SetBranchAddress("his2DRPhiROCN.",&his2DRPhiROCNC);
1261 treeHis->SetBranchAddress("his2DRPhiROCSum.",&his2DRPhiROCSumC);
1262 treeHis->SetBranchAddress("his2DRZROCN.",&his2DRZROCNC);
1263 treeHis->SetBranchAddress("his2DRZROCSum.",&his2DRZROCSumC);
1264 treeHis->SetBranchAddress("his3DROCN.",&his3DROCNC);
1265 treeHis->SetBranchAddress("his3DROCSum.",&his3DROCSumC);
1266 treeHis->SetBranchAddress("neventsAll",&neventsAllC);
1267 treeHis->SetBranchAddress("eventsUsed",&neventsChunkC);
1268 treeHis->GetEntry(chunk);
1269 neventsAll+=neventsAllC;
1270 neventsChunk+=neventsChunkC;
1271 //
1272 TH1 * histosC[8]={ his1DROCNC, his1DROCSumC, his2DRPhiROCNC, his2DRPhiROCSumC, his2DRZROCNC, his2DRZROCSumC, his3DROCNC, his3DROCSumC};
1273 if (ifile==0) for (Int_t ihis=0; ihis<8; ihis++) histos[ihis] = (TH1*)(histosC[ihis]->Clone());
1274 if (ifile>0) {
1275 for (Int_t ihis=0; ihis<8; ihis++) histos[ihis]->Add(histosC[ihis]);
1276 }
1277 }
1278 his1DROCN=(TH1D*)histos[0]; his1DROCSum=(TH1D*)histos[1];
1279 his2DRPhiROCN=(TH2D*)histos[2]; his2DRPhiROCSum=(TH2D*)histos[3]; his2DRZROCN=(TH2D*)histos[4]; his2DRZROCSum=(TH2D*)histos[5];
1280 his3DROCN=(TH3D*)histos[6]; his3DROCSum=(TH3D*)histos[7];
1281 //
1282 // Select input histogram
1283 //
1284 TH3D * hisInput= his3DROCSum;
1285 Int_t neventsCorr=0; // number of events used for the correction studies
1286 if (nfilesMerge>0){
1287 neventsCorr=neventsChunk;
1288 hisInput=his3DROCN;
1289 }else{
1290 neventsCorr=neventsAll;
1291 hisInput=his3DROCSum;
4b955de3 1292 hisInput->Scale(Double_t(fullNorm)/Double_t(neventsAll));
7212f8c3 1293 }
1294
1295 TObjArray *distortionArray = new TObjArray;
4b955de3 1296 TObjArray *histoArray = new TObjArray;
7212f8c3 1297 //
1298 // Make a reference - ideal distortion/correction
1299 //
1300 TH3D * his3DReference = NormalizeHistoQ(hisInput,kFALSE); // q normalized to the Q/m^3
1301 his3DReference->Scale(scale*0.000001/fgke0); //scale back to the C/cm^3/epsilon0
1302 AliTPCSpaceCharge3D *spaceChargeRef = new AliTPCSpaceCharge3D;
4b955de3 1303 spaceChargeRef->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
7212f8c3 1304 spaceChargeRef->SetInputSpaceCharge(his3DReference, his2DRPhiROCSum,his2DRPhiROCSum,1);
1305 spaceChargeRef->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
7212f8c3 1306 spaceChargeRef->AddVisualCorrection(spaceChargeRef,1);
1307 spaceChargeRef->SetName("DistRef");
4b955de3 1308 his3DReference->SetName("hisDistRef");
7212f8c3 1309 distortionArray->AddLast(spaceChargeRef);
4b955de3 1310 histoArray->AddLast(his3DReference);
1311 //
1312 // Draw histos
1313 TCanvas * canvasSC = new TCanvas("canvasSCDefault","canvasSCdefault",500,400);
1314 canvasSC->SetRightMargin(0.12);
1315 gStyle->SetTitleOffset(0.8,"z");
1316 canvasSC->SetRightMargin(0.13);
1317 spaceChargeRef->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1318 canvasSC->SaveAs(Form("canvasCreateHistoDRPhiinXY_Z10_%d_%d.pdf",nfilesMerge,sign));
1319 spaceChargeRef->CreateHistoDRinXY(10,250,250)->Draw("colz");
1320 canvasSC->SaveAs(Form("canvasCreateHistoDRinXY_Z10_%d_%d.pdf",nfilesMerge,sign));
1321 spaceChargeRef->CreateHistoSCinZR(0.05,250,250)->Draw("colz");
1322 canvasSC->SaveAs(Form("canvasCreateHistoSCinZR_Phi005_%d_%d.pdf",nfilesMerge,sign));
1323 spaceChargeRef->CreateHistoSCinXY(10.,250,250)->Draw("colz");
1324 canvasSC->SaveAs(Form("canvasCreateHistoSCinRPhi_Z10_%d_%d.pdf",nfilesMerge,sign));
1325
1326
7212f8c3 1327 //
1328 // Make Z scan corrections
1329 //
511fd208 1330 if (1){
4b955de3 1331 for (Int_t ihis=1; ihis<=9; ihis+=2){
1332 TH3 *his3DZ = PermutationHistoZ(his3DReference,16*(ihis));
7212f8c3 1333 AliTPCSpaceCharge3D *spaceChargeZ = new AliTPCSpaceCharge3D;
4b955de3 1334 spaceChargeZ->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
7212f8c3 1335 spaceChargeZ->SetInputSpaceCharge(his3DZ, his2DRPhiROCSum,his2DRPhiROCSum,1);
1336 spaceChargeZ->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
1337 spaceChargeZ->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1338 spaceChargeZ->AddVisualCorrection(spaceChargeZ,100+ihis);
4b955de3 1339 spaceChargeZ->SetName(Form("DistZ_%d", 16*(ihis)));
1340 his3DZ->SetName(Form("HisDistZ_%d", 16*(ihis)));
7212f8c3 1341 distortionArray->AddLast(spaceChargeZ);
4b955de3 1342 histoArray->AddLast(his3DZ);
7212f8c3 1343 }
1344 //
1345 // Make Sector scan corrections
1346 //
1347 for (Int_t ihis=1; ihis<=9; ihis+=2){
1348 TH3 *his3DSector = PermutationHistoPhi(his3DReference,TMath::Pi()*(ihis)/9.);
1349 AliTPCSpaceCharge3D *spaceChargeSector = new AliTPCSpaceCharge3D;
4b955de3 1350 spaceChargeSector->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
7212f8c3 1351 spaceChargeSector->SetInputSpaceCharge(his3DSector, his2DRPhiROCSum,his2DRPhiROCSum,1);
1352 spaceChargeSector->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
1353 spaceChargeSector->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1354 spaceChargeSector->AddVisualCorrection(spaceChargeSector,200+ihis);
1355 spaceChargeSector->SetName(Form("DistSector_%d", ihis));
4b955de3 1356 his3DSector->SetName(Form("DistSector_%d", ihis));
7212f8c3 1357 distortionArray->AddLast(spaceChargeSector);
4b955de3 1358 histoArray->AddLast(his3DSector);
1359 }
7212f8c3 1360 //
4b955de3 1361 // Make Local phi scan smear corrections
1362 //
1363 for (Int_t ihis=1; ihis<=8; ihis++){
1364 TH3 *his3DSector = PermutationHistoLocalPhi(his3DReference,ihis);
1365 AliTPCSpaceCharge3D *spaceChargeSector = new AliTPCSpaceCharge3D;
1366 spaceChargeSector->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
1367 spaceChargeSector->SetInputSpaceCharge(his3DSector, his2DRPhiROCSum,his2DRPhiROCSum,1);
1368 spaceChargeSector->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
1369 spaceChargeSector->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1370 spaceChargeSector->AddVisualCorrection(spaceChargeSector,300+ihis);
1371 spaceChargeSector->SetName(Form("DistPhi_%d", ihis));
1372 his3DSector->SetName(Form("HisDistPhi_%d", ihis));
1373 distortionArray->AddLast(spaceChargeSector);
1374 histoArray->AddLast(his3DSector);
1375 }
1376 // //
1377// // Rebin Z
1378// //
1379// for (Int_t ihis=2; ihis<=8; ihis+=2){
1380// TH3 *his3DSector = his3DReference->RebinZ(ihis,Form("RebinZ_%d",ihis));
1381// AliTPCSpaceCharge3D *spaceChargeSector = new AliTPCSpaceCharge3D;
1382// spaceChargeSector->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
1383// spaceChargeSector->SetInputSpaceCharge(his3DSector, his2DRPhiROCSum,his2DRPhiROCSum,1);
1384// spaceChargeSector->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
1385// spaceChargeSector->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1386// spaceChargeSector->AddVisualCorrection(spaceChargeSector,300+ihis);
1387// spaceChargeSector->SetName(Form("RebinZ_%d", ihis));
1388// his3DSector->SetName(Form("RebinZ_%d", ihis));
1389// distortionArray->AddLast(spaceChargeSector);
1390// histoArray->AddLast(his3DSector);
1391// }
1392// //
1393// // Rebin Phi
1394// //
1395// for (Int_t ihis=2; ihis<=5; ihis++){
1396// TH3 *his3DSector = his3DReference->RebinZ(ihis,Form("RebinPhi_%d",ihis));
1397// AliTPCSpaceCharge3D *spaceChargeSector = new AliTPCSpaceCharge3D;
1398// spaceChargeSector->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
1399// spaceChargeSector->SetInputSpaceCharge(his3DSector, his2DRPhiROCSum,his2DRPhiROCSum,1);
1400// spaceChargeSector->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
1401// spaceChargeSector->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1402// spaceChargeSector->AddVisualCorrection(spaceChargeSector,300+ihis);
1403// spaceChargeSector->SetName(Form("RebinZ_%d", ihis));
1404// his3DSector->SetName(Form("RebinZ_%d", ihis));
1405// distortionArray->AddLast(spaceChargeSector);
1406// histoArray->AddLast(his3DSector);
1407// }
1408 }
7212f8c3 1409 //
4b955de3 1410 // Space points scan
7212f8c3 1411 //
1412 Int_t nx = his3DROCN->GetXaxis()->GetNbins();
1413 Int_t ny = his3DROCN->GetYaxis()->GetNbins();
1414 Int_t nz = his3DROCN->GetZaxis()->GetNbins();
1415 Int_t nbins=nx*ny*nz;
1416 TVectorF vx(nbins), vy(nbins), vz(nbins), vq(nbins), vqall(nbins);
1417 //
1418 // charge in the ROC
1419 // for open gate data only fraction of ions enter to drift volume
1420 //
1421 const Int_t kbins=1000;
4b955de3 1422 Double_t deltaR[kbins], deltaZ[kbins],deltaRPhi[kbins], deltaQ[kbins];
7212f8c3 1423 Int_t ndist = distortionArray->GetEntries();
1424 for (Int_t ix=1; ix<=nx; ix+=2){ // phi bin loop
1425 for (Int_t iy=1; iy<=ny; iy+=2){ // r bin loop
1426 Double_t phi= his3DROCN->GetXaxis()->GetBinCenter(ix);
1427 Double_t r = his3DROCN->GetYaxis()->GetBinCenter(iy);
1428 Double_t x = r*TMath::Cos(phi);
1429 Double_t y = r*TMath::Sin(phi);
1430 //
1431 for (Int_t iz=1; iz<=nz; iz++){ // z bin loop
1432 Double_t z = his3DROCN->GetZaxis()->GetBinCenter(iz);
1433 Double_t qN= his3DROCN->GetBinContent(ix,iy,iz);
1434 Double_t qSum= his3DROCSum->GetBinContent(ix,iy,iz);
1435 // Double_t dV in cm = dphi * r * dz in cm**3
1436 Double_t dV= (his3DROCN->GetXaxis()->GetBinWidth(ix)*r)*his3DROCN->GetZaxis()->GetBinWidth(iz);
1437 Double_t norm= 1e6*ePerADC*TMath::Qe()/dV; //normalization factor to the Q/m^3 inside of the ROC;
1438 (*pcstream)<<"hisDump"<<
1439 "neventsAll="<<neventsAll<< // total number of events used for the Q reference
1440 "nfiles="<<nfiles<< // number of files to define properties
1441 "nfilesMerge="<<nfilesMerge<< // number of files to define propertiesneventsCorr
1442 "neventsCorr="<<neventsCorr<< // number of events used to define the corection
4b955de3 1443 "fullNorm="<<fullNorm<< // in case full statistic used this is the normalization coeficient
7212f8c3 1444
1445 "ix="<<ix<<
1446 "iy="<<iy<<
1447 "iz="<<iz<<
1448 // x,y,z
1449 "x="<<x<<
1450 "y="<<y<<
1451 "z="<<z<<
1452 // phi,r,z
1453 "phi="<<phi<<
1454 "r="<<r<<
1455 "z="<<z<<
1456 "norm="<<norm<<
1457 "qN="<<qN<<
1458 "qSum="<<qSum;
1459 for (Int_t idist=0; idist<ndist; idist++){
1460 AliTPCCorrection * corr = (AliTPCCorrection *)distortionArray->At(idist);
4b955de3 1461 TH3 * his = (TH3*)histoArray->At(idist);
7212f8c3 1462 Double_t phi0= TMath::ATan2(y,x);
1463 Int_t nsector=(z>=0) ? 0:18;
1464 Float_t distPoint[3]={x,y,z};
1465 corr->CorrectPoint(distPoint, nsector);
1466 Double_t r0=TMath::Sqrt(x*x+y*y);
1467 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1468 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1469 deltaR[idist] = r1-r0;
1470 deltaRPhi[idist] = (phi1-phi0)*r0;
1471 deltaZ[idist] = distPoint[2]-z;
4b955de3 1472 deltaQ[idist] = his->GetBinContent(ix,iy,iz);
7212f8c3 1473 //
1474 (*pcstream)<<"hisDump"<< //correct point - input point
4b955de3 1475 Form("%sQ=",corr->GetName())<<deltaQ[idist]<<
7212f8c3 1476 Form("%sDR=",corr->GetName())<<deltaR[idist]<<
1477 Form("%sDRPhi=",corr->GetName())<<deltaRPhi[idist]<<
1478 Form("%sDZ=",corr->GetName())<<deltaZ[idist];
1479 }
1480 (*pcstream)<<"hisDump"<<
1481 "\n";
1482 }
1483 }
1484 }
5a998ace 1485 pcstream->GetFile()->cd();
1486 for (Int_t idist=0; idist<ndist; idist++){
1487 AliTPCCorrection * corr = (AliTPCCorrection *)distortionArray->At(idist);
1488 corr->Write(corr->GetName());
1489 }
7212f8c3 1490 delete pcstream;
4b955de3 1491 //
1492 // generate track distortions
1493 //
633a08cb 1494 const Double_t xITSlayer[7]={2.2, 2.8 ,3.6 , 20, 22,41,43 }; // ITS layers R poition (http://arxiv.org/pdf/1304.1306v3.pdf)
1495 const Double_t resITSlayer[7]={0.0004, 0.0004 ,0.0004 , 0.0004, 0.0004, 0.0004, 0.0004 }; // ITS layers R poition (http://arxiv.org/pdf/1304.1306v3.pdf - pixel scenario)
1496 const Double_t kMaxSnp = 0.85;
1497 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
4b955de3 1498 if (nTracks>0){
1499 for(Int_t nt=1; nt<=nTracks; nt++){
1500 gRandom->SetSeed(nt);
1501 TObjArray trackArray(10000);
1502 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
1503 Double_t eta = gRandom->Uniform(-1, 1);
1504 Double_t pt = 1/(gRandom->Rndm()*5+0.00001); // momentum for f1
1505 Short_t psign=1;
1506 if(gRandom->Rndm() < 0.5){
1507 psign =1;
1508 }else{
1509 psign=-1;
1510 }
1511 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1512 Double_t pxyz[3];
1513 pxyz[0]=pt*TMath::Cos(phi);
1514 pxyz[1]=pt*TMath::Sin(phi);
1515 pxyz[2]=pt*TMath::Tan(theta);
1516 Double_t vertex[3]={0,0,0};
1517 Double_t cv[21]={0};
1518 AliExternalTrackParam *t= new AliExternalTrackParam(vertex, pxyz, cv, psign);
1519 Double_t refX0=85.;
1520 Double_t refX1=1.;
1521 Int_t dir=-1;
1522 (*pcstreamTrack)<<"trackFit"<<
5e4083cf 1523 "bsign="<<bsign<<
4b955de3 1524 "neventsAll="<<neventsAll<< // total number of events used for the Q reference
1525 "nfiles="<<nfiles<< // number of files to define properties
1526 "nfilesMerge="<<nfilesMerge<< // number of files to define propertiesneventsCorr
1527 "neventsCorr="<<neventsCorr<< // number of events used to define the corection
1528 "fullNorm="<<fullNorm<< // in case full statistic used this is the normalization coeficient
633a08cb 1529 "itrack="<<nt<< //
1530 "input.="<<t; //
5e4083cf 1531
1532 Bool_t isOK0=kTRUE;
1533 Bool_t isOK1=kTRUE;
1534 Bool_t itsOK=kTRUE;
1535 Bool_t itsUpdateOK=kTRUE;
1536
4b955de3 1537 for (Int_t idist=0; idist<ndist; idist++){
1538 AliTPCCorrection * corr = (AliTPCCorrection *)distortionArray->At(idist);
633a08cb 1539 // 0. TPC only information at the entrance
1540 // 1. TPC only information close to vertex ( refX=1 cm)
1541 // 2. TPC constrained information close to the primary vertex
1542 // 3. TPC +ITS
4b955de3 1543 AliExternalTrackParam *ot0= new AliExternalTrackParam(vertex, pxyz, cv, psign);
1544 AliExternalTrackParam *ot1= new AliExternalTrackParam(vertex, pxyz, cv, psign);
1545 AliExternalTrackParam *td0 = corr->FitDistortedTrack(*ot0, refX0, dir, 0);
1546 AliExternalTrackParam *td1 = corr->FitDistortedTrack(*ot1, refX1, dir, 0);
f0085660 1547 if (td0==0) { // if fit fail use dummy values
1548 ot0= new AliExternalTrackParam(vertex, pxyz, cv, psign);
1549 td0= new AliExternalTrackParam(vertex, pxyz, cv, psign);
5e4083cf 1550 printf("Propagation0 failed: track\t%d\tdistortion\t%d\n",nt,idist);
1551 isOK0=kFALSE;
f0085660 1552 }
1553 if (td1==0) {
1554 ot1= new AliExternalTrackParam(vertex, pxyz, cv, psign);
1555 td1= new AliExternalTrackParam(vertex, pxyz, cv, psign);
5e4083cf 1556 printf("Propagation1 failed: track\t%d\tdistortion\t%d\n",nt,idist);
1557 isOK1=kFALSE;
f0085660 1558 }
5e4083cf 1559 // 2. TPC constrained emulation
633a08cb 1560 AliExternalTrackParam *tdConstrained = new AliExternalTrackParam(*td1);
1561 tdConstrained->Rotate(ot1->GetAlpha());
1562 tdConstrained->PropagateTo(ot1->GetX(), AliTrackerBase::GetBz());
1563 Double_t pointPos[2]={ot1->GetY(),ot1->GetZ()}; // local y and local z of point
1564 Double_t pointCov[3]={0.0001,0,0.0001}; //
1565 tdConstrained->Update(pointPos,pointCov);
1566 // 3. TPC+ITS constrained umulation
1567 AliExternalTrackParam *tdITS = new AliExternalTrackParam(*td0);
1568 AliExternalTrackParam *tdITSOrig = new AliExternalTrackParam(*ot0);
633a08cb 1569 //
5e4083cf 1570 if ( isOK0 && isOK1 ) {
1571 for (Int_t ilayer=6; ilayer>=0; ilayer--){
1572 if (!AliTrackerBase::PropagateTrackTo(tdITSOrig,xITSlayer[ilayer],kMass,5.,kTRUE,kMaxSnp)) itsOK=kFALSE;
1573 if (!AliTrackerBase::PropagateTrackTo(tdITS,xITSlayer[ilayer],kMass,5.,kTRUE,kMaxSnp)) {
1574 itsOK=kFALSE;
1575 printf("PropagationITS failed: track\t%d\tdistortion\t%d\t%d\n",nt,idist,ilayer);
1576 }
1577 //
1578 tdITS->Rotate(tdITSOrig->GetAlpha());
1579 if (tdITS->PropagateTo(tdITSOrig->GetX(), AliTrackerBase::GetBz())){
1580 Double_t itspointPos[2]={tdITSOrig->GetY(),tdITSOrig->GetZ()}; // local y and local z of point
1581 Double_t itspointCov[3]={resITSlayer[ilayer]*resITSlayer[ilayer],0,resITSlayer[ilayer]*resITSlayer[ilayer]};
1582 if (!tdITS->Update(itspointPos,itspointCov)){
1583 itsUpdateOK=kFALSE;
1584 }
1585 }
1586 }
1587 }else{
1588 itsOK=kFALSE;
633a08cb 1589 }
1590 //
4b955de3 1591 trackArray.AddLast(td0);
1592 trackArray.AddLast(td1);
633a08cb 1593 trackArray.AddLast(tdConstrained);
1594 trackArray.AddLast(tdITS);
1595 trackArray.AddLast(tdITSOrig);
1596 //
4b955de3 1597 trackArray.AddLast(ot0);
1598 trackArray.AddLast(ot1);
633a08cb 1599 char name0[100], name1[1000], nameITS[1000];
1600 char oname0[100], oname1[1000], onameConstrained[1000], onameITS[1000];
4b955de3 1601 snprintf(name0, 100, "T_%s_0.=",corr->GetName());
1602 snprintf(name1, 100, "T_%s_1.=",corr->GetName());
1603 snprintf(oname0, 100, "OT_%s_0.=",corr->GetName());
f0085660 1604 snprintf(oname1, 100, "OT_%s_1.=",corr->GetName());
633a08cb 1605 snprintf(onameConstrained, 100, "OConst_%s_1.=",corr->GetName());
1606 //
1607 snprintf(nameITS, 100, "TPCITS_%s_1.=",corr->GetName());
1608 snprintf(onameITS, 100, "OTPCITS_%s_1.=",corr->GetName());
4b955de3 1609 (*pcstreamTrack)<<"trackFit"<<
633a08cb 1610 name0<<td0<< // distorted TPC track only at the refX=85
1611 name1<<td1<< // distorted TPC track only at the refX=1
1612 onameConstrained<<tdConstrained<< // distorted TPC constrained track only at the refX=1
1613 //
1614 onameITS<<tdITSOrig<< // original TPC+ITS track
1615 nameITS<<tdITS<< // distorted TPC+ (indistrted)ITS track fit
1616 //
1617 oname0<<ot0<< // original track at the entrance refX=85
1618 oname1<<ot1; // original track at the refX=1 cm (to be used for TPC only and also for the constrained
1619
4b955de3 1620 }
5e4083cf 1621 (*pcstreamTrack)<<"trackFit"<<
1622 "isOK0="<<isOK0<< // propagation good at the inner field cage
1623 "isOK1="<<isOK1<< // god at 1 cm (close to vertex)
1624 "itsOK="<<itsOK<< //
1625 "itsUpdateOK="<<itsOK<< //
1626 "\n";
4b955de3 1627 }
1628 }
1629 delete pcstreamTrack;
7212f8c3 1630 return;
1631
1632}
1633
1634
1635void MakePlotPoisson3D(const char *inputfile="fluctuation.root", const char *outputfile="SpaceCharge.root", Int_t event=0){
4e715689 1636 //
7212f8c3 1637 // draw "standard" plot to show radial and theta dependence of the space charge distortion
4e715689 1638 //
7212f8c3 1639 // const char *inputfile="fluctuation.root"; const char *outputfile="SpaceCharge.root"; Int_t event=0
4e715689 1640 //
1641 TFile *fhistos = TFile::Open(inputfile);
7212f8c3 1642 TH2D *his2DRPhiROCN=0, *his2DRPhiROCSum=0, *his2DRZROCN=0, *his2DRZROCSum=0;
1643 TH1D *his1DROCN=0, *his1DROCSum=0;
1644 TH3D *his3DROCN=0, *his3DROCSum=0;
1645 const Double_t ePerADC = 500.;
1646 const Double_t fgke0 = 8.854187817e-12;
4e715689 1647 TTree * treeHis = (TTree*)fhistos->Get("fluctuation");
7212f8c3 1648 treeHis->SetBranchAddress("his1DROCN.",&his1DROCN);
1649 treeHis->SetBranchAddress("his1DROCSum.",&his1DROCSum);
1650 treeHis->SetBranchAddress("his2DRPhiROCN.",&his2DRPhiROCN);
1651 treeHis->SetBranchAddress("his2DRPhiROCSum.",&his2DRPhiROCSum);
4e715689 1652 treeHis->SetBranchAddress("his2DRZROCN.",&his2DRZROCN);
1653 treeHis->SetBranchAddress("his2DRZROCSum.",&his2DRZROCSum);
7212f8c3 1654 treeHis->SetBranchAddress("his3DROCN.",&his3DROCN);
1655 treeHis->SetBranchAddress("his3DROCSum.",&his3DROCSum);
4e715689 1656 treeHis->GetEntry(event);
7212f8c3 1657
1658 his3DROCSum->Scale(ePerADC*TMath::Qe()/fgke0);
1659
1660 AliTPCSpaceCharge3D *spaceChargeOrig = new AliTPCSpaceCharge3D;
1661 spaceChargeOrig->SetOmegaTauT1T2(0.0,1,1); // Ne CO2
1662 spaceChargeOrig->SetInputSpaceCharge(his3DROCSum, his2DRPhiROCSum,his2DRPhiROCSum,10*ePerADC*TMath::Qe());
1663 spaceChargeOrig->InitSpaceCharge3DPoisson(129, 129, 144,100);
1664 spaceChargeOrig->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1665 spaceChargeOrig->AddVisualCorrection(spaceChargeOrig,1);
4e715689 1666 //
4b955de3 1667 //AliTPCSpaceCharge3D *spaceChargeRef= spaceChargeOrig;
1668
1669
7212f8c3 1670
4e715689 1671 //
7212f8c3 1672 Int_t nfuns=5;
1673 Double_t dmax=0.75, dmin=-0.75;
1674 Double_t phiRange=18;
1675 TCanvas *canvasDistortionP3D = new TCanvas("canvasdistortionP3D","canvasdistortionP3D",1000,700);
1676 canvasDistortionP3D->SetGrid(1,1);
1677 canvasDistortionP3D->Divide(1,2);
1678 canvasDistortionP3D->cd(1)->SetGrid(1,1);
1679 TLegend * legendR= new TLegend(0.11,0.11,0.45,0.35,"R scan (#Theta=0.1)");
1680 for (Int_t ifun1=0; ifun1<=nfuns; ifun1++){
1681 Double_t rfun= 85.+ifun1*(245.-85.)/nfuns;
1682 TF1 *pf1 = new TF1("f1",Form("AliTPCCorrection::GetCorrSector(x,%f,0.1,1,1)",rfun),0,phiRange);
1683 pf1->SetMinimum(dmin);
1684 pf1->SetMaximum(dmax);
1685 pf1->SetNpx(360);
1686 pf1->SetLineColor(1+ifun1);
1687 pf1->SetLineWidth(2);
1688 pf1->GetXaxis()->SetTitle("sector");
1689 pf1->GetXaxis()->SetNdivisions(530);
1690 pf1->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
1691 if (ifun1==0) pf1->Draw();
1692 pf1->Draw("same");
1693 legendR->AddEntry(pf1,Form("r=%1.0f",rfun));
1694 }
1695 legendR->Draw();
4e715689 1696 //
7212f8c3 1697 canvasDistortionP3D->cd(2)->SetGrid(1,1);
1698 TLegend * legendTheta= new TLegend(0.11,0.11,0.45,0.35,"#Theta scan (r=125 cm)");
1699 for (Int_t ifun1=0; ifun1<=nfuns; ifun1++){
1700 Double_t tfun= 0.1+ifun1*(0.8)/nfuns;
1701 TF1 *pf1 = new TF1("f1",Form("AliTPCCorrection::GetCorrSector(x,125,%f,1,1)",tfun),0,phiRange);
1702 pf1->SetMinimum(dmin);
1703 pf1->SetMaximum(dmax);
1704 pf1->SetNpx(360);
1705 pf1->SetLineColor(1+ifun1);
1706 pf1->SetLineWidth(2);
1707 pf1->GetXaxis()->SetTitle("sector");
1708 pf1->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
1709 pf1->GetXaxis()->SetNdivisions(530);
1710 if (ifun1==0) pf1->Draw();
1711 pf1->Draw("same");
1712 legendTheta->AddEntry(pf1,Form("#Theta=%1.2f",tfun));
1713 }
1714 legendTheta->Draw();
4e715689 1715
7212f8c3 1716}
1717
1718TH3D * NormalizeHistoQ(TH3D * hisInput, Bool_t normEpsilon){
1719 //
1720 // Renormalize the histogram to the Q/m^3
1721 // Input:
1722 // hisInput - input 3D histogram
1723 // normEpsilon - flag - normalize to epsilon0
1724 //
1725 const Double_t ePerADC = 500.;
1726 const Double_t fgkEpsilon0 = 8.854187817e-12;
1727 TH3D * hisOutput= new TH3D(*hisInput);
1728 Int_t nx = hisInput->GetXaxis()->GetNbins();
1729 Int_t ny = hisInput->GetYaxis()->GetNbins();
1730 Int_t nz = hisInput->GetZaxis()->GetNbins();
1731 for (Int_t ix=1; ix<=nx; ix++){
1732 for (Int_t iy=1; iy<=ny; iy++){
1733 for (Int_t iz=1; iz<=nz; iz++){
1734 // Double_t z = hisInput->GetZaxis()->GetBinCenter(iz);
1735 Double_t deltaRPhi = hisInput->GetXaxis()->GetBinWidth(ix)* hisInput->GetYaxis()->GetBinCenter(iy);
1736 Double_t deltaR= hisInput->GetYaxis()->GetBinWidth(iy);
1737 Double_t deltaZ= hisInput->GetYaxis()->GetBinWidth(iz);
1738 Double_t volume= (deltaRPhi*deltaR*deltaZ)/1000000.;
1739 Double_t q = hisInput->GetBinContent(ix,iy,iz)* ePerADC*TMath::Qe(); // Q in coulombs
1740 Double_t rho = q/volume; // rpho - density in Q/m^3
1741 if (normEpsilon) rho/=fgkEpsilon0;
1742 hisOutput->SetBinContent(ix,iy,iz,rho);
4e715689 1743 }
1744 }
1745 }
7212f8c3 1746 return hisOutput;
1747}
1748
1749
1750
1751TH3D * PermutationHistoZ(TH3D * hisInput, Double_t deltaZ){
1752 //
1753 // Used to estimate the effect of the imperfection of the lookup tables as function of update frequency
4e715689 1754 //
7212f8c3 1755 // Permute/rotate the conten of the histogram in z direction
1756 // Reshufle/shift content - Keeping the integral the same
1757 // Parameters:
1758 // hisInput - input 3D histogram (phi,r,z)
1759 // deltaZ - deltaZ -shift of the space charge
1760 Double_t zmax=250;
1761 TH3D * hisOutput= new TH3D(*hisInput);
1762 Int_t nx = hisInput->GetXaxis()->GetNbins();
1763 Int_t ny = hisInput->GetYaxis()->GetNbins();
1764 Int_t nz = hisInput->GetZaxis()->GetNbins();
4e715689 1765 //
7212f8c3 1766 //
1767 for (Int_t ix=1; ix<=nx; ix++){
1768 for (Int_t iy=1; iy<=ny; iy++){
1769 for (Int_t iz=1; iz<=nz; iz++){
1770 Double_t zold = hisInput->GetZaxis()->GetBinCenter(iz);
1771 Double_t z=zold;
1772 if (z>0){
1773 z+=deltaZ;
1774 if (z<0) z+=zmax;
1775 if (z>zmax) z-=zmax;
1776 }else{
1777 z-=deltaZ;
1778 if (z>0) z-=zmax;
1779 if (z<-zmax) z+=zmax; }
1780 Double_t kz= hisInput->GetZaxis()->FindBin(z);
1781 Double_t content = hisInput->GetBinContent(ix,iy,iz);
1782 hisOutput->SetBinContent(ix,iy,kz,content);
4e715689 1783 }
1784 }
1785 }
7212f8c3 1786 return hisOutput;
1787}
1788
4b955de3 1789
1790
1791
7212f8c3 1792TH3D * PermutationHistoPhi(TH3D * hisInput, Double_t deltaPhi){
1793 //
1794 // Used to estimate the effect of the imperfection of the lookup tables as function of update frequency
1795 //
1796 // Permute/rotate the conten of the histogram in phi
1797 // Reshufle/shift content - Keeping the integral the same
1798 // Parameters:
1799 // hisInput - input 3D histogram (phi,r,z)
1800 // deltaPhi - deltaPhi -shift of the space charge
1801 TH3D * hisOutput= new TH3D(*hisInput);
1802 Int_t nx = hisInput->GetXaxis()->GetNbins();
1803 Int_t ny = hisInput->GetYaxis()->GetNbins();
1804 Int_t nz = hisInput->GetZaxis()->GetNbins();
4e715689 1805 //
4e715689 1806 //
7212f8c3 1807 for (Int_t iy=1; iy<=ny; iy++){
1808 for (Int_t iz=1; iz<=nz; iz++){
1809 for (Int_t ix=1; ix<=nx; ix++){
1810 Double_t phiOld = hisInput->GetXaxis()->GetBinCenter(ix);
1811 Double_t phi=phiOld;
1812 phi+=deltaPhi;
1813 if (phi<0) phi+=TMath::TwoPi();
1814 if (phi>TMath::TwoPi()) phi-=TMath::TwoPi();
1815 Double_t kx= hisInput->GetXaxis()->FindBin(phi);
1816 Double_t content = hisInput->GetBinContent(ix,iy,iz);
1817 hisOutput->SetBinContent(kx,iy,iz,content);
1818 }
4e715689 1819 }
1820 }
7212f8c3 1821 return hisOutput;
1822}
4e715689 1823
7212f8c3 1824
4b955de3 1825TH3D * PermutationHistoLocalPhi(TH3D * hisInput, Int_t deltaPhi){
7212f8c3 1826 //
1827 // Used to estimate the effect of the imperfection of the lookup tables as function of update frequency
4b955de3 1828 // Use moving average of the content instead of the content
7212f8c3 1829 //
7212f8c3 1830 // Parameters:
1831 // hisInput - input 3D histogram (phi,r,z)
4b955de3 1832 // deltaPhi - moving average width
7212f8c3 1833 TH3D * hisOutput= new TH3D(*hisInput);
1834 Int_t nx = hisInput->GetXaxis()->GetNbins();
1835 Int_t ny = hisInput->GetYaxis()->GetNbins();
1836 Int_t nz = hisInput->GetZaxis()->GetNbins();
4b955de3 1837 Int_t binSector=nx/18;
7212f8c3 1838 //
1839 //
1840 for (Int_t iy=1; iy<=ny; iy++){
1841 for (Int_t iz=1; iz<=nz; iz++){
1842 for (Int_t ix=1; ix<=nx; ix++){
4b955de3 1843 Double_t sumRo=0,sumW=0;
1844 for (Int_t idx=-deltaPhi; idx<=deltaPhi; idx++){
1845 Int_t index=ix+idx;
1846 if (index<1) index+=nx+1; // underflow and overflow bins
1847 if (index>nx) index-=nx+1;
1848 Double_t content = hisInput->GetBinContent(index,iy,iz);
1849 sumRo+=content;
1850 sumW++;
1851 }
1852 Double_t meanCont= sumRo/sumW;
1853 hisOutput->SetBinContent(ix,iy,iz,meanCont);
1854 //printf("%d\t%f\n",ix,hisInput->GetBinContent(ix,iy,iz)/(hisInput->GetBinContent(ix,iy,iz)+meanCont));
1855 }
7212f8c3 1856 }
1857 }
1858 return hisOutput;
4e715689 1859}
1860
1861
4e715689 1862
7212f8c3 1863void ScanIterrationPrecision(TH3 * hisInput, Int_t offset){
1864 //
1865 //
1866 //
1867 for (Int_t iter=0; iter<=7; iter++){
1868 Int_t niter= 50.*TMath::Power(1.5,iter);
1869 AliTPCSpaceCharge3D *spaceChargeOrig = new AliTPCSpaceCharge3D;
1870 spaceChargeOrig->SetOmegaTauT1T2(0.0,1,1); // Ne CO2
1871 spaceChargeOrig->SetInputSpaceCharge(hisInput,0,0,1);
1872 spaceChargeOrig->InitSpaceCharge3DPoisson(129, 129, 144,niter);
1873 spaceChargeOrig->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1874 spaceChargeOrig->AddVisualCorrection(spaceChargeOrig,offset+iter+1);
1875 }
1876}
1877
1878
1879void DrawFluctuationSector(Int_t stat, Double_t norm){
4e715689 1880 //
7212f8c3 1881 // Draw correction - correction at rotated sector
1882 // The same set of events used
1883 // Int_t stat=0; Double_t norm=10000;
4b955de3 1884 //
1885 // Notes:
1886 // 1. (something wrong for the setup 2 pileups -problem with data 24.07
1887 //
1888 //
7212f8c3 1889 TFile *f0= TFile::Open(Form("SpaceChargeFluc%d.root",stat));
1890 TTree * tree0 = (TTree*)f0->Get("hisDump");
1891 tree0->SetCacheSize(1000000000);
1892 tree0->SetMarkerStyle(25);
1893 TObjArray * fitArray=new TObjArray(3);
1894 tree0->SetAlias("scNorm",Form("%f/neventsCorr",norm));
4e715689 1895 //
7212f8c3 1896 // Sector Scan
1897 //
1898 TH2 * hisSectorScan[5]={0};
1899 TH1 * hisSectorScanSigma[5]={0};
1900 for (Int_t ihis=0; ihis<5; ihis++){
1901 tree0->Draw(Form("(DistRefDR-DistSector_%dDR)*scNorm:r>>hisSec%d(50,84,245,100,-1,1)",ihis*2+1,ihis*2+1),"abs(z)<90","colzgoff");
1902 hisSectorScan[ihis]=(TH2*)tree0->GetHistogram();
1903 hisSectorScan[ihis]->FitSlicesY(0,0,-1,0,"QNR",fitArray);
1904 hisSectorScanSigma[ihis]=(TH1*)(fitArray->At(2)->Clone());
1905 hisSectorScanSigma[ihis]->SetMinimum(0);
1906 hisSectorScanSigma[ihis]->SetMaximum(0.2);
1907 }
1908 gStyle->SetOptStat(0);
1909 gStyle->SetOptTitle(1);
1910 TCanvas * canvasFlucSectorScan=new TCanvas("canvasFlucSectorScan","canvasFlucSectorScan",750,700);
1911 canvasFlucSectorScan->Divide(2,2,0,0);
1912 gStyle->SetPadBorderMode(0);
1913 for (Int_t ihis=0; ihis<4; ihis++){
1914 canvasFlucSectorScan->cd(ihis+1)->SetLogz(1);
1915 hisSectorScan[ihis]->GetXaxis()->SetTitle("r (cm)");
1916 hisSectorScan[ihis]->GetYaxis()->SetTitle("#Delta_{R} (cm)");
1917 hisSectorScan[ihis]->Draw("colz");
1918 TLegend * legendSec=new TLegend(0.5,0.7,0.89,0.89);
1919 legendSec->AddEntry(hisSectorScan[ihis],Form("Sector #Delta %d",(ihis*2+1)));
1920 legendSec->Draw();
1921 }
1922 canvasFlucSectorScan->SaveAs("canvasFlucSectorScan.pdf");
1923 canvasFlucSectorScan->SaveAs("canvasFlucSectorScan.png");
4e715689 1924 //
7212f8c3 1925 gStyle->SetOptTitle(0);
1926 TCanvas * canvasFlucSectorScanFit=new TCanvas("canvasFlucSectorScanFit","canvasFlucSectorScanFit",750,550);
1927 TLegend * legendSector = new TLegend(0.50,0.55,0.89,0.89,"Space charge: corr(sec)-corr(sec-#Delta_{sec})");
1928 for (Int_t ihis=0; ihis<5; ihis++){
1929 hisSectorScanSigma[ihis]->GetXaxis()->SetTitle("r (cm)");
1930 hisSectorScanSigma[ihis]->GetYaxis()->SetTitle("#sigma(#Delta_{R}) (cm)");
1931 hisSectorScanSigma[ihis]->SetMarkerStyle(21+ihis%5);
1932 hisSectorScanSigma[ihis]->SetMarkerColor(1+ihis%4);
1933 if (ihis==0) hisSectorScanSigma[ihis]->Draw("");
1934 hisSectorScanSigma[ihis]->Draw("same");
1935 legendSector->AddEntry(hisSectorScanSigma[ihis],Form("#Delta %d",(ihis*2+1)));
1936 }
1937 legendSector->Draw();
1938 canvasFlucSectorScanFit->SaveAs("canvasFlucSectorScanFit.pdf");
1939 canvasFlucSectorScanFit->SaveAs("canvasFlucSectorScanFit.png");
1940}
4e715689 1941
1942
4e715689 1943
7212f8c3 1944void DrawFluctuationdeltaZ(Int_t stat, Double_t norm){
1945 //
1946 // Draw correction - correction shifted z
1947 // The same set of events used
1948 //Int_t stat=0; Double_t norm=10000;
4b955de3 1949 Int_t deltaZ=16.;
7212f8c3 1950 TFile *f0= TFile::Open(Form("SpaceChargeFluc%d.root",stat));
1951 TTree * tree0 = 0;
1952 if (f0) tree0 = (TTree*)f0->Get("hisDump");
1953 if (!tree0){
4b955de3 1954 tree0 = AliXRDPROOFtoolkit::MakeChainRandom("space.list","hisDump",0,10);
7212f8c3 1955 }
1956 tree0->SetCacheSize(1000000000);
1957 tree0->SetMarkerStyle(25);
1958 TObjArray * fitArray=new TObjArray(3);
1959 tree0->SetAlias("scNorm",Form("%f/neventsCorr",norm));
1960 //
1961 // DeltaZ Scan
1962 //
1963 TH2 * hisDeltaZScan[6]={0};
1964 TH1 * hisDeltaZScanSigma[6]={0};
1965 for (Int_t ihis=0; ihis<6; ihis++){
4b955de3 1966 tree0->Draw(Form("(DistRefDR-DistZ_%dDR)*scNorm:r>>hisZ%d(50,84,245,100,-1,1)",(ihis+1)*deltaZ,(ihis+1)*deltaZ),"abs(z/r)<1","colzgoff");
7212f8c3 1967 hisDeltaZScan[ihis]=(TH2*)tree0->GetHistogram();
1968 hisDeltaZScan[ihis]->FitSlicesY(0,0,-1,0,"QNR",fitArray);
1969 hisDeltaZScanSigma[ihis]=(TH1*)(fitArray->At(2)->Clone());
1970 hisDeltaZScanSigma[ihis]->SetMinimum(0);
1971 hisDeltaZScanSigma[ihis]->SetMaximum(0.2);
1972 }
1973 gStyle->SetOptStat(0);
1974 gStyle->SetOptTitle(1);
1975 TCanvas * canvasFlucDeltaZScan=new TCanvas("canvasFlucDeltaZScan","canvasFlucDeltaZScan",700,700);
1976 canvasFlucDeltaZScan->Divide(3,2,0,0);
1977 gStyle->SetPadBorderMode(0);
1978 for (Int_t ihis=0; ihis<6; ihis++){
1979 canvasFlucDeltaZScan->cd(ihis+1)->SetLogz(1);
1980 hisDeltaZScan[ihis]->GetXaxis()->SetTitle("r (cm)");
1981 hisDeltaZScan[ihis]->GetYaxis()->SetTitle("#Delta_{R} (cm)");
1982 hisDeltaZScan[ihis]->Draw("colz");
1983 TLegend * legendSec=new TLegend(0.5,0.7,0.89,0.89);
4b955de3 1984 legendSec->AddEntry(hisDeltaZScan[ihis],Form("DeltaZ #Delta %d",(ihis+1)*deltaZ));
7212f8c3 1985 legendSec->Draw();
1986 }
1987 canvasFlucDeltaZScan->SaveAs(Form("canvasFlucDeltaZScan%d.pdf",stat));
1988 canvasFlucDeltaZScan->SaveAs(Form("canvasFlucDeltaZScan%d.png",stat));
1989
1990 //
1991 gStyle->SetOptTitle(0);
1992 TCanvas * canvasFlucDeltaZScanFit=new TCanvas("canvasFlucDeltaZScanFit","canvasFlucDeltaZScanFit");
1993 TLegend * legendDeltaZ = new TLegend(0.50,0.55,0.89,0.89,"Space charge: corr(z_{ref})-corr(z_{ref}-#Delta_{z})");
1994 for (Int_t ihis=0; ihis<5; ihis++){
1995 hisDeltaZScanSigma[ihis]->GetXaxis()->SetTitle("r (cm)");
1996 hisDeltaZScanSigma[ihis]->GetYaxis()->SetTitle("#sigma(#Delta_{R}) (cm)");
1997 hisDeltaZScanSigma[ihis]->SetMarkerStyle(21+ihis%5);
1998 hisDeltaZScanSigma[ihis]->SetMarkerColor(1+ihis%4);
1999 if (ihis==0) hisDeltaZScanSigma[ihis]->Draw("");
2000 hisDeltaZScanSigma[ihis]->Draw("same");
4b955de3 2001 legendDeltaZ->AddEntry(hisDeltaZScanSigma[ihis],Form("#Delta %d (cm)",(ihis+1)*deltaZ));
7212f8c3 2002 }
2003 legendDeltaZ->Draw();
2004 canvasFlucDeltaZScanFit->SaveAs(Form("canvasFlucDeltaZScanFit%d.pdf",stat));
2005 canvasFlucDeltaZScanFit->SaveAs(Form("canvasFlucDeltaZScanFit%d.png",stat));
4e715689 2006
7212f8c3 2007}
4e715689 2008
4e715689 2009
7212f8c3 2010void DrawDefault(Int_t stat){
2011 //
2012 // Draw correction - correction shifted z
2013 // The same set of events used
2014 // Int_t stat=0
2015 TFile *f0= TFile::Open(Form("SpaceChargeFluc%d.root",stat));
2016 TTree * tree0 = (TTree*)f0->Get("hisDump");
2017 tree0->SetCacheSize(1000000000);
2018 tree0->SetMarkerStyle(25);
2019 tree0->SetMarkerSize(0.4);
4b955de3 2020 // TObjArray * fitArray=new TObjArray(3);
7212f8c3 2021 tree0->Draw("10000*DistRefDR/neventsCorr:r:z/r","abs(z/r)<0.9&&z>0&&rndm>0.8","colz");
4b955de3 2022}
2023
2024
2025
2026
2027void DrawTrackFluctuation(){
2028 //
2029 // Function to make a fluctuation figures for differnt multiplicities of pileup space charge
2030 // it is assumed that the text files
2031 //
2032 //
2033 TObjArray arrayFit(3);
2034 const char *inputList;
2035 TH2F * hisCorrRef[5]={0};
2036 TH2F * hisCorrNo[5]={0};
2037 TH1 * hisCorrRefM[5], *hisCorrRefRMS[5];
2038 TH1 * hisCorrNoM[5], *hisCorrNoRMS[5];
2039 //
2040 // 1. Load chains for different statistic
2041 //
2042 TCut cutOut="abs(T_DistRef_0.fX-OT_DistRef_0.fX)<0.1&&T_DistRef_0.fX>1&&abs(OT_DistRef_0.fP[4])<4";
2043 TCut cutOutF="abs(R.T_DistRef_0.fX-R.OT_DistRef_0.fX)<0.1&&R.T_DistRef_0.fX>1&&abs(R.OT_DistRef_0.fP[4])<4";
2044 TChain * chains[5]={0};
5e4083cf 2045 TChain * chainR = AliXRDPROOFtoolkit::MakeChain("track0_1.list","trackFit",0,1000); // list of the reference data (full stat used)
4b955de3 2046 chainR->SetCacheSize(1000000000);
5e4083cf 2047 for (Int_t ichain=0; ichain<5; ichain++){ // create the chain for given mulitplicity bin
4b955de3 2048 chains[ichain] = AliXRDPROOFtoolkit::MakeChain(Form("track%d_1.list",2*(ichain+1)),"trackFit",0,1000);
2049 chains[ichain]->AddFriend(chainR,"R");
2050 chains[ichain]->SetCacheSize(1000000000);
2051 chains[ichain]->SetMarkerStyle(25);
2052 chains[ichain]->SetMarkerSize(0.5);
2053 chains[ichain]->SetAlias("meanNorm","(1+0.2*abs(neventsCorr/10000-1))"); // second order correction - renomalization of mean hardwired
2054 // to be fitted?
2055 }
2056 //
2057 // 2. fill histograms if not available in file
2058 //
2059 //
2060 TFile *ftrackFluctuation = TFile::Open("trackFluctuation.root","update");
2061 for (Int_t ihis=0; ihis<5; ihis++){
2062 ftrackFluctuation->cd();
2063 hisCorrRef[ihis] = (TH2F*)(ftrackFluctuation->Get(Form("DeltaRPhiCorr%d",(ihis+1)*2000)));
2064 hisCorrNo[ihis] = (TH2F*)(ftrackFluctuation->Get(Form("DeltaRPhi%d",(ihis+1)*2000)));
2065 if (hisCorrRef[ihis]==0) {
2066 chains[ihis]->Draw("(T_DistRef_0.fP[0]/meanNorm-neventsCorr*R.T_DistRef_0.fP[0]/10000):R.OT_DistRef_0.fP[4]>>his(10,-4,4,100,-0.25,0.25)",cutOut+cutOutF+"","colzgoff");
2067 hisCorrRef[ihis]=(TH2F*)(chains[ihis]->GetHistogram()->Clone());
2068 hisCorrRef[ihis]->SetName(Form("DeltaRPhiCorr%d",(ihis+1)*2000));
2069 hisCorrRef[ihis]->SetTitle(Form("Corrected #Delta r#phi - Pileup %d",(ihis+1)*2000));
2070 hisCorrRef[ihis]->GetXaxis()->SetTitle("1/p_{t} (1/GeV/c)");
2071 hisCorrRef[ihis]->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
2072 hisCorrRef[ihis]->Write();
2073 //
2074 chains[ihis]->Draw("(T_DistRef_0.fP[0]/meanNorm):R.OT_DistRef_0.fP[4]>>hisCorNo(10,-3,3,100,-4,4)",cutOut+cutOutF+"","colzgoff");
2075 hisCorrNo[ihis]=(TH2F*)(chains[ihis]->GetHistogram()->Clone());
2076 hisCorrNo[ihis]->SetName(Form("DeltaRPhi%d",(ihis+1)*2000));
2077 hisCorrNo[ihis]->SetTitle(Form("Delta r#phi = Pileup %d",(ihis+1)*2000));
2078 hisCorrNo[ihis]->GetXaxis()->SetTitle("1/p_{t} (1/GeV/c)");
2079 hisCorrNo[ihis]->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
2080 hisCorrNo[ihis]->Write();
2081 }
2082 }
2083 ftrackFluctuation->Flush();
2084 //
2085 //
2086 //
2087 for (Int_t ihis=0; ihis<5; ihis++){
2088 hisCorrRef[ihis]->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2089 hisCorrRefM[ihis] = (TH1*)arrayFit.At(1)->Clone();
2090 hisCorrRefRMS[ihis] = (TH1*)arrayFit.At(2)->Clone();
2091 hisCorrRefM[ihis]->GetXaxis()->SetTitle("1/p_{t} (1/GeV/c)");
2092 hisCorrRefM[ihis]->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
2093 hisCorrRefM[ihis]->SetMarkerStyle(20);
2094 hisCorrRefRMS[ihis]->SetMarkerStyle(21);
2095 hisCorrRefM[ihis]->SetMarkerColor(1);
2096 hisCorrRefRMS[ihis]->SetMarkerColor(2);
2097 hisCorrNo[ihis]->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2098 hisCorrNoM[ihis] = (TH1*)arrayFit.At(1)->Clone();
2099 hisCorrNoRMS[ihis] = (TH1*)arrayFit.At(2)->Clone();
2100 }
4e715689 2101
4b955de3 2102 //
2103 TCanvas *canvasMean = new TCanvas("canvasCorrectionMean","canvasCorrectionMean",900,1000);
2104 TCanvas *canvasMeanSummary = new TCanvas("canvasCorrectionMeanSummary","canvasCorrectionMeanSummary",700,600);
2105
2106 canvasMean->Divide(3,5);
2107 gStyle->SetOptStat(0);
2108 for (Int_t ihis=0; ihis<5; ihis++){
2109 TLegend * legend = new TLegend(0.11,0.11,0.5,0.3,Form("Pile up %d",(ihis+1)*2000));
2110 canvasMean->cd(3*ihis+1);
2111 hisCorrNo[ihis]->Draw("colz");
2112 canvasMean->cd(3*ihis+2);
2113 hisCorrRef[ihis]->Draw("colz");
2114 canvasMean->cd(3*ihis+3);
2115 hisCorrRefM[ihis]->SetMaximum(0.25);
2116 hisCorrRefM[ihis]->SetMinimum(-0.25);
2117 hisCorrRefM[ihis]->Draw("");
2118 hisCorrRefRMS[ihis]->Draw("same");
2119 legend->AddEntry(hisCorrRefM[ihis],"Mean");
2120 legend->AddEntry(hisCorrRefRMS[ihis],"RMS");
2121 legend->Draw();
2122 }
2123 canvasMeanSummary->cd();
2124 TLegend * legendMeanSummary = new TLegend(0.5,0.6,0.89,0.89,"Space charge correction fluctuation in r#phi");
2125 for (Int_t ihis=4; ihis>=0; ihis--){
2126 hisCorrRefRMS[ihis]->SetMarkerColor(1+ihis);
2127 hisCorrRefRMS[ihis]->SetMinimum(0);
2128 hisCorrRefRMS[ihis]->GetYaxis()->SetTitle("#sigma_{r#phi} (cm)");
2129 if (ihis==4) hisCorrRefRMS[ihis]->Draw("");
2130 hisCorrRefRMS[ihis]->Draw("same");
2131 legendMeanSummary->AddEntry(hisCorrRefRMS[ihis],Form("%d pile-up events",(ihis+1)*2000));
2132 }
2133 legendMeanSummary->Draw();
4e715689 2134
4b955de3 2135 canvasMean->SaveAs("canvasCorrectionMean.pdf");
2136 canvasMeanSummary->SaveAs("canvasCorrectionMeanSummary.pdf");
2137 //canvasMean->Write();
2138 //canvasMeanSummary->Write();
2139 ftrackFluctuation->Close();
7212f8c3 2140}
4b955de3 2141
2142void DrawTrackFluctuationZ(){
2143 //
2144 // Draw track fucutation dz
2145 //
2146 const Int_t kColors[6]={1,2,3,4,6,7};
2147 const Int_t kStyle[6]={20,21,24,25,24,25};
2148 TObjArray arrayFit(3);
2149 TCut cutOut="abs(T_DistRef_0.fX-OT_DistRef_0.fX)<0.1&&T_DistRef_0.fX>1&&abs(OT_DistRef_0.fP[4])<4";
2150 TCut cutOutF="abs(R.T_DistRef_0.fX-R.OT_DistRef_0.fX)<0.1&&R.T_DistRef_0.fX>1&&abs(R.OT_DistRef_0.fP[4])<4";
2151 TChain * chains[5]={0};
2152 TChain * chainR = AliXRDPROOFtoolkit::MakeChain("track0_1.list","trackFit",0,1000);
2153 chainR->SetCacheSize(1000000000);
2154 for (Int_t ichain=0; ichain<5; ichain++){
2155 chains[ichain] = AliXRDPROOFtoolkit::MakeChain(Form("track%d_1.list",2*(ichain+1)),"trackFit",0,1000);
2156 chains[ichain]->AddFriend(chainR,"R");
2157 chains[ichain]->SetCacheSize(1000000000);
2158 chains[ichain]->SetMarkerStyle(25);
2159 chains[ichain]->SetMarkerSize(0.5);
2160 }
2161 //
2162 // 2.) Create 2D histo or read from files
2163 //
2164 TObjArray * arrayHisto = new TObjArray(25);
2165 TFile *ftrackFluctuationZ = TFile::Open("trackFluctuationZ.root","update");
2166 for (Int_t ihis=0; ihis<5; ihis++){
2167 ftrackFluctuationZ->cd();
2168 for (Int_t idz=0; idz<5; idz++){
2169 Int_t z= 16+idz*32;
2170 TH2 *his= (TH2*)ftrackFluctuationZ->Get(Form("TrackDz%d_PileUp%d",z, (ihis+1)*2000));
2171 if (!his){
2172 chains[ihis]->Draw(Form("T_DistZ_%d_0.fP[0]-T_DistRef_0.fP[0]:T_DistRef_0.fP[4]>>his(10,-4,4,100,-0.25,0.25)",z),cutOut+"","colz");
2173 his = (TH2*)(chains[ihis]->GetHistogram()->Clone());
2174 his->SetName(Form("TrackDz%d_PileUp%d",z, (ihis+1)*2000));
2175 his->Write();
2176 }
2177 arrayHisto->AddAtAndExpand(his,ihis*5+idz);
2178 }
2179 }
2180 ftrackFluctuationZ->Flush();
2181
2182 //
2183 // 3.) Make fits
2184 //
2185 TCanvas *canvasDz = new TCanvas("canvasDz","canvasDz",800,800);
2186 canvasDz->Divide(2,2,0,0);
2187 for (Int_t ihis=3; ihis>=0; ihis--){
2188 canvasDz->cd(ihis+1)->SetTicks(3);
2189 TLegend * legend = new TLegend(0.31,0.51, 0.95,0.95,Form("Distortion due time/z delay (Pileup=%d)", (ihis+1)*2000));
2190 legend->SetBorderSize(0);
2191 for (Int_t idz=3; idz>=0; idz--){
2192 TH2 * his = (TH2*)arrayHisto->At(ihis*5+idz);
2193 his->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2194 TH1 * hisRMS = (TH1*)arrayFit.At(2)->Clone();
2195 hisRMS->SetMaximum(0.12);
2196 hisRMS->SetMinimum(0);
2197 hisRMS->GetXaxis()->SetTitle("1/p_{t} (GeV/c)");
2198 hisRMS->GetYaxis()->SetTitle("#sigma_{r#phi}(cm)");
2199 hisRMS->SetMarkerStyle(kStyle[idz]);
2200 hisRMS->SetMarkerColor(kColors[idz]);
2201 if (idz==3) hisRMS->Draw();
2202 legend->AddEntry(hisRMS,Form("#Delta_{z}=%d (cm)",16+idz*32));
2203 hisRMS->Draw("same");
2204 }
2205 legend->Draw();
2206 }
2207 canvasDz->SaveAs("spaceChargeDeltaZScan.pdf");
2208
2209}
2210
2211
2ff913b6 2212
2213
2214
2215void DrawTrackFluctuationFrame(){
2216 //
2217 // Function to make a fluctuation figures for differnt multiplicities of pileup space charge
2218 // it is assumed that the text files
2219 //
2220 //
2221 TObjArray arrayFit(3);
2222 const char *inputList;
2223 TH2F * hisCorrRef[10]={0};
2224 TH2F * hisCorrNo[10]={0};
2225 TH1 * hisCorrRefM[10], *hisCorrRefRMS[10];
2226 TH1 * hisCorrNoM[10], *hisCorrNoRMS[10];
2227 //
2228 // 1. Load chains for different statistic
2229 //
2230 TCut cutOut="abs(T_DistRef_0.fX-OT_DistRef_0.fX)<0.1&&T_DistRef_0.fX>1&&abs(OT_DistRef_0.fP[4])<4";
2231 TCut cutOutF="abs(R.T_DistRef_0.fX-R.OT_DistRef_0.fX)<0.1&&R.T_DistRef_0.fX>1&&abs(R.OT_DistRef_0.fP[4])<4";
2232 TCut cutFit="Entry$%4==0"; //use only subset of data for fit
2233
2234 TChain * chains[10]={0};
2235 TChain * chainR = AliXRDPROOFtoolkit::MakeChain("track0_1.list","trackFit",0,1000);
2236 chainR->SetCacheSize(1000000000);
2237 for (Int_t ichain=0; ichain<7; ichain++){
2238 chains[ichain] = AliXRDPROOFtoolkit::MakeChain(Form("track%d_1.list",2*(ichain+1)),"trackFit",0,1000);
2239 chains[ichain]->AddFriend(chainR,"R");
2240 chains[ichain]->SetCacheSize(1000000000);
2241 chains[ichain]->SetMarkerStyle(25);
2242 chains[ichain]->SetMarkerSize(0.5);
2243 chains[ichain]->SetAlias("meanNorm","(1+0.2*abs(neventsCorr/10000-1))"); // second order correction - renomalization of mean hardwired
2244 chains[ichain]->SetAlias("dMean0","(neventsCorr*R.T_DistRef_0.fP[0]/10000)");
2245 chains[ichain]->SetAlias("dMeas0","T_DistRef_0.fP[0]");
2246 chains[ichain]->SetAlias("dMean1","(neventsCorr*R.T_DistRef_1.fP[0]/10000)");
2247 chains[ichain]->SetAlias("dMeas1","T_DistRef_1.fP[0]");
2248 for (Int_t ig=0; ig<10;ig++) chains[ichain]->SetAlias(Form("FR%d",ig),Form("(abs(Entry$-%d)<1000)",ig*2000+1000));
2249 }
2250 //
2251 // 2. Get or Create histogram (do fit per frame)
2252 //
2253 TStatToolkit toolkit;
2254 Double_t chi2=0;
2255 Int_t npoints=0;
2256 TVectorD param;
2257 TMatrixD covar;
2258 TString fstringG=""; // global part
2259 fstringG+="dMean0++";
2260 TVectorD vec0,vec1;
2261 TString fstringF0=""; // global part
2262 for (Int_t ig=0; ig<10;ig++) fstringF0+=Form("FR%d++",ig);
2263 for (Int_t ig=0; ig<10;ig++) fstringF0+=Form("FR%d*dMean0++",ig);
2264 TString fstringF1=""; // global part
2265 for (Int_t ig=0; ig<10;ig++) fstringF1+=Form("FR%d++",ig);
2266 for (Int_t ig=0; ig<10;ig++) fstringF1+=Form("FR%d*dMean0++",ig);
2267 for (Int_t ig=0; ig<10;ig++) fstringF1+=Form("FR%d*dMean0*abs(T_DistRef_0.fP[3])++",ig);
2268 for (Int_t ig=0; ig<10;ig++) fstringF1+=Form("FR%d*dMean0*(T_DistRef_0.fP[3]^2)++",ig);
2269
2270 //
2271 //
2272 TH2F *hisA=0, *hisF0=0, *hisF1=0, *hisM=0;
2273 TObjArray * arrayHisto = new TObjArray(200);
2274 TFile *ftrackFit = TFile::Open("trackFluctuationFrame.root","update");
2275 for (Int_t ihis=0; ihis<7; ihis++){
2276 printf("\n\nProcessing frames\t%d\nnn",(ihis+1)*2000);
2277 hisM = (TH2F*)ftrackFit->Get(Form("hisMean_%d",(ihis+1)*2000));
2278 hisA = (TH2F*)ftrackFit->Get(Form("hisAll_%d",(ihis+1)*2000));
2279 hisF0 = (TH2F*)ftrackFit->Get(Form("hisFrame0_%d",(ihis+1)*2000));
2280 hisF1 = (TH2F*)ftrackFit->Get(Form("hisFrame1_%d",(ihis+1)*2000));
2281 if (!hisA){
2282 ftrackFit->cd();
2283 TString * fitResultAll = TStatToolkit::FitPlane(chains[ihis],"dMeas0", fstringG.Data(),cutOut+cutOutF+cutFit, chi2,npoints,param,covar,-1,0, 40*2000, kFALSE);
2284 chains[ihis]->SetAlias("fitAll",fitResultAll->Data());
2285 TString * fitResultF0 = TStatToolkit::FitPlane(chains[ihis],"dMeas0", fstringF0.Data(),cutOut+cutOutF+cutFit+"abs(dMeas0-fitAll)<0.3", chi2,npoints,vec0,covar,-1,0, 10*2000, kFALSE);
2286 chains[ihis]->SetAlias("fitF0",fitResultF0->Data());
2287 TString * fitResultF1 = TStatToolkit::FitPlane(chains[ihis],"dMeas0", fstringF1.Data(),cutOut+cutOutF+cutFit+"abs(dMeas0-fitAll)<0.3", chi2,npoints,vec1,covar,-1,0, 10*2000, kFALSE);
2288 chains[ihis]->SetAlias("fitF1",fitResultF1->Data());
2289 fitResultF0->Tokenize("++")->Print();
2290 chains[ihis]->Draw(Form("dMeas0-fitAll:T_DistRef_0.fP[4]>>hisAll_%d(20,-4,4,100,-0.25,0.25)",(ihis+1)*2000),cutOut+cutOutF,"colz",100000,0);
2291 hisA = (TH2F*)chains[ihis]->GetHistogram();
2292 chains[ihis]->Draw(Form("dMeas0-fitF0:T_DistRef_0.fP[4]>>hisFrame0_%d(20,-4,4,100,-0.10,0.10)",(ihis+1)*2000),cutOut+cutOutF,"colz",20000,0);
2293 hisF0 = (TH2F*)chains[ihis]->GetHistogram();
2294 chains[ihis]->Draw(Form("dMeas0-fitF1:T_DistRef_0.fP[4]>>hisFrame1_%d(20,-4,4,100,-0.10,0.10)",(ihis+1)*2000),cutOut+cutOutF,"colz",20000,0);
2295 hisF1 = (TH2F*)chains[ihis]->GetHistogram();
2296 chains[ihis]->Draw(Form("dMeas0-dMean0:T_DistRef_0.fP[4]>>hisMean_%d(20,-4,4,100,-0.25,0.25)",(ihis+1)*2000),cutOut+cutOutF,"colz",100000,0);
2297 hisM = (TH2F*)chains[ihis]->GetHistogram();
2298 hisM->Write(); hisA->Write();hisF0->Write(); hisF1->Write();
2299 ftrackFit->Flush();
2300 }
2301 }
5e4083cf 2302
2ff913b6 2303 for (Int_t ihis=0; ihis<7; ihis++){
2304 printf("\n\nProcessing frames\t%d\nnn",(ihis+1)*2000);
2305 hisM = (TH2F*)ftrackFit->Get(Form("hisMean_%d",(ihis+1)*2000));
2306 hisA = (TH2F*)ftrackFit->Get(Form("hisAll_%d",(ihis+1)*2000));
2307 hisF0 = (TH2F*)ftrackFit->Get(Form("hisFrame0_%d",(ihis+1)*2000));
2308 hisF1 = (TH2F*)ftrackFit->Get(Form("hisFrame1_%d",(ihis+1)*2000));
2309 arrayHisto->AddLast(hisA);
2310 arrayHisto->AddLast(hisF0);
2311 arrayHisto->AddLast(hisF1);
2312 arrayHisto->AddLast(hisM);
2313 }
2314 delete ftrackFit;
2315 //
2316 // 3. Draw figures
2317 //
2318 gStyle->SetOptStat(0);
2319 TCanvas *canvasFit = new TCanvas("canvasFitFrame","canvasFitframe",900,700);
2320 canvasFit->Divide(3,2,0,0);
2321 for (Int_t ihis=1; ihis<7; ihis++){
2322 //
2323 canvasFit->cd(ihis);
2324 char hname[10000];
2325 snprintf(hname,1000,"hisAll_%d",(ihis+1)*2000);
2326 hisA = (TH2F*)arrayHisto->FindObject(hname);
2327 snprintf(hname,1000,"hisFrame0_%d",(ihis+1)*2000);
2328 hisF0 = (TH2F*)arrayHisto->FindObject(hname);
2329 snprintf(hname,1000,"hisFrame1_%d",(ihis+1)*2000);
2330 hisF1 = (TH2F*)arrayHisto->FindObject(hname);
2331 snprintf(hname,1000,"hisMean_%d",(ihis+1)*2000);
2332 hisM = (TH2F*)arrayHisto->FindObject(hname);
2333 //
2334 //
2335 hisM->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2336 TH1 * hisRA= (TH1*)arrayFit.At(2)->Clone();
2337 hisF0->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2338 TH1 * hisRF0= (TH1*)arrayFit.At(2)->Clone();
2339 hisF1->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2340 TH1 * hisRF1= (TH1*)arrayFit.At(2)->Clone();
2341 //
2342 hisRA->SetMarkerStyle(20);
2343 hisRF0->SetMarkerStyle(21);
2344 hisRF1->SetMarkerStyle(21);
2345 hisRA->SetMarkerColor(1);
2346 hisRF0->SetMarkerColor(4);
2347 hisRF1->SetMarkerColor(2);
2348 TF1 * f1a= new TF1("f1a","pol1");
2349 TF1 * f1f0= new TF1("f1a0","pol1");
2350 TF1 * f1f1= new TF1("f1a1","pol1");
2351 f1a->SetLineColor(1);
2352 f1f0->SetLineColor(4);
2353 f1f1->SetLineColor(2);
2354 hisRA->Fit(f1a);
2355 hisRF0->Fit(f1f0);
2356 hisRF1->Fit(f1f1);
2357 hisRF1->SetMinimum(0);
2358 hisRF1->SetMaximum(0.05);
2359 // hisRA->Draw();
2360 hisRF1->GetXaxis()->SetTitle("q/p_{T} (1/GeV)");
2361 hisRF1->GetYaxis()->SetTitle("#sigma_{r#phi} (cm)");
2362 hisRF1->Draw("");
2363 hisRF0->Draw("same");
2364 TLegend * legend = new TLegend(0.11,0.11,0.65,0.25, Form("Track residual r#phi distortion: N_{ion}=%d",(ihis+1)*2000));
2365 legend->AddEntry(hisRF0,"a_{0}+a_{1}#rho");
2366 legend->AddEntry(hisRF1,"a_{0}+(a_{1}+a_{2}z+a_{3}z^2)#rho");
2367 legend->SetBorderSize(0);
2368 legend->Draw();
2369 }
2370 //
2371 canvasFit->SaveAs("canvasFrameFitRPhiVersion0.pdf");
2372 canvasFit->SaveAs("canvasFrameFitRPhiVersion0.png");
2373 //
2374}
5a998ace 2375
2376
2377
2378void MakeLocalDistortionPlotsGlobalFitPolDrift(Float_t xmin, Float_t xmax){
2379 //
2380 // Make local distortion plots correcting using global z fits of order 0,2,4,6,8,10,12
2381 // Currently polynomial correction as an multiplicative factor of the mean distortion map used
2382 // To be done - calculate numerical derivative of distortion maps
2383 // corresponding the local change of densities - after TDR?
2384 //
2385 // As input:
2386 // 1.) distortion file with the setup 10000 pileup events used
2387 // 2.) mean distortion file
2388 // distortions are fitted rescaling (z dependent) mean distortion file
2389 //
2390 TTreeSRedirector *pcstream = new TTreeSRedirector("localFit.root","update");
2391 TFile *fRef = TFile::Open("SpaceChargeFluc0_1.root");
2392 TFile *fCurrent = TFile::Open("SpaceChargeFluc10_1.root");
2393 TTree * treeRef = (TTree*)fRef->Get("hisDump");
2394 TTree * treeCurrent = (TTree*)fCurrent->Get("hisDump");
2395 treeCurrent->AddFriend(treeRef,"R");
2396 treeCurrent->SetAlias("refR","(neventsCorr*R.DistRefDR/10000)");
2397 treeCurrent->SetAlias("refRPhi","(neventsCorr*R.DistRefDRPhi/10000)");
2398 treeCurrent->SetCacheSize(1000000000);
2399 treeCurrent->SetAlias("drift","1.-abs(z/250.)");
2400 TStatToolkit toolkit;
2401 Double_t chi2=0;
2402 Int_t npoints=0;
2403 TVectorD param;
2404 TMatrixD covar;
2405 //
2406 TCut cut="z>0";
2407 TString fstringG0="refR++"; // global part - for z dependence we should use the Chebischev
2408 TString fstringG1="refR++"; // global part - for z dependence we should use the Chebischev
2409 TString fstringG2="refR++"; // global part - for z dependence we should use the Chebischev
2410 TString fstringG3="refR++"; // global part - for z dependence we should use the Chebischev
2411 TString fstringG4="refR++"; // global part - for z dependence we should use the Chebischev
2412 TString fstringG5="refR++"; // global part - for z dependence we should use the Chebischev
2413 TString fstringG6="refR++"; // global part - for z dependence we should use the Chebischev
2414 TString fstringG7="refR++"; // global part - for z dependence we should use the Chebischev
2415 TString fstringG8="refR++"; // global part - for z dependence we should use the Chebischev
2416 //
2417 for (Int_t i=1; i<=1; i++) fstringG1+=TString::Format("refR*pow(drift,%d)++",i);
2418 for (Int_t i=1; i<=2; i++) fstringG2+=TString::Format("refR*pow(drift,%d)++",i);
2419 for (Int_t i=1; i<=3; i++) fstringG3+=TString::Format("refR*pow(drift,%d)++",i);
2420 for (Int_t i=1; i<=4; i++) fstringG4+=TString::Format("refR*pow(drift,%d)++",i);
2421 for (Int_t i=1; i<=5; i++) fstringG5+=TString::Format("refR*pow(drift,%d)++",i);
2422 for (Int_t i=1; i<=6; i++) fstringG6+=TString::Format("refR*pow(drift,%d)++",i);
2423 for (Int_t i=1; i<=7; i++) fstringG7+=TString::Format("refR*pow(drift,%d)++",i);
2424 for (Int_t i=1; i<=8; i++) fstringG8+=TString::Format("refR*pow(drift,%d)++",i);
2425
2426
2427 TString * fitResultG0 = TStatToolkit::FitPlane(treeCurrent,"DistRefDR-refR", fstringG0.Data(),cut+"Entry$%7==0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
2428 TString * fitResultG1 = TStatToolkit::FitPlane(treeCurrent,"DistRefDR-refR", fstringG1.Data(),cut+"Entry$%7==0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
2429 TString * fitResultG2 = TStatToolkit::FitPlane(treeCurrent,"DistRefDR-refR", fstringG2.Data(),cut+"Entry$%7==0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
2430 TString * fitResultG3 = TStatToolkit::FitPlane(treeCurrent,"DistRefDR-refR", fstringG3.Data(),cut+"Entry$%7==0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
2431 TString * fitResultG4 = TStatToolkit::FitPlane(treeCurrent,"DistRefDR-refR", fstringG4.Data(),cut+"Entry$%7==0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
2432 TString * fitResultG5 = TStatToolkit::FitPlane(treeCurrent,"DistRefDR-refR", fstringG5.Data(),cut+"Entry$%7==0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
2433 TString * fitResultG6 = TStatToolkit::FitPlane(treeCurrent,"DistRefDR-refR", fstringG6.Data(),cut+"Entry$%7==0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
2434 TString * fitResultG7 = TStatToolkit::FitPlane(treeCurrent,"DistRefDR-refR", fstringG7.Data(),cut+"Entry$%7==0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
2435 TString * fitResultG8 = TStatToolkit::FitPlane(treeCurrent,"DistRefDR-refR", fstringG8.Data(),cut+"Entry$%7==0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
2436 //
2437 treeCurrent->SetAlias("fitG0",fitResultG0->Data());
2438 treeCurrent->SetAlias("fitG1",fitResultG1->Data());
2439 treeCurrent->SetAlias("fitG2",fitResultG2->Data());
2440 treeCurrent->SetAlias("fitG3",fitResultG3->Data());
2441 treeCurrent->SetAlias("fitG4",fitResultG4->Data());
2442 treeCurrent->SetAlias("fitG5",fitResultG5->Data());
2443 treeCurrent->SetAlias("fitG6",fitResultG6->Data());
2444 treeCurrent->SetAlias("fitG7",fitResultG7->Data());
2445 treeCurrent->SetAlias("fitG8",fitResultG8->Data());
2446 treeCurrent->SetMarkerStyle(25);
2447 treeCurrent->SetMarkerSize(0.2);
2448 //
2449 gStyle->SetOptFit(1);
2450 gStyle->SetOptStat(0);
2451 gStyle->SetOptTitle(1);
2452
2453 TCut cutDrawGraph="z>0&&abs(z-30)<5&&abs(r-90)<10";
2454 TCut cutDrawHisto="z>0&&abs(z)<125&&abs(r-90)<10";
2455 TH1F *hisFull=new TH1F("hisFull","hisFull",100,xmin,xmax);
2456 TH1F *hisG0=new TH1F("hisG0","hisG0",100,xmin,xmax);
2457 TH1F *hisG1=new TH1F("hisG1","hisG1",100,xmin,xmax);
2458 TH1F *hisG2=new TH1F("hisG2","hisG2",100,xmin,xmax);
2459 TH1F *hisG3=new TH1F("hisG3","hisG3",100,xmin,xmax);
2460 TH1F *hisG4=new TH1F("hisG4","hisG4",100,xmin,xmax);
2461 TH1F *hisG5=new TH1F("hisG5","hisG5",100,xmin,xmax);
2462 TH1F *hisG6=new TH1F("hisG6","hisG6",100,xmin,xmax);
2463 TH1F *hisG7=new TH1F("hisG7","hisG7",100,xmin,xmax);
2464 TH1F *hisG8=new TH1F("hisG8","hisG8",100,xmin,xmax);
2465 TH1F * hisResG[10]={hisFull, hisG0, hisG1,hisG2, hisG3,hisG4, hisG5,hisG6, hisG7,hisG8};
2466 treeCurrent->Draw("DistRefDR-refR>>hisFull",cutDrawHisto,"");
2467 for (Int_t ihis=0; ihis<9; ihis++){
2468 treeCurrent->Draw(TString::Format("DistRefDR-refR-fitG%d>>hisG%d",ihis,ihis),cutDrawHisto,"");
2469 }
2470 //
2471 TF1 *fg = new TF1("fg","gaus");
2472 TVectorD vecP(10), vecRes(10), vecResE(10);
2473 for (Int_t ihis=0; ihis<10; ihis++){
2474 hisResG[ihis]->Fit(fg);
2475 vecP[ihis]=ihis-1;
2476 vecRes[ihis]=fg->GetParameter(2);
2477 vecResE[ihis]=fg->GetParError(2);
2478 hisResG[ihis]->GetXaxis()->SetTitle("#Delta_{R} (cm)");
2479 pcstream->GetFile()->cd();
2480 hisResG[ihis]->Write();
2481 (*pcstream)<<"residuals"<<
2482 TString::Format("diffHis%d.=",ihis)<<hisResG[ihis];
2483 }
2484 TGraphErrors *gr = new TGraphErrors(10,vecP.GetMatrixArray(), vecRes.GetMatrixArray(),0, vecResE.GetMatrixArray());
2485 gr->SetMarkerStyle(25);
2486 gr->Draw("alp");
2487 (*pcstream)<<"residuals"<<
2488 "graph.="<<gr<<
2489 "\n";
2490
2491 TCanvas *canvasRes = new TCanvas("canvasRes","canvasRes",900,900);
2492 canvasRes->Divide(2,4);
2493 TH2* htemp=0;
2494 for (Int_t ihis=0; ihis<4; ihis++){
2495 canvasRes->cd(ihis*2+1);
2496 if (ihis==0) {
2497 treeCurrent->Draw("DistRefDR-refR:phi:r",cutDrawGraph,"colz");
2498 }
2499 if (ihis>0) {
2500 treeCurrent->Draw(TString::Format("DistRefDR-refR-fitG%d:phi:r",(ihis-1)*2),cutDrawGraph,"colz");
2501 }
2502 htemp = (TH2F*)gPad->GetPrimitive("htemp");
2503 htemp->GetXaxis()->SetTitle("#phi");
2504 htemp->GetYaxis()->SetTitle("#Delta_{R} (cm)");
2505 htemp->GetZaxis()->SetTitle("r (cm)");
2506 // htemp->SetTitle("1/pT difference");
2507 //htemp->GetYaxis()->SetRangeUser(xmin,xmax);
2508 canvasRes->cd(ihis*2+2);
2509 if (ihis>0) hisResG[(ihis-1)*2+1]->Draw();
2510 if (ihis==0) hisResG[0]->Draw();
2511 }
2512 delete pcstream;
2513
2514 canvasRes->SaveAs("locaFluctuationR.pdf");
2515 canvasRes->SaveAs("locaFluctuationR.png");
2516}
2517
2518void MakeLocalDistortionPlotsGlobalFitPolDriftSummary(Float_t xmin, Float_t xmax){
2519 //
2520 // Make local distortion plots correcting using global z fits of order 0,2,4,6,8,10,12
2521 // Currently polynomial correction as an multiplicative factor of the mean distortion map used
2522 // To be done - calculate numerical derivative of distortion maps
2523 // corresponding the local change of densities - after TDR?
2524 //
2525 // As input:
2526 // 1.) distortion file with the setup 10000 pileup events used
2527 // 2.) mean distortion file
2528 // distortions are fitted rescaling (z dependent) mean distortion file
2529 //
2530 gStyle->SetOptStat(0);
2531 gStyle->SetOptFit(1);
2532 gStyle->SetOptTitle(1);
2533 TTreeSRedirector *pcstream = new TTreeSRedirector("localFit.root","update");
2534 TH1 *hisResG[10] = {0};
2535 hisResG[0]=(TH1*)(pcstream->GetFile()->Get("hisFull"));
2536 TF1 *fg = new TF1("fg","gaus");
2537 TVectorD vecP(10), vecRes(10), vecResE(10);
2538 for (Int_t ihis=0; ihis<10; ihis++){
2539 hisResG[ihis+1]=(TH1*)(pcstream->GetFile()->Get(TString::Format("hisG%d",ihis)));
2540 hisResG[ihis]->Fit(fg);
2541 vecP[ihis]=ihis-1;
2542 vecRes[ihis]=fg->GetParameter(2);
2543 vecResE[ihis]=fg->GetParError(2);
2544 }
2545 //
2546 TCanvas *canvasRes = new TCanvas("canvasRes","canvasRes",800,800);
2547 canvasRes->Divide(2,3,0,0);
2548 for (Int_t ihis=0; ihis<6; ihis++){
2549 canvasRes->cd(ihis+1);
2550 hisResG[ihis]->GetXaxis()->SetTitle("#Delta_{R} (cm)");
2551 hisResG[ihis]->Draw();
2552 }
2553 canvasRes->SaveAs("fluctuationTableSummaryHist.pdf");
2554 canvasRes->SaveAs("fluctuationTableSummaryHist.png");
2555
2556 TCanvas *canvasFluctuationGraph = new TCanvas("canvasGraph","canvasGraph",600,500);
2557 TGraphErrors *gr = new TGraphErrors(10,vecP.GetMatrixArray(), vecRes.GetMatrixArray(),0, vecResE.GetMatrixArray());
2558 gr->SetMarkerStyle(25);
2559 gr->SetMinimum(0);
2560 gr->GetXaxis()->SetTitle("#Fit parameters");
2561 gr->GetYaxis()->SetTitle("#sigma_{R} (cm)");
2562 gr->Draw("alp");
2563 //
2564 canvasFluctuationGraph->SaveAs("canvasFluctuationGraphR.pdf");
2565 canvasFluctuationGraph->SaveAs("canvasFluctuationGraphR.png");
2566
2567}
2568
2569
2570
2571void MakeLocalDistortionPlots(Int_t npoints, Int_t npointsZ){
2572 //
2573 //
2574 TTreeSRedirector *pcstream = new TTreeSRedirector("localBins.root","update");
2575 TFile *fCurrent = TFile::Open("SpaceChargeFluc10_1.root");
2576 TFile *fRef = TFile::Open("SpaceChargeFluc0_1.root");
2577 //
2578 AliTPCSpaceCharge3D* distortion = ( AliTPCSpaceCharge3D*)fCurrent->Get("DistRef");
2579 AliTPCSpaceCharge3D* distortionRef = ( AliTPCSpaceCharge3D*)fRef->Get("DistRef");
2580 distortion->AddVisualCorrection(distortion,1);
2581 distortionRef->AddVisualCorrection(distortionRef,2);
2582 //
2583 //
2584 //
2585 TVectorD normZR(125), normZRPhi(125), normZZ(125), normZPos(125);
2586 TVectorD normZRChi2(125), normZRPhiChi2(125), normZZChi2(125);
2587 //
2588 for (Int_t iz =0; iz<125; iz++){
2589 Double_t z0 = -250+iz*4;
2590 TLinearFitter fitterR(2,"pol1");
2591 TLinearFitter fitterRPhi(2,"pol1");
2592 TLinearFitter fitterZ(2,"pol1");
2593 Double_t xvalue[10]={0};
2594 //fitterR.AddPoint(xvalue,0,0.001/npointsZ);
2595 //fitterRPhi.AddPoint(xvalue,0,0.000001/npointsZ);
2596 //fitterZ.AddPoint(xvalue,0,0.001/npointsZ);
2597 for (Int_t ipoint =0; ipoint<npointsZ; ipoint++){
2598 Double_t r0 = 85+gRandom->Rndm()*(245-85.);
2599 Double_t phi0 = gRandom->Rndm()*TMath::TwoPi();
2600 // some problematic parts to be skipped - investigated later
2601 if (TMath::Abs(distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,1))>50) continue;
2602 if (TMath::Abs(distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,2))>50) continue;
2603 if (TMath::Abs(distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1))>20) continue;
2604 if (TMath::Abs(distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,2))>20) continue;
2605 if (TMath::Abs(distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,1))>50) continue;
2606 if (TMath::Abs(distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,2))>50) continue;
2607 //
2608 //
2609 xvalue[0]=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,2);
2610 fitterR.AddPoint(xvalue,distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,1));
2611 xvalue[0]=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,2);
2612 fitterRPhi.AddPoint(xvalue,distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1));
2613 xvalue[0]=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,2);
2614 fitterZ.AddPoint(xvalue,distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,1));
2615 }
2616 fitterR.Eval();
2617 fitterRPhi.Eval();
2618 fitterZ.Eval();
2619 normZR[iz]=fitterR.GetParameter(1);
2620 normZRPhi[iz]=fitterRPhi.GetParameter(1);
2621 normZZ[iz]=fitterZ.GetParameter(1);
2622 normZRChi2[iz]=TMath::Sqrt(fitterR.GetChisquare()/fitterR.GetNpoints());
2623 normZRPhiChi2[iz]=TMath::Sqrt(fitterRPhi.GetChisquare()/fitterRPhi.GetNpoints());
2624 normZZChi2[iz]=TMath::Sqrt(fitterZ.GetChisquare()/fitterZ.GetNpoints());
2625
2626 normZPos[iz]=z0;
2627 }
2628
2629 {
2630 (*pcstream)<<"meanNormZ"<<
2631 "normZPos.="<<&normZPos<<
2632 //
2633 "normZR.="<<&normZR<< // mult. scaling to minimize R distortions
2634 "normZRPhi.="<<&normZRPhi<< // mult. scaling
2635 "normZZ.="<<&normZZ<<
2636 //
2637 "normZRChi2.="<<&normZRChi2<< // mult. scaling to minimize R distortions
2638 "normZRPhiChi2.="<<&normZRPhiChi2<< // mult. scaling
2639 "normZZChi2.="<<&normZZChi2<<
2640 "\n";
2641 }
2642 delete pcstream;
2643 pcstream = new TTreeSRedirector("localBins.root","update");
2644 //
2645 TTree * treeNormZ= (TTree*)pcstream->GetFile()->Get("meanNormZ");
2646 TGraphErrors * grZRfit= TStatToolkit::MakeGraphErrors( treeNormZ, "normZR.fElements:normZPos.fElements","",25,2,0.5);
2647 TGraphErrors * grZRPhifit= TStatToolkit::MakeGraphErrors( treeNormZ, "normZRPhi.fElements:normZPos.fElements","",25,4,0.5);
2648 grZRfit->Draw("alp");
2649 grZRPhifit->Draw("lp");
2650
2651 //
2652 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2653 //
2654 for (Int_t izNorm=0; izNorm<2; izNorm++){
2655 Double_t r0 = 85+gRandom->Rndm()*(245-85.);
2656 Double_t z0 = gRandom->Rndm()*250;
2657 Double_t phi0 = gRandom->Rndm()*TMath::TwoPi();
2658 Double_t fSector=18.*phi0/TMath::TwoPi();
2659 Double_t dSector=fSector-Int_t(fSector);
2660
2661
2662 Int_t iz0 = TMath::Nint(125.*(z0+250.)/500.);
2663 if (iz0<0) iz0=0;
2664 if (iz0>=125) iz0=124;
2665 Double_t rNorm=1,rphiNorm=1,zNorm=1;
2666 if (izNorm==1){
2667 rNorm=normZR[iz0];
2668 rphiNorm=normZRPhi[iz0];
2669 zNorm=normZZ[iz0];
2670 }
2671
2672 Double_t deltaR0 =distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,1);
2673 Double_t deltaRPhi0=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1);
2674 Double_t deltaZ0 =distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,1);
2675 //
2676 Double_t ddeltaR0 =distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,1)-rNorm*distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,2);
2677 Double_t ddeltaRPhi0=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1)-rphiNorm*distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,2);
2678 Double_t ddeltaZ0 =distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,1)-zNorm*distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,2);
2679 //
2680 //
2681 TVectorD vecDMeanRPhi(20), vecDMeanR(20), vecDMeanZ(20), vecDPhi(20);
2682 for (Int_t ideltaBin=0; ideltaBin<20; ideltaBin++){
2683 Double_t deltaPhi=ideltaBin*TMath::TwoPi()/360.;
2684 vecDPhi[ideltaBin]=deltaPhi;
2685 vecDMeanRPhi[ideltaBin]=0;
2686 vecDMeanR[ideltaBin]=0;
2687 vecDMeanZ[ideltaBin]=0;
2688
2689 for (Int_t i=-2; i<=2; i++){
2690 vecDMeanR[ideltaBin]+=0.2*(distortion->GetCorrXYZ(r0*TMath::Cos(phi0+deltaPhi*i), r0*TMath::Sin(phi0+deltaPhi*i), z0,0,1)-rNorm*distortion->GetCorrXYZ(r0*TMath::Cos(phi0+deltaPhi*i), r0*TMath::Sin(phi0+deltaPhi*i), z0,0,2));
2691 vecDMeanRPhi[ideltaBin]+=0.2*(distortion->GetCorrXYZ(r0*TMath::Cos(phi0+deltaPhi*i), r0*TMath::Sin(phi0+deltaPhi*i), z0,1,1)-rphiNorm*distortion->GetCorrXYZ(r0*TMath::Cos(phi0+deltaPhi*i), r0*TMath::Sin(phi0+deltaPhi*i), z0,1,2));
2692 vecDMeanZ[ideltaBin]+=0.2*(distortion->GetCorrXYZ(r0*TMath::Cos(phi0+deltaPhi*i), r0*TMath::Sin(phi0+deltaPhi*i), z0,2,1)-zNorm*distortion->GetCorrXYZ(r0*TMath::Cos(phi0+deltaPhi*i), r0*TMath::Sin(phi0+deltaPhi*i), z0,2,2));
2693 }
2694 }
2695 (*pcstream)<<"meanDistortion"<<
2696 "izNorm="<<izNorm<< // z normalizatio flag
2697 "fSector="<<fSector<< // sector position
2698 "dSector="<<dSector<< // distance to the sector boundary
2699 //
2700 "r0="<<r0<<
2701 "z0="<<z0<<
2702 "phi0="<<phi0<<
2703 //
2704 "rNorm="<<rNorm<<
2705 "rphiNorm="<<rphiNorm<<
2706 "zNorm="<<zNorm<<
2707 //
2708 "deltaR0="<<deltaR0<<
2709 "deltaZ0="<<deltaZ0<<
2710 "deltaRPhi0="<<deltaRPhi0<<
2711 //
2712 "ddeltaR0="<<ddeltaR0<<
2713 "ddeltaZ0="<<ddeltaZ0<<
2714 "ddeltaRPhi0="<<ddeltaRPhi0<<
2715 //
2716 "vecDMeanRPhi.="<<&vecDMeanRPhi<<
2717 "vecDMeanR.="<<&vecDMeanR<<
2718 "vecDMeanZ.="<<&vecDMeanZ<<
2719 "vecDPhi.="<<&vecDPhi<<
2720 "\n";
2721 }
2722 }
2723 delete pcstream;
2724 pcstream = new TTreeSRedirector("localBins.root","update");
2725 /*
2726 meanDistortion->Draw("vecDMeanR.fElements-ddeltaR0:vecDPhi.fElements","izNorm==1&&abs(phi0-pi)>0.2&&abs(ddeltaRPhi0)<10","")
2727
2728 */
2729}
2730
2731
2732void DrawLocalDistortionPlots(){
2733 //
2734 // Draw summary residuals plots after applying z dependent q normalization.
2735 // Plots used for the TPC TDR and internal note
2736 //
2737 // Two input trees are used:
2738 // meanNormZ - here we store the result of the q ze dependent fits for Radial, RPhi and Z distortions
2739 // meanDistortion - phi averaged residual distortion before and after applying q(z) dependent correction
2740 //
2741 // It is assumed that the correction table for randomly selected pile-up setup
2742 // can be approximated rescaling the mean correction table.
2743 // Rescaling coeficients are fitted separatelly in 125 z bins
2744 //
2745
2746 TTreeSRedirector *pcstream = new TTreeSRedirector("localBins.root","update");
2747 TTree * meanNormZ = (TTree*) pcstream->GetFile()->Get("meanNormZ");
2748 TTree * meanDistortion = (TTree*) pcstream->GetFile()->Get("meanDistortion");
2749 meanNormZ->SetMarkerStyle(25); meanNormZ->SetMarkerSize(0.5);
2750 Int_t colors[5]={1,2,3,4,6};
2751 Int_t markers[5]={21,20,23,24,25};
2752 TH2 * his2D;
2753 TObjArray arrFit(3);
2754 //
2755 // 1. Z dependence of the normalization is smooth function - about 10 bins to represent
2756 //
2757 TGraphErrors *graphZRnorm[100]= {0};
2758 TGraphErrors *graphZRRPhinorm[100]= {0};
2759 TCanvas * canvasRFit = new TCanvas("canvasRFit","canvasRFit",600,500);
2760 TLegend * legendRFit = new TLegend(0.12,0.12,0.88,0.4,"Q normalization fit. #DeltaR=c(z)#DeltaR_{ref}");
2761 legendRFit->SetBorderSize(0);
2762 for (Int_t igraph=0; igraph<5; igraph++){
2763 Int_t color = colors[igraph];
2764 Int_t marker = markers[igraph];
2765 Int_t event=gRandom->Rndm()*meanNormZ->GetEntries();
2766 graphZRnorm[igraph] = TStatToolkit::MakeGraphErrors( meanNormZ, "normZR.fElements:normZPos.fElements",TString::Format("Entry$==%d&&abs(normZPos.fElements)<220",event),marker,color,0.7);
2767 graphZRnorm[igraph]->SetTitle(TString::Format("Pile-up setup=%d",igraph));
2768 graphZRnorm[igraph]->SetMinimum(0.60);
2769 graphZRnorm[igraph]->SetMaximum(1.2);
2770 graphZRnorm[igraph]->GetXaxis()->SetTitle("z (cm)");
2771 graphZRnorm[igraph]->GetYaxis()->SetTitle("c(z)");
2772 if (igraph==0) graphZRnorm[igraph]->Draw("alp");
2773 graphZRnorm[igraph]->Draw("lp");
2774 legendRFit->AddEntry( graphZRnorm[igraph],TString::Format("Pile-up setup=%d",event),"p");
2775 }
2776 legendRFit->Draw();
2777 canvasRFit->SaveAs("canvasZRFit5Random.pdf"); // ~/hera/alice/miranov/SpaceCharge/Fluctuations/PbPbWithGain/dirmergeAll/dEpsilon20/canvasZRRPhiFit5Random.png
2778 canvasRFit->SaveAs("canvasZRFit5Random.png");
2779 //
2780 // 2.
2781 //
2782 TCanvas * canvasRRPhiFit = new TCanvas("canvasRRPhiFit","canvasRRPhiFit",600,500);
2783 TLegend * legendRRPhiFit = new TLegend(0.12,0.12,0.88,0.4,"Q normalization fit. #DeltaR=c_{R}(z)#DeltaR_{ref} #DeltaR#phi=c_{R#phi}(z)#Delta_{R#phi}_{ref}");
2784 legendRRPhiFit->SetBorderSize(0);
2785 for (Int_t igraph=0; igraph<5; igraph++){
2786 Int_t color = colors[igraph];
2787 Int_t marker = markers[igraph];
2788 Int_t event=gRandom->Rndm()*meanNormZ->GetEntries();
2789 graphZRRPhinorm[igraph] = TStatToolkit::MakeGraphErrors( meanNormZ, "normZR.fElements:normZRPhi.fElements",TString::Format("Entry$==%d&&abs(normZPos.fElements)<220",event),marker,color,0.7);
2790 graphZRRPhinorm[igraph]->GetXaxis()->SetLimits(0.6,1.2);
2791 graphZRRPhinorm[igraph]->SetTitle(TString::Format("Pile-up setup=%d",igraph));
2792 graphZRRPhinorm[igraph]->SetMinimum(0.6);
2793 graphZRRPhinorm[igraph]->SetMaximum(1.2);
2794 graphZRRPhinorm[igraph]->GetXaxis()->SetTitle("c_{R#phi}");
2795 graphZRRPhinorm[igraph]->GetYaxis()->SetTitle("c_{R}");
2796 if (igraph==0) graphZRRPhinorm[igraph]->Draw("alp");
2797 graphZRRPhinorm[igraph]->Draw("lp");
2798 legendRRPhiFit->AddEntry( graphZRRPhinorm[igraph],TString::Format("Pile-up setup=%d",event),"p");
2799 }
2800 legendRRPhiFit->Draw();
2801 canvasRRPhiFit->SaveAs("canvasZRRPhiFit5Random.pdf"); // ~/hera/alice/miranov/SpaceCharge/Fluctuations/PbPbWithGain/dirmergeAll/dEpsilon20/canvasZRRPhiFit5Random.png
2802 canvasRRPhiFit->SaveAs("canvasZRRPhiFit5Random.png");
2803 //
2804
2805 //
2806 // 3. Residual distortion after z dependent q distortion
2807 //
2808 gStyle->SetOptStat(0);
2809 TH1D * hisRRes, *hisRRphiRes=0;
2810 meanDistortion->Draw("ddeltaRPhi0:r0>>hisdRPhi0(40,85,245,100,-0.25,0.25)","abs(z0)<85.&&izNorm==1","colz");
2811 his2D = (TH2D*)meanDistortion->GetHistogram();
2812 his2D->FitSlicesY(0,0,-1,0,"QNR",&arrFit);
2813 hisRRphiRes=(TH1D*)arrFit.At(2)->Clone();
2814 meanDistortion->Draw("ddeltaR0:r0>>hisdR(40,85,245,100,-0.25,0.25)","abs(z0)<85.&&izNorm==1","colz");
2815 his2D = (TH2D*)meanDistortion->GetHistogram();
2816 his2D->FitSlicesY(0,0,-1,0,"QNR",&arrFit);
2817 hisRRes=(TH1D*)arrFit.At(2)->Clone();
2818 hisRRphiRes->SetMarkerStyle(25);
2819 hisRRes->SetMarkerStyle(21);
2820 hisRRphiRes->SetMarkerColor(2);
2821 hisRRes->SetMarkerColor(4);
2822
2823 hisRRes->GetXaxis()->SetTitle("R (cm)");
2824 hisRRes->GetYaxis()->SetTitle("#sigma_{res} (cm)");
2825 hisRRes->SetMinimum(0);
2826 TCanvas * canvasResidualsFit = new TCanvas("canvasResidualsFit","canvasResidualsFit",600,500);
2827 TLegend * legendRRPhiFitSigma = new TLegend(0.2,0.6,0.88,0.88,"Distortion residuals. RMS(#Delta-c(z)#Delta_{ref}) (|z|<85)");
2828 legendRRPhiFit->SetBorderSize(0);
2829 hisRRes->Draw("");
2830 hisRRphiRes->Draw("same");
2831 legendRRPhiFitSigma->AddEntry(hisRRes,"Radial");
2832 legendRRPhiFitSigma->AddEntry(hisRRphiRes,"R#phi");
2833 legendRRPhiFitSigma->Draw();
2834 canvasResidualsFit->SaveAs("canvasResidualsFit.pdf"); //~/hera/alice/miranov/SpaceCharge/Fluctuations/PbPbWithGain/dirmergeAll/dEpsilon20/canvasResidualsFit.png
2835 canvasResidualsFit->SaveAs("canvasResidualsFit.png");
2836 //
2837 // 4. Residual distortion after z dependent q distortion - local phi average
2838 //
2839 TH1D * hisRResPhi, *hisRRphiResPhi=0;
2840 meanDistortion->Draw("ddeltaR0-vecDMeanR.fElements:vecDPhi.fElements>>hisRResPhi(18,0.0,0.32,100,-0.3,0.3)","abs(z0)<85.&&izNorm==1&&r0<120","colz",100000);
2841 his2D = (TH2D*)meanDistortion->GetHistogram()->Clone();
2842 his2D->FitSlicesY(0,0,-1,0,"QNR",&arrFit);
2843 hisRResPhi=(TH1D*)arrFit.At(2)->Clone();
2844 //
2845 meanDistortion->Draw("ddeltaRPhi0-vecDMeanRPhi.fElements:vecDPhi.fElements>>hisRRphResPhi(18,0.0,0.32,100,-0.3,0.3)","abs(z0)<85.&&izNorm==1","colz",100000);
2846 his2D = (TH2D*)meanDistortion->GetHistogram()->Clone();
2847 his2D->FitSlicesY(0,0,-1,0,"QNR",&arrFit);
2848 hisRRphiResPhi=(TH1D*)arrFit.At(2)->Clone();
2849 //
2850 hisRRphiResPhi->SetMarkerStyle(25);
2851 hisRResPhi->SetMarkerStyle(21);
2852 hisRRphiResPhi->SetMarkerColor(2);
2853 hisRResPhi->SetMarkerColor(4);
2854 hisRResPhi->GetXaxis()->SetTitle("#Delta#phi bin width");
2855 hisRResPhi->GetYaxis()->SetTitle("#sigma_{res} (cm)");
2856 hisRResPhi->SetMinimum(0);
2857 hisRResPhi->SetMaximum(1.5*hisRResPhi->GetMaximum());
2858 gStyle->SetOptStat(0);
2859 TCanvas * canvasResidualsFitPhi = new TCanvas("canvasResidualsFitPhi","canvasResidualsFitPhi",600,500);
2860 TLegend * legendRRPhiFitSigmaPhi = new TLegend(0.2,0.6,0.88,0.88,"Distortion residuals-mean in bin. RMS((#Delta-c(z)#Delta_{ref})) (|z|<85)");
2861 {
2862 hisRResPhi->Draw("");
2863 hisRRphiResPhi->Draw("same");
2864 legendRRPhiFitSigmaPhi->AddEntry(hisRResPhi,"Radial");
2865 legendRRPhiFitSigmaPhi->SetBorderSize(0);
2866 legendRRPhiFitSigmaPhi->AddEntry(hisRRphiResPhi,"R#phi");
2867 legendRRPhiFitSigmaPhi->Draw();
2868 }
2869
2870 canvasResidualsFitPhi->SaveAs("canvasResidualsFitPhi.pdf"); //~/hera/alice/miranov/SpaceCharge/Fluctuations/PbPbWithGain/dirmergeAll/dEpsilon20/canvasResidualsFit.png
2871 canvasResidualsFitPhi->SaveAs("canvasResidualsFitPhi.png");
2872
2873
2874}
2875
2876