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