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