]> 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
o Add Qmax information
[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"
52//
53// constants
54//
55Double_t omegaTau=0.325;
7212f8c3 56//
57// Function declaration
58//
59// TOY MC
60void spaceChargeFluctuationToyMC(Int_t nframes, Double_t interactionRate);
61void spaceChargeFluctuationToyDraw();
62void spaceChargeFluctuationToyDrawSummary();
c0172f82 63
7212f8c3 64//
65// RAW data analysis
66//
76dbc35e 67TH1 * GenerateMapRawIons(Int_t useGain,const char *fileName="raw.root", const char *outputName="histo.root", Int_t maxEvents=25);
0eff8c53 68void DoMerge();
7212f8c3 69void AnalyzeMaps1D(); // make nice plots
70void MakeFluctuationStudy3D(Int_t nhistos, Int_t nevents, Int_t niter);
71TH3D * NormalizeHistoQ(TH3D * hisInput, Bool_t normEpsilon);
4b955de3 72//
7212f8c3 73TH3D * PermutationHistoZ(TH3D * hisInput, Double_t deltaZ);
74TH3D * PermutationHistoPhi(TH3D * hisInput, Double_t deltaPhi);
4b955de3 75TH3D * PermutationHistoLocalPhi(TH3D * hisInput, Int_t deltaPhi);
76void MakeSpaceChargeFluctuationScan(Double_t scale, Int_t nfiles, Int_t sign);
7212f8c3 77void DrawFluctuationdeltaZ(Int_t stat=0, Double_t norm=10000);
78void DrawFluctuationSector(Int_t stat=0, Double_t norm=10000);
c0172f82 79
4b955de3 80void spaceChargeFluctuation(Int_t mode=0, Float_t arg0=0, Float_t arg1=0, Float_t arg2=0){
c0172f82 81 //
7212f8c3 82 // function called from the shell script
c0172f82 83 //
7212f8c3 84 gRandom->SetSeed(0);
76dbc35e 85 if (mode==0) GenerateMapRawIons(arg0);
0eff8c53 86 if (mode==1) DoMerge();
7212f8c3 87 if (mode==2) spaceChargeFluctuationToyMC(arg0,arg1);
88 if (mode==3) MakeFluctuationStudy3D(10000, arg0, arg1);
4b955de3 89 if (mode==4) MakeSpaceChargeFluctuationScan(arg0,arg1,arg2); // param: scale, nfiles, sign Bz
7212f8c3 90 if (mode==5) {
91 DrawFluctuationdeltaZ(arg0,arg1);
92 DrawFluctuationSector(arg0,arg1);
93 }
c0172f82 94}
06ab15a8 95
96
97Double_t RndmdNchdY(Double_t s){
98 //
99 // dNch/deta - last 2 points inventeted (to find it somewhere ?)
100 //
101 // http://arxiv.org/pdf/1012.1657v2.pdf - table 1. ALICE PbPb
102 // Scaled according s^0.15
103 // http://arxiv.org/pdf/1210.3615v2.pdf
104 // This we can cite.
105 // Usage example::
106 /*
107 TH1F his550("his550","his550",1000,0,3000)
108 for (Int_t i=0; i<300000; i++) his550->Fill(RndmdNchdY(5.5));
109 his550->Draw();
110 TF1 f1("f1","[0]*x^(-(0.00001+abs([1])))",1,2000)
111 f1.SetParameters(1,-1)
112 his550->Fit("f1","","",10,3000);
113 TH1F his276("his276","his276",1000,0,3000)
2ff913b6 114 for (Int_t i=0; i<300000; i++) his276->Fill(RndmdNchdY(2.76));
06ab15a8 115 his276->Draw();
116
117 */
118 static TSpline3 * spline276=0;
2ff913b6 119 const Double_t sref=2.76; // reference s
120
06ab15a8 121 if (!spline276){
2ff913b6 122 // Refence multiplicities for 2.76 TeV
06ab15a8 123 // multplicity from archive except of the last point was set to 0
2ff913b6 124 //
125 const Double_t mult[20]={1601, 1294, 966, 649, 426, 261, 149, 76, 35, 0.001};
126 const Double_t cent[20]={2.5, 7.5, 15, 25, 35, 45, 55, 65, 75, 100.};
06ab15a8 127 TGraphErrors * gr = new TGraphErrors(10,cent,mult);
128 spline276 = new TSpline3("spline276",gr);
129 }
2ff913b6 130 Double_t norm = TMath::Power((s*s)/(sref*sref),0.15);
06ab15a8 131 spline276->Eval(gRandom->Rndm()*100.);
132 return spline276->Eval(gRandom->Rndm()*100.)*norm;
133}
134
135
136
137
138
4b955de3 139void pileUpToyMC(Int_t nframes){
140 //
141 //
142 //
143 /*
144 Int)t nframes=1000;
145 */
146 TTreeSRedirector *pcstream = new TTreeSRedirector("pileup.root","recreate");
4b955de3 147 Double_t central = 2350;
148 Double_t pmean=5;
149 TVectorD vectorT(nframes);
150 //
151 for (Int_t irate=1; irate<10; irate++){
152 printf("rate\t%d\n",irate);
153 for (Int_t iframe=0; iframe<nframes; iframe++){
154 if (iframe%100000==0)printf("iframe=%d\n",iframe);
155 Int_t ntracksAll=0;
156 Int_t nevents=gRandom->Poisson(irate);
157 Int_t ntracks=0; // to be taken from the MB primary distribution
158 Bool_t hasCentral=0;
159 for (Int_t ievent=0; ievent<nevents; ievent++){
06ab15a8 160 ntracks=RndmdNchdY(5.5);
4b955de3 161 ntracksAll+=ntracks;
162 if (ntracks>central) hasCentral = kTRUE;
163 }
164 (*pcstream)<<"pileupFrame"<<
165 "rate="<<irate<<
166 "nevents="<<nevents<<
167 "ntracks="<<ntracks<<
168 "ntracksAll="<<ntracksAll<<
169 "hasCentral"<<hasCentral<<
170 "\n";
171 vectorT[iframe]=ntracksAll;
172 }
173 Double_t mean = TMath::Mean(nframes, vectorT.GetMatrixArray());
174 Double_t rms = TMath::RMS(nframes, vectorT.GetMatrixArray());
175 Double_t median = TMath::Median(nframes, vectorT.GetMatrixArray());
176 Double_t ord90 = TMath::KOrdStat(nframes,vectorT.GetMatrixArray() , Int_t(nframes*0.90));
177 Double_t ord95 = TMath::KOrdStat(nframes,vectorT.GetMatrixArray() , Int_t(nframes*0.95));
178 Double_t ord99 = TMath::KOrdStat(nframes,vectorT.GetMatrixArray() , Int_t(nframes*0.99));
179 Double_t ord999 = TMath::KOrdStat(nframes,vectorT.GetMatrixArray() , Int_t(nframes*0.999));
180 Double_t ord9999 = TMath::KOrdStat(nframes,vectorT.GetMatrixArray() , Int_t(nframes*0.9999));
181 (*pcstream)<<"pileup"<<
182 "rate="<<irate<<
183 "mean="<<mean<<
184 "rms="<<rms<<
185 "median="<<median<<
186 "ord90="<<ord90<<
187 "ord95="<<ord95<<
188 "ord99="<<ord99<<
189 "ord999="<<ord999<<
190 "ord9999="<<ord9999<<
191 "\n";
192 }
193 delete pcstream;
194 // Draw
195 pcstream = new TTreeSRedirector("pileup.root","update");
196 TTree * treeStat = (TTree*)(pcstream->GetFile()->Get("pileup"));
197 TTree * treeFrame = (TTree*)(pcstream->GetFile()->Get("pileupFrame"));
198 Int_t mentries = treeStat->Draw("ord999","1","goff");
199 Double_t maximum = TMath::MaxElement(mentries, treeStat->GetV1());
200 const char * names[6]={"mean","median","ord90","ord95","ord99","ord999"};
201 const char * titles[6]={"Mean","Median","90 %","95 %","99 %","99.9 %"};
202 const Int_t mcolors[6]={1,2,3,4,6,7};
203 //
204 //
205 TF1 * f1 = new TF1("f1","[0]*x+[1]*sqrt(x)");
206 Double_t par0=0;
207 //
208 TCanvas * canvasMult = new TCanvas("canvasCumul","canvasCumul");
209 canvasMult->SetLeftMargin(0.13);
210 TLegend * legend= new TLegend(0.14,0.6,0.45,0.89, "Effective dN_{ch}/d#eta");
211 TGraphErrors *graphs[6]={0};
212 for (Int_t igr=0; igr<6; igr++){
213 graphs[igr] = TStatToolkit::MakeGraphErrors(treeStat,Form("%s:rate",names[igr]),"1",21+(igr%5),mcolors[igr],0);
214 graphs[igr]->SetMinimum(0);
215 graphs[igr]->GetYaxis()->SetTitleOffset(1.3);
216 graphs[igr]->SetMaximum(maximum*1.1);
217 graphs[igr]->GetXaxis()->SetTitle("<N_{ev}>");
218 graphs[igr]->GetYaxis()->SetTitle("dN_{ch}/d#eta");
219 TF1 * f2 = new TF1("f2","[0]*x+[1]*sqrt(x)");
220 f2->SetLineColor(mcolors[igr]);
221 f2->SetLineWidth(0.5);
222 if (igr>0) f2->FixParameter(0,par0);
223 graphs[igr]->Fit(f2,"","");
224 if (igr==0) par0=f2->GetParameter(0);
225 if (igr==0) graphs[igr]->Draw("ap");
226 graphs[igr]->Draw("p");
227 legend->AddEntry(graphs[igr], titles[igr],"p");
228 }
229 legend->SetBorderSize(0);
230 legend->Draw();
231
232 canvasMult->SaveAs("effectiveMult.pdf");
233 canvasMult->SaveAs("effectiveMult.png");
234 gStyle->SetOptStat(0);
235 TH2F * hisMult = new TH2F("ntracksNevent","ntracksnevents",9,1,10,100,0,2*maximum);
236 {
237 treeFrame->Draw("ntracksAll:rate>>ntracksNevent","","colz");
238 hisMult->GetXaxis()->SetTitle("<N_{ev}>");
239 hisMult->GetYaxis()->SetTitle("dN_{ch}/d#eta");
240 hisMult->GetYaxis()->SetTitleOffset(1.3);
241 hisMult->Draw("colz");
242 }
243 canvasMult->SaveAs("effectiveMultColz.pdf");
244 canvasMult->SaveAs("effectiveMultColz.png");
245 //
246 //
247 //
248 TH2F * hisMult5 = new TH2F("ntracksNevent5","ntracksnEvents5",9,1,10,100,0,maximum);
249 {
250 treeFrame->Draw("ntracksAll:nevents>>ntracksNevent5","abs(rate-5)<0.5","colz");
251 hisMult5->GetXaxis()->SetTitle("N_{ev}");
252 hisMult5->GetYaxis()->SetTitle("dN_{ch}/d#eta");
253 hisMult5->GetYaxis()->SetTitleOffset(1.3);
254 hisMult5->Draw("colz");
255 }
256 canvasMult->SaveAs("effectiveMultF5.pdf");
257 canvasMult->SaveAs("effectiveMultF5.png");
2ff913b6 258
259 {
260 gStyle->SetOptFit(1);
261 gStyle->SetOptStat(1);
262 gStyle->SetOptTitle(1);
263 TCanvas * canvasMultH = new TCanvas("canvasCumulH","canvasCumulH",700,700);
264 canvasMultH->Divide(1,2);
265 canvasMultH->cd(1);
266 TH1F his550("his550","his550",1000,0,3000);
267 TH1F his276("his276","his276",1000,0,3000);
268 for (Int_t i=0; i<300000; i++) his550.Fill(RndmdNchdY(5.5));
269 for (Int_t i=0; i<300000; i++) his276.Fill(RndmdNchdY(2.76));
270 TF1 f1("f1","[0]*x^(-(0.00001+abs([1])))",1,2000);
271 f1.SetParameters(1,-1);
272 his550.GetXaxis()->SetTitle("dN_{ch}/d#eta");
273 his276.GetXaxis()->SetTitle("dN_{ch}/d#eta");
274 his550.Fit("f1","","",10,3000);
275 his276.Fit("f1","","",10,3000);
276 canvasMultH->cd(1)->SetLogx(1);
277 canvasMultH->cd(1)->SetLogy(1);
278 his550.Draw();
279 canvasMultH->cd(2)->SetLogx(1);
280 canvasMultH->cd(2)->SetLogy(1);
281 his276.Draw("");
282 canvasMultH->SaveAs("dNchdEta.pdf");
283 }
4b955de3 284 delete pcstream;
285}
c0172f82 286
7212f8c3 287void spaceChargeFluctuationToyMC(Int_t nframes, Double_t interactionRate){
c0172f82 288 //
7212f8c3 289 // Toy MC to generate space charge fluctuation, to estimate the fluctuation of the integral space charge in part of the
290 // TPC
291 // Parameters:
292 // nframes - number of frames to simulate
293 // 1. Make a toy simulation part for given setup
294 // 2. Make a summary plots for given setups - see function spaceChargeFluctuationToyMCDraw()
c0172f82 295 //
296 TTreeSRedirector *pcstream = new TTreeSRedirector("spaceChargeFluctuation.root","recreate");
4e715689 297 Double_t driftTime=0.1;
4e715689 298 Double_t eventMean=interactionRate*driftTime;
299 Double_t trackMean=500;
300 Double_t FPOT=1.0, EEND=3000;
301 Double_t EEXPO=0.8567;
302 const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
c0172f82 303
304 for (Int_t iframe=0; iframe<nframes; iframe++){
4e715689 305 printf("iframe=%d\n",iframe);
c0172f82 306 Int_t nevents=gRandom->Poisson(interactionRate*driftTime);
307 Int_t ntracksAll=0;
4e715689 308 TVectorD vecTracksPhi180(180);
309 TVectorD vecTracksPhi36(36);
310 TVectorD vecEPhi180(180);
311 TVectorD vecEPhi36(36);
312 Double_t dESum=0;
c0172f82 313 for (Int_t ievent=0; ievent<nevents; ievent++){
4e715689 314 Int_t ntracks=gRandom->Exp(trackMean); // to be taken from the MB primary distribution
315 Float_t RAN = gRandom->Rndm();
316 ntracks=TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/2.;
c0172f82 317 ntracksAll+=ntracks;
4e715689 318 for (Int_t itrack=0; itrack<ntracks; itrack++){
319 Double_t phi = gRandom->Rndm();
320 vecTracksPhi180(Int_t(phi*180))+=1;
321 vecTracksPhi36(Int_t(phi*36)) +=1;
322 // simplified MC to get track length including loopers
323 Double_t theta= gRandom->Rndm();
324 Double_t pt = gRandom->Exp(0.5)+0.05;
325 Double_t crv = TMath::Abs(5*kB2C/pt); //GetC(b); // bz*kB2C/pt;
326 Double_t deltaPhi=0;
327 if (TMath::Abs(2*crv*(245-85)/2.) <1.) deltaPhi=TMath::ASin(crv*(245-85)/2.);
328 else
329 deltaPhi=TMath::Pi();
330 Double_t dE=deltaPhi/crv;
331 Double_t xloop=1;
332 if (1./crv<250) {
333 xloop = TMath::Min(1./(TMath::Abs(theta)+0.0001),10.);
334 if (xloop<1) xloop=1;
335 }
336 dESum+=xloop*dE;
337 if (itrack==0) (*pcstream)<<"track"<<
338 "pt="<<pt<<
339 "crv="<<crv<<
340 "theta="<<theta<<
341 "dE="<<dE<<
342 "xloop="<<xloop<<
343 "\n";
344
345 vecEPhi180(Int_t(phi*180)) +=dE*xloop;
346 vecEPhi36(Int_t(phi*36)) +=dE*xloop;
c0172f82 347 }
4e715689 348 (*pcstream)<<"event"<<
349 "ntracks="<<ntracks<<
350 "nevents="<<nevents<<
351 "\n";
c0172f82 352 }
353 (*pcstream)<<"ntracks"<<
7212f8c3 354 "rate="<<interactionRate<< // interaction rate
355 "eventMean="<<eventMean<< // mean number of events per frame
356 "trackMean="<<trackMean<< // assumed mean of the tracks per event
357 //
4e715689 358 "nevents="<<nevents<< // number of events withing time frame
359 "ntracksAll="<<ntracksAll<< // number of tracks within time frame
360 "dESum="<<dESum<< // sum of the energy loss
361 "vecTracksPhi36.="<<&vecTracksPhi36<< // number of tracks in phi bin (36 bins) within time frame
362 "vecTracksPhi180.="<<&vecTracksPhi180<< // number of tracks in phi bin (180 bins) within time frame
363 "vecEPhi36.="<<&vecEPhi36<< // number of tracks in phi bin (36 bins) within time frame
364 "vecEPhi180.="<<&vecEPhi180<< // number of tracks in phi bin (180 bins) within time frame
c0172f82 365 "\n";
366 }
367 delete pcstream;
7212f8c3 368 spaceChargeFluctuationToyDraw();
369}
370
371
372void spaceChargeFluctuationToyDraw(){
4e715689 373 //
374 // Toy MC to simulate the space charge integral fluctuation
7212f8c3 375 // Draw function for given setup
376 // for MC generation part see : void spaceChargeFluctuationToyMC
377 TTreeSRedirector *pcstream = new TTreeSRedirector("spaceChargeFluctuation.root","update");
378 TFile * f = pcstream->GetFile();
4e715689 379 TTree * treeStat = (TTree*)f->Get("ntracks");
380 TTree * treedE = (TTree*)f->Get("track");
381 TTree * treeEv = (TTree*)f->Get("event");
7212f8c3 382
4e715689 383 Int_t nentries=treedE->Draw("dE*xloop","1","",1000000);
7212f8c3 384
4e715689 385 Double_t meandE=TMath::Mean(nentries,treedE->GetV1());
386 Double_t rmsdE=TMath::RMS(nentries,treedE->GetV1());
387 treeStat->SetAlias("meandE",Form("(%f+0)",meandE));
388 treeStat->SetAlias("rmsdE",Form("(%f+0)",rmsdE));
389 nentries=treeEv->Draw("ntracks","1","",1000000);
390 Double_t meanT=TMath::Mean(nentries,treeEv->GetV1());
391 Double_t rmsT=TMath::RMS(nentries,treeEv->GetV1());
392 treeStat->SetAlias("tracksMean",Form("(%f+0)",meanT));
393 treeStat->SetAlias("tracksRMS",Form("(%f+0)",rmsT));
7212f8c3 394 nentries = treeStat->Draw("eventMean","","");
395 Double_t meanEvents =TMath::Mean(nentries,treeStat->GetV1());
396 treeStat->SetMarkerStyle(21);
397 treeStat->SetMarkerSize(0.4);
4e715689 398 //
7212f8c3 399 const Int_t kColors[6]={1,2,3,4,6,7};
400 const Int_t kStyle[6]={20,21,24,25,24,25};
401 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)"};
402 const char * hnames[6]={"Events","Tracks","TracksPhi180","QPhi180", "TracksPhi36","QPhi36"};
403
404 TH1* hisFluc[6]={0};
405 TH1* hisPull[6]={0};
406 TVectorD *vecFitFluc[6]={0};
407 TVectorD *vecFitFlucPull[6]={0};
408 //
409 // histograms
410 //
411 treeStat->Draw("nevents/eventMean>>hisEv(100,0.85,1.15)","");
412 hisFluc[0]=(TH1*)treeStat->GetHistogram()->Clone();
413 treeStat->Draw("ntracksAll/(eventMean*tracksMean)>>hisTrackAll(100,0.85,1.1)","","same");
414 hisFluc[1]=(TH1*)treeStat->GetHistogram()->Clone();
415 treeStat->Draw("vecTracksPhi180.fElements/(eventMean*tracksMean/180)>>hisTrackSector(100,0.85,1.1)","1/180","same");
416 hisFluc[2]=(TH1*)treeStat->GetHistogram()->Clone();
417 treeStat->Draw("vecEPhi180.fElements/(eventMean*tracksMean*meandE/180)>>hisdESector(100,0.85,1.1)","1/180","same");
418 hisFluc[3]=(TH1*)treeStat->GetHistogram()->Clone();
419 treeStat->Draw("vecTracksPhi36.fElements/(eventMean*tracksMean/36)>>hisTrackSector36(100,0.85,1.1)","1/36","same");
420 hisFluc[4]=(TH1*)treeStat->GetHistogram()->Clone();
421 treeStat->Draw("vecEPhi36.fElements/(eventMean*tracksMean*meandE/36)>>hisdESector36(100,0.85,1.1)","1/36","same");
422 hisFluc[5]=(TH1*)treeStat->GetHistogram()->Clone();
423 //
424 // pulls
4e715689 425 //
7212f8c3 426 treeStat->Draw("((nevents/eventMean)-1)/sqrt(1/eventMean)>>pullEvent(100,-6,6)","","err"); //tracks All pull
427 hisPull[0]=(TH1*)treeStat->GetHistogram()->Clone();
428 treeStat->Draw("(ntracksAll/(eventMean*tracksMean)-1)/sqrt(1/eventMean+(tracksRMS/tracksMean)**2/eventMean)>>pullTrackAll(100,-6,6)","","err"); //tracks All pull
429 hisPull[1]=(TH1*)treeStat->GetHistogram()->Clone();
430 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
431 hisPull[2]=(TH1*)treeStat->GetHistogram()->Clone();
432 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
433 hisPull[3]=(TH1*)treeStat->GetHistogram()->Clone();
434 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
435 hisPull[4]=(TH1*)treeStat->GetHistogram()->Clone();
436 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
437 hisPull[5]=(TH1*)treeStat->GetHistogram()->Clone();
4e715689 438 //
7212f8c3 439
440 for (Int_t ihis=0; ihis<6; ihis++) {
441 vecFitFluc[ihis] = new TVectorD(3);
442 vecFitFlucPull[ihis] = new TVectorD(3);
443 TF1 * fg = new TF1(Form("fg%d",ihis),"gaus");
444 fg->SetLineWidth(0.5);
445 fg->SetLineColor(kColors[ihis]);
446 hisFluc[ihis]->Fit(fg,"","");
447 fg->GetParameters( vecFitFluc[ihis]->GetMatrixArray());
448 hisPull[ihis]->Fit(fg,"","");
449 fg->GetParameters( vecFitFlucPull[ihis]->GetMatrixArray());
450 hisFluc[ihis]->SetName(Form("Fluctuation%s",hnames[ihis]));
451 hisFluc[ihis]->SetTitle(Form("Fluctuation%s",htitles[ihis]));
452 hisPull[ihis]->SetName(Form("Pull%s",hnames[ihis]));
453 hisPull[ihis]->SetTitle(Form("Pull%s",htitles[ihis]));
454 }
455
456 gStyle->SetOptStat(0);
457 TCanvas * canvasQFluc= new TCanvas("SpaceChargeFluc","SpaceChargeFluc",600,700);
458 canvasQFluc->Divide(1,2);
459 canvasQFluc->cd(1);
460 TLegend * legendFluc = new TLegend(0.11,0.55,0.45,0.89,"Relative fluctuation");
461 TLegend * legendPull = new TLegend(0.11,0.55,0.45,0.89,"Fluctuation pulls");
462 for (Int_t ihis=0; ihis<6; ihis++){
463 hisFluc[ihis]->SetMarkerStyle(kStyle[ihis]);
464 hisFluc[ihis]->SetMarkerColor(kColors[ihis]);
465 hisFluc[ihis]->SetMarkerSize(0.8);
466 if (ihis==0) hisFluc[ihis]->Draw("err");
467 hisFluc[ihis]->Draw("errsame");
468 legendFluc->AddEntry(hisFluc[ihis],htitles[ihis]);
469 }
470 legendFluc->Draw();
471
472 canvasQFluc->cd(2);
473 for (Int_t ihis=0; ihis<6; ihis++){
474 hisPull[ihis]->SetMarkerStyle(kStyle[ihis]);
475 hisPull[ihis]->SetMarkerColor(kColors[ihis]);
476 hisPull[ihis]->SetMarkerSize(0.8);
477 if (ihis==0) hisPull[ihis]->Draw("err");
478 hisPull[ihis]->Draw("errsame");
479 legendPull->AddEntry(hisPull[ihis],htitles[ihis]);
480 }
481 legendPull->Draw();
4e715689 482 //
7212f8c3 483 for (Int_t ihis=0; ihis<6; ihis++){
484 hisFluc[ihis]->Write();
485 hisPull[ihis]->Write();
486 }
487 (*pcstream)<<"summary"<< // summary information for given setup
488 "meanEvents="<<meanEvents<< // mean number of events in the frame
489 "meandE="<<meandE<< // mean "energy loss" of track
490 "rmsdE="<<rmsdE<< // rms
491 "meanT="<<meanT<< // mean number of tracks per MB event
492 "rmsT="<<rmsT<< // rms of onumber of tracks
493 // // fit of the relative fluctuation
494 "vflucE.="<<vecFitFluc[0]<< // in events
495 "vflucEP.="<<vecFitFlucPull[0]<< // in events pull
496 "vflucTr.="<<vecFitFluc[1]<< // in tracks
497 "vflucTrP.="<<vecFitFlucPull[1]<<
498 //
499 "vflucTr180.="<<vecFitFluc[2]<<
500 "vflucTr180P.="<<vecFitFlucPull[2]<<
501 "vflucE180.="<<vecFitFluc[3]<<
502 "vflucE180P.="<<vecFitFlucPull[3]<<
503 //
504 "vflucTr36.="<<vecFitFluc[4]<<
505 "vflucTr36P.="<<vecFitFlucPull[4]<<
506 "vflucE36.="<<vecFitFluc[5]<<
507 "vflucE36P.="<<vecFitFlucPull[5]<<
508 "\n";
509 canvasQFluc->SaveAs("CanvasFluctuation.pdf");
510 canvasQFluc->SaveAs("CanvasFluctuation.png");
511 delete pcstream;
c0172f82 512
4e715689 513}
7212f8c3 514
515void spaceChargeFluctuationToyDrawSummary(){
4e715689 516 //
7212f8c3 517 // make a summary information plots using several runs with differnt mean IR setting
518 // Input:
519 // space.list - list of root files produced by spaceChargeFluctuationToyDraw
520 // Output:
521 // canvas saved in current directory
4e715689 522 //
7212f8c3 523 TChain * chain = AliXRDPROOFtoolkit::MakeChain("space.list","summary",0,100);
524 chain->SetMarkerStyle(21);
525 const Int_t kColors[6]={1,2,3,4,6,7};
526 const Int_t kStyle[6]={20,21,24,25,24,25};
527 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)"};
528 // const char * hnames[6]={"Events","Tracks","TracksPhi180","QPhi180", "TracksPhi36","QPhi36"};
529 //
530 Double_t meanT,rmsT=0;
531 Double_t meandE,rmsdE=0;
532 Int_t entries = chain->Draw("meanT:rmsT:meandE:rmsdE","1","goff");
533 meanT =TMath::Mean(entries, chain->GetV1());
534 rmsT =TMath::Mean(entries, chain->GetV2());
535 meandE =TMath::Mean(entries, chain->GetV3());
536 rmsdE =TMath::Mean(entries, chain->GetV4());
537 //
538 //
539 //
540 TGraphErrors * graphs[6]={0};
541 TF1 * functions[6]={0};
542
543 graphs[5]=TStatToolkit::MakeGraphErrors(chain,"vflucE36.fElements[2]:meanEvents:0.025*vflucE36.fElements[2]","1",kStyle[5],kColors[5],1);
544 graphs[4]=TStatToolkit::MakeGraphErrors(chain,"vflucTr36.fElements[2]:meanEvents:0.025*vflucTr36.fElements[2]","1",kStyle[4],kColors[4],1);
545 graphs[3]=TStatToolkit::MakeGraphErrors(chain,"vflucE180.fElements[2]:meanEvents:0.025*vflucE180.fElements[2]","1",kStyle[3],kColors[3],1);
546 graphs[2]=TStatToolkit::MakeGraphErrors(chain,"vflucTr180.fElements[2]:meanEvents:0.025*vflucTr180.fElements[2]","1",kStyle[2],kColors[2],1);
547 graphs[1]=TStatToolkit::MakeGraphErrors(chain,"vflucTr.fElements[2]:meanEvents:0.025*vflucTr.fElements[2]","1",kStyle[1],kColors[1],1);
548 graphs[0]=TStatToolkit::MakeGraphErrors(chain,"vflucE.fElements[2]:meanEvents:0.025*vflucE.fElements[2]","1",kStyle[0],kColors[0],1);
549
550 functions[5]=new TF1("fe","sqrt(1+[0]**2+[1]/[2]+[1]*[3]**2/[2])/sqrt(x)",2000,200000);
551 functions[5]->SetParameters(rmsT/meanT,36.,meanT,rmsdE/meandE);
552 functions[4]=new TF1("fe","sqrt(1+[0]**2+[1]/[2])/sqrt(x)",2000,200000);
553 functions[4]->SetParameters(rmsT/meanT,36.,meanT,0);
554 functions[3]=new TF1("fe","sqrt(1+[0]**2+[1]/[2]+[1]*[3]**2/[2])/sqrt(x)",2000,200000);
555 functions[3]->SetParameters(rmsT/meanT,180.,meanT,rmsdE/meandE);
556 functions[2]=new TF1("fe","sqrt(1+[0]**2+[1]/[2])/sqrt(x)",2000,200000);
557 functions[2]->SetParameters(rmsT/meanT,180.,meanT,0);
558 functions[1]=new TF1("fe","sqrt(1+[0]**2)/sqrt(x)",2000,200000);
559 functions[1]->SetParameters(rmsT/meanT,0);
560 functions[0]=new TF1("fe","sqrt(1)/sqrt(x)",2000,200000);
561
562
563 TCanvas *canvasF= new TCanvas("fluc","fluc",600,500);
564 // TLegend *legend = new TLegend(0.5,0.65,0.89,0.89,"Relative fluctuation #sigma=#sqrt{1+#frac{#sigma_{T}^{2}}{#mu_{T}^{2}}}");
565 TLegend *legendF = new TLegend(0.45,0.5,0.89,0.89,"Relative fluctuation of charge");
566 for (Int_t ihis=0; ihis<4; ihis++){
567 graphs[ihis]->SetMinimum(0.00);
568 graphs[ihis]->SetMaximum(0.05);
569 if (ihis==0) graphs[ihis]->Draw("ap");
570 graphs[ihis]->GetXaxis()->SetTitle("events");
571 graphs[ihis]->GetXaxis()->SetNdivisions(507);
572 graphs[ihis]->GetYaxis()->SetTitle("#frac{#sigma}{#mu}");
573 graphs[ihis]->Draw("p");
574 legendF->AddEntry(graphs[ihis],htitles[ihis],"p");
575 if (functions[ihis]){
576 functions[ihis]->SetLineColor(kColors[ihis]);
577 functions[ihis]->SetLineWidth(0.5);
578 functions[ihis]->Draw("same");
579 }
580 }
581 legendF->Draw();
582 canvasF->SaveAs("spaceChargeFlucScan.pdf");
583 canvasF->SaveAs("spaceChargeFlucScan.png");
584
585 TCanvas *canvasF36= new TCanvas("fluc36","fluc36",600,500);
586 // TLegend *legend = new TLegend(0.5,0.65,0.89,0.89,"Relative fluctuation #sigma=#sqrt{1+#frac{#sigma_{T}^{2}}{#mu_{T}^{2}}}");
587 TLegend *legendF36 = new TLegend(0.45,0.5,0.89,0.89,"Relative fluctuation of charge");
588 for (Int_t ihis=0; ihis<6; ihis++){
589 if (ihis==2 || ihis==3) continue;
590 graphs[ihis]->SetMinimum(0.00);
591 graphs[ihis]->SetMaximum(0.05);
592 if (ihis==0) graphs[ihis]->Draw("ap");
593 graphs[ihis]->GetXaxis()->SetTitle("events");
594 graphs[ihis]->GetXaxis()->SetNdivisions(507);
595 graphs[ihis]->GetYaxis()->SetTitle("#frac{#sigma}{#mu}");
596 graphs[ihis]->Draw("p");
597 legendF36->AddEntry(graphs[ihis],htitles[ihis],"p");
598 if (functions[ihis]){
599 functions[ihis]->SetLineColor(kColors[ihis]);
600 functions[ihis]->SetLineWidth(0.5);
601 functions[ihis]->Draw("same");
602 }
603 }
604 legendF36->Draw();
605 canvasF36->SaveAs("spaceChargeFlucScan36.pdf");
606 canvasF36->SaveAs("spaceChargeFlucScan36.png");
607
608
609}
610
611
612
613void FitMultiplicity(const char * fname="mult_dist_pbpb.root"){
614 //
615 // Fit multiplicity distribution using as a power law in limited range
616 // const char * fname="mult_dist_pbpb.root"
617 TFile *fmult=TFile::Open(fname);
4e715689 618 TF1 f1("f1","[0]*(x+abs([2]))**(-abs([1]))",1,3000);
619 TH1* his = (TH1*) fmult->Get("mult_dist_PbPb_normalizedbywidth");
620 f1.SetParameters(his->GetEntries(),1,1);
621 his->Fit(&f1,"","",2,3000);
622
623 Double_t FPOT=1.0, EEND=3000, EEXPO= TMath::Abs(f1.GetParameter(1));
624 EEXPO=0.8567;
625 const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
626 TH1F *hisr= new TH1F("aaa","aaa",4000,0,4000);
627 for (Int_t i=0; i<400000; i++){
628 Float_t RAN = gRandom->Rndm();
629 hisr->Fill(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO));
630 }
4e715689 631}
c0172f82 632
c0172f82 633
634
635
76dbc35e 636TH1 * GenerateMapRawIons(Int_t useGainMap, const char *fileName, const char *outputName, Int_t maxEvents){
c0172f82 637 //
638 // Generate 3D maps of the space charge for the rad data maps
639 // different threshold considered
640 // Paramaters:
0eff8c53 641 // useGainMap - switch usage of the gain map
c0172f82 642 // fileName - name of input raw file
643 // outputName - name of output file with the space charge histograms
644 // maxEvents - grouping of the events
645 //
646 //
7212f8c3 647 gRandom->SetSeed(0); //set initial seed to be independent for different jobs
648
c0172f82 649 TTreeSRedirector * pcstream = new TTreeSRedirector(outputName, "recreate");
650 const char * ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";
651 AliCDBManager * man = AliCDBManager::Instance();
652 man->SetDefaultStorage(ocdbpath);
653 man->SetRun(0);
654 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 2.76/2.));
655 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
656 AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
76dbc35e 657 AliTPCCalPad * gain = AliTPCcalibDB::Instance()->GetDedxGainFactor();
658 AliTPCCalPad * noise = AliTPCcalibDB::Instance()->GetPadNoise();
659
c0172f82 660 TStopwatch timer;
661 timer.Start();
662 // arrays of space charges - different elements corresponds to different threshold to accumulate charge
76dbc35e 663 TH1D * hisQ1D[3]={0};
76dbc35e 664 TH1D * hisQ1DROC[3]={0};
7212f8c3 665 TH2D * hisQ2DRPhi[3]={0};
666 TH2D * hisQ2DRZ[3]={0};
667 TH2D * hisQ2DRPhiROC[3]={0};
668 TH2D * hisQ2DRZROC[3]={0};
669 TH3D * hisQ3D[3]={0}; // 3D maps space charge from drift volume
670 TH3D * hisQ3DROC[3]={0}; // 3D maps space charge from ROC
671
5c70ec78 672 Int_t nbinsRow=param->GetNRowLow()+param->GetNRowUp();
673 Double_t *xbins = new Double_t[nbinsRow+1];
674 xbins[0]=param->GetPadRowRadiiLow(0)-1; //underflow bin
675 for (Int_t ibin=0; ibin<param->GetNRowLow();ibin++) xbins[1+ibin]=param->GetPadRowRadiiLow(ibin);
676 for (Int_t ibin=0; ibin<param->GetNRowUp();ibin++) xbins[1+ibin+param->GetNRowLow()]=param->GetPadRowRadiiUp(ibin);
c0172f82 677 //
678 for (Int_t ith=0; ith<3; ith++){
679 char chname[100];
7212f8c3 680 // 1D
681 snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
682 hisQ1D[ith] = new TH1D(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
683 snprintf(chname,100,"hisQ1DROC_Th%d",2*ith+2);
684 hisQ1DROC[ith] = new TH1D(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
685 // 3D
c0172f82 686 snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
7212f8c3 687 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);
688 snprintf(chname,100,"hisQ3DROC_Th%d",2*ith+2);
689 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);
690 // 2D
691 snprintf(chname,100,"hisQ2DRPhi_Th%d",2*ith+2);
692 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 693 snprintf(chname,100,"hisQ2DRZ_Th%d",2*ith+2);
76dbc35e 694 hisQ2DRZ[ith] = new TH2D(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36), 125,-250,250);
c0172f82 695 //
7212f8c3 696 snprintf(chname,100,"hisQ2DRPhiROC_Th%d",2*ith+2);
697 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 698 snprintf(chname,100,"hisQ2DRZROC_Th%d",2*ith+2);
699 hisQ2DRZROC[ith] = new TH2D(chname,chname,param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36), 125,-250,250);
7212f8c3 700 //
5c70ec78 701 hisQ1D[ith]->GetXaxis()->Set(nbinsRow,xbins);
5c70ec78 702 hisQ1DROC[ith]->GetXaxis()->Set(nbinsRow,xbins);
7212f8c3 703 hisQ3D[ith]->GetYaxis()->Set(nbinsRow,xbins);
5c70ec78 704 hisQ3DROC[ith]->GetYaxis()->Set(nbinsRow,xbins);
c0172f82 705 //
7212f8c3 706 hisQ2DRPhi[ith]->GetYaxis()->Set(nbinsRow,xbins);
707 hisQ2DRZ[ith]->GetXaxis()->Set(nbinsRow,xbins);
708 hisQ2DRPhiROC[ith]->GetYaxis()->Set(nbinsRow,xbins);
709 hisQ2DRZROC[ith]->GetXaxis()->Set(nbinsRow,xbins);
710 //
c0172f82 711 hisQ1D[ith]->SetDirectory(0);
7212f8c3 712 hisQ1DROC[ith]->SetDirectory(0);
713 hisQ3D[ith]->SetDirectory(0);
c0172f82 714 hisQ3DROC[ith]->SetDirectory(0);
7212f8c3 715 //
716 hisQ2DRPhi[ith]->SetDirectory(0);
717 hisQ2DRZ[ith]->SetDirectory(0);
76dbc35e 718 hisQ2DRZROC[ith]->SetDirectory(0);
7212f8c3 719 hisQ2DRPhiROC[ith]->SetDirectory(0);
c0172f82 720 }
721 //
722 //
723 AliRawReader *reader = new AliRawReaderRoot(fileName);
724 reader->Reset();
725 AliAltroRawStream* stream = new AliAltroRawStream(reader);
726 stream->SelectRawData("TPC");
727 Int_t evtnr=0;
728 Int_t chunkNr=0;
729 //
730
731 while (reader->NextEvent()) {
0eff8c53 732 Double_t shiftZ= gRandom->Rndm()*250.;
c0172f82 733 //
734 if(evtnr>=maxEvents) {
735 chunkNr++;
c0172f82 736 pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
737 pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
738 for (Int_t ith=0; ith<3; ith++){
739 hisQ1D[ith]->Write(Form("His1DDrift_%d",ith));
7212f8c3 740 hisQ2DRPhi[ith]->Write(Form("His2DRPhiDrift_%d",ith));
741 hisQ2DRZ[ith]->Write(Form("His2DRZDrift_%d",ith));
c0172f82 742 hisQ3D[ith]->Write(Form("His3DDrift_%d",ith));
743 hisQ1DROC[ith]->Write(Form("His1DROC_%d",ith));
7212f8c3 744 hisQ2DRPhiROC[ith]->Write(Form("His2DRPhiROC_%d",ith));
76dbc35e 745 hisQ2DRZROC[ith]->Write(Form("His2DRZROC_%d",ith));
c0172f82 746 hisQ3DROC[ith]->Write(Form("His3DROC_%d",ith));
5c70ec78 747 (*pcstream)<<"histo"<<
76dbc35e 748 "events="<<evtnr<<
749 "useGain="<<useGainMap<<
5c70ec78 750 Form("hist1D_%d.=",ith*2+2)<<hisQ1D[ith]<<
7212f8c3 751 Form("hist2DRPhi_%d.=",ith*2+2)<<hisQ2DRPhi[ith]<<
76dbc35e 752 Form("hist2DRZ_%d.=",ith*2+2)<<hisQ2DRZ[ith]<<
753 Form("hist3D_%d.=",ith*2+2)<<hisQ3D[ith]<<
5c70ec78 754 Form("hist1DROC_%d.=",ith*2+2)<<hisQ1DROC[ith]<<
7212f8c3 755 Form("hist2DRPhiROC_%d.=",ith*2+2)<<hisQ2DRPhiROC[ith]<<
76dbc35e 756 Form("hist2DRZROC_%d.=",ith*2+2)<<hisQ2DRZROC[ith]<<
757 Form("hist3DROC_%d.=",ith*2+2)<<hisQ3DROC[ith];
5c70ec78 758 }
759 (*pcstream)<<"histo"<<"\n";
760 for (Int_t ith=0; ith<3; ith++){
c0172f82 761 hisQ1D[ith]->Reset();
7212f8c3 762 hisQ2DRPhi[ith]->Reset();
76dbc35e 763 hisQ2DRZ[ith]->Reset();
c0172f82 764 hisQ3D[ith]->Reset();
765 hisQ1DROC[ith]->Reset();
7212f8c3 766 hisQ2DRPhiROC[ith]->Reset();
76dbc35e 767 hisQ2DRZROC[ith]->Reset();
c0172f82 768 hisQ3DROC[ith]->Reset();
769 }
76dbc35e 770 evtnr=0;
c0172f82 771 }
772 cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
773 evtnr++;
774 AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
775 AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
776 //
777 while (input.NextDDL()){
76dbc35e 778 Int_t sector = input.GetSector();
779 AliTPCCalROC * gainROC =gain->GetCalROC(sector);
780 AliTPCCalROC * noiseROC =noise->GetCalROC(sector);
c0172f82 781 while ( input.NextChannel() ) {
782 Int_t row = input.GetRow();
c0172f82 783 Int_t pad = input.GetPad();
784 Int_t nPads = param->GetNPads(sector,row);
785 Double_t localX = param->GetPadRowRadii(sector,row);
7212f8c3 786 Double_t localY = (pad-nPads/2)*param->GetPadPitchWidth(sector);
c0172f82 787 Double_t localPhi= TMath::ATan2(localY,localX);
788 Double_t phi = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
c0172f82 789 Double_t padLength=param->GetPadPitchLength(sector,row);
76dbc35e 790 Double_t gainPad = gainROC->GetValue(row,pad);
791 Double_t noisePad = noiseROC->GetValue(row,pad);
c0172f82 792 //
793 while ( input.NextBunch() ){
794 Int_t startTbin = (Int_t)input.GetStartTimeBin();
795 Int_t bunchlength = (Int_t)input.GetBunchLength();
796 const UShort_t *sig = input.GetSignals();
797 Int_t aboveTh[3]={0};
798 for (Int_t i=0; i<bunchlength; i++){
76dbc35e 799 if (sig[i]<4*noisePad) continue;
c0172f82 800 for (Int_t ith=0; ith<3; ith++){
801 if (sig[i]>(ith*2)+2) aboveTh[ith]++;
802 }
803 }
804 for (Int_t ith=0; ith<3; ith++){
805 if (aboveTh[ith%3]>1){
806 for (Int_t i=0; i<bunchlength; i++){
807 //
808 // normalization
809 //
76dbc35e 810 Double_t zIonDrift =(param->GetZLength()-startTbin*param->GetZWidth());
c0172f82 811 zIonDrift+=shiftZ;
76dbc35e 812 Double_t signal=sig[i];
813 if (useGainMap) signal/=gainPad;
7212f8c3 814 Double_t shiftPhi = ((sector%36)<18) ? 0: TMath::TwoPi();
c0172f82 815 if (TMath::Abs(zIonDrift)<param->GetZLength()){
816 if ((sector%36)>=18) zIonDrift*=-1; // c side has opposite sign
76dbc35e 817 if (sector%36<18) hisQ1D[ith]->Fill(localX, signal/padLength);
7212f8c3 818 hisQ2DRPhi[ith]->Fill(phi+shiftPhi,localX, signal/padLength);
76dbc35e 819 hisQ2DRZ[ith]->Fill(localX, zIonDrift, signal/padLength);
820 hisQ3D[ith]->Fill(phi,localX,zIonDrift,signal/padLength);
c0172f82 821 }
822 //
823 Double_t zIonROC = ((sector%36)<18)? shiftZ: -shiftZ; // z position of the "ion disc" - A side C side opposite sign
76dbc35e 824 if (sector%36<18) hisQ1DROC[ith]->Fill(localX, signal/padLength);
7212f8c3 825 hisQ2DRPhiROC[ith]->Fill(phi+shiftPhi,localX, signal/padLength);
76dbc35e 826 hisQ2DRZROC[ith]->Fill(localX, zIonROC, signal/padLength);
827 hisQ3DROC[ith]->Fill(phi,localX,zIonROC,signal/padLength);
c0172f82 828 }
829 }
830 }
831 }
832 }
833 }
834 }
835 timer.Print();
836 delete pcstream;
837 return 0;
838}
76dbc35e 839
840
7212f8c3 841void DoMerge(){
76dbc35e 842 //
7212f8c3 843 // Merge results to the tree
76dbc35e 844 //
7212f8c3 845 TFile * fhisto = new TFile("histo.root","recreate");
846 TTree * tree = 0;
847 TChain *chain = AliXRDPROOFtoolkit::MakeChainRandom("histo.list","histo",0,100,1);
848 chain->SetBranchStatus("hist3DROC_6*",kFALSE);
849 chain->SetBranchStatus("hist3DROC_4*",kFALSE);
850 tree = chain->CopyTree("1");
851 tree->Write("histo");
852 delete fhisto;
853}
854
855
856
857
858void AnalyzeMaps1D(){
76dbc35e 859 //
7212f8c3 860 // Analyze space charge maps stored as s hitograms in trees
76dbc35e 861 //
7212f8c3 862 TFile * fhisto = new TFile("histo.root");
863 TTree * tree = (TTree*)fhisto->Get("histo");
76dbc35e 864 //
865 TH1 *his1Th[3]={0,0,0};
866 TF1 *fq1DStep= new TF1("fq1DStep","([0]+[1]*(x>134))/x**min(abs([2]),3)",85,245);
867 fq1DStep->SetParameters(1,-0.5,1);
868 tree->Draw("hist1DROC_2.fArray:hist1D_2.fXaxis.fXbins.fArray>>his(40,85,245)","","prof");
869 tree->GetHistogram()->Fit(fq1DStep);
870 // normalize step between the IROC-OROC
871 tree->SetAlias("normQ",Form("(1+%f*(hist1D_2.fXaxis.fXbins.fArray>136))",fq1DStep->GetParameter(1)/fq1DStep->GetParameter(0)));
872 //
873 {
874 Int_t entries= tree->Draw("hist1DROC_2.fArray/(events*normQ)","1","goff");
875 Double_t median=TMath::Median(entries,tree->GetV1());
876 TCut cut10Median = Form("hist1DROC_2.fArray/(events*normQ)<%f",10*median);
877 //
878 tree->Draw("hist1DROC_2.fArray/(events*normQ):hist1D_2.fXaxis.fXbins.fArray>>his1Th0(40,86,245)",cut10Median+"","prof");
879 his1Th[0] = tree->GetHistogram();
880 tree->Draw("hist1DROC_4.fArray/(events*normQ):hist1D_2.fXaxis.fXbins.fArray>>his1Th1(40,86,245)",cut10Median+"","prof");
881 his1Th[1] = tree->GetHistogram();
882 tree->Draw("hist1DROC_6.fArray/(events*normQ):hist1D_2.fXaxis.fXbins.fArray>>his1Th2(40,86,245)",cut10Median+"","prof");
883 his1Th[2]=tree->GetHistogram();
884 }
885 //
886 TCanvas *canvasR = new TCanvas("canvasR","canvasR",600,500);
887 canvasR->cd();
888 for (Int_t i=0; i<3; i++){
889 his1Th[i]->SetMarkerStyle(21);
890 his1Th[i]->SetMarkerColor(i+2);
891 fq1DStep->SetLineColor(i+2);
892 his1Th[i]->Fit(fq1DStep,"","");
893 his1Th[i]->GetXaxis()->SetTitle("r (cm)");
894 his1Th[i]->GetYaxis()->SetTitle("#frac{N_{el}}{N_{ev}}(ADC/cm)");
895 }
896 TLegend * legend = new TLegend(0.11,0.11,0.7,0.39,"1D space Charge map (ROC part) (z,phi integrated)");
897 for (Int_t i=0; i<3; i++){
898 his1Th[i]->SetMinimum(0);fq1DStep->SetLineColor(i+2);
899 his1Th[i]->Fit(fq1DStep,"qnr","qnr");
900 if (i==0) his1Th[i]->Draw("");
901 his1Th[i]->Draw("same");
902 legend->AddEntry(his1Th[i],Form("Thr=%d Slope=%2.2f",2*i+2,fq1DStep->GetParameter(2)));
903 }
904 legend->Draw();
7212f8c3 905 canvasR->SaveAs("spaceCharge1d.png");
906 canvasR->SaveAs("spaceCharge1d.eps");
76dbc35e 907 //
908 //
909 //
910}
0eff8c53 911void MakeFluctuationStudy3D(Int_t nhistos, Int_t nevents, Int_t niter){
912 //
913 //
914 //
0eff8c53 915 //
7212f8c3 916 // Input:
917 // nhistos - maximal number of histograms to be used for sum
918 // nevents - number of events to make a fluctuation studies
919 // niter - number of itterations
920 // Algortihm:
0eff8c53 921 // 1. Make a summary integral 3D/2D/1D maps
922 // 2. Create several maps with niter events - Poisson flucturation in n
923 // 3. Store results 3D maps in the tree (and also as histogram) current and mean
7212f8c3 924 //
925
926 TFile * fhisto = TFile::Open("histo.root");
76dbc35e 927 TTree * tree = (TTree*)fhisto->Get("histo");
0eff8c53 928 tree->SetCacheSize(10000000000);
76dbc35e 929
7212f8c3 930 TTreeSRedirector * pcstream = new TTreeSRedirector("fluctuation.root", "update");
931
932
76dbc35e 933 TH1D * his1DROC=0, * his1DROCSum=0, * his1DROCN=0;
934 TH1D * his1DDrift=0, * his1DDriftSum=0, * his1DDriftN=0 ;
7212f8c3 935 TH2D * his2DRPhiROC=0, * his2DRPhiROCSum=0, * his2DRPhiROCN=0;
0eff8c53 936 TH2D * his2DRZROC=0, * his2DRZROCSum=0, * his2DRZROCN=0;
7212f8c3 937 TH2D * his2DRPhiDrift=0, * his2DRPhiDriftSum=0, * his2DRPhiDriftN=0;
0eff8c53 938 TH2D * his2DRZDrift=0, * his2DRZDriftSum=0, * his2DRZDriftN=0;
76dbc35e 939 TH3D * his3DROC=0, * his3DROCSum=0, * his3DROCN=0;
940 TH3D * his3DDrift=0, * his3DDriftSum=0, * his3DDriftN=0;
941 //
0eff8c53 942 if (nhistos<0 || nhistos> tree->GetEntries()) nhistos = tree->GetEntries();
943 Int_t eventsPerChunk=0;
76dbc35e 944 tree->SetBranchAddress("hist1D_2.",&his1DDrift);
945 tree->SetBranchAddress("hist1DROC_2.",&his1DROC);
7212f8c3 946 tree->SetBranchAddress("hist2DRPhi_2.",&his2DRPhiDrift);
0eff8c53 947 tree->SetBranchAddress("hist2DRZ_2.",&his2DRZDrift);
7212f8c3 948 tree->SetBranchAddress("hist2DRPhiROC_2.",&his2DRPhiROC);
76dbc35e 949 tree->SetBranchAddress("hist3D_2.",&his3DDrift);
950 tree->SetBranchAddress("hist3DROC_2.",&his3DROC);
0eff8c53 951 tree->SetBranchAddress("hist2DRZROC_2.",&his2DRZROC);
952 tree->SetBranchAddress("events",&eventsPerChunk);
953 //
954 // 1. Make a summary integral 3D/2D/1D maps
76dbc35e 955 //
0eff8c53 956 Int_t neventsAll=0;
76dbc35e 957 for (Int_t i=0; i<nhistos; i++){
958 tree->GetEntry(i);
959 if (i%25==0) printf("%d\n",i);
960 if (his1DROCSum==0) his1DROCSum=new TH1D(*his1DROC);
961 if (his1DDriftSum==0) his1DDriftSum=new TH1D(*his1DDrift);
7212f8c3 962 if (his2DRPhiROCSum==0) his2DRPhiROCSum=new TH2D(*his2DRPhiROC);
0eff8c53 963 if (his2DRZROCSum==0) his2DRZROCSum=new TH2D(*his2DRZROC);
7212f8c3 964 if (his2DRPhiDriftSum==0) his2DRPhiDriftSum=new TH2D(*his2DRPhiDrift);
0eff8c53 965 if (his2DRZDriftSum==0) his2DRZDriftSum=new TH2D(*his2DRZDrift);
76dbc35e 966 if (his3DROCSum==0) his3DROCSum=new TH3D(*his3DROC);
967 if (his3DDriftSum==0) his3DDriftSum=new TH3D(*his3DDrift);
968 his1DROCSum->Add(his1DROC);
969 his1DDriftSum->Add(his1DDrift);
7212f8c3 970 his2DRPhiROCSum->Add(his2DRPhiROC);
0eff8c53 971 his2DRZROCSum->Add(his2DRZROC);
7212f8c3 972 his2DRPhiDriftSum->Add(his2DRPhiDrift);
0eff8c53 973 his2DRZDriftSum->Add(his2DRZDrift);
76dbc35e 974 his3DROCSum->Add(his3DROC);
975 his3DDriftSum->Add(his3DDrift);
0eff8c53 976 neventsAll+=eventsPerChunk;
76dbc35e 977 }
978 //
0eff8c53 979 // 2. Create several maps with niter events - Poisson flucturation in n
980 //
981 for (Int_t iter=0; iter<niter; iter++){
982 printf("Itteration=\t%d\n",iter);
0eff8c53 983 Int_t nchunks=gRandom->Poisson(nevents)/eventsPerChunk; // chunks with n typically 25 events
984 for (Int_t i=0; i<nchunks; i++){
0eff8c53 985 tree->GetEntry(gRandom->Rndm()*nhistos);
986 if (i%10==0) printf("%d\t%d\n",iter, i);
987 if (his1DROCN==0) his1DROCN=new TH1D(*his1DROC);
988 if (his1DDriftN==0) his1DDriftN=new TH1D(*his1DDrift);
7212f8c3 989 if (his2DRPhiROCN==0) his2DRPhiROCN=new TH2D(*his2DRPhiROC);
990 if (his2DRPhiDriftN==0) his2DRPhiDriftN=new TH2D(*his2DRPhiDrift);
0eff8c53 991 if (his2DRZROCN==0) his2DRZROCN=new TH2D(*his2DRZROC);
992 if (his2DRZDriftN==0) his2DRZDriftN=new TH2D(*his2DRZDrift);
993 if (his3DROCN==0) his3DROCN=new TH3D(*his3DROC);
994 if (his3DDriftN==0) his3DDriftN=new TH3D(*his3DDrift);
995 his1DROCN->Add(his1DROC);
996 his1DDriftN->Add(his1DDrift);
7212f8c3 997 his2DRPhiROCN->Add(his2DRPhiROC);
0eff8c53 998 his2DRZDriftN->Add(his2DRZDrift);
999 his2DRZROCN->Add(his2DRZROC);
7212f8c3 1000 his2DRPhiDriftN->Add(his2DRPhiDrift);
0eff8c53 1001 his3DROCN->Add(his3DROC);
1002 his3DDriftN->Add(his3DDrift);
1003 }
1004 //
1005 // 3. Store results 3D maps in the tree (and also as histogram) current and mea
1006 //
1007 Int_t eventsUsed= nchunks*eventsPerChunk;
1008 (*pcstream)<<"fluctuation"<<
1009 "neventsAll="<<neventsAll<< // total number of event to define mean
1010 "nmean="<<nevents<< // mean number of events used
1011 "eventsUsed="<<eventsUsed<< // number of chunks used for one fluct. study
1012 //
1013 // 1,2,3D histogram per group and total
1014 "his1DROCN.="<<his1DROCN<<
1015 "his1DROCSum.="<<his1DROCSum<<
1016 "his1DDriftN.="<<his1DDriftN<<
1017 "his1DDriftSum.="<<his1DDriftSum<<
7212f8c3 1018 "his2DRPhiROCN.="<<his2DRPhiROCN<<
1019 "his2DRPhiROCSum.="<<his2DRPhiROCSum<<
1020 "his2DRPhiDriftN.="<<his2DRPhiDriftN<<
1021 "his2DRPhiDriftSum.="<<his2DRPhiDriftSum<<
0eff8c53 1022 "his2DRZROCN.="<<his2DRZROCN<<
1023 "his2DRZROCSum.="<<his2DRZROCSum<<
1024 "his2DRZDriftN.="<<his2DRZDriftN<<
1025 "his2DRZDriftSum.="<<his2DRZDriftSum<<
1026 "his3DROCN.="<<his3DROCN<<
1027 "his3DROCSum.="<<his3DROCSum<<
1028 "his3DDriftN.="<<his3DDriftN<<
1029 "his3DDriftSum.="<<his3DDriftSum<<
1030 "\n";
1031 pcstream->GetFile()->mkdir(Form("Fluc%d",iter));
1032 pcstream->GetFile()->cd(Form("Fluc%d",iter));
4e715689 1033 //
7212f8c3 1034 his2DRPhiROCN->Write("his2DRPhiROCN");
0eff8c53 1035 his2DRZROCN->Write("his2DRZROCN");
4e715689 1036 //
7212f8c3 1037 his2DRPhiROCSum->Write("his2DRPhiROCSum");
0eff8c53 1038 his2DRZROCSum->Write("his2DRZROCSum");
4e715689 1039 //
7212f8c3 1040 his2DRPhiDriftN->Write("his2DRPhiDriftN");
0eff8c53 1041 his2DRZDriftN->Write("his2DRZDriftN");
4e715689 1042 //
7212f8c3 1043 his2DRPhiDriftSum->Write("his2DRPhiDriftSum");
0eff8c53 1044 his2DRZDriftSum->Write("his2DRZDriftSum");
1045 //
1046 his3DROCN->Write("his3DROCN");
1047 his3DROCSum->Write("his3DROCSum");
1048 his3DDriftN->Write("his3DDriftN");
1049 his3DDriftSum->Write("his3DDriftSum");
1050
1051 his1DROCN->Reset();
1052 his1DDriftN->Reset();
7212f8c3 1053 his2DRPhiROCN->Reset();
0eff8c53 1054 his2DRZDriftN->Reset();
1055 his2DRZROCN->Reset();
7212f8c3 1056 his2DRPhiDriftN->Reset();
0eff8c53 1057 his3DROCN->Reset();
1058 his3DDriftN->Reset();
76dbc35e 1059 }
0eff8c53 1060
1061 delete pcstream;
76dbc35e 1062}
1063
1064
1065void DrawDCARPhiTrendTime(){
1066 //
0eff8c53 1067 // Macros to draw the DCA correlation with the luminosity (estimated from the occupancy)
76dbc35e 1068 //
0eff8c53 1069 // A side and c side 0 differnt behaviour -
1070 // A side - space charge effect
1071 // C side - space charge effect+ FC charging:
1072 // Variables to query from the QA/calibration DB - tree:
1073 // QA.TPC.CPass1.dcar_posA_0 -dca rphi in cm - offset
1074 // Calib.TPC.occQA.Sum() - luminosity is estimated using the mean occupancy per run
1075 //
76dbc35e 1076 TFile *fdb = TFile::Open("outAll.root");
0eff8c53 1077 if (!fdb) fdb = TFile::Open("http://www-alice.gsi.de/TPC/CPassMonitor/outAll.root");
76dbc35e 1078 TTree * tree = (TTree*)fdb->Get("joinAll");
1079 tree->SetCacheSize(100000000);
1080 tree->SetMarkerStyle(25);
1081
1082 //QA.TPC.CPass1.dcar_posA_0 QA.TPC.CPass1.dcar_posA_0_Err QA.TPC.CPass1.meanMult Calib.TPC.occQA. DAQ.L3_magnetCurrent
1083
1084 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);
1085 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);
1086 Double_t mean,rms;
1087 TStatToolkit::EvaluateUni(grA->GetN(),grA->GetY(), mean,rms,grA->GetN()*0.8);
1088 grA->SetMinimum(mean-5*rms);
1089 grA->SetMaximum(mean+3*rms);
1090
1091
1092 grA->GetXaxis()->SetTitle("occ*sign(bz)");
1093 grA->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
1094 grA->Draw("ap");
1095 grC->Draw("p");
1096 TLegend* legend = new TLegend(0.11,0.11,0.5,0.3,"DCA_{rphi} as function of IR (2013)" );
1097 legend->AddEntry(grA,"A side","p");
1098 legend->AddEntry(grC,"C side","p");
1099 legend->Draw();
76dbc35e 1100}
1101
1102
1103
1104void DrawOpenGate(){
1105 //
0eff8c53 1106 // Make nice plot to demonstrate the space charge effect in run with the open gating grid
1107 // For the moment the inmput is harwired - the CPass0 calibration data used
1108 // Make nice drawing (with axis labels):
1109 // To fix (longer term)
1110 // the distortion map to be recalculated - using gaussian fit (currently we use mean)
1111 // the histogram should be extended
76dbc35e 1112 TFile f("/hera/alice/alien/alice/data/2013/LHC13g/000197470/cpass0/OCDB/root_archive.zip#meanITSVertex.root");
1113 TFile fref("/hera/alice/alien/alice/data/2013/LHC13g/000197584/cpass0/OCDB/root_archive.zip#meanITSVertex.root");
1114 //
1115 TTree * treeTOFdy=(TTree*)f.Get("TOFdy");
1116 TTree * treeTOFdyRef=(TTree*)fref.Get("TOFdy");
1117 treeTOFdy->AddFriend(treeTOFdyRef,"R");
1118 treeTOFdy->SetMarkerStyle(25);
1119 TTree * treeITSdy=(TTree*)f.Get("ITSdy");
1120 TTree * treeITSdyRef=(TTree*)fref.Get("ITSdy");
1121 treeITSdy->AddFriend(treeITSdyRef,"R");
1122 treeITSdy->SetMarkerStyle(25);
1123 TTree * treeVertexdy=(TTree*)f.Get("Vertexdy");
1124 TTree * treeVertexdyRef=(TTree*)fref.Get("Vertexdy");
1125 treeVertexdy->AddFriend(treeVertexdyRef,"R");
1126 treeVertexdy->SetMarkerStyle(25);
0eff8c53 1127
1128 // treeITSdy->Draw("mean-R.mean:sector:abs(theta)","entries>50&&abs(snp)<0.1&&theta<0","colz")
76dbc35e 1129
0eff8c53 1130 treeITSdy->Draw("mean-R.mean:sector:abs(theta)","entries>50&&abs(snp)<0.1&&theta>0","colz");
c634db45 1131}
1132
0eff8c53 1133
ae4cea45 1134void DrawCurrent(const char * ocdb="/cvmfs/alice.gsi.de/alice/data/2013/OCDB/TPC/Calib/HighVoltage", Int_t run0=100000, Int_t run1=110000){
c634db45 1135 //
1136 //
1137 /*
1138 const char * ocdb="/cvmfs/alice.gsi.de/alice/data/2013/OCDB/TPC/Calib/HighVoltage";
1139 Int_t run0=197460;
1140 Int_t run1=197480;
1141 */
1142 const Int_t knpoints=100000;
1143 TVectorD vecTime(knpoints);
1144 TVectorD vecI(knpoints);
1145 Int_t npoints=0;
1146 for (Int_t irun=run0; irun<run1; irun++){
1147 TFile * f = TFile::Open(Form("%s/Run%d_%d_v1_s0.root",ocdb,irun,irun));
1148 if (!f) continue;
1149 AliCDBEntry * entry = (AliCDBEntry *)f->Get("AliCDBEntry");
1150 if (!entry) continue;
1151 AliDCSSensorArray * array = (AliDCSSensorArray *)entry->GetObject();
1152 if (!array) continue;
1153 AliDCSSensor * sensor = array->GetSensor("TPC_VHV_D_I_MON");
1154 //sensor->Draw(Form("%d",irun));
1155 TGraph *graph = sensor->GetGraph();
1156 for (Int_t ipoint=0; ipoint<graph->GetN(); ipoint++){
1157 vecTime[npoints]=sensor->GetStartTime()+graph->GetX()[ipoint]*3600;
1158 vecI[npoints]=graph->GetY()[ipoint];
1159 npoints++;
1160 }
1161 }
1162 TGraph * graph = new TGraph(npoints, vecTime.GetMatrixArray(), vecI.GetMatrixArray());
1163 graph->Draw("alp");
1164
76dbc35e 1165
1166}
4e715689 1167
1168
4b955de3 1169void MakeSpaceChargeFluctuationScan(Double_t scale, Int_t nfilesMerge, Int_t sign){
7212f8c3 1170 //
7212f8c3 1171 //
1172 // Input:
1173 // scale - scaling of the space charge
1174 // nfilesMerge - amount of chunks to merge
1175 // - =0 all chunks used
1176 // <0 subset form full statistic
1177 // >0 subset from the limited (1000 mean) stistic
4b955de3 1178 // Make fluctuation scan:
1179 // 1.) Shift of z disk - to show which granularity in time needed
1180 // 2.) Shift in sector - to show influence of the gass gain and epsilon
1181 // 3.) Smearing in phi - to define phi granularity needed
1182 // 4.) Rebin z
1183 // 5.) Rebin phi
1184 // For given SC setups the distortion on the space point and track level characterezed
1185 // SpaceChargeFluc%d_%d.root
1186 // SpaceChargeTrackFluc%d%d.root
1187 //
1188
1189
1190 //
1191 // Some constant definition
1192 //
1193 Int_t nitteration=100; // number of itteration in the lookup
1194 Int_t fullNorm =10000; // normalization fro the full statistic
1195
1196 //
1197 // Init magnetic field and OCDB
1198 //
1199
1200 Double_t bsign= sign;
1201 if (bsign>1) bsign=-1;
2ff913b6 1202 const Int_t nTracks=2000;
4b955de3 1203 const char *ocdb="local://$ALICE_ROOT/OCDB/";
1204 AliCDBManager::Instance()->SetDefaultStorage(ocdb);
1205 AliCDBManager::Instance()->SetRun(0);
1206 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", bsign, bsign, AliMagF::k5kG));
1207 //
1208
1209
1210 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("SpaceChargeFluc%d_%d.root",nfilesMerge,sign),"recreate");
1211 TTreeSRedirector *pcstreamTrack = new TTreeSRedirector(Form("SpaceChargeTrackFluc%d_%d.root",nfilesMerge,sign),"recreate");
7212f8c3 1212 TH1D *his1DROCN=0, *his1DROCSum=0;
1213 TH2D *his2DRPhiROCN=0, *his2DRPhiROCSum=0, *his2DRZROCN=0, *his2DRZROCSum=0;
1214 TH3D *his3DROCN=0, *his3DROCSum=0;
1215 TH1D *his1DROCNC=0, *his1DROCSumC=0;
1216 TH2D *his2DRPhiROCNC=0, *his2DRPhiROCSumC=0, *his2DRZROCNC=0, *his2DRZROCSumC=0;
1217 TH3D *his3DROCNC=0, *his3DROCSumC=0;
1218 TH1 * histos[8]={his1DROCN, his1DROCSum, his2DRPhiROCN, his2DRPhiROCSum, his2DRZROCN, his2DRZROCSum, his3DROCN, his3DROCSum};
1219 Int_t neventsAll=0, neventsAllC=0;
1220 Int_t neventsChunk=0, neventsChunkC=0;
1221 const Double_t ePerADC = 500.;
1222 const Double_t fgke0 = 8.854187817e-12;
1223 //
1224 //
1225 //
1226 const char *inputFile="fluctuation.root";
1227 TObjArray * fileList = (gSystem->GetFromPipe("cat fluctuation.list")).Tokenize("\n");
1228 if (fileList->GetEntries()==0) fileList->AddLast(new TObjString(inputFile));
1229 Int_t nfiles = fileList->GetEntries();
1230 Int_t indexPer[1000];
1231 Double_t numbersPer[10000];
1232 for (Int_t i=0; i<nfiles; i++) numbersPer[i]=gRandom->Rndm();
1233 TMath::Sort(nfiles, numbersPer,indexPer);
1234
1235 for (Int_t ifile=0; ifile<nfiles; ifile++){
1236 if (nfilesMerge>0 && ifile>=nfilesMerge) continue; // merge only limited amount if specified by argument
1237 TFile *fhistos = TFile::Open(fileList->At(indexPer[ifile])->GetName());
1238 if (!fhistos) continue;
1239 TTree * treeHis = (TTree*)fhistos->Get("fluctuation");
1240 if (!treeHis) { printf("file %s does not exist or tree does not exist\n",fileList->At(ifile)->GetName()); continue;}
1241 Int_t nchunks=treeHis->GetEntries();
1242 Int_t chunk=nchunks*gRandom->Rndm();
1243 treeHis->SetBranchAddress("his1DROCN.",&his1DROCNC);
1244 treeHis->SetBranchAddress("his1DROCSum.",&his1DROCSumC);
1245 treeHis->SetBranchAddress("his2DRPhiROCN.",&his2DRPhiROCNC);
1246 treeHis->SetBranchAddress("his2DRPhiROCSum.",&his2DRPhiROCSumC);
1247 treeHis->SetBranchAddress("his2DRZROCN.",&his2DRZROCNC);
1248 treeHis->SetBranchAddress("his2DRZROCSum.",&his2DRZROCSumC);
1249 treeHis->SetBranchAddress("his3DROCN.",&his3DROCNC);
1250 treeHis->SetBranchAddress("his3DROCSum.",&his3DROCSumC);
1251 treeHis->SetBranchAddress("neventsAll",&neventsAllC);
1252 treeHis->SetBranchAddress("eventsUsed",&neventsChunkC);
1253 treeHis->GetEntry(chunk);
1254 neventsAll+=neventsAllC;
1255 neventsChunk+=neventsChunkC;
1256 //
1257 TH1 * histosC[8]={ his1DROCNC, his1DROCSumC, his2DRPhiROCNC, his2DRPhiROCSumC, his2DRZROCNC, his2DRZROCSumC, his3DROCNC, his3DROCSumC};
1258 if (ifile==0) for (Int_t ihis=0; ihis<8; ihis++) histos[ihis] = (TH1*)(histosC[ihis]->Clone());
1259 if (ifile>0) {
1260 for (Int_t ihis=0; ihis<8; ihis++) histos[ihis]->Add(histosC[ihis]);
1261 }
1262 }
1263 his1DROCN=(TH1D*)histos[0]; his1DROCSum=(TH1D*)histos[1];
1264 his2DRPhiROCN=(TH2D*)histos[2]; his2DRPhiROCSum=(TH2D*)histos[3]; his2DRZROCN=(TH2D*)histos[4]; his2DRZROCSum=(TH2D*)histos[5];
1265 his3DROCN=(TH3D*)histos[6]; his3DROCSum=(TH3D*)histos[7];
1266 //
1267 // Select input histogram
1268 //
1269 TH3D * hisInput= his3DROCSum;
1270 Int_t neventsCorr=0; // number of events used for the correction studies
1271 if (nfilesMerge>0){
1272 neventsCorr=neventsChunk;
1273 hisInput=his3DROCN;
1274 }else{
1275 neventsCorr=neventsAll;
1276 hisInput=his3DROCSum;
4b955de3 1277 hisInput->Scale(Double_t(fullNorm)/Double_t(neventsAll));
7212f8c3 1278 }
1279
1280 TObjArray *distortionArray = new TObjArray;
4b955de3 1281 TObjArray *histoArray = new TObjArray;
7212f8c3 1282 //
1283 // Make a reference - ideal distortion/correction
1284 //
1285 TH3D * his3DReference = NormalizeHistoQ(hisInput,kFALSE); // q normalized to the Q/m^3
1286 his3DReference->Scale(scale*0.000001/fgke0); //scale back to the C/cm^3/epsilon0
1287 AliTPCSpaceCharge3D *spaceChargeRef = new AliTPCSpaceCharge3D;
4b955de3 1288 spaceChargeRef->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
7212f8c3 1289 spaceChargeRef->SetInputSpaceCharge(his3DReference, his2DRPhiROCSum,his2DRPhiROCSum,1);
1290 spaceChargeRef->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
7212f8c3 1291 spaceChargeRef->AddVisualCorrection(spaceChargeRef,1);
1292 spaceChargeRef->SetName("DistRef");
4b955de3 1293 his3DReference->SetName("hisDistRef");
7212f8c3 1294 distortionArray->AddLast(spaceChargeRef);
4b955de3 1295 histoArray->AddLast(his3DReference);
1296 //
1297 // Draw histos
1298 TCanvas * canvasSC = new TCanvas("canvasSCDefault","canvasSCdefault",500,400);
1299 canvasSC->SetRightMargin(0.12);
1300 gStyle->SetTitleOffset(0.8,"z");
1301 canvasSC->SetRightMargin(0.13);
1302 spaceChargeRef->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1303 canvasSC->SaveAs(Form("canvasCreateHistoDRPhiinXY_Z10_%d_%d.pdf",nfilesMerge,sign));
1304 spaceChargeRef->CreateHistoDRinXY(10,250,250)->Draw("colz");
1305 canvasSC->SaveAs(Form("canvasCreateHistoDRinXY_Z10_%d_%d.pdf",nfilesMerge,sign));
1306 spaceChargeRef->CreateHistoSCinZR(0.05,250,250)->Draw("colz");
1307 canvasSC->SaveAs(Form("canvasCreateHistoSCinZR_Phi005_%d_%d.pdf",nfilesMerge,sign));
1308 spaceChargeRef->CreateHistoSCinXY(10.,250,250)->Draw("colz");
1309 canvasSC->SaveAs(Form("canvasCreateHistoSCinRPhi_Z10_%d_%d.pdf",nfilesMerge,sign));
1310
1311
7212f8c3 1312 //
1313 // Make Z scan corrections
1314 //
4b955de3 1315 if (1){
1316 for (Int_t ihis=1; ihis<=9; ihis+=2){
1317 TH3 *his3DZ = PermutationHistoZ(his3DReference,16*(ihis));
7212f8c3 1318 AliTPCSpaceCharge3D *spaceChargeZ = new AliTPCSpaceCharge3D;
4b955de3 1319 spaceChargeZ->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
7212f8c3 1320 spaceChargeZ->SetInputSpaceCharge(his3DZ, his2DRPhiROCSum,his2DRPhiROCSum,1);
1321 spaceChargeZ->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
1322 spaceChargeZ->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1323 spaceChargeZ->AddVisualCorrection(spaceChargeZ,100+ihis);
4b955de3 1324 spaceChargeZ->SetName(Form("DistZ_%d", 16*(ihis)));
1325 his3DZ->SetName(Form("HisDistZ_%d", 16*(ihis)));
7212f8c3 1326 distortionArray->AddLast(spaceChargeZ);
4b955de3 1327 histoArray->AddLast(his3DZ);
7212f8c3 1328 }
1329 //
1330 // Make Sector scan corrections
1331 //
1332 for (Int_t ihis=1; ihis<=9; ihis+=2){
1333 TH3 *his3DSector = PermutationHistoPhi(his3DReference,TMath::Pi()*(ihis)/9.);
1334 AliTPCSpaceCharge3D *spaceChargeSector = new AliTPCSpaceCharge3D;
4b955de3 1335 spaceChargeSector->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
7212f8c3 1336 spaceChargeSector->SetInputSpaceCharge(his3DSector, his2DRPhiROCSum,his2DRPhiROCSum,1);
1337 spaceChargeSector->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
1338 spaceChargeSector->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1339 spaceChargeSector->AddVisualCorrection(spaceChargeSector,200+ihis);
1340 spaceChargeSector->SetName(Form("DistSector_%d", ihis));
4b955de3 1341 his3DSector->SetName(Form("DistSector_%d", ihis));
7212f8c3 1342 distortionArray->AddLast(spaceChargeSector);
4b955de3 1343 histoArray->AddLast(his3DSector);
1344 }
7212f8c3 1345 //
4b955de3 1346 // Make Local phi scan smear corrections
1347 //
1348 for (Int_t ihis=1; ihis<=8; ihis++){
1349 TH3 *his3DSector = PermutationHistoLocalPhi(his3DReference,ihis);
1350 AliTPCSpaceCharge3D *spaceChargeSector = new AliTPCSpaceCharge3D;
1351 spaceChargeSector->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
1352 spaceChargeSector->SetInputSpaceCharge(his3DSector, his2DRPhiROCSum,his2DRPhiROCSum,1);
1353 spaceChargeSector->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
1354 spaceChargeSector->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1355 spaceChargeSector->AddVisualCorrection(spaceChargeSector,300+ihis);
1356 spaceChargeSector->SetName(Form("DistPhi_%d", ihis));
1357 his3DSector->SetName(Form("HisDistPhi_%d", ihis));
1358 distortionArray->AddLast(spaceChargeSector);
1359 histoArray->AddLast(his3DSector);
1360 }
1361 // //
1362// // Rebin Z
1363// //
1364// for (Int_t ihis=2; ihis<=8; ihis+=2){
1365// TH3 *his3DSector = his3DReference->RebinZ(ihis,Form("RebinZ_%d",ihis));
1366// AliTPCSpaceCharge3D *spaceChargeSector = new AliTPCSpaceCharge3D;
1367// spaceChargeSector->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
1368// spaceChargeSector->SetInputSpaceCharge(his3DSector, his2DRPhiROCSum,his2DRPhiROCSum,1);
1369// spaceChargeSector->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
1370// spaceChargeSector->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1371// spaceChargeSector->AddVisualCorrection(spaceChargeSector,300+ihis);
1372// spaceChargeSector->SetName(Form("RebinZ_%d", ihis));
1373// his3DSector->SetName(Form("RebinZ_%d", ihis));
1374// distortionArray->AddLast(spaceChargeSector);
1375// histoArray->AddLast(his3DSector);
1376// }
1377// //
1378// // Rebin Phi
1379// //
1380// for (Int_t ihis=2; ihis<=5; ihis++){
1381// TH3 *his3DSector = his3DReference->RebinZ(ihis,Form("RebinPhi_%d",ihis));
1382// AliTPCSpaceCharge3D *spaceChargeSector = new AliTPCSpaceCharge3D;
1383// spaceChargeSector->SetOmegaTauT1T2(omegaTau*bsign,1,1); // Ne CO2
1384// spaceChargeSector->SetInputSpaceCharge(his3DSector, his2DRPhiROCSum,his2DRPhiROCSum,1);
1385// spaceChargeSector->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
1386// spaceChargeSector->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1387// spaceChargeSector->AddVisualCorrection(spaceChargeSector,300+ihis);
1388// spaceChargeSector->SetName(Form("RebinZ_%d", ihis));
1389// his3DSector->SetName(Form("RebinZ_%d", ihis));
1390// distortionArray->AddLast(spaceChargeSector);
1391// histoArray->AddLast(his3DSector);
1392// }
1393 }
7212f8c3 1394 //
4b955de3 1395 // Space points scan
7212f8c3 1396 //
1397 Int_t nx = his3DROCN->GetXaxis()->GetNbins();
1398 Int_t ny = his3DROCN->GetYaxis()->GetNbins();
1399 Int_t nz = his3DROCN->GetZaxis()->GetNbins();
1400 Int_t nbins=nx*ny*nz;
1401 TVectorF vx(nbins), vy(nbins), vz(nbins), vq(nbins), vqall(nbins);
1402 //
1403 // charge in the ROC
1404 // for open gate data only fraction of ions enter to drift volume
1405 //
1406 const Int_t kbins=1000;
4b955de3 1407 Double_t deltaR[kbins], deltaZ[kbins],deltaRPhi[kbins], deltaQ[kbins];
7212f8c3 1408 Int_t ndist = distortionArray->GetEntries();
1409 for (Int_t ix=1; ix<=nx; ix+=2){ // phi bin loop
1410 for (Int_t iy=1; iy<=ny; iy+=2){ // r bin loop
1411 Double_t phi= his3DROCN->GetXaxis()->GetBinCenter(ix);
1412 Double_t r = his3DROCN->GetYaxis()->GetBinCenter(iy);
1413 Double_t x = r*TMath::Cos(phi);
1414 Double_t y = r*TMath::Sin(phi);
1415 //
1416 for (Int_t iz=1; iz<=nz; iz++){ // z bin loop
1417 Double_t z = his3DROCN->GetZaxis()->GetBinCenter(iz);
1418 Double_t qN= his3DROCN->GetBinContent(ix,iy,iz);
1419 Double_t qSum= his3DROCSum->GetBinContent(ix,iy,iz);
1420 // Double_t dV in cm = dphi * r * dz in cm**3
1421 Double_t dV= (his3DROCN->GetXaxis()->GetBinWidth(ix)*r)*his3DROCN->GetZaxis()->GetBinWidth(iz);
1422 Double_t norm= 1e6*ePerADC*TMath::Qe()/dV; //normalization factor to the Q/m^3 inside of the ROC;
1423 (*pcstream)<<"hisDump"<<
1424 "neventsAll="<<neventsAll<< // total number of events used for the Q reference
1425 "nfiles="<<nfiles<< // number of files to define properties
1426 "nfilesMerge="<<nfilesMerge<< // number of files to define propertiesneventsCorr
1427 "neventsCorr="<<neventsCorr<< // number of events used to define the corection
4b955de3 1428 "fullNorm="<<fullNorm<< // in case full statistic used this is the normalization coeficient
7212f8c3 1429
1430 "ix="<<ix<<
1431 "iy="<<iy<<
1432 "iz="<<iz<<
1433 // x,y,z
1434 "x="<<x<<
1435 "y="<<y<<
1436 "z="<<z<<
1437 // phi,r,z
1438 "phi="<<phi<<
1439 "r="<<r<<
1440 "z="<<z<<
1441 "norm="<<norm<<
1442 "qN="<<qN<<
1443 "qSum="<<qSum;
1444 for (Int_t idist=0; idist<ndist; idist++){
1445 AliTPCCorrection * corr = (AliTPCCorrection *)distortionArray->At(idist);
4b955de3 1446 TH3 * his = (TH3*)histoArray->At(idist);
7212f8c3 1447 Double_t phi0= TMath::ATan2(y,x);
1448 Int_t nsector=(z>=0) ? 0:18;
1449 Float_t distPoint[3]={x,y,z};
1450 corr->CorrectPoint(distPoint, nsector);
1451 Double_t r0=TMath::Sqrt(x*x+y*y);
1452 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1453 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1454 deltaR[idist] = r1-r0;
1455 deltaRPhi[idist] = (phi1-phi0)*r0;
1456 deltaZ[idist] = distPoint[2]-z;
4b955de3 1457 deltaQ[idist] = his->GetBinContent(ix,iy,iz);
7212f8c3 1458 //
1459 (*pcstream)<<"hisDump"<< //correct point - input point
4b955de3 1460 Form("%sQ=",corr->GetName())<<deltaQ[idist]<<
7212f8c3 1461 Form("%sDR=",corr->GetName())<<deltaR[idist]<<
1462 Form("%sDRPhi=",corr->GetName())<<deltaRPhi[idist]<<
1463 Form("%sDZ=",corr->GetName())<<deltaZ[idist];
1464 }
1465 (*pcstream)<<"hisDump"<<
1466 "\n";
1467 }
1468 }
1469 }
1470 delete pcstream;
4b955de3 1471 //
1472 // generate track distortions
1473 //
1474 if (nTracks>0){
1475 for(Int_t nt=1; nt<=nTracks; nt++){
1476 gRandom->SetSeed(nt);
1477 TObjArray trackArray(10000);
1478 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
1479 Double_t eta = gRandom->Uniform(-1, 1);
1480 Double_t pt = 1/(gRandom->Rndm()*5+0.00001); // momentum for f1
1481 Short_t psign=1;
1482 if(gRandom->Rndm() < 0.5){
1483 psign =1;
1484 }else{
1485 psign=-1;
1486 }
1487 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1488 Double_t pxyz[3];
1489 pxyz[0]=pt*TMath::Cos(phi);
1490 pxyz[1]=pt*TMath::Sin(phi);
1491 pxyz[2]=pt*TMath::Tan(theta);
1492 Double_t vertex[3]={0,0,0};
1493 Double_t cv[21]={0};
1494 AliExternalTrackParam *t= new AliExternalTrackParam(vertex, pxyz, cv, psign);
1495 Double_t refX0=85.;
1496 Double_t refX1=1.;
1497 Int_t dir=-1;
1498 (*pcstreamTrack)<<"trackFit"<<
1499 "neventsAll="<<neventsAll<< // total number of events used for the Q reference
1500 "nfiles="<<nfiles<< // number of files to define properties
1501 "nfilesMerge="<<nfilesMerge<< // number of files to define propertiesneventsCorr
1502 "neventsCorr="<<neventsCorr<< // number of events used to define the corection
1503 "fullNorm="<<fullNorm<< // in case full statistic used this is the normalization coeficient
1504 "itrack="<<nt<<
1505 "input.="<<t;
1506 for (Int_t idist=0; idist<ndist; idist++){
1507 AliTPCCorrection * corr = (AliTPCCorrection *)distortionArray->At(idist);
1508 AliExternalTrackParam *ot0= new AliExternalTrackParam(vertex, pxyz, cv, psign);
1509 AliExternalTrackParam *ot1= new AliExternalTrackParam(vertex, pxyz, cv, psign);
1510 AliExternalTrackParam *td0 = corr->FitDistortedTrack(*ot0, refX0, dir, 0);
1511 AliExternalTrackParam *td1 = corr->FitDistortedTrack(*ot1, refX1, dir, 0);
1512 trackArray.AddLast(td0);
1513 trackArray.AddLast(td1);
1514 trackArray.AddLast(ot0);
1515 trackArray.AddLast(ot1);
1516 char name0[100], name1[1000];
1517 char oname0[100], oname1[1000];
1518 snprintf(name0, 100, "T_%s_0.=",corr->GetName());
1519 snprintf(name1, 100, "T_%s_1.=",corr->GetName());
1520 snprintf(oname0, 100, "OT_%s_0.=",corr->GetName());
1521 snprintf(oname1, 100, "OT_%s_1.=",corr->GetName());
1522 (*pcstreamTrack)<<"trackFit"<<
1523 name0<<td0<<
1524 name1<<td1<<
1525 oname0<<ot0<<
1526 oname1<<ot1;
1527 }
1528 (*pcstreamTrack)<<"trackFit"<<"\n";
1529 }
1530 }
1531 delete pcstreamTrack;
7212f8c3 1532 return;
1533
1534}
1535
1536
1537void MakePlotPoisson3D(const char *inputfile="fluctuation.root", const char *outputfile="SpaceCharge.root", Int_t event=0){
4e715689 1538 //
7212f8c3 1539 // draw "standard" plot to show radial and theta dependence of the space charge distortion
4e715689 1540 //
7212f8c3 1541 // const char *inputfile="fluctuation.root"; const char *outputfile="SpaceCharge.root"; Int_t event=0
4e715689 1542 //
1543 TFile *fhistos = TFile::Open(inputfile);
7212f8c3 1544 TH2D *his2DRPhiROCN=0, *his2DRPhiROCSum=0, *his2DRZROCN=0, *his2DRZROCSum=0;
1545 TH1D *his1DROCN=0, *his1DROCSum=0;
1546 TH3D *his3DROCN=0, *his3DROCSum=0;
1547 const Double_t ePerADC = 500.;
1548 const Double_t fgke0 = 8.854187817e-12;
4e715689 1549 TTree * treeHis = (TTree*)fhistos->Get("fluctuation");
7212f8c3 1550 treeHis->SetBranchAddress("his1DROCN.",&his1DROCN);
1551 treeHis->SetBranchAddress("his1DROCSum.",&his1DROCSum);
1552 treeHis->SetBranchAddress("his2DRPhiROCN.",&his2DRPhiROCN);
1553 treeHis->SetBranchAddress("his2DRPhiROCSum.",&his2DRPhiROCSum);
4e715689 1554 treeHis->SetBranchAddress("his2DRZROCN.",&his2DRZROCN);
1555 treeHis->SetBranchAddress("his2DRZROCSum.",&his2DRZROCSum);
7212f8c3 1556 treeHis->SetBranchAddress("his3DROCN.",&his3DROCN);
1557 treeHis->SetBranchAddress("his3DROCSum.",&his3DROCSum);
4e715689 1558 treeHis->GetEntry(event);
7212f8c3 1559
1560 his3DROCSum->Scale(ePerADC*TMath::Qe()/fgke0);
1561
1562 AliTPCSpaceCharge3D *spaceChargeOrig = new AliTPCSpaceCharge3D;
1563 spaceChargeOrig->SetOmegaTauT1T2(0.0,1,1); // Ne CO2
1564 spaceChargeOrig->SetInputSpaceCharge(his3DROCSum, his2DRPhiROCSum,his2DRPhiROCSum,10*ePerADC*TMath::Qe());
1565 spaceChargeOrig->InitSpaceCharge3DPoisson(129, 129, 144,100);
1566 spaceChargeOrig->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1567 spaceChargeOrig->AddVisualCorrection(spaceChargeOrig,1);
4e715689 1568 //
4b955de3 1569 //AliTPCSpaceCharge3D *spaceChargeRef= spaceChargeOrig;
1570
1571
7212f8c3 1572
4e715689 1573 //
7212f8c3 1574 Int_t nfuns=5;
1575 Double_t dmax=0.75, dmin=-0.75;
1576 Double_t phiRange=18;
1577 TCanvas *canvasDistortionP3D = new TCanvas("canvasdistortionP3D","canvasdistortionP3D",1000,700);
1578 canvasDistortionP3D->SetGrid(1,1);
1579 canvasDistortionP3D->Divide(1,2);
1580 canvasDistortionP3D->cd(1)->SetGrid(1,1);
1581 TLegend * legendR= new TLegend(0.11,0.11,0.45,0.35,"R scan (#Theta=0.1)");
1582 for (Int_t ifun1=0; ifun1<=nfuns; ifun1++){
1583 Double_t rfun= 85.+ifun1*(245.-85.)/nfuns;
1584 TF1 *pf1 = new TF1("f1",Form("AliTPCCorrection::GetCorrSector(x,%f,0.1,1,1)",rfun),0,phiRange);
1585 pf1->SetMinimum(dmin);
1586 pf1->SetMaximum(dmax);
1587 pf1->SetNpx(360);
1588 pf1->SetLineColor(1+ifun1);
1589 pf1->SetLineWidth(2);
1590 pf1->GetXaxis()->SetTitle("sector");
1591 pf1->GetXaxis()->SetNdivisions(530);
1592 pf1->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
1593 if (ifun1==0) pf1->Draw();
1594 pf1->Draw("same");
1595 legendR->AddEntry(pf1,Form("r=%1.0f",rfun));
1596 }
1597 legendR->Draw();
4e715689 1598 //
7212f8c3 1599 canvasDistortionP3D->cd(2)->SetGrid(1,1);
1600 TLegend * legendTheta= new TLegend(0.11,0.11,0.45,0.35,"#Theta scan (r=125 cm)");
1601 for (Int_t ifun1=0; ifun1<=nfuns; ifun1++){
1602 Double_t tfun= 0.1+ifun1*(0.8)/nfuns;
1603 TF1 *pf1 = new TF1("f1",Form("AliTPCCorrection::GetCorrSector(x,125,%f,1,1)",tfun),0,phiRange);
1604 pf1->SetMinimum(dmin);
1605 pf1->SetMaximum(dmax);
1606 pf1->SetNpx(360);
1607 pf1->SetLineColor(1+ifun1);
1608 pf1->SetLineWidth(2);
1609 pf1->GetXaxis()->SetTitle("sector");
1610 pf1->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
1611 pf1->GetXaxis()->SetNdivisions(530);
1612 if (ifun1==0) pf1->Draw();
1613 pf1->Draw("same");
1614 legendTheta->AddEntry(pf1,Form("#Theta=%1.2f",tfun));
1615 }
1616 legendTheta->Draw();
4e715689 1617
7212f8c3 1618}
1619
1620TH3D * NormalizeHistoQ(TH3D * hisInput, Bool_t normEpsilon){
1621 //
1622 // Renormalize the histogram to the Q/m^3
1623 // Input:
1624 // hisInput - input 3D histogram
1625 // normEpsilon - flag - normalize to epsilon0
1626 //
1627 const Double_t ePerADC = 500.;
1628 const Double_t fgkEpsilon0 = 8.854187817e-12;
1629 TH3D * hisOutput= new TH3D(*hisInput);
1630 Int_t nx = hisInput->GetXaxis()->GetNbins();
1631 Int_t ny = hisInput->GetYaxis()->GetNbins();
1632 Int_t nz = hisInput->GetZaxis()->GetNbins();
1633 for (Int_t ix=1; ix<=nx; ix++){
1634 for (Int_t iy=1; iy<=ny; iy++){
1635 for (Int_t iz=1; iz<=nz; iz++){
1636 // Double_t z = hisInput->GetZaxis()->GetBinCenter(iz);
1637 Double_t deltaRPhi = hisInput->GetXaxis()->GetBinWidth(ix)* hisInput->GetYaxis()->GetBinCenter(iy);
1638 Double_t deltaR= hisInput->GetYaxis()->GetBinWidth(iy);
1639 Double_t deltaZ= hisInput->GetYaxis()->GetBinWidth(iz);
1640 Double_t volume= (deltaRPhi*deltaR*deltaZ)/1000000.;
1641 Double_t q = hisInput->GetBinContent(ix,iy,iz)* ePerADC*TMath::Qe(); // Q in coulombs
1642 Double_t rho = q/volume; // rpho - density in Q/m^3
1643 if (normEpsilon) rho/=fgkEpsilon0;
1644 hisOutput->SetBinContent(ix,iy,iz,rho);
4e715689 1645 }
1646 }
1647 }
7212f8c3 1648 return hisOutput;
1649}
1650
1651
1652
1653TH3D * PermutationHistoZ(TH3D * hisInput, Double_t deltaZ){
1654 //
1655 // Used to estimate the effect of the imperfection of the lookup tables as function of update frequency
4e715689 1656 //
7212f8c3 1657 // Permute/rotate the conten of the histogram in z direction
1658 // Reshufle/shift content - Keeping the integral the same
1659 // Parameters:
1660 // hisInput - input 3D histogram (phi,r,z)
1661 // deltaZ - deltaZ -shift of the space charge
1662 Double_t zmax=250;
1663 TH3D * hisOutput= new TH3D(*hisInput);
1664 Int_t nx = hisInput->GetXaxis()->GetNbins();
1665 Int_t ny = hisInput->GetYaxis()->GetNbins();
1666 Int_t nz = hisInput->GetZaxis()->GetNbins();
4e715689 1667 //
7212f8c3 1668 //
1669 for (Int_t ix=1; ix<=nx; ix++){
1670 for (Int_t iy=1; iy<=ny; iy++){
1671 for (Int_t iz=1; iz<=nz; iz++){
1672 Double_t zold = hisInput->GetZaxis()->GetBinCenter(iz);
1673 Double_t z=zold;
1674 if (z>0){
1675 z+=deltaZ;
1676 if (z<0) z+=zmax;
1677 if (z>zmax) z-=zmax;
1678 }else{
1679 z-=deltaZ;
1680 if (z>0) z-=zmax;
1681 if (z<-zmax) z+=zmax; }
1682 Double_t kz= hisInput->GetZaxis()->FindBin(z);
1683 Double_t content = hisInput->GetBinContent(ix,iy,iz);
1684 hisOutput->SetBinContent(ix,iy,kz,content);
4e715689 1685 }
1686 }
1687 }
7212f8c3 1688 return hisOutput;
1689}
1690
4b955de3 1691
1692
1693
7212f8c3 1694TH3D * PermutationHistoPhi(TH3D * hisInput, Double_t deltaPhi){
1695 //
1696 // Used to estimate the effect of the imperfection of the lookup tables as function of update frequency
1697 //
1698 // Permute/rotate the conten of the histogram in phi
1699 // Reshufle/shift content - Keeping the integral the same
1700 // Parameters:
1701 // hisInput - input 3D histogram (phi,r,z)
1702 // deltaPhi - deltaPhi -shift of the space charge
1703 TH3D * hisOutput= new TH3D(*hisInput);
1704 Int_t nx = hisInput->GetXaxis()->GetNbins();
1705 Int_t ny = hisInput->GetYaxis()->GetNbins();
1706 Int_t nz = hisInput->GetZaxis()->GetNbins();
4e715689 1707 //
4e715689 1708 //
7212f8c3 1709 for (Int_t iy=1; iy<=ny; iy++){
1710 for (Int_t iz=1; iz<=nz; iz++){
1711 for (Int_t ix=1; ix<=nx; ix++){
1712 Double_t phiOld = hisInput->GetXaxis()->GetBinCenter(ix);
1713 Double_t phi=phiOld;
1714 phi+=deltaPhi;
1715 if (phi<0) phi+=TMath::TwoPi();
1716 if (phi>TMath::TwoPi()) phi-=TMath::TwoPi();
1717 Double_t kx= hisInput->GetXaxis()->FindBin(phi);
1718 Double_t content = hisInput->GetBinContent(ix,iy,iz);
1719 hisOutput->SetBinContent(kx,iy,iz,content);
1720 }
4e715689 1721 }
1722 }
7212f8c3 1723 return hisOutput;
1724}
4e715689 1725
7212f8c3 1726
4b955de3 1727TH3D * PermutationHistoLocalPhi(TH3D * hisInput, Int_t deltaPhi){
7212f8c3 1728 //
1729 // Used to estimate the effect of the imperfection of the lookup tables as function of update frequency
4b955de3 1730 // Use moving average of the content instead of the content
7212f8c3 1731 //
7212f8c3 1732 // Parameters:
1733 // hisInput - input 3D histogram (phi,r,z)
4b955de3 1734 // deltaPhi - moving average width
7212f8c3 1735 TH3D * hisOutput= new TH3D(*hisInput);
1736 Int_t nx = hisInput->GetXaxis()->GetNbins();
1737 Int_t ny = hisInput->GetYaxis()->GetNbins();
1738 Int_t nz = hisInput->GetZaxis()->GetNbins();
4b955de3 1739 Int_t binSector=nx/18;
7212f8c3 1740 //
1741 //
1742 for (Int_t iy=1; iy<=ny; iy++){
1743 for (Int_t iz=1; iz<=nz; iz++){
1744 for (Int_t ix=1; ix<=nx; ix++){
4b955de3 1745 Double_t sumRo=0,sumW=0;
1746 for (Int_t idx=-deltaPhi; idx<=deltaPhi; idx++){
1747 Int_t index=ix+idx;
1748 if (index<1) index+=nx+1; // underflow and overflow bins
1749 if (index>nx) index-=nx+1;
1750 Double_t content = hisInput->GetBinContent(index,iy,iz);
1751 sumRo+=content;
1752 sumW++;
1753 }
1754 Double_t meanCont= sumRo/sumW;
1755 hisOutput->SetBinContent(ix,iy,iz,meanCont);
1756 //printf("%d\t%f\n",ix,hisInput->GetBinContent(ix,iy,iz)/(hisInput->GetBinContent(ix,iy,iz)+meanCont));
1757 }
7212f8c3 1758 }
1759 }
1760 return hisOutput;
4e715689 1761}
1762
1763
4e715689 1764
7212f8c3 1765void ScanIterrationPrecision(TH3 * hisInput, Int_t offset){
1766 //
1767 //
1768 //
1769 for (Int_t iter=0; iter<=7; iter++){
1770 Int_t niter= 50.*TMath::Power(1.5,iter);
1771 AliTPCSpaceCharge3D *spaceChargeOrig = new AliTPCSpaceCharge3D;
1772 spaceChargeOrig->SetOmegaTauT1T2(0.0,1,1); // Ne CO2
1773 spaceChargeOrig->SetInputSpaceCharge(hisInput,0,0,1);
1774 spaceChargeOrig->InitSpaceCharge3DPoisson(129, 129, 144,niter);
1775 spaceChargeOrig->CreateHistoDRPhiinXY(10,250,250)->Draw("colz");
1776 spaceChargeOrig->AddVisualCorrection(spaceChargeOrig,offset+iter+1);
1777 }
1778}
1779
1780
1781void DrawFluctuationSector(Int_t stat, Double_t norm){
4e715689 1782 //
7212f8c3 1783 // Draw correction - correction at rotated sector
1784 // The same set of events used
1785 // Int_t stat=0; Double_t norm=10000;
4b955de3 1786 //
1787 // Notes:
1788 // 1. (something wrong for the setup 2 pileups -problem with data 24.07
1789 //
1790 //
7212f8c3 1791 TFile *f0= TFile::Open(Form("SpaceChargeFluc%d.root",stat));
1792 TTree * tree0 = (TTree*)f0->Get("hisDump");
1793 tree0->SetCacheSize(1000000000);
1794 tree0->SetMarkerStyle(25);
1795 TObjArray * fitArray=new TObjArray(3);
1796 tree0->SetAlias("scNorm",Form("%f/neventsCorr",norm));
4e715689 1797 //
7212f8c3 1798 // Sector Scan
1799 //
1800 TH2 * hisSectorScan[5]={0};
1801 TH1 * hisSectorScanSigma[5]={0};
1802 for (Int_t ihis=0; ihis<5; ihis++){
1803 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");
1804 hisSectorScan[ihis]=(TH2*)tree0->GetHistogram();
1805 hisSectorScan[ihis]->FitSlicesY(0,0,-1,0,"QNR",fitArray);
1806 hisSectorScanSigma[ihis]=(TH1*)(fitArray->At(2)->Clone());
1807 hisSectorScanSigma[ihis]->SetMinimum(0);
1808 hisSectorScanSigma[ihis]->SetMaximum(0.2);
1809 }
1810 gStyle->SetOptStat(0);
1811 gStyle->SetOptTitle(1);
1812 TCanvas * canvasFlucSectorScan=new TCanvas("canvasFlucSectorScan","canvasFlucSectorScan",750,700);
1813 canvasFlucSectorScan->Divide(2,2,0,0);
1814 gStyle->SetPadBorderMode(0);
1815 for (Int_t ihis=0; ihis<4; ihis++){
1816 canvasFlucSectorScan->cd(ihis+1)->SetLogz(1);
1817 hisSectorScan[ihis]->GetXaxis()->SetTitle("r (cm)");
1818 hisSectorScan[ihis]->GetYaxis()->SetTitle("#Delta_{R} (cm)");
1819 hisSectorScan[ihis]->Draw("colz");
1820 TLegend * legendSec=new TLegend(0.5,0.7,0.89,0.89);
1821 legendSec->AddEntry(hisSectorScan[ihis],Form("Sector #Delta %d",(ihis*2+1)));
1822 legendSec->Draw();
1823 }
1824 canvasFlucSectorScan->SaveAs("canvasFlucSectorScan.pdf");
1825 canvasFlucSectorScan->SaveAs("canvasFlucSectorScan.png");
4e715689 1826 //
7212f8c3 1827 gStyle->SetOptTitle(0);
1828 TCanvas * canvasFlucSectorScanFit=new TCanvas("canvasFlucSectorScanFit","canvasFlucSectorScanFit",750,550);
1829 TLegend * legendSector = new TLegend(0.50,0.55,0.89,0.89,"Space charge: corr(sec)-corr(sec-#Delta_{sec})");
1830 for (Int_t ihis=0; ihis<5; ihis++){
1831 hisSectorScanSigma[ihis]->GetXaxis()->SetTitle("r (cm)");
1832 hisSectorScanSigma[ihis]->GetYaxis()->SetTitle("#sigma(#Delta_{R}) (cm)");
1833 hisSectorScanSigma[ihis]->SetMarkerStyle(21+ihis%5);
1834 hisSectorScanSigma[ihis]->SetMarkerColor(1+ihis%4);
1835 if (ihis==0) hisSectorScanSigma[ihis]->Draw("");
1836 hisSectorScanSigma[ihis]->Draw("same");
1837 legendSector->AddEntry(hisSectorScanSigma[ihis],Form("#Delta %d",(ihis*2+1)));
1838 }
1839 legendSector->Draw();
1840 canvasFlucSectorScanFit->SaveAs("canvasFlucSectorScanFit.pdf");
1841 canvasFlucSectorScanFit->SaveAs("canvasFlucSectorScanFit.png");
1842}
4e715689 1843
1844
4e715689 1845
7212f8c3 1846void DrawFluctuationdeltaZ(Int_t stat, Double_t norm){
1847 //
1848 // Draw correction - correction shifted z
1849 // The same set of events used
1850 //Int_t stat=0; Double_t norm=10000;
4b955de3 1851 Int_t deltaZ=16.;
7212f8c3 1852 TFile *f0= TFile::Open(Form("SpaceChargeFluc%d.root",stat));
1853 TTree * tree0 = 0;
1854 if (f0) tree0 = (TTree*)f0->Get("hisDump");
1855 if (!tree0){
4b955de3 1856 tree0 = AliXRDPROOFtoolkit::MakeChainRandom("space.list","hisDump",0,10);
7212f8c3 1857 }
1858 tree0->SetCacheSize(1000000000);
1859 tree0->SetMarkerStyle(25);
1860 TObjArray * fitArray=new TObjArray(3);
1861 tree0->SetAlias("scNorm",Form("%f/neventsCorr",norm));
1862 //
1863 // DeltaZ Scan
1864 //
1865 TH2 * hisDeltaZScan[6]={0};
1866 TH1 * hisDeltaZScanSigma[6]={0};
1867 for (Int_t ihis=0; ihis<6; ihis++){
4b955de3 1868 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 1869 hisDeltaZScan[ihis]=(TH2*)tree0->GetHistogram();
1870 hisDeltaZScan[ihis]->FitSlicesY(0,0,-1,0,"QNR",fitArray);
1871 hisDeltaZScanSigma[ihis]=(TH1*)(fitArray->At(2)->Clone());
1872 hisDeltaZScanSigma[ihis]->SetMinimum(0);
1873 hisDeltaZScanSigma[ihis]->SetMaximum(0.2);
1874 }
1875 gStyle->SetOptStat(0);
1876 gStyle->SetOptTitle(1);
1877 TCanvas * canvasFlucDeltaZScan=new TCanvas("canvasFlucDeltaZScan","canvasFlucDeltaZScan",700,700);
1878 canvasFlucDeltaZScan->Divide(3,2,0,0);
1879 gStyle->SetPadBorderMode(0);
1880 for (Int_t ihis=0; ihis<6; ihis++){
1881 canvasFlucDeltaZScan->cd(ihis+1)->SetLogz(1);
1882 hisDeltaZScan[ihis]->GetXaxis()->SetTitle("r (cm)");
1883 hisDeltaZScan[ihis]->GetYaxis()->SetTitle("#Delta_{R} (cm)");
1884 hisDeltaZScan[ihis]->Draw("colz");
1885 TLegend * legendSec=new TLegend(0.5,0.7,0.89,0.89);
4b955de3 1886 legendSec->AddEntry(hisDeltaZScan[ihis],Form("DeltaZ #Delta %d",(ihis+1)*deltaZ));
7212f8c3 1887 legendSec->Draw();
1888 }
1889 canvasFlucDeltaZScan->SaveAs(Form("canvasFlucDeltaZScan%d.pdf",stat));
1890 canvasFlucDeltaZScan->SaveAs(Form("canvasFlucDeltaZScan%d.png",stat));
1891
1892 //
1893 gStyle->SetOptTitle(0);
1894 TCanvas * canvasFlucDeltaZScanFit=new TCanvas("canvasFlucDeltaZScanFit","canvasFlucDeltaZScanFit");
1895 TLegend * legendDeltaZ = new TLegend(0.50,0.55,0.89,0.89,"Space charge: corr(z_{ref})-corr(z_{ref}-#Delta_{z})");
1896 for (Int_t ihis=0; ihis<5; ihis++){
1897 hisDeltaZScanSigma[ihis]->GetXaxis()->SetTitle("r (cm)");
1898 hisDeltaZScanSigma[ihis]->GetYaxis()->SetTitle("#sigma(#Delta_{R}) (cm)");
1899 hisDeltaZScanSigma[ihis]->SetMarkerStyle(21+ihis%5);
1900 hisDeltaZScanSigma[ihis]->SetMarkerColor(1+ihis%4);
1901 if (ihis==0) hisDeltaZScanSigma[ihis]->Draw("");
1902 hisDeltaZScanSigma[ihis]->Draw("same");
4b955de3 1903 legendDeltaZ->AddEntry(hisDeltaZScanSigma[ihis],Form("#Delta %d (cm)",(ihis+1)*deltaZ));
7212f8c3 1904 }
1905 legendDeltaZ->Draw();
1906 canvasFlucDeltaZScanFit->SaveAs(Form("canvasFlucDeltaZScanFit%d.pdf",stat));
1907 canvasFlucDeltaZScanFit->SaveAs(Form("canvasFlucDeltaZScanFit%d.png",stat));
4e715689 1908
7212f8c3 1909}
4e715689 1910
4e715689 1911
7212f8c3 1912void DrawDefault(Int_t stat){
1913 //
1914 // Draw correction - correction shifted z
1915 // The same set of events used
1916 // Int_t stat=0
1917 TFile *f0= TFile::Open(Form("SpaceChargeFluc%d.root",stat));
1918 TTree * tree0 = (TTree*)f0->Get("hisDump");
1919 tree0->SetCacheSize(1000000000);
1920 tree0->SetMarkerStyle(25);
1921 tree0->SetMarkerSize(0.4);
4b955de3 1922 // TObjArray * fitArray=new TObjArray(3);
7212f8c3 1923 tree0->Draw("10000*DistRefDR/neventsCorr:r:z/r","abs(z/r)<0.9&&z>0&&rndm>0.8","colz");
4b955de3 1924}
1925
1926
1927
1928
1929void DrawTrackFluctuation(){
1930 //
1931 // Function to make a fluctuation figures for differnt multiplicities of pileup space charge
1932 // it is assumed that the text files
1933 //
1934 //
1935 TObjArray arrayFit(3);
1936 const char *inputList;
1937 TH2F * hisCorrRef[5]={0};
1938 TH2F * hisCorrNo[5]={0};
1939 TH1 * hisCorrRefM[5], *hisCorrRefRMS[5];
1940 TH1 * hisCorrNoM[5], *hisCorrNoRMS[5];
1941 //
1942 // 1. Load chains for different statistic
1943 //
1944 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";
1945 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";
1946 TChain * chains[5]={0};
1947 TChain * chainR = AliXRDPROOFtoolkit::MakeChain("track0_1.list","trackFit",0,1000);
1948 chainR->SetCacheSize(1000000000);
1949 for (Int_t ichain=0; ichain<5; ichain++){
1950 chains[ichain] = AliXRDPROOFtoolkit::MakeChain(Form("track%d_1.list",2*(ichain+1)),"trackFit",0,1000);
1951 chains[ichain]->AddFriend(chainR,"R");
1952 chains[ichain]->SetCacheSize(1000000000);
1953 chains[ichain]->SetMarkerStyle(25);
1954 chains[ichain]->SetMarkerSize(0.5);
1955 chains[ichain]->SetAlias("meanNorm","(1+0.2*abs(neventsCorr/10000-1))"); // second order correction - renomalization of mean hardwired
1956 // to be fitted?
1957 }
1958 //
1959 // 2. fill histograms if not available in file
1960 //
1961 //
1962 TFile *ftrackFluctuation = TFile::Open("trackFluctuation.root","update");
1963 for (Int_t ihis=0; ihis<5; ihis++){
1964 ftrackFluctuation->cd();
1965 hisCorrRef[ihis] = (TH2F*)(ftrackFluctuation->Get(Form("DeltaRPhiCorr%d",(ihis+1)*2000)));
1966 hisCorrNo[ihis] = (TH2F*)(ftrackFluctuation->Get(Form("DeltaRPhi%d",(ihis+1)*2000)));
1967 if (hisCorrRef[ihis]==0) {
1968 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");
1969 hisCorrRef[ihis]=(TH2F*)(chains[ihis]->GetHistogram()->Clone());
1970 hisCorrRef[ihis]->SetName(Form("DeltaRPhiCorr%d",(ihis+1)*2000));
1971 hisCorrRef[ihis]->SetTitle(Form("Corrected #Delta r#phi - Pileup %d",(ihis+1)*2000));
1972 hisCorrRef[ihis]->GetXaxis()->SetTitle("1/p_{t} (1/GeV/c)");
1973 hisCorrRef[ihis]->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
1974 hisCorrRef[ihis]->Write();
1975 //
1976 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");
1977 hisCorrNo[ihis]=(TH2F*)(chains[ihis]->GetHistogram()->Clone());
1978 hisCorrNo[ihis]->SetName(Form("DeltaRPhi%d",(ihis+1)*2000));
1979 hisCorrNo[ihis]->SetTitle(Form("Delta r#phi = Pileup %d",(ihis+1)*2000));
1980 hisCorrNo[ihis]->GetXaxis()->SetTitle("1/p_{t} (1/GeV/c)");
1981 hisCorrNo[ihis]->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
1982 hisCorrNo[ihis]->Write();
1983 }
1984 }
1985 ftrackFluctuation->Flush();
1986 //
1987 //
1988 //
1989 for (Int_t ihis=0; ihis<5; ihis++){
1990 hisCorrRef[ihis]->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
1991 hisCorrRefM[ihis] = (TH1*)arrayFit.At(1)->Clone();
1992 hisCorrRefRMS[ihis] = (TH1*)arrayFit.At(2)->Clone();
1993 hisCorrRefM[ihis]->GetXaxis()->SetTitle("1/p_{t} (1/GeV/c)");
1994 hisCorrRefM[ihis]->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
1995 hisCorrRefM[ihis]->SetMarkerStyle(20);
1996 hisCorrRefRMS[ihis]->SetMarkerStyle(21);
1997 hisCorrRefM[ihis]->SetMarkerColor(1);
1998 hisCorrRefRMS[ihis]->SetMarkerColor(2);
1999 hisCorrNo[ihis]->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2000 hisCorrNoM[ihis] = (TH1*)arrayFit.At(1)->Clone();
2001 hisCorrNoRMS[ihis] = (TH1*)arrayFit.At(2)->Clone();
2002 }
4e715689 2003
4b955de3 2004 //
2005 TCanvas *canvasMean = new TCanvas("canvasCorrectionMean","canvasCorrectionMean",900,1000);
2006 TCanvas *canvasMeanSummary = new TCanvas("canvasCorrectionMeanSummary","canvasCorrectionMeanSummary",700,600);
2007
2008 canvasMean->Divide(3,5);
2009 gStyle->SetOptStat(0);
2010 for (Int_t ihis=0; ihis<5; ihis++){
2011 TLegend * legend = new TLegend(0.11,0.11,0.5,0.3,Form("Pile up %d",(ihis+1)*2000));
2012 canvasMean->cd(3*ihis+1);
2013 hisCorrNo[ihis]->Draw("colz");
2014 canvasMean->cd(3*ihis+2);
2015 hisCorrRef[ihis]->Draw("colz");
2016 canvasMean->cd(3*ihis+3);
2017 hisCorrRefM[ihis]->SetMaximum(0.25);
2018 hisCorrRefM[ihis]->SetMinimum(-0.25);
2019 hisCorrRefM[ihis]->Draw("");
2020 hisCorrRefRMS[ihis]->Draw("same");
2021 legend->AddEntry(hisCorrRefM[ihis],"Mean");
2022 legend->AddEntry(hisCorrRefRMS[ihis],"RMS");
2023 legend->Draw();
2024 }
2025 canvasMeanSummary->cd();
2026 TLegend * legendMeanSummary = new TLegend(0.5,0.6,0.89,0.89,"Space charge correction fluctuation in r#phi");
2027 for (Int_t ihis=4; ihis>=0; ihis--){
2028 hisCorrRefRMS[ihis]->SetMarkerColor(1+ihis);
2029 hisCorrRefRMS[ihis]->SetMinimum(0);
2030 hisCorrRefRMS[ihis]->GetYaxis()->SetTitle("#sigma_{r#phi} (cm)");
2031 if (ihis==4) hisCorrRefRMS[ihis]->Draw("");
2032 hisCorrRefRMS[ihis]->Draw("same");
2033 legendMeanSummary->AddEntry(hisCorrRefRMS[ihis],Form("%d pile-up events",(ihis+1)*2000));
2034 }
2035 legendMeanSummary->Draw();
4e715689 2036
4b955de3 2037 canvasMean->SaveAs("canvasCorrectionMean.pdf");
2038 canvasMeanSummary->SaveAs("canvasCorrectionMeanSummary.pdf");
2039 //canvasMean->Write();
2040 //canvasMeanSummary->Write();
2041 ftrackFluctuation->Close();
7212f8c3 2042}
4b955de3 2043
2044void DrawTrackFluctuationZ(){
2045 //
2046 // Draw track fucutation dz
2047 //
2048 const Int_t kColors[6]={1,2,3,4,6,7};
2049 const Int_t kStyle[6]={20,21,24,25,24,25};
2050 TObjArray arrayFit(3);
2051 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";
2052 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";
2053 TChain * chains[5]={0};
2054 TChain * chainR = AliXRDPROOFtoolkit::MakeChain("track0_1.list","trackFit",0,1000);
2055 chainR->SetCacheSize(1000000000);
2056 for (Int_t ichain=0; ichain<5; ichain++){
2057 chains[ichain] = AliXRDPROOFtoolkit::MakeChain(Form("track%d_1.list",2*(ichain+1)),"trackFit",0,1000);
2058 chains[ichain]->AddFriend(chainR,"R");
2059 chains[ichain]->SetCacheSize(1000000000);
2060 chains[ichain]->SetMarkerStyle(25);
2061 chains[ichain]->SetMarkerSize(0.5);
2062 }
2063 //
2064 // 2.) Create 2D histo or read from files
2065 //
2066 TObjArray * arrayHisto = new TObjArray(25);
2067 TFile *ftrackFluctuationZ = TFile::Open("trackFluctuationZ.root","update");
2068 for (Int_t ihis=0; ihis<5; ihis++){
2069 ftrackFluctuationZ->cd();
2070 for (Int_t idz=0; idz<5; idz++){
2071 Int_t z= 16+idz*32;
2072 TH2 *his= (TH2*)ftrackFluctuationZ->Get(Form("TrackDz%d_PileUp%d",z, (ihis+1)*2000));
2073 if (!his){
2074 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");
2075 his = (TH2*)(chains[ihis]->GetHistogram()->Clone());
2076 his->SetName(Form("TrackDz%d_PileUp%d",z, (ihis+1)*2000));
2077 his->Write();
2078 }
2079 arrayHisto->AddAtAndExpand(his,ihis*5+idz);
2080 }
2081 }
2082 ftrackFluctuationZ->Flush();
2083
2084 //
2085 // 3.) Make fits
2086 //
2087 TCanvas *canvasDz = new TCanvas("canvasDz","canvasDz",800,800);
2088 canvasDz->Divide(2,2,0,0);
2089 for (Int_t ihis=3; ihis>=0; ihis--){
2090 canvasDz->cd(ihis+1)->SetTicks(3);
2091 TLegend * legend = new TLegend(0.31,0.51, 0.95,0.95,Form("Distortion due time/z delay (Pileup=%d)", (ihis+1)*2000));
2092 legend->SetBorderSize(0);
2093 for (Int_t idz=3; idz>=0; idz--){
2094 TH2 * his = (TH2*)arrayHisto->At(ihis*5+idz);
2095 his->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2096 TH1 * hisRMS = (TH1*)arrayFit.At(2)->Clone();
2097 hisRMS->SetMaximum(0.12);
2098 hisRMS->SetMinimum(0);
2099 hisRMS->GetXaxis()->SetTitle("1/p_{t} (GeV/c)");
2100 hisRMS->GetYaxis()->SetTitle("#sigma_{r#phi}(cm)");
2101 hisRMS->SetMarkerStyle(kStyle[idz]);
2102 hisRMS->SetMarkerColor(kColors[idz]);
2103 if (idz==3) hisRMS->Draw();
2104 legend->AddEntry(hisRMS,Form("#Delta_{z}=%d (cm)",16+idz*32));
2105 hisRMS->Draw("same");
2106 }
2107 legend->Draw();
2108 }
2109 canvasDz->SaveAs("spaceChargeDeltaZScan.pdf");
2110
2111}
2112
2113
2ff913b6 2114
2115
2116
2117void DrawTrackFluctuationFrame(){
2118 //
2119 // Function to make a fluctuation figures for differnt multiplicities of pileup space charge
2120 // it is assumed that the text files
2121 //
2122 //
2123 TObjArray arrayFit(3);
2124 const char *inputList;
2125 TH2F * hisCorrRef[10]={0};
2126 TH2F * hisCorrNo[10]={0};
2127 TH1 * hisCorrRefM[10], *hisCorrRefRMS[10];
2128 TH1 * hisCorrNoM[10], *hisCorrNoRMS[10];
2129 //
2130 // 1. Load chains for different statistic
2131 //
2132 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";
2133 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";
2134 TCut cutFit="Entry$%4==0"; //use only subset of data for fit
2135
2136 TChain * chains[10]={0};
2137 TChain * chainR = AliXRDPROOFtoolkit::MakeChain("track0_1.list","trackFit",0,1000);
2138 chainR->SetCacheSize(1000000000);
2139 for (Int_t ichain=0; ichain<7; ichain++){
2140 chains[ichain] = AliXRDPROOFtoolkit::MakeChain(Form("track%d_1.list",2*(ichain+1)),"trackFit",0,1000);
2141 chains[ichain]->AddFriend(chainR,"R");
2142 chains[ichain]->SetCacheSize(1000000000);
2143 chains[ichain]->SetMarkerStyle(25);
2144 chains[ichain]->SetMarkerSize(0.5);
2145 chains[ichain]->SetAlias("meanNorm","(1+0.2*abs(neventsCorr/10000-1))"); // second order correction - renomalization of mean hardwired
2146 chains[ichain]->SetAlias("dMean0","(neventsCorr*R.T_DistRef_0.fP[0]/10000)");
2147 chains[ichain]->SetAlias("dMeas0","T_DistRef_0.fP[0]");
2148 chains[ichain]->SetAlias("dMean1","(neventsCorr*R.T_DistRef_1.fP[0]/10000)");
2149 chains[ichain]->SetAlias("dMeas1","T_DistRef_1.fP[0]");
2150 for (Int_t ig=0; ig<10;ig++) chains[ichain]->SetAlias(Form("FR%d",ig),Form("(abs(Entry$-%d)<1000)",ig*2000+1000));
2151 }
2152 //
2153 // 2. Get or Create histogram (do fit per frame)
2154 //
2155 TStatToolkit toolkit;
2156 Double_t chi2=0;
2157 Int_t npoints=0;
2158 TVectorD param;
2159 TMatrixD covar;
2160 TString fstringG=""; // global part
2161 fstringG+="dMean0++";
2162 TVectorD vec0,vec1;
2163 TString fstringF0=""; // global part
2164 for (Int_t ig=0; ig<10;ig++) fstringF0+=Form("FR%d++",ig);
2165 for (Int_t ig=0; ig<10;ig++) fstringF0+=Form("FR%d*dMean0++",ig);
2166 TString fstringF1=""; // global part
2167 for (Int_t ig=0; ig<10;ig++) fstringF1+=Form("FR%d++",ig);
2168 for (Int_t ig=0; ig<10;ig++) fstringF1+=Form("FR%d*dMean0++",ig);
2169 for (Int_t ig=0; ig<10;ig++) fstringF1+=Form("FR%d*dMean0*abs(T_DistRef_0.fP[3])++",ig);
2170 for (Int_t ig=0; ig<10;ig++) fstringF1+=Form("FR%d*dMean0*(T_DistRef_0.fP[3]^2)++",ig);
2171
2172 //
2173 //
2174 TH2F *hisA=0, *hisF0=0, *hisF1=0, *hisM=0;
2175 TObjArray * arrayHisto = new TObjArray(200);
2176 TFile *ftrackFit = TFile::Open("trackFluctuationFrame.root","update");
2177 for (Int_t ihis=0; ihis<7; ihis++){
2178 printf("\n\nProcessing frames\t%d\nnn",(ihis+1)*2000);
2179 hisM = (TH2F*)ftrackFit->Get(Form("hisMean_%d",(ihis+1)*2000));
2180 hisA = (TH2F*)ftrackFit->Get(Form("hisAll_%d",(ihis+1)*2000));
2181 hisF0 = (TH2F*)ftrackFit->Get(Form("hisFrame0_%d",(ihis+1)*2000));
2182 hisF1 = (TH2F*)ftrackFit->Get(Form("hisFrame1_%d",(ihis+1)*2000));
2183 if (!hisA){
2184 ftrackFit->cd();
2185 TString * fitResultAll = TStatToolkit::FitPlane(chains[ihis],"dMeas0", fstringG.Data(),cutOut+cutOutF+cutFit, chi2,npoints,param,covar,-1,0, 40*2000, kFALSE);
2186 chains[ihis]->SetAlias("fitAll",fitResultAll->Data());
2187 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);
2188 chains[ihis]->SetAlias("fitF0",fitResultF0->Data());
2189 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);
2190 chains[ihis]->SetAlias("fitF1",fitResultF1->Data());
2191 fitResultF0->Tokenize("++")->Print();
2192 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);
2193 hisA = (TH2F*)chains[ihis]->GetHistogram();
2194 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);
2195 hisF0 = (TH2F*)chains[ihis]->GetHistogram();
2196 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);
2197 hisF1 = (TH2F*)chains[ihis]->GetHistogram();
2198 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);
2199 hisM = (TH2F*)chains[ihis]->GetHistogram();
2200 hisM->Write(); hisA->Write();hisF0->Write(); hisF1->Write();
2201 ftrackFit->Flush();
2202 }
2203 }
2204
2205 for (Int_t ihis=0; ihis<7; ihis++){
2206 printf("\n\nProcessing frames\t%d\nnn",(ihis+1)*2000);
2207 hisM = (TH2F*)ftrackFit->Get(Form("hisMean_%d",(ihis+1)*2000));
2208 hisA = (TH2F*)ftrackFit->Get(Form("hisAll_%d",(ihis+1)*2000));
2209 hisF0 = (TH2F*)ftrackFit->Get(Form("hisFrame0_%d",(ihis+1)*2000));
2210 hisF1 = (TH2F*)ftrackFit->Get(Form("hisFrame1_%d",(ihis+1)*2000));
2211 arrayHisto->AddLast(hisA);
2212 arrayHisto->AddLast(hisF0);
2213 arrayHisto->AddLast(hisF1);
2214 arrayHisto->AddLast(hisM);
2215 }
2216 delete ftrackFit;
2217 //
2218 // 3. Draw figures
2219 //
2220 gStyle->SetOptStat(0);
2221 TCanvas *canvasFit = new TCanvas("canvasFitFrame","canvasFitframe",900,700);
2222 canvasFit->Divide(3,2,0,0);
2223 for (Int_t ihis=1; ihis<7; ihis++){
2224 //
2225 canvasFit->cd(ihis);
2226 char hname[10000];
2227 snprintf(hname,1000,"hisAll_%d",(ihis+1)*2000);
2228 hisA = (TH2F*)arrayHisto->FindObject(hname);
2229 snprintf(hname,1000,"hisFrame0_%d",(ihis+1)*2000);
2230 hisF0 = (TH2F*)arrayHisto->FindObject(hname);
2231 snprintf(hname,1000,"hisFrame1_%d",(ihis+1)*2000);
2232 hisF1 = (TH2F*)arrayHisto->FindObject(hname);
2233 snprintf(hname,1000,"hisMean_%d",(ihis+1)*2000);
2234 hisM = (TH2F*)arrayHisto->FindObject(hname);
2235 //
2236 //
2237 hisM->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2238 TH1 * hisRA= (TH1*)arrayFit.At(2)->Clone();
2239 hisF0->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2240 TH1 * hisRF0= (TH1*)arrayFit.At(2)->Clone();
2241 hisF1->FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2242 TH1 * hisRF1= (TH1*)arrayFit.At(2)->Clone();
2243 //
2244 hisRA->SetMarkerStyle(20);
2245 hisRF0->SetMarkerStyle(21);
2246 hisRF1->SetMarkerStyle(21);
2247 hisRA->SetMarkerColor(1);
2248 hisRF0->SetMarkerColor(4);
2249 hisRF1->SetMarkerColor(2);
2250 TF1 * f1a= new TF1("f1a","pol1");
2251 TF1 * f1f0= new TF1("f1a0","pol1");
2252 TF1 * f1f1= new TF1("f1a1","pol1");
2253 f1a->SetLineColor(1);
2254 f1f0->SetLineColor(4);
2255 f1f1->SetLineColor(2);
2256 hisRA->Fit(f1a);
2257 hisRF0->Fit(f1f0);
2258 hisRF1->Fit(f1f1);
2259 hisRF1->SetMinimum(0);
2260 hisRF1->SetMaximum(0.05);
2261 // hisRA->Draw();
2262 hisRF1->GetXaxis()->SetTitle("q/p_{T} (1/GeV)");
2263 hisRF1->GetYaxis()->SetTitle("#sigma_{r#phi} (cm)");
2264 hisRF1->Draw("");
2265 hisRF0->Draw("same");
2266 TLegend * legend = new TLegend(0.11,0.11,0.65,0.25, Form("Track residual r#phi distortion: N_{ion}=%d",(ihis+1)*2000));
2267 legend->AddEntry(hisRF0,"a_{0}+a_{1}#rho");
2268 legend->AddEntry(hisRF1,"a_{0}+(a_{1}+a_{2}z+a_{3}z^2)#rho");
2269 legend->SetBorderSize(0);
2270 legend->Draw();
2271 }
2272 //
2273 canvasFit->SaveAs("canvasFrameFitRPhiVersion0.pdf");
2274 canvasFit->SaveAs("canvasFrameFitRPhiVersion0.png");
2275 //
2276}