EMCAL/jetDst.C

   1 //*-- Author: Aleksei Pavlinov (WSU)
   2 //*--
   3 #include "jetDst.h"
   4 #if !defined(__CINT__) || defined(__MAKECINT__)
   5 // These files are needed for compilation only
   6 // (see directory /auto/u/pavlinov/macros on PDSF)
   7 // Using with caution if you are not PAI
   8 #include "macroIncludePai.h"
   9 #include "macroIncludeAlice.h"
  10 #endif
  11
  12 // For convinience only
  13 class AliEMCALJetFinder; AliEMCALJetFinder* jetFinder=0;
  14 TH2F *hLego, *hLegoW, *hTrig;
  15 TH1F *hETrigPatch, *hJetPartPt;
  16 TCanvas *c1=0;
  17 // directory with file and number of files
  18 TString dirIn("/auto/alice/sahal/PPR3/11/"), nFileIn;
  19 Int_t file1=1, file2=20, nevMax=0;
  20 // change dirOut if you want.
  21 TString dirOut("RES/FILE/"), nFileOut;
  22 TFile *file=0;
  23 // Default value for jet fineer - 11-feb-2202
  24 Int_t   debugJet=1;
  25 Float_t coneRadius = 0.3;
  26 Float_t etSeed     = 4.;
  27 Float_t minJetEt   = 40.;
  28 Float_t minCellEt  = 1.; // be carefull
  29 Float_t minPtCut   = 2.; // for including charge particles to file
  30 // For BG subtraction
  31 Int_t modeSubtraction = 1;
  32 Float_t minMove = 0.05;
  33 Float_t maxMove = 0.15;
  34 Float_t precBg  = 0.035;
  35 // key 13-feb-2002
  36 Bool_t hadronCorrection = kTRUE;
  37 //
  38 class AliEMCALJetMicroDst; AliEMCALJetMicroDst *dst=0;
  39 Int_t nevSum = 0, nevF;
  40 TArrayI nevInFiles(100);
  41
  42 extern "C++" {void loadlibs();}
  43
  44 void jetDst(Int_t mode, Int_t var, Int_t nmax)
  45 {
  46 //
  47 // mode - variable for selection input mode (see function defineInput());
  48 // var  - variable for selection jetfinder mode (see function defineJetFinder()).
  49 //
  50 // Dynamically link some shared libs
  51     if (gClassTable->GetID("AliRun") < 0) {
  52         gROOT->LoadMacro("../macros/loadlibs.C");
  53         loadlibs();
  54     }
  55
  56     if(!defineInput(mode)) return;
  57
  58     dst = new AliEMCALJetMicroDst;
  59     dst->SetDebug(1);
  60     if(!nFileOut.Contains(".root")){
  61       nFileOut += var;
  62       nFileOut += ".root";
  63     }
  64     dst->Create(nFileOut.Data());
  65     gROOT->GetListOfBrowsables()->Add(dst);
  66
  67     for(Int_t j=file1; j<=file2; j++){
  68       nevF = 0;
  69       if(file2>1) {// multiple file
  70         nFileIn = dirIn;
  71         if       (mode>=11 && mode<=11) {
  72           nFileIn += j;
  73           nFileIn += ("/galice.root");
  74         } else if(mode>=21 && mode<=51){
  75           nFileIn += "galice";
  76           if(j>0) nFileIn += j;
  77           nFileIn += ".root";
  78         }
  79       }
  80 // Connect the Root Galice file containing Geometry, Kine and Hits
  81       file = new TFile(nFileIn.Data(), "READ");
  82       printf("\n FILE %i -> %s\n", j, nFileIn.Data());
  83
  84       if (file == 0) continue;
  85
  86 // Get AliRun object from file or create it if not on file
  87       gAlice = 0;
  88       gAlice = (AliRun*)(file->Get("gAlice"));
  89       if(!gAlice)  return;
  90       printf("AliRun object found on file\n");
  91
  92       defineJetFinder(mode,var);
  93 //
  94 //   Loop over events in file
  95 //
  96
  97 //      if(!c1) c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
  98       Int_t nhit=0;
  99       TString sw;
 100       if(nmax>0) nevMax = nmax;
 101       for (Int_t nev = 0; nev<= nevMax; nev++) {
 102         jetFinder->ResetJets();
 103         Int_t nparticles = gAlice->GetEvent(nev);
 104         if (nparticles <= 0) break; // last event
 105         nevSum++;
 106         nevF++;
 107 // Partons jet
 108 /*
 109         jetFinder->FillFromPartons();
 110         jetFinder->Find();
 111         for(Int_t ij=0; ij<jetFinder->Njets(); ij++){
 112           hJetPartPt->Fill(jetFinder->JetEnergy(ij));
 113         }
 114         jetFinder->ResetJets();
 115 */
 116 // ECAL information
 117         jetFinder->FillFromHits();
 118
 119         //        gROOT->cd(); be carefull with this staff
 120         hLegoW = new TH2F((*jetFinder->GetLego())); // for saving
 121         sw = hLegoW->GetName();
 122         sw += nev;
 123         hLegoW->SetName(sw.Data());
 124         sw = "Energy in ECAL:Event";
 125         sw += nev;
 126         hLegoW->SetTitle(sw.Data());
 127
 128         if(hTrig==0) {
 129            hLego = jetFinder->GetLego();
 130            hTrig = bookTrigHist(hLego);
 131         }
 132         fillTriggerPatch(hLegoW, hTrig, hETrigPatch);
 133
 134 //      jetFinder->FillFromDigits();
 135 // TPC  information
 136         jetFinder->FillFromTracks(1, 0);
 137 //                                ^ preserves info from hit
 138
 139 // TPC  information from Hits associated to the EMCAL
 140 //      jetFinder->FillFromHitFlaggedTracks(1);
 141 //
 142         jetFinder->Find();
 143         dst->Fill(gAlice, jetFinder);
 144 //
 145 // Look at results
 146         printf("\n Events : %d  Jets : %d \n", nev, jetFinder->Njets());
 147         Int_t njet = jetFinder->Njets();
 148         for (Int_t nj=0; nj<njet; nj++)
 149         {
 150             printf("\n Jet Energy %5d %8.2f \n",
 151                    nj, jetFinder->JetEnergy(nj));
 152             printf("Jet Phi    %5d %8.2f %8.2f \n",
 153                    nj, jetFinder->JetPhiL(nj), jetFinder->JetPhiW(nj));
 154             printf("Jet Eta    %5d %8.2f %8.2f \n",
 155                    nj, jetFinder->JetEtaL(nj), jetFinder->JetEtaW(nj));
 156         }
 157
 158         Int_t in;
 159         if (c1 && nev <-10 && njet==1 && jetFinder->JetEnergy(0)>150) {
 160            c1->cd();
 161            cout << "Event "<< nev <<endl;
 162            hLego->Draw("lego");
 163            c1->Update();
 164            cin  >> in;
 165            if(in<=0) return;
 166         }
 167       } // event loop
 168       nevInFiles[j-1] = nevF;
 169       if (gAlice) {
 170          delete gAlice;  // don't delete after closing file
 171          gAlice = 0;
 172       }
 173       file->Close();
 174       if(c1) hLego->Draw("lego");
 175     } // file loop
 176     if(nmax==0) dst->Close();
 177     nevInFiles.Set(file2);
 178     printf("\n Summ.events %6i -> %6i\n", nevSum, (Int_t) nevInFiles.GetSum());
 179
 180     printf("Events    ");
 181     for(Int_t j=file1; j<=file2; j++){
 182       printf(" %i4 ",nevInFiles[j-1]);
 183       if((j-file1)%10==0 && j!=file1) printf("\n          ");
 184     }
 185 }
 186
 187 // ###################################################################
 188 TString* defineInput(Int_t mode)
 189 {// 7-feb-2002
 190   Char_t ctmp[200];
 191   switch(mode){
 192 // case -5 needing for testing after merging
 193   case -5:
 194   case -4:
 195     hadronCorrection = kFALSE; // for checking Tom's idea
 196   case -3:
 197   case -2:
 198   case -1:
 199     if     (mode==-2) minPtCut = 1.;
 200     else if(mode<=-3 || mode>=-5) {
 201       minCellEt = 0.0; // 13-feb-2002
 202       minPtCut  = 0.0;
 203     }
 204
 205     coneRadius = 0.5;
 206     minJetEt   = 20.0;
 207     file1=1;
 208     file2=1;
 209     modeSubtraction = 0;
 210     nevMax = 1;
 211     nevMax = 999;
 212     dirIn    = ".";
 213     nFileIn  = "galice.root";
 214     sprintf(ctmp,"PythiaMode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
 215             mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
 216     nFileOut = dirOut + ctmp;
 217     if(!hadronCorrection) nFileOut.ReplaceAll(".root","_NoHadrCorr.root");
 218     break;
 219   case 1:
 220     nFileOut = dirOut + "dstTest.root";
 221     break;
 222
 223   case 11:
 224     // jj100 41 events per file - see Log file
 225     dirIn = "/auto/alice/sahal/PPR4/1/";
 226     file1 = 1;
 227     file2 = 30; // max 75
 228     nevMax = 41;
 229     sprintf(ctmp,"HijingMode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
 230             mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
 231     nFileOut = dirOut + ctmp;
 232     break;
 233     //=============================================================
 234     // central events with jet trigger - 15 feb 2002
 235     //=============================================================
 236   case 21:
 237     dirIn  = "/auto/alice/sahal/kHijing_jj50/";
 238     file1  = 0;
 239     file2  = 9;
 240     nevMax = 41;
 241     sprintf(ctmp,"Hijing_jj50Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
 242             mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
 243     nFileOut = dirOut + ctmp;
 244     break;
 245   case 22:
 246     dirIn  = "/auto/alice/sahal/kHijing_jj75/";
 247     file1  = 0;
 248     file2  = 9;
 249     nevMax = 41;
 250     sprintf(ctmp,"Hijing_jj75Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
 251             mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
 252     nFileOut = dirOut + ctmp;
 253     break;
 254   case 23:
 255     dirIn  = "/auto/alice/sahal/kHijing_jj100/";
 256     file1  = 0;
 257     file2  = 4;
 258     nevMax = 41;
 259     sprintf(ctmp,"Hijing_jj100Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
 260             mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
 261     nFileOut = dirOut + ctmp;
 262     break;
 263   case 24:
 264     dirIn  = "/auto/alice/sahal/kHijing_jj125/";
 265     file1  = 0;
 266     file2  = 4;
 267     nevMax = 41;
 268     sprintf(ctmp,"Hijing_jj125Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
 269             mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
 270     nFileOut = dirOut + ctmp;
 271     break;
 272
 273   case 51:
 274     // BG - HIJING events
 275     dirIn = "/auto/alice/sahal/background/";
 276     file1 = 1;
 277     file2 = 9;
 278     nevMax = 100; // ??
 279     sprintf(ctmp,"HijingBGMode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
 280             mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
 281     nFileOut = dirOut + ctmp;
 282     break;
 283   default:
 284      printf("\n<I> Wrong input mode => %i\ncout", mode);
 285      return 0;
 286   }
 287   printf("<I> dir %s file1 %3i file2 %3i\n", dirIn.Data(), file1, file2);
 288   return &nFileOut;
 289 }
 290
 291 void  defineJetFinder(Int_t mode, Int_t var)
 292 {// 11-feb-2002
 293 //  Create and configure jet finder
 294 //  gAlice must be defined already !!!
 295   if(!jetFinder) {
 296     jetFinder =  new AliEMCALJetFinder("UA1 Jet Finder", "Pre Production");
 297     jetFinder->Init();
 298
 299     jetFinder->SetDebug(debugJet);
 300     jetFinder->SetConeRadius(coneRadius);
 301     jetFinder->SetEtSeed(etSeed);
 302     jetFinder->SetMinJetEt(minJetEt);
 303     jetFinder->SetMinCellEt(minCellEt);
 304     jetFinder->SetPtCut(minPtCut);
 305
 306     if(var>=2) jetFinder->SetMomentumSmearing(1);
 307     if(var>=3) jetFinder->SetEfficiencySim(1);
 308
 309     jetFinder->SetHadronCorrection(0);
 310     if(var>=4 && hadronCorrection) {
 311       jetFinder->SetHadronCorrection(1);
 312       jetFinder->SetHadronCorrector(AliEMCALHadronCorrectionv0::Instance());
 313     }
 314     if(var>=5) { // 11-feb-2002
 315       jetFinder->SetParametersForBgSubtraction
 316       (modeSubtraction,minMove,maxMove,precBg);
 317     }
 318
 319     if(mode<0) jetFinder->SetSamplingFraction(1.);//must use only for PYTHIA file
 320     else       jetFinder->SetSamplingFraction(12.9);
 321     gROOT->GetListOfBrowsables()->Add(jetFinder);
 322   }
 323 }
 324
 325 void   testDST(char* fname)
 326 {// 8-feb-2002
 327   dst = new AliEMCALJetMicroDst;
 328   dst->Open(fname);
 329   dst->Test();
 330   gROOT->GetListOfBrowsables()->Add(dst);
 331 }
 332
 333 //----------------- 4-feb-2002 --------------
 334 TH2F* bookTrigHist(TH2F *hid)
 335 {
 336   if(hid==0) return 0;
 337   TAxis *xax  = hid->GetXaxis();
 338   Double_t xmin = xax->GetXmin();
 339   Double_t xmax = xax->GetXmax();
 340   Int_t nx = 6;
 341   TAxis *yax  = hid->GetYaxis();
 342   Double_t ymin = yax->GetXmin();
 343   Double_t ymax = yax->GetXmax();
 344   Int_t ny = 6;
 345   TH2F *hout = new TH2F("hTrig","Energy grid for trigger patch",
 346   nx,xmin,xmax, ny,ymin,ymax);
 347   hout->AddDirectory(0);
 348
 349   // time solution
 350   hETrigPatch = new TH1F("hETrigPatch","Energy in trigger patch",300, 0.0, 300.);
 351   hJetPartPt  = new TH1F("hJetPartPt", "E_{T} for partons jet", 300, 0.0, 300.);
 352   hETrigPatch->AddDirectory(0);
 353   return hout;
 354 }
 355
 356 void fillTriggerPatch(TH2F *hid, TH2F *htrig, TH1F *hE)
 357 {
 358   TAxis *xax  = hid->GetXaxis();
 359   TAxis *yax  = hid->GetYaxis();
 360   Int_t nx = xax->GetNbins();
 361   Int_t ny = yax->GetNbins();
 362
 363   htrig->Reset();
 364   Double_t x, y, e;
 365   for(Int_t ix=1; ix<=nx; ix++){
 366      x = xax->GetBinCenter(ix);
 367      for(Int_t iy=1; iy<=ny; iy++){
 368         y = yax->GetBinCenter(iy);
 369         e = hid->GetBinContent(ix, iy);
 370         htrig->Fill(x,y,e);
 371      }
 372   }
 373
 374   xax  = htrig->GetXaxis();
 375   yax  = htrig->GetYaxis();
 376   nx = xax->GetNbins();
 377   ny = yax->GetNbins();
 378   for(Int_t ix=1; ix<=nx; ix++){
 379      for(Int_t iy=1; iy<=ny; iy++){
 380        hE->Fill(htrig->GetBinContent(ix, iy));
 381      }
 382   }
 383 }