PWG1/macros/VtxInfoLHC.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <Riostream.h>
   3 #include <TCanvas.h>
   4 #include <TLegend.h>
   5 #include <TLegendEntry.h>
   6 #include <TFile.h>
   7 #include <TLine.h>
   8 #include <TNtuple.h>
   9 #include <TStyle.h>
  10 #include <TString.h>
  11 #include <TSystem.h>
  12 #include <TH1F.h>
  13 #include <TGraphErrors.h>
  14 #include <TF1.h>
  15 #include <TH2F.h>
  16 #include <fstream.h>
  17 #include <TMath.h>
  18 #include <TProfile.h>
  19
  20 #endif
  21
  22 /* Macro to produce LHC file for the luminous region.
  23    Author: caffarri@pd.infn.it
  24 */
  25
  26
  27 void VtxInfoLHC(TString vtxtype="TRKnc",
  28                 TString fname="Vertex.Performance.root",
  29                 TString ntname="fNtupleBeamSpot"){
  30
  31   // const Int_t nbinsmult=11;
  32   // Float_t limmult[nbinsmult+1]={0.,1.5, 2.5, 5.5,10.5,14.5,24.5,30.5,44.5,60.5,74.5, 99999.};
  33
  34   Float_t totev=0,totevtriggered=0;
  35
  36   Float_t resolVtx=510.;
  37   Float_t errResolVtx = 10;
  38   Float_t p2= 1.3;
  39   Float_t errp2 = 0.1;
  40   Float_t meanMult=38;
  41   Float_t errTotResol =(errResolVtx*errResolVtx)+(0.3*(errp2*errp2)/TMath::Power(meanMult,p2+1));
  42
  43   Float_t timeBinning = 120.;
  44
  45
  46   Char_t hisnam[20];
  47   Double_t xRange = 1.;
  48   Int_t xBins = 125.;
  49   Double_t zRange = 40.;
  50   Int_t zBins = 40.;
  51
  52   /*
  53     Int_t binsHRes=250*0.5;
  54     Float_t rangeHRes=10000.*0.5;
  55     Float_t rangeGrRms = 5000; // micron
  56     Float_t rangeGrPull = 15;
  57     Float_t rangeHPull=10.;
  58     Float_t rangeHResS = 1.;
  59     Float_t rangeHResS = 10000*0.5.;
  60   */
  61
  62   TFile *f=new TFile(fname.Data());
  63   TList *cOutput = (TList*)f->Get("cOutputVtxESD");
  64   TNtuple *nt=(TNtuple*)cOutput->FindObject(ntname.Data());
  65   // TNtuple *nt=(TNtuple*)f->FindObject(ntname.Data());
  66   Int_t nnnev=nt->GetEntries();
  67   printf("Events = %d\n",nnnev);
  68
  69   Float_t xVtx,xVtxSc,xerrVtx;
  70   Float_t yVtx,yVtxSc,yerrVtx;
  71   Float_t zVtx,zVtxSc,zerrVtx;
  72   Float_t ntrklets,ncontrVtx,triggered;
  73   Float_t run, tstamp;
  74
  75   TString sxx="x"; sxx.Append(vtxtype.Data());
  76   nt->SetBranchAddress(sxx.Data(),&xVtx);
  77   TString syy="y"; syy.Append(vtxtype.Data());
  78   nt->SetBranchAddress(syy.Data(),&yVtx);
  79   TString szz="z"; szz.Append(vtxtype.Data());
  80   nt->SetBranchAddress(szz.Data(),&zVtx);
  81
  82   TString xerr="xerr"; xerr.Append(vtxtype.Data());
  83   nt->SetBranchAddress(xerr.Data(),&xerrVtx);
  84   TString yerr="yerr"; yerr.Append(vtxtype.Data());
  85   nt->SetBranchAddress(yerr.Data(),&yerrVtx);
  86   TString zerr="zerr"; zerr.Append(vtxtype.Data());
  87   nt->SetBranchAddress(zerr.Data(),&zerrVtx);
  88
  89
  90   TCanvas *stampCanvas = new TCanvas ("stampCanvas", "stampCanvas");
  91   stampCanvas->Divide(3);
  92   TCanvas *stampCanvasSig = new TCanvas ("stampCanvasSig", "stampCanvasSig");
  93   stampCanvasSig->Divide(3);
  94   TCanvas *diamCanvas = new TCanvas ("diamCanvas", "diamCanvas");
  95   diamCanvas->Divide(2);
  96
  97   TString trkstitle;
  98   nt->SetBranchAddress("ntrklets",&ntrklets);
  99   trkstitle="Tracklets Multiplicity";
 100
 101   TString ntrks="ntrks"; ntrks.Append(vtxtype.Data());
 102   nt->SetBranchAddress(ntrks.Data(),&ncontrVtx);
 103   // nt->SetBranchAddress(ntrks.Data(),&ntrklets);
 104   nt->SetBranchAddress("triggered",&triggered);
 105   nt->SetBranchAddress("run", &run);
 106   nt->SetBranchAddress("tstamp",&tstamp);
 107
 108   Float_t minTstamp=10E+13;
 109   Float_t maxTstamp=-1;
 110   Int_t subevents;
 111
 112   Float_t minRun = 10E+7;
 113   Float_t maxRun = 1;
 114
 115   for (Int_t ientries=0; ientries<nt->GetEntriesFast(); ientries++){
 116     nt->GetEntry(ientries);
 117
 118     if (tstamp<minTstamp) minTstamp=tstamp;
 119     if (tstamp>maxTstamp) maxTstamp=tstamp;
 120
 121     if (run<minRun) minRun=run-10;
 122     if (run>maxRun) maxRun=run+10;
 123   }
 124
 125   Int_t nTimeBins = (Int_t)(maxTstamp-minTstamp) / timeBinning;
 126
 127   printf ("number of time bins %d\n", nTimeBins);
 128   subevents=0;
 129
 130   TH1F **hxVtx = new TH1F*[nTimeBins+1]; //"hxVtx", "hxVtx", 125,-1,1);
 131   TH1F **hyVtx = new TH1F*[nTimeBins+1]; //("hyVtx", "hyVtx", 125,-1,1);
 132   TH1F **hzVtx = new TH1F*[nTimeBins+1]; //("hzVtx", "hzVtx", 40,-40,40);
 133
 134   TH1F **hxVtxMult = new TH1F*[nTimeBins+1]; //
 135   TH1F **hyVtxMult = new TH1F*[nTimeBins+1]; //
 136
 137   for(Int_t kTimeBin=0;kTimeBin<nTimeBins+1;kTimeBin++){
 138     printf ("here we are, kTimeBin = %d\n", kTimeBin);
 139
 140     sprintf(hisnam,"hxVtx%d",kTimeBin);
 141     hxVtx[kTimeBin]=new TH1F(hisnam,"",xBins,-xRange,xRange);
 142     sprintf(hisnam,"hyVtx%d",kTimeBin);
 143     hyVtx[kTimeBin]=new TH1F(hisnam,"",xBins,-xRange,xRange);
 144     sprintf(hisnam,"hzVtx%d",kTimeBin);
 145     hzVtx[kTimeBin]=new TH1F(hisnam,"",zBins,-zRange,zRange);
 146     sprintf(hisnam,"hxVtxMult%d",kTimeBin);
 147     hxVtxMult[kTimeBin]=new TH1F(hisnam,"",xBins,-xRange,xRange);
 148     sprintf(hisnam,"hyVtxMult%d",kTimeBin);
 149     hyVtxMult[kTimeBin]=new TH1F(hisnam,"",xBins,-xRange,xRange);
 150
 151   }
 152
 153   TH1F *htstampX = new TH1F("tstamp Vx","tstamp Vx", nTimeBins+1, minTstamp, maxTstamp);
 154   TH1F *htstampY = new TH1F("tstamp Vy","tstamp Vy", nTimeBins+1, minTstamp, maxTstamp);
 155   TH1F *htstampZ = new TH1F("tstamp Vz","tstamp Vz", nTimeBins+1, minTstamp, maxTstamp);
 156
 157   TH1F *htstampSigX = new TH1F("tstamp Sig x","tstamp Sig x", nTimeBins+1, minTstamp, maxTstamp);
 158   TH1F *htstampSigY = new TH1F("tstamp Sig y","tstamp Sig y", nTimeBins+1, minTstamp, maxTstamp);
 159   TH1F *htstampSigZ = new TH1F("tstamp Sig z","tstamp Sig z", nTimeBins+1, minTstamp, maxTstamp);
 160
 161   TH1F *htstampXdiam = new TH1F("tstamp Diam X","tstamp Diam X", nTimeBins+1, minTstamp, maxTstamp);
 162   TH1F *htstampYdiam = new TH1F("tstamp Diam Y","tstamp Diam Y", nTimeBins+1, minTstamp, maxTstamp);
 163
 164   TCanvas *cVtx = new TCanvas("cVtx", "Vtx Canvas", 1000,700);
 165   cVtx->Divide(3);
 166
 167   Int_t *fill = new Int_t [nTimeBins+1];
 168
 169   for(Int_t iev=0;iev<nnnev;iev++){
 170     nt->GetEvent(iev);
 171
 172     Int_t rightTimeBin =(Int_t)(tstamp-minTstamp)/ 120;
 173     printf ("right Time Bin %d \n", rightTimeBin);
 174
 175     if ((run>=114783) && (run<=114798)) fill[rightTimeBin]=1005;
 176     if ((run>=114916) && (run<=114931)) fill[rightTimeBin]=1013;
 177     if ((run>=115186) && (run<=115193)) fill[rightTimeBin]=1019;
 178     if ((run>=115310) && (run<=115345)) fill[rightTimeBin]=1022;
 179     if ((run>=115393) && (run<=115414)) fill[rightTimeBin]=1023;
 180     if ((run>=115514) && (run<=115521)) fill[rightTimeBin]=1026;
 181     if ((run>=115880) && (run<=115892)) fill[rightTimeBin]=1031;
 182     if ((run>=116079) && (run<=116081)) fill[rightTimeBin]=1033;
 183     if ((run>=116102) && (run<=116134)) fill[rightTimeBin]=1034;
 184     if ((run>=116197) && (run<=116204)) fill[rightTimeBin]=1035;
 185     if ((run>=116287) && (run<=116288)) fill[rightTimeBin]=1038;
 186     if ((run>=116351) && (run<=116354)) fill[rightTimeBin]=1042;
 187     if ((run>=116401) && (run<=116403)) fill[rightTimeBin]=1044;
 188     if ((run>=116559) && (run<=116574)) fill[rightTimeBin]=1045;
 189     if ((run>=116609) && (run<=116611)) fill[rightTimeBin]=1046;
 190     if ((run>=116640) && (run<=116645)) fill[rightTimeBin]=1047;
 191     if ((run>=116681) && (run<=116684)) fill[rightTimeBin]=1049;
 192
 193
 194     //      if((tstamp<stopTime)&&(tstamp>=iTstamp)){
 195     subevents++;
 196     // ncontrVtx = ntrklets;
 197
 198     //printf ("ntrksTRKnc = %d ncontr = %d \n ", ntrklets, ncontrVtx);
 199
 200     xVtxSc = 10000.*xVtx;
 201     yVtxSc = 10000.*yVtx;
 202     zVtxSc = 10000.*zVtx;
 203
 204     if (ncontrVtx>0){
 205       hxVtx[rightTimeBin]->Fill(xVtx);
 206       hyVtx[rightTimeBin]->Fill(yVtx);
 207       hzVtx[rightTimeBin]->Fill(zVtx);
 208
 209       if((ntrklets>30)&&(ntrklets<45)){
 210         hxVtxMult[rightTimeBin]->Fill(xVtx);
 211         hyVtxMult[rightTimeBin]->Fill(yVtx);
 212       }
 213     }
 214
 215   }
 216
 217   Double_t xMeanVtx, xMeanVtxErr, xSigmaVtx, xSigmaVtxErr;
 218   Double_t yMeanVtx, yMeanVtxErr, ySigmaVtx, ySigmaVtxErr;
 219   Double_t zMeanVtx, zMeanVtxErr, zSigmaVtx, zSigmaVtxErr;
 220   Double_t xVtxRms, yVtxRms, zVtxRms;
 221   Double_t xVtxMult, yVtxMult;
 222   Double_t xDiam,  xDiamErr, yDiam,  yDiamErr, zDiam, zDiamErr;
 223
 224   /*
 225     TString outFileName ="aliceDataVertex";
 226     TString outFileExt = ".txt";
 227     TString outFileFill ;
 228   */
 229
 230   ofstream outFile("1045_lumireg_ALICE.txt");
 231
 232   for (Int_t iTimeBin=0; iTimeBin<nTimeBins+1; iTimeBin++){
 233
 234     if (fill[iTimeBin] != 1045) continue;
 235
 236     Int_t corrTstamp = iTimeBin*120+minTstamp;
 237
 238     if (iTimeBin == nTimeBins+1) {
 239       cVtx->cd(1);
 240       hxVtxMult[iTimeBin]->Draw();
 241       cVtx->cd(2);
 242       hyVtxMult[iTimeBin]->Draw();
 243
 244     }
 245
 246     TF1 *fitVtx, *fitVtxMult;
 247     if((hxVtx[iTimeBin]->GetEffectiveEntries()<500)){// && (hxVtx[iTimeBin]->GetEffectiveEntries()<500)){
 248       //outFile << "0." <<" " << "0."<<" ";
 249       xMeanVtx =0.;
 250       xMeanVtxErr = 0.;
 251       xSigmaVtx = 0.;
 252       xSigmaVtxErr = 0.;
 253       xDiam = 0.;
 254       xDiamErr =0.;
 255       xVtxRms =0.;
 256     }
 257     else{
 258       hxVtx[iTimeBin]->Fit("gaus", "M");
 259       fitVtx = hxVtx[iTimeBin]->GetFunction("gaus");
 260       xMeanVtx = fitVtx->GetParameter(1);
 261       xMeanVtxErr = fitVtx->GetParError(1);
 262       xSigmaVtx = fitVtx->GetParameter(2);
 263       xSigmaVtxErr = fitVtx->GetParError(2);
 264       xVtxRms = hxVtx[iTimeBin]->GetRMS();
 265
 266       hxVtxMult[iTimeBin]->Fit("gaus", "M");
 267       fitVtxMult= hxVtxMult[iTimeBin]->GetFunction("gaus");
 268       xVtxMult = fitVtxMult->GetParameter(2)*10000;
 269       xDiam=TMath::Sqrt(xVtxMult*xVtxMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
 270       if(xDiam >= 500 ) xDiam=0.;
 271       xDiamErr = TMath::Sqrt(fitVtxMult->GetParError(2)*10000+ errTotResol);
 272       printf ("xDiam = %f errDiam = %f \n", xDiam, xDiamErr);
 273
 274       }
 275
 276
 277     if((hyVtx[iTimeBin] ->GetEffectiveEntries()<500)){// && (hyVtxMult[iTimeBin]->GetEffectiveEntries()<500)){
 278       //outFile << "0." <<" " << "0."<<" ";
 279       yMeanVtx=0.;
 280       yMeanVtxErr = 0.;
 281       ySigmaVtx = 0.;
 282       ySigmaVtxErr = 0.;
 283       yDiam = 0;
 284       yDiamErr =0;
 285       yVtxRms=0.;
 286     }
 287     else{
 288       hyVtx[iTimeBin]->Fit("gaus", "M");
 289       fitVtx = hyVtx[iTimeBin]->GetFunction("gaus");
 290       yMeanVtx = fitVtx->GetParameter(1);
 291       yMeanVtxErr = fitVtx->GetParError(1);
 292       ySigmaVtx = fitVtx->GetParameter(2);
 293       ySigmaVtxErr = fitVtx->GetParError(2);
 294       yVtxRms = hyVtx[iTimeBin]->GetRMS();
 295
 296       hyVtxMult[iTimeBin]->Fit("gaus", "M");
 297       fitVtxMult= hyVtxMult[iTimeBin]->GetFunction("gaus");
 298       yVtxMult = fitVtxMult->GetParameter(2)*10000;
 299       yDiam=TMath::Sqrt(yVtxMult*yVtxMult - ((resolVtx*resolVtx)/TMath::Power(meanMult, p2)));
 300       if(yDiam >= 500 ) yDiam=0.;
 301       yDiamErr = TMath::Sqrt(fitVtxMult->GetParError(2)*10000 + errTotResol);
 302       printf("yDiam = %f  errDiam= %f \n", yDiam, yDiamErr);
 303     }
 304
 305     if((hzVtx[iTimeBin]->GetEffectiveEntries()<500)){
 306       //  outFile << "0." <<" " << "0."<<" ";
 307       zMeanVtx=0.;
 308       zMeanVtxErr = 0.;
 309       zSigmaVtx = 0.;
 310       zSigmaVtxErr = 0.;
 311       zDiam = 0.;
 312       zDiamErr =0.;
 313       zVtxRms =0.;
 314     }
 315     else{
 316       hzVtx[iTimeBin]->Fit("gaus", "M");
 317       fitVtx = hzVtx[iTimeBin]->GetFunction("gaus");
 318       zMeanVtx = fitVtx->GetParameter(1);
 319       zMeanVtxErr = fitVtx->GetParError(1);
 320       zSigmaVtx = fitVtx->GetParameter(2);
 321       zSigmaVtxErr = fitVtx->GetParError(2);
 322       zVtxRms = hzVtx[iTimeBin]->GetRMS();
 323
 324       zDiam = zSigmaVtx;
 325       zDiamErr = zSigmaVtxErr;
 326
 327
 328     }
 329
 330     //outfile<<
 331
 332     Int_t iBin = htstampX->FindBin(corrTstamp);
 333     htstampX->SetBinContent(iBin, xMeanVtx);
 334     htstampX->SetBinError(iBin, xMeanVtxErr);
 335
 336     iBin = htstampY->FindBin(corrTstamp);
 337     htstampY->SetBinContent(iBin, yMeanVtx);
 338     htstampY->SetBinError(iBin, yMeanVtxErr);
 339
 340     iBin = htstampZ->FindBin(corrTstamp);
 341     htstampZ->SetBinContent(iBin, zMeanVtx);
 342     htstampZ->SetBinError(iBin, zMeanVtxErr);
 343
 344     iBin = htstampSigX->FindBin(corrTstamp);
 345     htstampSigX->SetBinContent(iBin, xSigmaVtx);
 346     htstampSigX->SetBinError(iBin, xSigmaVtxErr);
 347
 348     iBin = htstampSigY->FindBin(corrTstamp);
 349     htstampSigY->SetBinContent(iBin, ySigmaVtx);
 350     htstampSigY->SetBinError(iBin, ySigmaVtxErr);
 351
 352     iBin = htstampSigZ->FindBin(corrTstamp);
 353     htstampSigZ->SetBinContent(iBin, zSigmaVtx);
 354     htstampSigZ->SetBinError(iBin, zSigmaVtxErr);
 355
 356     iBin = htstampXdiam->FindBin(corrTstamp);
 357     htstampXdiam->SetBinContent(iBin, xDiam);
 358     htstampXdiam->SetBinError(iBin, xDiamErr);
 359
 360     iBin = htstampYdiam->FindBin(corrTstamp);
 361     htstampYdiam->SetBinContent(iBin,yDiam);
 362     htstampYdiam->SetBinError(iBin, yDiamErr);
 363
 364     outFile <<"2"<<" "<<fill[iTimeBin]<<" "<<corrTstamp<<" ";
 365     outFile <<xMeanVtx*10<<" "<<xMeanVtxErr*10<<" ";
 366     outFile <<yMeanVtx*10<<" "<<yMeanVtxErr*10<<" ";
 367     outFile <<zMeanVtx*10<<" "<<zMeanVtxErr*10<<" ";
 368     outFile <<xSigmaVtx*10<<" "<<xSigmaVtxErr*10<<" ";
 369     outFile <<ySigmaVtx*10<<" "<<ySigmaVtxErr*10<<" ";
 370     outFile <<zSigmaVtx*10<<" "<<zSigmaVtxErr*10<<" ";
 371     outFile <<xDiam*0.001<<" "<<xDiamErr*0.001<<" ";
 372     outFile <<yDiam*0.001<<" "<<yDiamErr*0.001<<" ";
 373     outFile <<zDiam*10<<" "<<zDiamErr*10<<" ";
 374     outFile <<xVtxRms*10 <<" "<<yVtxRms*10<<" "<<zVtxRms*10<<endl;
 375
 376
 377
 378
 379   }
 380   printf ("Err Resol = %f \n", errTotResol);
 381   stampCanvas->cd(1);
 382   htstampX->Draw();
 383   stampCanvas->cd(2);
 384   htstampY->Draw();
 385   stampCanvas->cd(3);
 386   htstampZ->Draw();
 387
 388   stampCanvasSig->cd(1);
 389   htstampSigX->Draw();
 390   stampCanvasSig->cd(2);
 391   htstampSigY->Draw();
 392   stampCanvasSig->cd(3);
 393   htstampSigZ->Draw();
 394
 395   diamCanvas->cd(1);
 396   htstampXdiam->Draw();
 397   diamCanvas->cd(2);
 398   htstampYdiam->Draw();
 399 }