TPC/Upgrade/macros/spaceChargeFluctuation.C

   1 /*
   2   .x ~/rootlogon.C
   3   .L $ALICE_ROOT/TPC/Upgrade/macros/spaceChargeFluctuation.C+
   4 //  GenerateMapRaw("/hera/alice/local/filtered/alice/data/2010/LHC10e/000129527/raw/merged_highpt_1.root","output.root", 50);
   5   GenerateMapRawIons( kFALSE, "/hera/alice/local/filtered/alice/data/2010/LHC10e/000129527/raw/merged_highpt_1.root","output.root", 25);
   6
   7  */
   8 #include "TMath.h"
   9 #include "TRandom.h"
  10 #include "TTreeStream.h"
  11 #include "TVectorD.h"
  12 #include "TCanvas.h"
  13 #include "TStopwatch.h"
  14 #include "AliTPCParam.h"
  15 #include "AliTPCcalibDB.h"
  16 #include "AliTPCAltroMapping.h"
  17 #include "AliAltroRawStream.h"
  18 #include "AliSysInfo.h"
  19 #include "AliTPCRawStreamV3.h"
  20 #include "AliCDBManager.h"
  21 #include "TGeoGlobalMagField.h"
  22 #include "AliMagF.h"
  23 #include "AliRawReaderRoot.h"
  24 #include "AliRawReader.h"
  25 #include "TH3.h"
  26 #include "TH2.h"
  27 #include "AliTPCCalPad.h"
  28 #include "AliTPCCalROC.h"
  29 #include "TChain.h"
  30 #include "AliXRDPROOFtoolkit.h"
  31 #include "TLegend.h"
  32 #include "TCut.h"
  33 #include "TGraphErrors.h"
  34 #include "TStatToolkit.h"
  35
  36
  37 TH1 * GenerateMapRawIons(Int_t useGain,const char *fileName="raw.root", const char *outputName="histo.root", Int_t maxEvents=25);
  38 void DoMerge();
  39
  40 void spaceChargeFluctuation(Int_t mode=0, Int_t arg0=0){
  41   //
  42   //
  43   if (mode==0) GenerateMapRawIons(arg0);
  44   if (mode==1) DoMerge();
  45 }
  46
  47 void spaceChargeFluctuationMC(Int_t nframes){
  48   //
  49   // Toy MC to generate space charge fluctuation
  50   //
  51   //
  52   TTreeSRedirector *pcstream = new TTreeSRedirector("spaceChargeFluctuation.root","recreate");
  53   Double_t interactionRate=100000;
  54   Double_t driftTime=0.1;
  55
  56   for (Int_t iframe=0; iframe<nframes; iframe++){
  57     Int_t nevents=gRandom->Poisson(interactionRate*driftTime);
  58     Int_t ntracksAll=0;
  59     TVectorD vecAllPhi128(128);
  60     TVectorD vecAllPhi36(36);
  61     for (Int_t ievent=0; ievent<nevents; ievent++){
  62       //
  63       //
  64       Int_t ntracks=gRandom->Exp(500)+gRandom->Exp(500); // to be taken form the distribution
  65       ntracksAll+=ntracks;
  66       for (Int_t isec=0; isec<128; isec++){
  67         vecAllPhi128(isec)+=ntracks*gRandom->Rndm()/128;
  68       }
  69       for (Int_t isec=0; isec<36; isec++){
  70         vecAllPhi36(isec)+=ntracks*gRandom->Rndm()/36;
  71       }
  72     }
  73     (*pcstream)<<"ntracks"<<
  74       "nevents="<<nevents<<                  // number of events withing time frame
  75       "ntracksAll="<<ntracksAll<<            // number of tracks within  time frame
  76       "vecAllPhi36="<<&vecAllPhi36<<         // number of tracks in phi bin (36 bins)    within time frame
  77       "vecAllPhi128="<<&vecAllPhi128<<       // number of tracks in phi bin (128 bins)   within time frame
  78       "\n";
  79   }
  80   delete pcstream;
  81 }
  82
  83
  84 // void spaceChargeFluctuationDraw(){
  85 //   //
  86 //   //
  87 //   //
  88 //   TFile *f = TFile::Open("spaceChargeFluctuation.root");
  89 //   TTree * tree = (TTree*)f->Get("ntracks");
  90 //   TCanvas * canvasFluc = new TCanvas("cavasFluc","canvasFluc");
  91 // }
  92
  93
  94
  95
  96
  97
  98 TH1 * GenerateMapRawIons(Int_t useGainMap, const char *fileName, const char *outputName, Int_t maxEvents){
  99   //
 100   // Generate 3D maps of the space charge for the rad data maps
 101   // different threshold considered
 102   // Paramaters:
 103   //    useGainMap    - switch usage of the gain map
 104   //    fileName      - name of input raw file
 105   //    outputName    - name of output file with the space charge histograms
 106   //    maxEvents     - grouping of the events
 107   //
 108   //
 109   TTreeSRedirector * pcstream  = new TTreeSRedirector(outputName, "recreate");
 110   const char *  ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";
 111   AliCDBManager * man = AliCDBManager::Instance();
 112   man->SetDefaultStorage(ocdbpath);
 113   man->SetRun(0);
 114   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG,       AliMagF::kBeamTypepp, 2.76/2.));
 115   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
 116   AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
 117   AliTPCCalPad * gain = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 118   AliTPCCalPad * noise = AliTPCcalibDB::Instance()->GetPadNoise();
 119
 120   TStopwatch timer;
 121   timer.Start();
 122   //   arrays of space charges - different elements corresponds to different threshold to accumulate charge
 123   TH3D * hisQ3D[3]={0};                // 3D maps space charge from drift volume
 124   TH2D * hisQ2D[3]={0};
 125   TH2D * hisQ2DRZ[3]={0};
 126   TH1D * hisQ1D[3]={0};
 127   TH3D * hisQ3DROC[3]={0};             // 3D maps space charge from ROC
 128   TH2D * hisQ2DROC[3]={0};
 129   TH2D * hisQ2DRZROC[3]={0};
 130   TH1D * hisQ1DROC[3]={0};
 131   Int_t nbinsRow=param->GetNRowLow()+param->GetNRowUp();
 132   Double_t *xbins = new Double_t[nbinsRow+1];
 133   xbins[0]=param->GetPadRowRadiiLow(0)-1;   //underflow bin
 134   for (Int_t ibin=0; ibin<param->GetNRowLow();ibin++) xbins[1+ibin]=param->GetPadRowRadiiLow(ibin);
 135   for (Int_t ibin=0; ibin<param->GetNRowUp();ibin++)  xbins[1+ibin+param->GetNRowLow()]=param->GetPadRowRadiiUp(ibin);
 136
 137
 138   //
 139   for (Int_t ith=0; ith<3; ith++){
 140     char chname[100];
 141     snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
 142     hisQ3D[ith] = new TH3D(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
 143     snprintf(chname,100,"hisQ2D_Th%d",2*ith+2);
 144     hisQ2D[ith] = new TH2D(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
 145     snprintf(chname,100,"hisQ2DRZ_Th%d",2*ith+2);
 146
 147     hisQ2DRZ[ith] = new TH2D(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36),  125,-250,250);
 148
 149     snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
 150     hisQ1D[ith] = new TH1D(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
 151     //
 152     snprintf(chname,100,"hisQ3DROC_Th%d",2*ith+2);
 153     hisQ3DROC[ith] = new TH3D(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
 154     snprintf(chname,100,"hisQ2DROC_Th%d",2*ith+2);
 155     hisQ2DROC[ith] = new TH2D(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
 156     snprintf(chname,100,"hisQ2DRZROC_Th%d",2*ith+2);
 157     hisQ2DRZROC[ith] = new TH2D(chname,chname,param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36), 125,-250,250);
 158     snprintf(chname,100,"hisQ1DROC_Th%d",2*ith+2);
 159     hisQ1DROC[ith] = new TH1D(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
 160     hisQ1D[ith]->GetXaxis()->Set(nbinsRow,xbins);
 161     hisQ2D[ith]->GetYaxis()->Set(nbinsRow,xbins);
 162     hisQ2DRZ[ith]->GetXaxis()->Set(nbinsRow,xbins);
 163     hisQ3D[ith]->GetYaxis()->Set(nbinsRow,xbins);
 164     hisQ1DROC[ith]->GetXaxis()->Set(nbinsRow,xbins);
 165     hisQ2DROC[ith]->GetYaxis()->Set(nbinsRow,xbins);
 166     hisQ2DRZROC[ith]->GetXaxis()->Set(nbinsRow,xbins);
 167     hisQ3DROC[ith]->GetYaxis()->Set(nbinsRow,xbins);
 168     //
 169     hisQ3D[ith]->SetDirectory(0);
 170     hisQ2D[ith]->SetDirectory(0);
 171     hisQ2DRZ[ith]->SetDirectory(0);
 172     hisQ1D[ith]->SetDirectory(0);
 173     hisQ3DROC[ith]->SetDirectory(0);
 174     hisQ2DRZROC[ith]->SetDirectory(0);
 175     hisQ2DROC[ith]->SetDirectory(0);
 176     hisQ1DROC[ith]->SetDirectory(0);
 177   }
 178   //
 179   //
 180   AliRawReader *reader = new AliRawReaderRoot(fileName);
 181   reader->Reset();
 182   AliAltroRawStream* stream = new AliAltroRawStream(reader);
 183   stream->SelectRawData("TPC");
 184   Int_t evtnr=0;
 185   Int_t chunkNr=0;
 186   //
 187
 188   while (reader->NextEvent()) {
 189     Double_t shiftZ= gRandom->Rndm()*250.;
 190     //
 191     if(evtnr>=maxEvents) {
 192       chunkNr++;
 193       pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
 194       pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
 195       for (Int_t ith=0; ith<3; ith++){
 196         hisQ1D[ith]->Write(Form("His1DDrift_%d",ith));
 197         hisQ2D[ith]->Write(Form("His2DDrift_%d",ith));
 198         hisQ2DRZ[ith]->Write(Form("His2DDriftRZ_%d",ith));
 199         hisQ3D[ith]->Write(Form("His3DDrift_%d",ith));
 200         hisQ1DROC[ith]->Write(Form("His1DROC_%d",ith));
 201         hisQ2DROC[ith]->Write(Form("His2DROC_%d",ith));
 202         hisQ2DRZROC[ith]->Write(Form("His2DRZROC_%d",ith));
 203         hisQ3DROC[ith]->Write(Form("His3DROC_%d",ith));
 204         (*pcstream)<<"histo"<<
 205           "events="<<evtnr<<
 206           "useGain="<<useGainMap<<
 207           Form("hist1D_%d.=",ith*2+2)<<hisQ1D[ith]<<
 208           Form("hist2D_%d.=",ith*2+2)<<hisQ2D[ith]<<
 209           Form("hist2DRZ_%d.=",ith*2+2)<<hisQ2DRZ[ith]<<
 210           Form("hist3D_%d.=",ith*2+2)<<hisQ3D[ith]<<
 211           Form("hist1DROC_%d.=",ith*2+2)<<hisQ1DROC[ith]<<
 212           Form("hist2DROC_%d.=",ith*2+2)<<hisQ2DROC[ith]<<
 213           Form("hist2DRZROC_%d.=",ith*2+2)<<hisQ2DRZROC[ith]<<
 214           Form("hist3DROC_%d.=",ith*2+2)<<hisQ3DROC[ith];
 215       }
 216       (*pcstream)<<"histo"<<"\n";
 217       for (Int_t ith=0; ith<3; ith++){
 218         hisQ1D[ith]->Reset();
 219         hisQ2D[ith]->Reset();
 220         hisQ2DRZ[ith]->Reset();
 221         hisQ3D[ith]->Reset();
 222         hisQ1DROC[ith]->Reset();
 223         hisQ2DROC[ith]->Reset();
 224         hisQ2DRZROC[ith]->Reset();
 225         hisQ3DROC[ith]->Reset();
 226       }
 227       evtnr=0;
 228     }
 229     cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
 230     evtnr++;
 231     AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
 232     AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
 233     //
 234     while (input.NextDDL()){
 235       Int_t sector = input.GetSector();
 236       AliTPCCalROC * gainROC =gain->GetCalROC(sector);
 237       AliTPCCalROC * noiseROC =noise->GetCalROC(sector);
 238       while ( input.NextChannel() ) {
 239         Int_t    row    = input.GetRow();
 240         Int_t rowAbs    = row+ ((sector<36) ? 0: param->GetNRow(0));
 241         Int_t    pad    = input.GetPad();
 242         Int_t    nPads   = param->GetNPads(sector,row);
 243         Double_t localX  = param->GetPadRowRadii(sector,row);
 244         Double_t localY  = (pad-nPads)*param->GetPadPitchWidth(sector);
 245         Double_t localPhi= TMath::ATan2(localY,localX);
 246         Double_t phi     = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
 247         if (sector%36>=18) phi+=TMath::TwoPi();
 248         Double_t padLength=param->GetPadPitchLength(sector,row);
 249         //  Double_t padWidth=param->GetPadPitchWidth(sector);
 250         //  Double_t zLength=param->GetZLength();
 251         Double_t deltaPhi=hisQ3D[0]->GetXaxis()->GetBinWidth(0);
 252         Double_t deltaZ= hisQ3D[0]->GetZaxis()->GetBinWidth(0);
 253         Double_t gainPad = gainROC->GetValue(row,pad);
 254         Double_t noisePad = noiseROC->GetValue(row,pad);
 255         //
 256         while ( input.NextBunch() ){
 257           Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
 258           Int_t  bunchlength  = (Int_t)input.GetBunchLength();
 259           const UShort_t *sig = input.GetSignals();
 260           Int_t aboveTh[3]={0};
 261           for (Int_t i=0; i<bunchlength; i++){
 262             if (sig[i]<4*noisePad) continue;
 263             for (Int_t ith=0; ith<3; ith++){
 264               if (sig[i]>(ith*2)+2) aboveTh[ith]++;
 265             }
 266           }
 267           for (Int_t ith=0; ith<3; ith++){
 268             if (aboveTh[ith%3]>1){
 269               for (Int_t i=0; i<bunchlength; i++){
 270                 //
 271                 // normalization
 272                 //
 273                 Double_t zIonDrift   =(param->GetZLength()-startTbin*param->GetZWidth());
 274                 zIonDrift+=shiftZ;
 275                 Double_t signal=sig[i];
 276                 if (useGainMap) signal/=gainPad;
 277                 if (TMath::Abs(zIonDrift)<param->GetZLength()){
 278                   if ((sector%36)>=18) zIonDrift*=-1;   // c side has opposite sign
 279                   if (sector%36<18) hisQ1D[ith]->Fill(localX, signal/padLength);
 280                   hisQ2D[ith]->Fill(phi,localX, signal/padLength);
 281                   hisQ2DRZ[ith]->Fill(localX, zIonDrift, signal/padLength);
 282                   hisQ3D[ith]->Fill(phi,localX,zIonDrift,signal/padLength);
 283                 }
 284                 //
 285                 Double_t zIonROC = ((sector%36)<18)? shiftZ: -shiftZ;  // z position of the "ion disc" -  A side C side opposite sign
 286                 if (sector%36<18) hisQ1DROC[ith]->Fill(localX, signal/padLength);
 287                 hisQ2DROC[ith]->Fill(phi,localX, signal/padLength);
 288                 hisQ2DRZROC[ith]->Fill(localX, zIonROC, signal/padLength);
 289                 hisQ3DROC[ith]->Fill(phi,localX,zIonROC,signal/padLength);
 290               }
 291             }
 292           }
 293         }
 294       }
 295     }
 296   }
 297   timer.Print();
 298   delete pcstream;
 299   return 0;
 300 }
 301
 302
 303 void AnalyzeMaps1D(){
 304   //
 305   // Analyze space charge maps stored as shitograms in trees
 306   //
 307   TFile *  fhisto = new TFile("histo.root","update");
 308   TTree * tree = (TTree*)fhisto->Get("histo");
 309   if (!tree){
 310     TChain *chain = AliXRDPROOFtoolkit::MakeChainRandom("histo.list","histo",0,1000,1);
 311     tree = chain->CopyTree("1");
 312     tree->Write("histo");
 313   }
 314   //
 315   //
 316   //
 317   TH1 *his1Th[3]={0,0,0};
 318   TF1 *fq1DStep= new TF1("fq1DStep","([0]+[1]*(x>134))/x**min(abs([2]),3)",85,245);
 319   fq1DStep->SetParameters(1,-0.5,1);
 320   tree->Draw("hist1DROC_2.fArray:hist1D_2.fXaxis.fXbins.fArray>>his(40,85,245)","","prof");
 321   tree->GetHistogram()->Fit(fq1DStep);
 322   // normalize step between the IROC-OROC
 323   tree->SetAlias("normQ",Form("(1+%f*(hist1D_2.fXaxis.fXbins.fArray>136))",fq1DStep->GetParameter(1)/fq1DStep->GetParameter(0)));
 324   //
 325   {
 326     Int_t entries= tree->Draw("hist1DROC_2.fArray/(events*normQ)","1","goff");
 327     Double_t median=TMath::Median(entries,tree->GetV1());
 328     TCut cut10Median = Form("hist1DROC_2.fArray/(events*normQ)<%f",10*median);
 329     //
 330     tree->Draw("hist1DROC_2.fArray/(events*normQ):hist1D_2.fXaxis.fXbins.fArray>>his1Th0(40,86,245)",cut10Median+"","prof");
 331     his1Th[0] = tree->GetHistogram();
 332     tree->Draw("hist1DROC_4.fArray/(events*normQ):hist1D_2.fXaxis.fXbins.fArray>>his1Th1(40,86,245)",cut10Median+"","prof");
 333     his1Th[1] = tree->GetHistogram();
 334     tree->Draw("hist1DROC_6.fArray/(events*normQ):hist1D_2.fXaxis.fXbins.fArray>>his1Th2(40,86,245)",cut10Median+"","prof");
 335     his1Th[2]=tree->GetHistogram();
 336   }
 337   //
 338   TCanvas *canvasR = new TCanvas("canvasR","canvasR",600,500);
 339   canvasR->cd();
 340   for (Int_t i=0; i<3; i++){
 341     his1Th[i]->SetMarkerStyle(21);
 342     his1Th[i]->SetMarkerColor(i+2);
 343     fq1DStep->SetLineColor(i+2);
 344     his1Th[i]->Fit(fq1DStep,"","");
 345     his1Th[i]->GetXaxis()->SetTitle("r (cm)");
 346     his1Th[i]->GetYaxis()->SetTitle("#frac{N_{el}}{N_{ev}}(ADC/cm)");
 347   }
 348   TLegend * legend  = new TLegend(0.11,0.11,0.7,0.39,"1D space Charge map (ROC part) (z,phi integrated)");
 349   for (Int_t i=0; i<3; i++){
 350     his1Th[i]->SetMinimum(0);fq1DStep->SetLineColor(i+2);
 351     his1Th[i]->Fit(fq1DStep,"qnr","qnr");
 352     if (i==0) his1Th[i]->Draw("");
 353     his1Th[i]->Draw("same");
 354     legend->AddEntry(his1Th[i],Form("Thr=%d Slope=%2.2f",2*i+2,fq1DStep->GetParameter(2)));
 355   }
 356   legend->Draw();
 357   legend->SaveAs("spaceCharge1d.png");
 358   //
 359   //
 360   //
 361 }
 362
 363 void DoMerge(){
 364   //
 365   //
 366   //
 367   TFile *  fhisto = new TFile("histo.root","recreate");
 368   TTree * tree = 0;
 369   TChain *chain = AliXRDPROOFtoolkit::MakeChainRandom("histo.list","histo",0,50,1);
 370   tree = chain->CopyTree("1");
 371   tree->Write("histo");
 372   delte fhist0;
 373 }
 374
 375
 376
 377 void MakeFluctuationStudy3D(Int_t nhistos, Int_t nevents, Int_t niter){
 378   //
 379   //
 380   //
 381   // Fix z binning before creating of official plot (PbPb  data)
 382   //
 383   // 1. Make a summary integral   3D/2D/1D maps
 384   // 2. Create several maps with niter events  - Poisson flucturation in n
 385   // 3. Store results 3D maps in the tree (and also as histogram)  current and mean
 386   //
 387   TTreeSRedirector * pcstream  = new TTreeSRedirector("histo.root", "update");
 388   TFile *  fhisto = pcstream->GetFile();
 389
 390   TTree * tree = (TTree*)fhisto->Get("histo");
 391   tree->SetCacheSize(10000000000);
 392
 393   TH1D * his1DROC=0,    * his1DROCSum=0,  * his1DROCN=0;
 394   TH1D * his1DDrift=0,  * his1DDriftSum=0, * his1DDriftN=0 ;
 395   TH2D * his2DROC=0,    * his2DROCSum=0,  * his2DROCN=0;
 396   TH2D * his2DRZROC=0,    * his2DRZROCSum=0,  * his2DRZROCN=0;
 397   TH2D * his2DDrift=0,  * his2DDriftSum=0, * his2DDriftN=0;
 398   TH2D * his2DRZDrift=0,  * his2DRZDriftSum=0, * his2DRZDriftN=0;
 399   TH3D * his3DROC=0,    * his3DROCSum=0,  * his3DROCN=0;
 400   TH3D * his3DDrift=0,  * his3DDriftSum=0, * his3DDriftN=0;
 401   //
 402   if (nhistos<0 || nhistos> tree->GetEntries()) nhistos = tree->GetEntries();
 403   Int_t  eventsPerChunk=0;
 404   tree->SetBranchAddress("hist1D_2.",&his1DDrift);
 405   tree->SetBranchAddress("hist1DROC_2.",&his1DROC);
 406   tree->SetBranchAddress("hist2D_2.",&his2DDrift);
 407   tree->SetBranchAddress("hist2DRZ_2.",&his2DRZDrift);
 408   tree->SetBranchAddress("hist2DROC_2.",&his2DROC);
 409   tree->SetBranchAddress("hist3D_2.",&his3DDrift);
 410   tree->SetBranchAddress("hist3DROC_2.",&his3DROC);
 411   tree->SetBranchAddress("hist2DRZROC_2.",&his2DRZROC);
 412   tree->SetBranchAddress("events",&eventsPerChunk);
 413   //
 414   // 1. Make a summary integral   3D/2D/1D maps
 415   //
 416   Int_t neventsAll=0;
 417   for (Int_t i=0; i<nhistos; i++){
 418     tree->GetEntry(i);
 419     if (i%25==0) printf("%d\n",i);
 420     if (his1DROCSum==0)     his1DROCSum=new TH1D(*his1DROC);
 421     if (his1DDriftSum==0)   his1DDriftSum=new TH1D(*his1DDrift);
 422     if (his2DROCSum==0)     his2DROCSum=new TH2D(*his2DROC);
 423     if (his2DRZROCSum==0)     his2DRZROCSum=new TH2D(*his2DRZROC);
 424     if (his2DDriftSum==0)   his2DDriftSum=new TH2D(*his2DDrift);
 425     if (his2DRZDriftSum==0)   his2DRZDriftSum=new TH2D(*his2DRZDrift);
 426     if (his3DROCSum==0)     his3DROCSum=new TH3D(*his3DROC);
 427     if (his3DDriftSum==0)   his3DDriftSum=new TH3D(*his3DDrift);
 428     his1DROCSum->Add(his1DROC);
 429     his1DDriftSum->Add(his1DDrift);
 430     his2DROCSum->Add(his2DROC);
 431     his2DRZROCSum->Add(his2DRZROC);
 432     his2DDriftSum->Add(his2DDrift);
 433     his2DRZDriftSum->Add(his2DRZDrift);
 434     his3DROCSum->Add(his3DROC);
 435     his3DDriftSum->Add(his3DDrift);
 436     neventsAll+=eventsPerChunk;
 437   }
 438   //
 439   // 2. Create several maps with niter events  - Poisson flucturation in n
 440   //
 441   for (Int_t iter=0; iter<niter; iter++){
 442     printf("Itteration=\t%d\n",iter);
 443     Int_t ilast=gRandom->Rndm()*nhistos;
 444     Int_t nchunks=gRandom->Poisson(nevents)/eventsPerChunk;  // chunks with n typically 25 events
 445     for (Int_t i=0; i<nchunks; i++){
 446       //    ilast+=gRandom->Rndm()
 447       tree->GetEntry(gRandom->Rndm()*nhistos);
 448       if (i%10==0) printf("%d\t%d\n",iter, i);
 449       if (his1DROCN==0)     his1DROCN=new TH1D(*his1DROC);
 450       if (his1DDriftN==0)   his1DDriftN=new TH1D(*his1DDrift);
 451       if (his2DROCN==0)     his2DROCN=new TH2D(*his2DROC);
 452       if (his2DDriftN==0)   his2DDriftN=new TH2D(*his2DDrift);
 453       if (his2DRZROCN==0)     his2DRZROCN=new TH2D(*his2DRZROC);
 454       if (his2DRZDriftN==0)   his2DRZDriftN=new TH2D(*his2DRZDrift);
 455       if (his3DROCN==0)     his3DROCN=new TH3D(*his3DROC);
 456       if (his3DDriftN==0)   his3DDriftN=new TH3D(*his3DDrift);
 457       his1DROCN->Add(his1DROC);
 458       his1DDriftN->Add(his1DDrift);
 459       his2DROCN->Add(his2DROC);
 460       his2DRZDriftN->Add(his2DRZDrift);
 461       his2DRZROCN->Add(his2DRZROC);
 462       his2DDriftN->Add(his2DDrift);
 463       his3DROCN->Add(his3DROC);
 464       his3DDriftN->Add(his3DDrift);
 465     }
 466     //
 467     // 3. Store results 3D maps in the tree (and also as histogram)  current and mea
 468     //
 469     Int_t eventsUsed=  nchunks*eventsPerChunk;
 470     (*pcstream)<<"fluctuation"<<
 471       "neventsAll="<<neventsAll<<   // total number of event to define mean
 472       "nmean="<<nevents<<         // mean number of events used
 473       "eventsUsed="<<eventsUsed<<         // number of chunks used for one fluct. study
 474       //
 475       // 1,2,3D histogram per group and total
 476       "his1DROCN.="<<his1DROCN<<
 477       "his1DROCSum.="<<his1DROCSum<<
 478       "his1DDriftN.="<<his1DDriftN<<
 479       "his1DDriftSum.="<<his1DDriftSum<<
 480       "his2DROCN.="<<his2DROCN<<
 481       "his2DROCSum.="<<his2DROCSum<<
 482       "his2DDriftN.="<<his2DDriftN<<
 483       "his2DDriftSum.="<<his2DDriftSum<<
 484       "his2DRZROCN.="<<his2DRZROCN<<
 485       "his2DRZROCSum.="<<his2DRZROCSum<<
 486       "his2DRZDriftN.="<<his2DRZDriftN<<
 487       "his2DRZDriftSum.="<<his2DRZDriftSum<<
 488       "his3DROCN.="<<his3DROCN<<
 489       "his3DROCSum.="<<his3DROCSum<<
 490       "his3DDriftN.="<<his3DDriftN<<
 491       "his3DDriftSum.="<<his3DDriftSum<<
 492       "\n";
 493     pcstream->GetFile()->mkdir(Form("Fluc%d",iter));
 494     pcstream->GetFile()->cd(Form("Fluc%d",iter));
 495     his2DRZROCN->Write("his2DRZROCN");
 496     his2DRZROCSum->Write("his2DRZROCSum");
 497     his2DRZDriftN->Write("his2DRZDriftN");
 498     his2DRZDriftSum->Write("his2DRZDriftSum");
 499     //
 500     his3DROCN->Write("his3DROCN");
 501     his3DROCSum->Write("his3DROCSum");
 502     his3DDriftN->Write("his3DDriftN");
 503     his3DDriftSum->Write("his3DDriftSum");
 504
 505     his1DROCN->Reset();
 506     his1DDriftN->Reset();
 507     his2DROCN->Reset();
 508     his2DRZDriftN->Reset();
 509     his2DRZROCN->Reset();
 510     his2DDriftN->Reset();
 511     his3DROCN->Reset();
 512     his3DDriftN->Reset();
 513   }
 514
 515   delete pcstream;
 516
 517
 518   //
 519   Double_t mean,rms;
 520   //
 521   mean=TMath::Median(his2DROCSum->GetNbinsX()*his2DROCSum->GetNbinsY(),his2DROCSum->GetArray());
 522   rms=TMath::Median(his2DROCSum->GetNbinsX()*his2DROCSum->GetNbinsY(),his2DROCSum->GetArray());
 523   his2DROCSum->SetMaximum(mean+4*rms);
 524   his2DROCSum->SetMinimum(mean-4*rms);
 525
 526 }
 527
 528
 529 void DrawDCARPhiTrendTime(){
 530   //
 531   // Macros to draw the DCA correlation with the luminosity (estimated from the occupancy)
 532   //
 533   // A side and c side  0 differnt behaviour -
 534   // A side - space charge effect
 535   // C side - space charge effect+ FC charging:
 536   //   Variables  to query from the QA/calibration DB - tree:
 537   //   QA.TPC.CPass1.dcar_posA_0   -dca rphi in cm - offset
 538   //   Calib.TPC.occQA.Sum()       - luminosity is estimated using the mean occupancy per run
 539   //
 540   TFile *fdb = TFile::Open("outAll.root");
 541   if (!fdb)  fdb = TFile::Open("http://www-alice.gsi.de/TPC/CPassMonitor/outAll.root");
 542   TTree * tree = (TTree*)fdb->Get("joinAll");
 543   tree->SetCacheSize(100000000);
 544   tree->SetMarkerStyle(25);
 545
 546   //QA.TPC.CPass1.dcar_posA_0 QA.TPC.CPass1.dcar_posA_0_Err  QA.TPC.CPass1.meanMult  Calib.TPC.occQA.  DAQ.L3_magnetCurrent
 547
 548   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);
 549   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);
 550   Double_t mean,rms;
 551   TStatToolkit::EvaluateUni(grA->GetN(),grA->GetY(), mean,rms,grA->GetN()*0.8);
 552   grA->SetMinimum(mean-5*rms);
 553   grA->SetMaximum(mean+3*rms);
 554
 555
 556   grA->GetXaxis()->SetTitle("occ*sign(bz)");
 557   grA->GetYaxis()->SetTitle("#Delta_{r#phi} (cm)");
 558   grA->Draw("ap");
 559   grC->Draw("p");
 560   TLegend* legend = new TLegend(0.11,0.11,0.5,0.3,"DCA_{rphi} as function of IR (2013)" );
 561   legend->AddEntry(grA,"A side","p");
 562   legend->AddEntry(grC,"C side","p");
 563   legend->Draw();
 564 }
 565
 566
 567
 568 void DrawOpenGate(){
 569   //
 570   //  Make nice plot to demonstrate the space charge effect in run with the open gating grid
 571   //  For the moment the inmput is harwired - the CPass0 calibration data used
 572   //  Make nice drawing (with axis labels):
 573   //  To fix (longer term)
 574   //     the distortion map to be recalculated - using gaussian fit (currently we use mean)
 575   //     the histogram should be extended
 576   TFile f("/hera/alice/alien/alice/data/2013/LHC13g/000197470/cpass0/OCDB/root_archive.zip#meanITSVertex.root");
 577   TFile fref("/hera/alice/alien/alice/data/2013/LHC13g/000197584/cpass0/OCDB/root_archive.zip#meanITSVertex.root");
 578   //
 579   TTree * treeTOFdy=(TTree*)f.Get("TOFdy");
 580   TTree * treeTOFdyRef=(TTree*)fref.Get("TOFdy");
 581   treeTOFdy->AddFriend(treeTOFdyRef,"R");
 582   treeTOFdy->SetMarkerStyle(25);
 583   TTree * treeITSdy=(TTree*)f.Get("ITSdy");
 584   TTree * treeITSdyRef=(TTree*)fref.Get("ITSdy");
 585   treeITSdy->AddFriend(treeITSdyRef,"R");
 586   treeITSdy->SetMarkerStyle(25);
 587   TTree * treeVertexdy=(TTree*)f.Get("Vertexdy");
 588   TTree * treeVertexdyRef=(TTree*)fref.Get("Vertexdy");
 589   treeVertexdy->AddFriend(treeVertexdyRef,"R");
 590   treeVertexdy->SetMarkerStyle(25);
 591
 592   //  treeITSdy->Draw("mean-R.mean:sector:abs(theta)","entries>50&&abs(snp)<0.1&&theta<0","colz")
 593
 594   treeITSdy->Draw("mean-R.mean:sector:abs(theta)","entries>50&&abs(snp)<0.1&&theta>0","colz");
 595
 596
 597 }