TPC/Upgrade/macros/spaceChargeFluctuation.C

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