EMCAL/anaDst.C

   1 //*-- Author: Aleksei Pavlinov(WSU)
   2 //            8-feb-2002 new version with AliEMCALJetMicroDst
   3 #include "anaDst.h"
   4 #include "jetDst.h"
   5 #if !defined(__CINT__) || defined(__MAKECINT__)
   6 // These files are needed for compilation only
   7 // (see directory /auto/u/pavlinov/macros on PDSF)
   8 // Using with caution if you are not PAI
   9 #include "macroIncludePai.h"
  10 #include "macroIncludeAlice.h"
  11 #endif
  12
  13 class TH1F;
  14 TH1F* hJetEt, *hJetEta, *hJetPhi, *hNtracksInJet, *hPtTracksInJet, *hTrackDR, *hNjet;
  15
  16 class TH2F; TH2F* hPartPt, *hPartEta, *hPartPhi;
  17 TH1F* hPartDiffPt, *hPartDiffEta, * hPartDiffPhi;
  18
  19 TH1F *hJetPartAng, *hDeta, *hDphi, *hDpt; // Comparing between jet and nearest parton
  20
  21 class TCanvas; TCanvas *c1, *c2, *c3;
  22 class TPad;    TPad  *p1;
  23 class TList;   TList   *lPartons, *lJets, *lCompPJ;
  24
  25 class TVector3;
  26
  27 Int_t nparticles, nJet;
  28 Double_t phiMax = TMath::Pi()*120./180.;
  29 Double_t rJet = 0.5; // this default
  30
  31 TFile *file;
  32 class AliEMCALJetMicroDst; AliEMCALJetMicroDst *jDst=0;
  33 Int_t nevMax = 1000;
  34 TString fname;
  35 // fitting
  36 class TF1 *g1;// gauss
  37 extern "C++" {void loadlibs();}
  38
  39 void
  40 anaDst(Int_t mode)
  41 {
  42 // Dynamically link some shared libs
  43   if (gClassTable->GetID("AliRun") < 0) {
  44       gROOT->LoadMacro("../macros/loadlibs.C");
  45       loadlibs();
  46   }
  47
  48 //partons after HSC and jet in acceptance.
  49   TVector3 *vJet;
  50   TVector3 *vPart1, *vPart2;
  51   vJet   = new TVector3;
  52   vPart1 = new TVector3;
  53   vPart2 = new TVector3;
  54
  55 // for getting information about nFileOut
  56 //#if !defined(__CINT__)
  57 //  gROOT->ProcessLine(".L jetDst.C+");
  58 //#else
  59 //  gROOT->LoadMacro("jetDst.C");
  60 //#endif
  61   if(!defineMicroDst(mode)) return;
  62
  63   bookHist1();
  64   bookHistPartonsAfterHardSc();
  65
  66     Float_t phiT[50], etaT[50], ptT[50];
  67
  68     Int_t nbytes=0, nb=0, nentries = Int_t(jDst->GetTree()->GetEntries());
  69     if(nevMax>nentries) nevMax = nentries;
  70     for (int nev=0; nev< nevMax; nev++) {
  71         nb = jDst->GetEntry(nev);
  72         nbytes += nb;
  73         printf("\n Event .............%d bytes %d\n", nev, nb);
  74         //        jDst->Print();
  75         if(nb<0) break; // last event or something wrong
  76         //        return;
  77
  78         Bool_t inEmc = fillInfoForPartons(vPart1, vPart2);
  79         //        continue;
  80
  81         nJet = jDst->GetNjet();
  82         Double_t et, eta, phi;
  83
  84         for(Int_t ij=0; ij < nJet; ij++) {
  85
  86            jDst->GetJet(ij, 1, (*vJet)); // W
  87
  88            et  = vJet->Pt();
  89            eta = vJet->Eta();
  90            phi = vJet->Phi();
  91
  92            printf("\n Jet:%d E %f phi %f eta %f tracks %d\n", ij,
  93                        et, phi, eta, 0);
  94
  95                 //              hNtracksInJet->Fill((Float_t)lJet->NTracks());
  96                 // Safe zone - rJet - 13-feb-2002
  97            if(nJet==1 && TMath::Abs(eta)<(0.7-rJet)&&
  98                phi>rJet && phi<(phiMax-rJet)) {
  99               hJetEt->Fill(et);
 100               hJetEta->Fill(eta);
 101               hJetPhi->Fill(phi);
 102               if(inEmc) compareJetPartons(vJet, vPart1, vPart2);
 103            }
 104
 105         } // jet
 106         cleanUpEvent();
 107    } // event
 108
 109     gStyle->SetOptStat(1111111);
 110
 111     //drawPartons();
 112
 113 #if !defined(__CINT__)
 114     if(hJetPartAng&&hJetPartAng->GetEntries()>2.){
 115       _drawListOfHist(lCompPJ, 2, 2);
 116       fitGauss(hDpt, 0);
 117     }
 118 #endif
 119
 120     gROOT->cd();
 121     //    printf("\n Variant %1i : Alice Root file => %s \n",
 122     //mode, nf.Data());
 123 }
 124
 125 Bool_t defineMicroDst(Int_t mode)
 126 {
 127   TString dir("/auto/u/pavlinov/ALICE/RES/FILE/");
 128   //  TString  *sTmp = defineInput(mode); // see jetDst.h
 129   switch(mode){
 130     //  case -4:
 131     //if(sTmp) fname = sTmp->Data();
 132     //else     return kFALSE;
 133     //break;
 134   case -3:
 135     fname = dir + "PythiaMode-3R0.5Seed4.0JEt20.0CEt0.0PtCut0.0.root";
 136     break;
 137   case -2:
 138     fname = dir + "PythiaMode-2R0.5Seed4.0JEt20.0CEt1.0PtCut1.0.root";
 139     break;
 140   case -1:
 141     fname = dir + "PythiaMode-1R0.5Seed4.0JEt20.0CEt1.0PtCut2.0.root";
 142     break;
 143
 144   case 11:
 145     fname = dir + "HijingMode11R0.3Seed4.0JEt40.0CEt1.0PtCut2.0.root";
 146     break;
 147   default:
 148      printf("\n<I> Wrong input mode => %i\n", mode);
 149      return kFALSE;
 150   }
 151   jDst = new AliEMCALJetMicroDst;
 152   jDst->Open(fname.Data());
 153   gROOT->GetListOfBrowsables()->Add(jDst);
 154   return kTRUE;
 155 }
 156
 157 Bool_t fillInfoForPartons(TVector3 *vec1, TVector3 *vec2)
 158 {
 159   // npart=4 for HIJING and 2 for PYTHIA
 160   Int_t npart = jDst->GetNpart();
 161   if(!jDst->GetParton(npart-2, (*vec1))||!jDst->GetParton(npart-1, (*vec2)) )
 162   return kFALSE;
 163
 164   static Double_t dphi;
 165
 166   Bool_t inEmc = kFALSE;
 167   TVector3 *vec = vec1;
 168   for (Int_t i=0; i<2; i++) {
 169     Double_t ri = Double_t(i);
 170     hPartPt->Fill(vec->Pt(), ri);
 171     hPartEta->Fill(vec->Eta(), ri);
 172     hPartPhi->Fill(vec->Phi(), ri);
 173     if(fabs(vec->Eta())<0.7 && vec->Phi()>0 && vec->Phi()< phiMax) inEmc=kTRUE;
 174     vec = vec2;
 175   }
 176   hPartDiffPt->Fill(vec1->Pt()-vec2->Pt());
 177   hPartDiffEta->Fill(vec1->Eta()-vec2->Eta());
 178   dphi = vec1->Phi() - vec2->Phi();
 179   dphi = fabs(TVector2::Phi_mpi_pi(dphi));
 180   hPartDiffPhi->Fill(dphi);
 181   return inEmc;
 182 }
 183
 184 void compareJetPartons(TVector3 *vecJet, TVector3 *vec1, TVector3 *vec2)
 185 {// 30-jan-2002
 186   TVector3 *v, *vw;
 187   Double_t ang, angw;
 188
 189   v    = vec1;
 190   ang  = vecJet->Angle((*v));
 191   vw   = vec2;
 192   angw = vecJet->Angle((*vw));
 193   if(ang>angw) {
 194     ang = angw;
 195     v   = vw;
 196   }
 197   hJetPartAng->Fill(ang);
 198   hDeta->Fill(vecJet->Eta() - v->Eta());
 199   hDphi->Fill(vecJet->Phi() - v->Phi());
 200   hDpt->Fill(vecJet->Pt() - v->Pt());
 201 }
 202
 203 void bookHist1()
 204 {
 205     gROOT->cd();
 206 // Book histos
 207     hNjet= new TH1F("hNjet","# jets",    11,  -0.5, 10.5);
 208     hJetEt   = new TH1F("hJetEt","Energy of jet",    150,  0.0, 450.);
 209     hJetEta  = new TH1F("hJetEta","#eta of jet",     180, -0.9,  0.9);
 210     hJetPhi = new TH1F("hJetPhi","#phi of jet",      62, -3.1,  3.1);
 211     hNtracksInJet   = new TH1F("hNtracksInJet","n tracks",   30,  0.5, 29.5);
 212     hPtTracksInJet  = new TH1F("hPtTracksInJet","p_{T} of tracks in jets cone",
 213     100., 0., 100.);
 214     hTrackDR  = new TH1F("hTrackDR","Track dR", 120., 0.,   6.);
 215 #if !defined(__CINT__)
 216     lJets = _moveHistsToList("Hist for jets");
 217 #endif
 218 // QA
 219     hJetPartAng =  new TH1F("hJetPartAng","angle between jet and nearest parton", 100,0.0,0.2);
 220     hDeta =  new TH1F("hDeta","#eta (jets) - #eta (partons)", 100,  -.1, +.1);
 221     hDphi =  new TH1F("hDphi","#phi (jets} - #phi (partons)", 100,  -.1, +.1);
 222     hDpt  =  new TH1F("hDpt","E_{t} (jets) - E_{t} (partons)",100, -200., +200.);
 223 #if !defined(__CINT__)
 224     lCompPJ = _moveHistsToList("Hist for jet vs parton");
 225 #endif
 226 }
 227
 228 void bookHistPartonsAfterHardSc()
 229 {
 230   hPartPt  = new TH2F("hPartPt","Partons p_{t}", 500,  30., 1030., 2, -0.5, 1.5);
 231   hPartEta = new TH2F("hPartEta","Partons #eta", 100,  -5., 5., 2,-0.5,1.5);
 232   hPartPhi = new TH2F("hPartPhi","Partons #phi", 120,-0.2,2.2, 2,-0.5,1.5);
 233   hPartDiffPt  = new TH1F("hPartDiffPt","p_{t} difference for partons",100, -100., 100.);
 234   hPartDiffEta  = new TH1F("hPartDiffEta","#eta difference for partons",100, -5., 5.);
 235   hPartDiffPhi  = new TH1F("hPartDiffPhi","#phi angle for partons",320, 0.0, 3.2);
 236 #if !defined(__CINT__)
 237   lPartons = _moveHistsToList("Hist for partons");
 238 #endif
 239 }
 240
 241 void drawPartons()
 242 {
 243   gStyle->SetOptStat(1111111);
 244   c3 = new TCanvas("c3","Partons canvas",400,10,600,700);
 245   c3->Divide(2,3);
 246   TH1D *hid;
 247
 248   c3->cd(1);
 249   hid = hPartPt->ProjectionX("_pt1",1,1); _drawHist(hid, 2);
 250   hid = hPartPt->ProjectionX("_pt2",2,2); _drawHist(hid, 1, 2, "same");
 251   c3->cd(2); _drawHist(hPartDiffPt,2);
 252
 253   c3->cd(3);
 254   hid = hPartEta->ProjectionX("_eta1",1,1); _drawHist(hid, 2);
 255   hid = hPartEta->ProjectionX("_eta2",2,2); _drawHist(hid, 1, 2, "same");
 256   c3->cd(4); _drawHist(hPartDiffEta,2);
 257
 258   c3->cd(5);
 259   hid = hPartPhi->ProjectionX("_phi1",1,1); _drawHist(hid, 2);
 260   hid = hPartPhi->ProjectionX("_phi2",2,2); _drawHist(hid, 1, 2, "same");
 261   c3->cd(6); _drawHist(hPartDiffPhi,2);
 262 }
 263
 264 void fitGauss(TH1F* hid, Int_t opt)
 265 { // see /home/pavlinov/cosmic/pythia/geant/lin.C
 266   Double_t c = 2.5;
 267   Double_t xmax = hid->GetBinCenter(hid->GetMaximumBin());
 268   Double_t rms  = hid->GetRMS();
 269   Double_t xmi=xmax -c*rms, xma=xmax+c*rms;
 270   _fit(111);
 271   if(!g1) {
 272     g1 = new TF1("g1","gaus", xmi, xma);
 273   }
 274   g1->SetLineColor(2); // red
 275   g1->SetLineWidth(3);
 276
 277   TString stmp = fname(0,fname.Length()-5); // discard .root
 278   if(!c1) {
 279     c1 = new TCanvas("c1","Jet resolution",400,10,600,700);
 280     p1 = _pad(stmp.Data(),c1);
 281   }
 282   p1->cd(); p1->Clear();
 283   hid->Fit("g1","R+");
 284
 285   if(opt) {
 286     _drawHist(hid, 2);
 287     if(opt>=10) {
 288       stmp += "_res.ps";
 289       stmp.ReplaceAll("/FILE","");
 290       p1->Print(stmp.Data());
 291     }
 292   }
 293 }
 294
 295 void cleanUpEvent()
 296 {
 297   /* ???
 298   if(lPart[0]) {
 299     lPart[0]->Delete();
 300     lPart[0]=0;
 301     lPart[1]->Delete();
 302     lPart[1]=0;
 303   }
 304   */
 305 }