JETAN/AliJetAnalysis.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //---------------------------------------------------------------------
  17 // JetAnalysis class
  18 // Analyse Jets
  19 // Author: andreas.morsch@cern.ch
  20 //---------------------------------------------------------------------
  21
  22 #include "AliJetAnalysis.h"
  23 ClassImp(AliJetAnalysis)
  24
  25
  26 ////////////////////////////////////////////////////////////////////////
  27 #include <TH1F.h>
  28 #include <TProfile.h>
  29 #include <TCanvas.h>
  30 #include <TFile.h>
  31 #include <TTree.h>
  32 #include <TStyle.h>
  33 #include <TSystem.h>
  34 #include <TLorentzVector.h>
  35
  36 #include "AliJetProductionDataPDC2004.h"
  37 #include "AliJet.h"
  38 #include "AliJetESDReaderHeader.h"
  39 #include "AliUA1JetHeader.h"
  40 #include "AliLeading.h"
  41
  42     AliJetAnalysis::AliJetAnalysis()
  43 {
  44     // Constructor
  45     fDirectory    = 0x0;
  46     fEventMin     =  0;
  47     fEventMax     = -1;
  48     fRunMin       =  0;
  49     fRunMax       = 11;
  50 }
  51
  52 void AliJetAnalysis::Analyse()
  53 {
  54     //
  55     // Some histos
  56     //
  57     TH1F::AddDirectory(0);
  58     TProfile::AddDirectory(0);
  59
  60     TH1F* e0H    = new TH1F("e0H"  ,"Jet Energy (reconstructed)",      40,  0., 200.);
  61     TH1F* e1H    = new TH1F("e1H" , "Jet Energy (generated)",          40,  0., 200.);
  62     TH1F* e2H    = new TH1F("e2H" , "Jet Energy (generated, nrec = 0", 40,  0., 200.);
  63     TH1F* e3H    = new TH1F("e3H" , "Jet Energy (leading)",            40,  0., 200.);
  64     TH1F* e4H    = new TH1F("e4H" , "Jet Energy (reconstructed: 105 < Egen < 125", 40,  0., 200.);
  65
  66     TH1F* e5H    = new TH1F("e5H" , "Jet Energy (generated)", 40,  0., 200.);
  67     TH1F* e6H    = new TH1F("e6H" , "Jet Energy (generated)", 40,  0., 200.);
  68     TH1F* e7H    = new TH1F("e7H" , "Jet Energy (generated)", 40,  0., 200.);
  69     TH1F* e8H    = new TH1F("e8H" , "Jet Energy (generated)", 40,  0., 200.);
  70
  71     TProfile* r5H    = new TProfile("r5H" , "rec/generated", 20,  0., 200, 0., 1., "S");
  72     TProfile* r6H    = new TProfile("r6H" , "rec/generated", 20,  0., 200, 0., 1., "S");
  73
  74     TProfile* r7H    = new TProfile("r7H" , "rec/generated", 20,  0., 200, 0., 1., "S");
  75     TProfile* r8H    = new TProfile("r8H" , "rec/generated", 20,  0., 200, 0., 1., "S");
  76
  77
  78     TH1F* dr1H = new TH1F("dr1H", "delta R",  160, 0.,   2.);
  79     TH1F* dr2H = new TH1F("dr2H", "delta R",  160, 0.,   2.);
  80     TH1F* dr3H = new TH1F("dr4H", "delta R",  160, 0.,   2.);
  81
  82     TH1F* etaH  = new TH1F("etaH",  "eta",  160, -2.,   2.);
  83     TH1F* eta1H = new TH1F("eta1H", "eta",  160, -2.,   2.);
  84     TH1F* eta2H = new TH1F("eta2H", "eta",  160, -2.,   2.);
  85
  86     TH1F* phiH   = new TH1F("phiH",   "phi",  160, -3.,   3.);
  87     TH1F* dphiH  = new TH1F("dphiH",  "phi",  160,  0.,   3.1415);
  88     TH1F* phi1H  = new TH1F("phi1H", "phi",   160,  0.,   6.28);
  89     TH1F* phi2H  = new TH1F("phi2H", "phi",   160,  0.,   6.28);
  90
  91
  92     TProfile* drP1    = new TProfile("drP1" , "Delta_R", 20,  0., 200, -1., 1., "S");
  93     TProfile* drP2    = new TProfile("drP2" , "Delta_R", 20,  0., 200, -1., 1., "S");
  94
  95   // Run data
  96     AliJetProductionDataPDC2004* runData = new AliJetProductionDataPDC2004();
  97
  98
  99     // Loop over runs
 100     TFile* jFile = 0x0;
 101
 102     for (Int_t iRun = fRunMin; iRun <= fRunMax; iRun++)
 103     {
 104
 105         // Open file
 106         char fn[20];
 107         sprintf(fn, "%s/%s.root", fDirectory, (runData->GetRunTitle(iRun)).Data());
 108
 109
 110         jFile = new TFile(fn);
 111
 112         printf("  Analyzing run: %d %s\n", iRun,fn);
 113         // Get jet header and display parameters
 114         AliUA1JetHeader* jHeader =
 115             (AliUA1JetHeader*) (jFile->Get("AliUA1JetHeader"));
 116         // jHeader->PrintParameters();
 117
 118         // Get reader header and events to be looped over
 119         AliJetESDReaderHeader *jReaderH =
 120             (AliJetESDReaderHeader*)(jFile->Get("AliJetKineReaderHeader"));
 121
 122         if (fEventMin == -1) fEventMin =  jReaderH->GetFirstEvent();
 123         if (fEventMax == -1) {
 124             fEventMax =  jReaderH->GetLastEvent();
 125         } else {
 126             fEventMax = TMath::Min(fEventMax, jReaderH->GetLastEvent());
 127         }
 128
 129
 130         // Calculate weight
 131         Float_t wgt = runData->GetWeight(iRun) / Float_t(fEventMax - fEventMin + 1);
 132         Float_t ptmin, ptmax;
 133         runData->GetPtHardLimits(iRun, ptmin, ptmax);
 134
 135
 136         // Loop over events
 137         AliJet *jets  = 0x0;
 138         AliJet *gjets = 0x0;
 139         AliLeading *leading = 0x0;
 140         Float_t egen, etag, econ, erec;
 141
 142
 143         for (Int_t i = fEventMin; i < fEventMax; i++) {
 144             printf("  Analyzing run: %d  Event %d / %d \n",
 145                    iRun, i, fEventMax);
 146
 147             // Het next tree with AliJet
 148             char nameT[100];
 149             sprintf(nameT, "TreeJ%d",i);
 150             TTree *jetT =(TTree *)(jFile->Get(nameT));
 151             jetT->SetBranchAddress("FoundJet",    &jets);
 152             jetT->SetBranchAddress("GenJet",      &gjets);
 153             jetT->SetBranchAddress("LeadingPart", &leading);
 154             jetT->GetEntry(0);
 155
 156 //
 157 //    Find the jet with the highest E_T within fiducial region
 158 //
 159             Int_t njets = jets->GetNJets();
 160             Int_t imax = -1;
 161             Int_t jmax = -1;
 162             Float_t emax = 0.;
 163
 164             for (Int_t ij = 0; ij < njets; ij++) {
 165                 if (jets->GetPt(ij) > emax &&
 166                     TMath::Abs(jets->GetEta(ij)) < 0.60) {
 167                     emax = jets->GetPt(ij);
 168                     jmax = imax;
 169                     imax = ij;
 170                 }
 171             }
 172
 173
 174             if (imax == -1) {
 175                 Int_t ngen = gjets->GetNJets();
 176                 if(ngen>0) e2H->Fill(gjets->GetPt(0), wgt);
 177             } else {
 178 //            printf("Reconstructed Jet %5d %13.3f %13.3f %13.3f\n",
 179 //                    imax, emax, jets->GetEta(imax), jets->GetPhi(imax));
 180                 econ = jets->GetPt(imax);
 181                 erec = jets->GetPt(imax) / 0.65;
 182                 dr1H->Fill(jets->GetEta(imax) - gjets->GetEta(0));
 183 //
 184 //    Find the generated jet closest to the reconstructed
 185 //
 186                 Float_t rmin;
 187                 Float_t etamin = 1.e6;
 188
 189                 Int_t   igen;
 190                 Float_t etaj = jets->GetEta(imax);
 191                 Float_t phij = jets->GetPhi(imax);
 192
 193                 Int_t ngen = gjets->GetNJets();
 194
 195                 if (ngen != 0) {
 196                     rmin = 1.e6;
 197                     igen = -1;
 198                     for (Int_t j = 0; j < ngen; j++) {
 199                         etag = gjets->GetEta(j);
 200                         Float_t phig = gjets->GetPhi(j);
 201                         Float_t deta = etag - etaj;
 202                         Float_t dphi = TMath::Abs(phig - phij);
 203                         if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
 204                         Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
 205                         if (r  < rmin) {
 206                             rmin = r;
 207                             etamin = deta;
 208                             igen = j;
 209                         }
 210                     }
 211
 212                     egen = gjets->GetPt(igen);
 213                     e1H->Fill(gjets->GetPt(igen), wgt);
 214                     etag = gjets->GetEta(igen);
 215                     Float_t phig = gjets->GetPhi(igen);
 216                     Float_t dphi = phig - phij;
 217
 218 //                if (econ < ptmax) {
 219                     e0H->Fill(erec, wgt);
 220 //                } else {
 221 //                    e0H->Fill(erec, 6.7e-6);
 222 //                }
 223
 224
 225
 226                     if (egen > 20. && egen < 40.) {
 227                         phiH->Fill((dphi));
 228                         etaH->Fill(etag - etaj, wgt);
 229                         phi1H->Fill(phig);
 230                         phi2H->Fill(phij);
 231                         e4H->Fill(jets->GetPt(imax));
 232                     }
 233
 234                     if (erec > 90. && erec < 110. && rmin < 0.1) {
 235                         e5H->Fill(egen, wgt);
 236                         dr2H->Fill(rmin);
 237                         if (egen < 30.) {
 238                             printf("Strange jet %6d %13.3f %13.3f %13.3f \n",
 239                                    imax, etaj, phij, erec);
 240                             for (Int_t j = 0; j < ngen; j++) {
 241                                 printf("Generated %6d %13.3f %13.3f %13.3f\n",
 242                                        j, gjets->GetEta(j),
 243                                        gjets->GetPhi(j), gjets->GetPt(j));
 244
 245                             }
 246                         }
 247                     }
 248
 249                     if (rmin < 0.1) {
 250                         r5H->Fill(egen, jets->GetPt(imax)/egen, wgt);
 251                         r6H->Fill(jets->GetPt(imax)
 252                                   / 0.4, jets->GetPt(imax)/egen, wgt);
 253                         e8H->Fill(erec, wgt);
 254                     }
 255
 256                     if (rmin < 0.1) {
 257                         drP1->Fill(egen, etamin, wgt);
 258                     }
 259                 } // ngen !=0
 260             } // has reconstructed jet
 261
 262 //
 263 //   Leading particle
 264 //
 265             if (leading->LeadingFound()) {
 266                 Float_t etal = leading->GetLeading()->Eta();
 267                 Float_t phil = leading->GetLeading()->Phi();
 268                 Float_t el   = leading->GetLeading()->E();
 269 //        printf("Leading %f %f %f \n", etal, phil, el);
 270
 271                 e3H->Fill(el, wgt);
 272
 273                 Float_t rmin = 1.e6;
 274                 Float_t etamin = 1.e6;
 275
 276                 Int_t igen = -1;
 277                 Int_t ngen = gjets->GetNJets();
 278                 for (Int_t j = 0; j < ngen; j++) {
 279                     etag = gjets->GetEta(j);
 280                     Float_t phig = gjets->GetPhi(j);
 281                     Float_t deta = etag-etal;
 282                     Float_t dphi = TMath::Abs(phig - phil);
 283                     if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
 284
 285                     Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
 286
 287                     if (r  < rmin) {
 288                         rmin = r;
 289                         igen = j;
 290                         etamin = deta;
 291                     }
 292                 }
 293                 if (egen > 20. && egen < 40.) {
 294                     dr3H->Fill(rmin);
 295                     eta1H->Fill(etag-etal, wgt);
 296
 297                 }
 298
 299                 if (el > 54. && el  < 66.) e6H->Fill(egen, wgt);
 300                 if (rmin < 0.3) {
 301                     r7H->Fill(egen, el/egen, wgt);
 302                     r8H->Fill(el / 0.2, el/egen, wgt);
 303                 }
 304
 305                 if (rmin < 0.1) {
 306                     drP2->Fill(egen, etamin, wgt);
 307                 }
 308             } // if leading particle
 309
 310
 311
 312 //
 313 // Generated Jet
 314 //
 315             Int_t ngen = gjets->GetNJets();
 316             emax = 0.;
 317             imax = -1;
 318
 319             for (Int_t j = 0; j < ngen; j++) {
 320                 if (gjets->GetPt(j) >
 321                     emax && TMath::Abs(gjets->GetEta(j)) < 0.5) {
 322                     emax = gjets->GetPt(j);
 323                     imax = j;
 324                 }
 325             }
 326             if (imax != -1) e7H->Fill(emax, wgt);
 327
 328
 329             delete jetT;
 330
 331         } // events
 332         if (jFile) jFile->Close();
 333         delete jFile;
 334
 335     } // runs
 336
 337 //  delete jFile;
 338 //  if (jFile) jFile->Close();
 339 /*
 340   TFile* f = new TFile("j.root", "recreate");
 341   e0H->Write();
 342   e1H->Write();
 343   e2H->Write();
 344   e3H->Write();
 345   e4H->Write();
 346   e7H->Write();
 347   e8H->Write();
 348   f->Close();
 349 */
 350
 351     // Get Control Plots
 352 //  gStyle->SetOptStat(0);
 353
 354     TCanvas* c1 = new TCanvas("c1");
 355   e0H->Draw();
 356   e1H->SetLineColor(2);
 357   e2H->SetLineColor(4);
 358   e3H->SetLineColor(5);
 359   e1H->Draw("same");
 360   e3H->Draw("same");
 361
 362
 363   TCanvas* c2 = new TCanvas("c2");
 364 //  dr1H->Draw();
 365   dr2H->SetLineColor(2);
 366   dr2H->Draw("");
 367
 368   TCanvas* c3 = new TCanvas("c3");
 369   dr2H->Draw();
 370   dr3H->Draw("same");
 371
 372   TCanvas* c4 = new TCanvas("c4");
 373   e0H->Draw();
 374
 375   TCanvas* c5 = new TCanvas("c5");
 376   etaH->SetXTitle("#eta_{rec} - #eta_{gen}");
 377
 378   etaH->Draw();
 379   eta1H->SetLineColor(2);
 380
 381   eta1H->Draw("same");
 382
 383   TCanvas* c5a = new TCanvas("c5a");
 384   eta1H->Draw();
 385
 386   TCanvas* c5b = new TCanvas("c5b");
 387   eta2H->Draw();
 388
 389   TCanvas* c6 = new TCanvas("c6");
 390   e4H->Draw();
 391   TCanvas* c7 = new TCanvas("c7");
 392   phiH->Draw();
 393
 394   TCanvas* c7a = new TCanvas("c7a");
 395   phi1H->Draw();
 396   TCanvas* c7b = new TCanvas("c7b");
 397   phi2H->Draw();
 398
 399   TCanvas* c8 = new TCanvas("c8");
 400   e5H->SetXTitle("E_{gen} (GeV)");
 401   e5H->Draw();
 402   e6H->SetLineColor(2);
 403   e6H->Draw("same");
 404
 405   TCanvas* c9 = new TCanvas("c9");
 406   e6H->Draw();
 407
 408
 409
 410   TCanvas* c10 = new TCanvas("c10");
 411
 412
 413   r5H->SetMaximum(1.);
 414   r5H->Draw();
 415   r5H->SetXTitle("E_{gen} (GeV)");
 416   r5H->SetYTitle("E_{leading} / E_{gen}");
 417   r6H->SetLineColor(2);
 418   r6H->Draw("same");
 419
 420   TCanvas* c11 = new TCanvas("c11");
 421 //
 422 // Leading particle rec/gen
 423 //
 424   r7H->SetMaximum(1.);
 425
 426   r7H->Draw();
 427   r7H->SetXTitle("E_{gen} (GeV)");
 428   r7H->SetYTitle("E_{leading} / E_{gen}");
 429
 430   r8H->SetLineColor(2);
 431   r8H->Draw("same");
 432
 433   TCanvas* c12 = new TCanvas("c12");
 434   drP1->SetXTitle("E_{gen} (GeV)");
 435   drP1->SetYTitle("#eta_{rec} - #eta_{gen}");
 436   drP1->Draw();
 437
 438   TCanvas* c13 = new TCanvas("c13");
 439   drP2->SetXTitle("E_{gen} (GeV)");
 440   drP2->SetYTitle("#eta_{leading} - #eta_{gen}");
 441
 442   drP2->Draw();
 443
 444   TCanvas* c14 = new TCanvas("c14");
 445   dphiH->Draw();
 446
 447 /*
 448   e1H->Write();
 449   e2H->Write();
 450   e3H->Write();
 451   e4H->Write();
 452 */
 453 }
 454