HLT/exa/PlotPythiaEvent.C

   1 void PlotPythiaEvent(char* path="/data1/AliRoot/pp/pythia/14000")
   2 {
   3   //gSystem->Load("$ROOTSYS/lib/libEG.so");                   // Root Event Generator interface
   4
   5   TH1F *hVr = new TH1F("hVr","Secondary vertex distribution (r direction)",100,0,100);
   6   TH1F *hVr2 = new TH1F("hVr2","Secondary vertex distribution near primary vertex",100,0,50);
   7
   8   TH1F *hVz = new TH1F("hVz","Secondary vertex distribution",100,-100,100);
   9   TH1F *hVz10 = new TH1F("hVz10","Secondary vertex distribution",10,-10,10);
  10   TH1F *hVz20 = new TH1F("hVz20","Secondary vertex distribution",200,-20,20);
  11   TH1F *hVz30 = new TH1F("hVz30","Secondary vertex distribution",300,-30,30);
  12   TH1F *hVz40 = new TH1F("hVz40","Secondary vertex distribution",400,-40,40);
  13   TH1F *hVz50 = new TH1F("hVz50","Secondary vertex distribution",500,-50,50);
  14   TH1F *hVz60 = new TH1F("hVz60","Secondary vertex distribution",600,-60,60);
  15   TH1F *hVz70 = new TH1F("hVz70","Secondary vertex distribution",700,-70,70);
  16   TH1F *hVz80 = new TH1F("hVz80","Secondary vertex distribution",800,-80,80);
  17   TH1F *hVz90 = new TH1F("hVz90","Secondary vertex distribution",900,-90,90);
  18   TH1F *hVz100 = new TH1F("hVz100","Secondary vertex distribution",100,-100,100);
  19
  20   TH1F *hPt = new TH1F("hPt","Transversal momentum",50,0,10);
  21   TH1F *hPt10 = new TH1F("hPt10","Transversal momentum for particles within the vertex cut",50,0,10);
  22   TH1F *hPt20 = new TH1F("hPt20","Transversal momentum for particles within the vertex cut",50,0,10);
  23   TH1F *hPt30 = new TH1F("hPt30","Transversal momentum for particles within the vertex cut",50,0,10);
  24   TH1F *hPt40 = new TH1F("hPt40","Transversal momentum for particles within the vertex cut",50,0,10);
  25   TH1F *hPt50 = new TH1F("hPt50","Transversal momentum for particles within the vertex cut",50,0,10);
  26   TH1F *hPt60 = new TH1F("hPt60","Transversal momentum for particles within the vertex cut",50,0,10);
  27   TH1F *hPt70 = new TH1F("hPt70","Transversal momentum for particles within the vertex cut",50,0,10);
  28   TH1F *hPt80 = new TH1F("hPt80","Transversal momentum for particles within the vertex cut",50,0,10);
  29   TH1F *hPt90 = new TH1F("hPt90","Transversal momentum for particles within the vertex cut",50,0,10);
  30   TH1F *hPt100 = new TH1F("hPt100","Transversal momentum for particles within the vertex cut",50,0,10);
  31
  32 /*TH1F *hPt3 = new TH1F("hPt3","Transversal momentum for particles outside the vertex cut",100,0,10);
  33   hPt3->SetXTitle("GeV/c");
  34   hPt3->SetYTitle("N");
  35
  36   TH1F *hPt4 = new TH1F("hPt4","hPt2/hPt",100,0,10);
  37   hPt4->SetXTitle("GeV/c");
  38   hPt4->SetYTitle("N");
  39
  40   TH1F *hPt5 = new TH1F("hPt5","hPt3/hPt",100,0,10);
  41   hPt5->SetXTitle("GeV/c");
  42   hPt5->SetYTitle("N");*/
  43
  44   Char_t fname[256];
  45   Int_t nFiles = 4;
  46   for (Int_t iFile = 0; iFile < nFiles; iFile++) {
  47
  48     sprintf(fname,"%s/pythia6_%d.root",path,iFile);
  49     TFile *f = new TFile(fname);
  50     cout << "Opened file " << fname << endl;
  51
  52     TClonesArray* arr = new TClonesArray("TParticle");
  53
  54     TTree *tree = (TTree *) f->Get("pythia");
  55     tree->GetBranch("particles")->SetAddress(&arr);
  56     Int_t nEntries = tree->GetEntries();
  57
  58     TParticle *p;
  59     TParticle d1;
  60     //TParticle d2;
  61
  62     for (Int_t i = 0; i < nEntries; i++) {
  63       tree->GetEvent(i);
  64       TClonesArray &part = *arr;
  65       Int_t nPart = part.GetEntriesFast();
  66       for (Int_t iPart = 0; iPart < nPart; iPart++) {
  67         p = (TParticle*)part[iPart];
  68         if (p->GetPdgCode()==3122){
  69
  70           Int_t mother = p->GetFirstMother();
  71
  72           Int_t fdaughter = p->GetFirstDaughter();
  73           Int_t ldaughter = p->GetLastDaughter();
  74           d1 = (TParticle*)part[fdaughter-1];
  75           //if (!d1) continue;
  76           d2 = (TParticle*)part[ldaughter-1];
  77           //cout << fdaughter << " " << ldaughter << endl;
  78           //cout << d1->GetPdgCode() << " " << d2->GetPdgCode() << endl;
  79
  80           Double_t vx = d2->Vx();
  81           Double_t vy = d2->Vy();
  82           Double_t vz = d2->Vz();
  83           Double_t vt = TMath::Sqrt((vx*vx)+(vy*vy));
  84           Double_t vr = TMath::Sqrt((vx*vx)+(vy*vy)+(vz*vz));
  85           Double_t pt = p->Pt();
  86
  87           hVr->Fill(vr);
  88           if (vr < 10) hVr2->Fill(vr);
  89           hVz->Fill(vz);
  90           hPt->Fill(pt);
  91           if (TMath::Abs(vz) < 10)
  92             {
  93               hVz10->Fill(vz);
  94               hPt10->Fill(pt);
  95             }
  96           if (TMath::Abs(vz) < 20)
  97             {
  98               hVz20->Fill(vz);
  99               hPt20->Fill(pt);
 100             }
 101           if (TMath::Abs(vz) < 30)
 102             {
 103               hVz30->Fill(vz);
 104               hPt30->Fill(pt);
 105             }
 106           if (TMath::Abs(vz) < 40)
 107             {
 108               hVz40->Fill(vz);
 109               hPt40->Fill(pt);
 110             }
 111           if (TMath::Abs(vz) < 50)
 112             {
 113               hVz50->Fill(vz);
 114               hPt50->Fill(pt);
 115             }
 116           if (TMath::Abs(vz) < 60)
 117             {
 118               hVz60->Fill(vz);
 119               hPt60->Fill(pt);
 120             }
 121           if (TMath::Abs(vz) < 70)
 122             {
 123               hVz70->Fill(vz);
 124               hPt70->Fill(pt);
 125             }
 126           if (TMath::Abs(vz) < 80)
 127             {
 128               hVz80->Fill(vz);
 129               hPt80->Fill(pt);
 130             }
 131           if (TMath::Abs(vz) < 90)
 132             {
 133               hVz90->Fill(vz);
 134               hPt90->Fill(pt);
 135             }
 136           if (TMath::Abs(vz) < 100)
 137             {
 138               hVz100->Fill(vz);
 139               hPt100->Fill(pt);
 140             }
 141
 142         }
 143
 144       }
 145
 146     }
 147
 148   }
 149
 150   //TCanvas *c = new TCanvas("c","c",50,50,600,800);
 151   //c->Divide(1,2);
 152   //c->cd(1);
 153   //gPad->SetLogy();
 154   //hVz->Draw();
 155   //hPt->Draw();
 156   //c->cd(2);
 157   //gPad->SetLogy();
 158   //hVz2->Draw();
 159   //hPt3->Draw();
 160
 161   //hPt4->Divide(hPt2,hPt);
 162   //hPt5->Divide(hPt3,hPt);
 163   //TCanvas *c2 = new TCanvas("c2","c2",50,50,600,800);
 164   //c2->Divide(1,2);
 165   //c2->cd(1);
 166   //hPt4->Draw();
 167   //c2->cd(2);
 168   //hPt5->Draw();
 169
 170   TFile *outfile = new TFile("/tmp/results.root","recreate");
 171   hVr->Write();
 172   hVr2->Write();
 173
 174   hVz->Write();
 175   hVz10->Write();
 176   hVz20->Write();
 177   hVz30->Write();
 178   hVz40->Write();
 179   hVz50->Write();
 180   hVz60->Write();
 181   hVz70->Write();
 182   hVz80->Write();
 183   hVz90->Write();
 184   hVz100->Write();
 185
 186   hPt->Write();
 187   hPt10->Write();
 188   hPt20->Write();
 189   hPt30->Write();
 190   hPt40->Write();
 191   hPt50->Write();
 192   hPt60->Write();
 193   hPt70->Write();
 194   hPt80->Write();
 195   hPt90->Write();
 196   hPt100->Write();
 197
 198   outfile->Close();
 199   f->Close();
 200
 201 }
 202