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