A macro for checking a few basic pileup properties in MC
authorbelikov <Iouri.Belikov@cern.ch>
Wed, 1 Oct 2014 10:24:21 +0000 (12:24 +0200)
committerbelikov <Iouri.Belikov@cern.ch>
Wed, 1 Oct 2014 10:24:21 +0000 (12:24 +0200)
ITS/UPGRADE/macros/QA/CheckPileupVerticesMC.C [new file with mode: 0644]

diff --git a/ITS/UPGRADE/macros/QA/CheckPileupVerticesMC.C b/ITS/UPGRADE/macros/QA/CheckPileupVerticesMC.C
new file mode 100644 (file)
index 0000000..bf31a70
--- /dev/null
@@ -0,0 +1,116 @@
+//
+// This macro checks a few basic properties of MC pileup vertices
+//
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+   #include <Riostream.h>
+   #include <TFile.h>
+   #include <TTree.h>
+   #include <TH1F.h>
+   #include <TCanvas.h>
+   #include <TMath.h>
+   #include <TClonesArray.h>
+   #include <TParticle.h>
+   #include <TParticlePDG.h>
+
+   #include "AliRunLoader.h"
+   #include "AliHeader.h"
+   #include "AliStack.h"
+   #include "AliGenCocktailEventHeader.h"
+#endif
+
+
+
+Int_t FindContributors(Float_t t, AliStack *stack) {
+  Int_t ntrk=0;
+  Int_t np=stack->GetNtrack();
+  for (Int_t i=0; i<np; i++) {
+      if (!stack->IsPhysicalPrimary(i)) continue;
+      TParticle *part=stack->Particle(i);
+      if (!part) continue;
+      TParticlePDG *partPDG = part->GetPDG();
+      if (!partPDG) continue;
+      if (TMath::Abs(partPDG->Charge())<1e-10) continue;
+      if (part->Pt()<0.5) continue;
+      Float_t dt=0.5*(t - part->T())/(t + part->T());
+      if (TMath::Abs(dt)>1e-5) continue;
+      ntrk++;
+  }
+  return ntrk;
+}
+
+
+void CheckPileupVerticesMC() {
+    TH1F *hnv=new TH1F("hnv","Number of vertices",100,0,300);
+    TH1F *hnt=new TH1F("hnt","Number of tracks per vertex",100,0,200);
+    TH1F *ht=new TH1F("ht","Interaction time;t (sec);", 100,-0.05e-3,0.05e-3);
+    TH1F *hdt=new TH1F("hdt","Difference in time between interactions; #Deltat (sec)",100,0,0.01e-3);
+    TH1F *hz=new TH1F("hz","Z position of vertices; Z (cm)",100,-20,20);
+    TH1F *hdz=new TH1F("hdz","Difference in Z position between vertices; #DeltaZ (cm)",100,0,2.5);
+
+    AliRunLoader *rl = AliRunLoader::Open("galice.root","something");
+    rl->LoadHeader();
+    rl->LoadKinematics();
+
+    Int_t nEvents=rl->GetNumberOfEvents();
+    for (Int_t i=0; i<nEvents; i++) {
+        cout<<"Event "<<i<<" out of "<<nEvents<<endl;
+        rl->GetEvent(i);
+        AliStack *stack=rl->Stack();    
+
+        AliHeader *ah=rl->GetHeader();
+        AliGenCocktailEventHeader *cock=
+            (AliGenCocktailEventHeader*)ah->GenEventHeader();
+        TList *headers=cock->GetHeaders();
+        const Int_t nvtx=headers->GetEntries();
+
+        hnv->Fill(nvtx);
+        Double_t *z=new Double_t[nvtx];
+        Double_t *t=new Double_t[nvtx];
+        Int_t *idx=new Int_t[nvtx];
+        for (Int_t v=0; v<nvtx; v++) {
+            AliGenEventHeader *h=(AliGenEventHeader *)headers->At(v);
+            t[v]=h->InteractionTime();
+            ht->Fill(t[v]);
+            Int_t nt=FindContributors(t[v],stack);
+            hnt->Fill(nt);
+
+            TArrayF vtx(3); h->PrimaryVertex(vtx);
+            z[v]=vtx[2];
+            hz->Fill(z[v]);
+        }
+
+       TMath::Sort(nvtx,t,idx);
+        for (Int_t v=0; v<nvtx-1; v++) {
+         Int_t i1=idx[v];
+         Int_t i2=idx[v+1];
+          Float_t dt=t[i1]-t[i2];
+          hdt->Fill(dt);
+       }        
+       TMath::Sort(nvtx,z,idx);
+        for (Int_t v=0; v<nvtx-1; v++) {
+         Int_t i1=idx[v];
+         Int_t i2=idx[v+1];
+          Float_t dz=z[i1]-z[i2];
+          hdz->Fill(dz);
+       }        
+
+        delete[] t;
+        delete[] z;
+        delete[] idx;
+
+    }
+
+    TCanvas *c1=new TCanvas("c1","",0,0,750,1000); c1->Divide(1,2); 
+    c1->cd(1); hnv->Draw();
+    c1->cd(2); gPad->SetLogy(); hnt->Draw();
+
+    TCanvas *c2=new TCanvas("c2","",0,0,750,1000); c2->Divide(1,2); 
+    c2->cd(1); ht->Draw();
+    c2->cd(2); hdt->Fit("expo");
+
+    TCanvas *c3=new TCanvas("c3","",0,0,750,1000); c3->Divide(1,2); 
+    c3->cd(1); hz->Fit("gaus");
+    c3->cd(2); hdz->Draw(); //hdz->Fit("expo");
+
+}