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