PWG3/vertexingHF/TestMultiVector.C

   1 void TestMultiVector(){
   2
   3   // Example of usage of AliMultiDimVector and AliSignificanceCalculator classes
   4
   5   gSystem->Load("libANALYSIS");
   6   gSystem->Load("libANALYSISalice");
   7   gSystem->Load("libAOD");
   8   gSystem->Load("libCORRFW");
   9   gSystem->Load("libPWG3");
  10   gSystem->Load("libPWG3vertexingHF");
  11
  12   Int_t nptbins=4;
  13   const Int_t npars=3;
  14   Int_t nofcells[npars]={20,20,20};
  15   Float_t looses[npars]={0.,700.,0.8};
  16   Float_t tights[npars]={1500.,0.,1.};
  17   TString AxisTitle[npars];
  18   AxisTitle[0]="DistPS (micron)";
  19   AxisTitle[1]="TrackDispersion";
  20   AxisTitle[2]="CosPoint";
  21
  22   AliMultiDimVector* mvsig=new AliMultiDimVector("Signal","Signal",npars,nptbins,nofcells,looses,tights,AxisTitle);
  23   AliMultiDimVector* mvsig2=new AliMultiDimVector("Signal","Signal",npars,nptbins,nofcells,looses,tights,AxisTitle);
  24   AliMultiDimVector* mvbkg=new AliMultiDimVector("Background","Background",npars,nptbins,nofcells,looses,tights,AxisTitle);
  25
  26
  27   TF1* dsig=new TF1("dsig","exp(-x/310.)",0.,5000.);
  28   TF1* dbak=new TF1("dbak","exp(-x/50)",0.,5000.);
  29   gStyle->SetOptStat(2);
  30   for(Int_t isignev=0;isignev<5000;isignev++){
  31     Int_t ptbin=gRandom->Integer(nptbins);
  32     Float_t dists=dsig->GetRandom();
  33     Float_t distb=dbak->GetRandom();
  34     Float_t cpas=1-TMath::Abs(gRandom->Gaus(0.,0.02));
  35     Float_t cpab=0.8+gRandom->Rndm()*0.2;
  36     Float_t sigverts=gRandom->Gaus(200.,25.);
  37     Float_t sigvertb=gRandom->Gaus(250.,25.);
  38     Float_t vals[npars]={dists,sigverts,cpas};
  39     Float_t valb[npars]={distb,sigvertb,cpab};
  40     mvsig->FillAndIntegrate(vals,ptbin);
  41     mvsig2->Fill(vals,ptbin);  // alternative way of filling
  42     mvbkg->FillAndIntegrate(valb,ptbin);
  43   }
  44   mvsig2->Integrate(); // mvsig2 is now equal to mvsig
  45
  46   // Merge Pt bins
  47   AliMultiDimVector* mvsigallpt=mvsig->ShrinkPtBins(0,3);
  48   AliMultiDimVector* mvbkgallpt=mvbkg->ShrinkPtBins(0,3);
  49
  50   //caluclate significance
  51   AliSignificanceCalculator* cal=new AliSignificanceCalculator(mvsigallpt,mvbkgallpt,1.,5.);
  52   AliMultiDimVector* mvsts=cal->GetSignificance();
  53   AliMultiDimVector* mvess=cal->GetSignificanceError();
  54   AliMultiDimVector* mvsob=cal->CalculateSOverB();
  55   AliMultiDimVector* mvesob=cal->CalculateSOverBError();
  56   AliMultiDimVector* mvpur=cal->CalculatePurity();
  57   AliMultiDimVector* mvepur=cal->CalculatePurityError();
  58   Int_t fixed[3]={0,0,0};
  59   gStyle->SetPalette(1);
  60
  61   Int_t maxInd[3];
  62   Float_t sigMax=cal->GetMaxSignificance(maxInd,0);
  63   Float_t cut0=mvsigallpt->GetCutValue(0,maxInd[0]);
  64   Float_t cut1=mvsigallpt->GetCutValue(1,maxInd[1]);
  65   Float_t cut2=mvsigallpt->GetCutValue(2,maxInd[2]);
  66
  67   printf("=========== Pt Integrated ==============\n");
  68   printf("Maximum of significance found in bin %d %d %d\n",maxInd[0],maxInd[1],maxInd[2]);
  69   printf("                                 cuts= %f %f %f\n",cut0,cut1,cut2);
  70   printf("Significance = %f +- %f\n",mvsts->GetElement(maxInd,0),mvess->GetElement(maxInd,0));
  71   printf("S/B          = %f +- %f\n",mvsob->GetElement(maxInd,0),mvesob->GetElement(maxInd,0));
  72   printf("Purity       = %f +- %f\n",mvpur->GetElement(maxInd,0),mvepur->GetElement(maxInd,0));
  73
  74
  75   // 2D projections
  76   TH2F* hsig1 = mvsigallpt->Project(0,2,fixed,0);
  77   TH2F* hbkg1 = mvbkgallpt->Project(0,2,fixed,0);
  78   TH2F* hsts1 = mvsts->Project(0,2,fixed,0);
  79   TH2F* hess1 = mvess->Project(0,2,fixed,0);
  80   TH2F* hpur1 = mvpur->Project(0,2,fixed,0);
  81   TH2F* hepur1 = mvepur->Project(0,2,fixed,0);
  82   TH2F* hsob1 = mvsob->Project(0,2,fixed,0);
  83   TH2F* hesob1 = mvesob->Project(0,2,fixed,0);
  84
  85   TH2F* hsig2 = mvsigallpt->Project(1,2,fixed,0);
  86   TH2F* hbkg2 = mvbkgallpt->Project(1,2,fixed,0);
  87   TH2F* hsts2 = mvsts->Project(1,2,fixed,0);
  88   TH2F* hess2 = mvess->Project(1,2,fixed,0);
  89   TH2F* hpur2 = mvpur->Project(1,2,fixed,0);
  90   TH2F* hepur2 = mvepur->Project(1,2,fixed,0);
  91   TH2F* hsob2 = mvsob->Project(1,2,fixed,0);
  92   TH2F* hesob2 = mvesob->Project(1,2,fixed,0);
  93
  94   TCanvas* c1=new TCanvas("c1","Var 0 vs. 2",1000,800);
  95   c1->Divide(4,2);
  96   c1->cd(1);
  97   hsig1->Draw("colz");
  98   c1->cd(2);
  99   hbkg1->Draw("colz");
 100   c1->cd(3);
 101   hsts1->Draw("colz");
 102   c1->cd(4);
 103   hess1->Draw("colz");
 104   c1->cd(5);
 105   hsob1->Draw("colz");
 106   c1->cd(6);
 107   hesob1->Draw("colz");
 108   c1->cd(7);
 109   hpur1->Draw("colz");
 110   c1->cd(8);
 111   hepur1->Draw("colz");
 112
 113   TCanvas* c2=new TCanvas("c2","Var 1 vs. 2",1000,800);
 114   c2->Divide(4,2);
 115   c2->cd(1);
 116   hsig2->Draw("colz");
 117   c2->cd(2);
 118   hbkg2->Draw("colz");
 119   c2->cd(3);
 120   hsts2->Draw("colz");
 121   c2->cd(4);
 122   hess2->Draw("colz");
 123   c2->cd(5);
 124   hsob2->Draw("colz");
 125   c2->cd(6);
 126   hesob2->Draw("colz");
 127   c2->cd(7);
 128   hpur2->Draw("colz");
 129   c2->cd(8);
 130   hepur2->Draw("colz");
 131
 132 }