1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //---------------------------------------------------------------------
18 // Analyse Jets (already found jets)
19 // Authors: andreas.morsch@cern.ch, jgcn@mail.cern.ch
20 // mercedes.lopez.noriega@cern.ch
21 //---------------------------------------------------------------------
23 #include <Riostream.h>
24 #include "AliJetAnalysis.h"
25 ClassImp(AliJetAnalysis)
28 #include <Riostream.h>
36 #include <TLorentzVector.h>
39 #include "AliJetProductionDataPDC2004.h"
41 #include "AliUA1JetHeader.h"
42 #include "AliLeading.h"
43 #include "AliJetReaderHeader.h"
45 AliJetAnalysis::AliJetAnalysis():
146 // Default constructor
147 fFile = "anaJets.root";
149 // initialize weight for dE/dr histo
153 AliJetAnalysis::~AliJetAnalysis()
159 ////////////////////////////////////////////////////////////////////////
160 // define histogrames
162 void AliJetAnalysis::DefineHistograms()
164 // Define the histograms to be filled
165 if (fDoKine) DefineKineH();
166 if (fDoCorr) DefineCorrH();
167 if (fDoCorr50) DefineCorr50H();
168 if (fDoShap) DefineShapH();
169 if (fDoFrag) DefineFragH();
170 if (fDoTrig) DefineTrigH();
171 if (fDoJt) DefineJtH();
172 if (fDodNdxi) DefinedNdxiH();
175 void AliJetAnalysis::DefineKineH()
179 fPKineEneH = new TH1F("PKineEne","Energy of leading particle",50,0,200);
180 SetProperties(fPKineEneH,"Energy (GeV)","Entries");
181 fPKinePtH = new TH1F("PKinePt","Pt of leading particle",50,0,200);
182 SetProperties(fPKinePtH,"P_{T} (GeV)","Entries");
183 fPKinePhiH = new TH1F("PKinePhiH","Azimuthal angle of leading particle",
184 90,0.,2.0*TMath::Pi());
185 SetProperties(fPKinePhiH,"#phi","Entries");
186 fPKineEtaH = new TH1F("PKineEtaH","Pseudorapidity of leading particle",
188 SetProperties(fPKineEtaH,"#eta","Entries");
190 // leading generated jet
192 fGKineEneH = new TH1F("GKineEne","Energy of generated jet",50,0,200);
193 SetProperties(fGKineEneH,"Energy (GeV)","Entries");
194 fGKinePtH = new TH1F("GKinePt","Pt of generated jet",50,0,200);
195 SetProperties(fGKinePtH,"P_{T} (GeV)","Entries");
196 fGKinePhiH = new TH1F("GKinePhiH","Azimuthal angle of generated jet",
197 90,0.,2.0*TMath::Pi());
198 SetProperties(fGKinePhiH,"#phi","Entries");
199 fGKineEtaH = new TH1F("GKineEtaH","Pseudorapidity of generated jet",
201 SetProperties(fGKineEtaH,"#eta","Entries");
203 // leading reconstructed jet
205 fRKineEneH = new TH1F("RKineEne","Energy of reconstructed jet",50,0,200);
206 SetProperties(fRKineEneH,"Energy (GeV)","Entries");
207 fRKinePtH = new TH1F("RKinePt","Pt of reconstructed jet",50,0,200);
208 SetProperties(fRKinePtH,"P_{T} (GeV)","Entries");
209 fRKinePhiH = new TH1F("RKinePhiH","Azimuthal angle of reconstructed jet",
210 90,0.,2.0*TMath::Pi());
211 SetProperties(fRKinePhiH,"#phi","Entries");
212 fRKineEtaH = new TH1F("RKineEtaH","Pseudorapidity of reconstructed jet",
214 SetProperties(fRKineEtaH,"#eta","Entries");
218 void AliJetAnalysis::DefineCorrH()
221 if (fDoPart && fDoGenJ) {
222 fPGCorrEneH = new TH2F("PGCorrEne","Energy correlation part-gen jet",
224 SetProperties(fPGCorrEneH,"Part Energy (GeV)","Gen Jet Energy (GeV)");
225 fPGCorrPtH = new TH2F("PGCorrPt","Pt correlation part-gen jet",
227 SetProperties(fPGCorrPtH,"Part P_{T} (GeV)","Gen Jet P_{T} (GeV)");
228 fPGCorrEtaH = new TH2F("PGCorrEta","Pseudorapidity correlation part-gen jet",
229 40,-1.0,1.0,40,-1.0,1.0);
230 SetProperties(fPGCorrEtaH,"Part #eta","Gen Jet #eta");
231 fPGCorrPhiH = new TH2F("PGCorrPhi","Azimuthal angle correlation part-gen jet",
232 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
233 SetProperties(fPGCorrPhiH,"Part #phi","Gen Jet #phi");
235 if (fDoPart && fDoRecJ) {
236 fPRCorrEneH = new TH2F("PRCorrEne","Energy correlation part-rec jet",
238 SetProperties(fPRCorrEneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
239 fPRCorrPtH = new TH2F("PRCorrPt","Pt correlation part-rec jet",
241 SetProperties(fPRCorrPtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
242 fPRCorrEtaH = new TH2F("PRCorrEta","Pseudorapidity correlation part-rec jet",
243 40,-1.0,1.0,40,-1.0,1.0);
244 SetProperties(fPRCorrEtaH,"part #eta","Rec Jet #eta");
245 fPRCorrPhiH = new TH2F("PRCorrPhi","Azimuthal angle correlation part-rec jet",
246 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
247 SetProperties(fPRCorrPhiH,"Part #phi","Rec Jet #phi");
249 if (fDoGenJ && fDoRecJ) {
250 fRGCorrEneH = new TH2F("RGCorrEne","Energy correlation rec jet-gen jet",
252 SetProperties(fRGCorrEneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
253 fRGCorrPtH = new TH2F("RGCorrPt","Pt correlation rec jet-gen jet",
255 SetProperties(fRGCorrPtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
256 fRGCorrEtaH = new TH2F("RGCorrEta","Pseudorapidity correlation rec jet-gen jet",
257 40,-1.0,1.0,40,-1.0,1.0);
258 SetProperties(fRGCorrEtaH,"Rec Jet #eta","Gen Jet #eta");
259 fRGCorrPhiH = new TH2F("RGCorrPhi","Azimuthal angle correlation rec jet-gen jet",
260 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
261 SetProperties(fRGCorrPhiH,"Rec Jet #phi","Gen Jet #phi");
265 void AliJetAnalysis::DefineCorr50H()
268 if (fDoPart && fDoRecJ) {
269 fPRCorr50EneH = new TH2F("PRCorr50Ene","Energy correlation part-rec jet",
271 SetProperties(fPRCorr50EneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
272 fPRCorr50PtH = new TH2F("PRCorr50Pt","Pt correlation part-rec jet",
274 SetProperties(fPRCorr50PtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
275 fPRCorr50EtaH = new TH2F("PRCorr50Eta","Pseudorapidity correlation part-rec jet",
276 40,-1.0,1.0,40,-1.0,1.0);
277 SetProperties(fPRCorr50EtaH,"part #eta","Rec Jet #eta");
278 fPRCorr50PhiH = new TH2F("PRCorr50Phi","Azimuthal angle correlation part-rec jet",
279 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
280 SetProperties(fPRCorr50PhiH,"Part #phi","Rec Jet #phi");
283 if (fDoGenJ && fDoRecJ) {
284 fRGCorr50EneH = new TH2F("RGCorr50Ene","Energy correlation rec jet-gen jet",
286 SetProperties(fRGCorr50EneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
287 fRGCorr50PtH = new TH2F("RGCorr50Pt","Pt correlation rec jet-gen jet",
289 SetProperties(fRGCorr50PtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
290 fRGCorr50EtaH = new TH2F("RGCorr50Eta","Pseudorapidity correlation rec jet-gen jet",
291 40,-1.0,1.0,40,-1.0,1.0);
292 SetProperties(fRGCorr50EtaH,"Rec Jet #eta","Gen Jet #eta");
293 fRGCorr50PhiH = new TH2F("RGCorr50Phi","Azimuthal angle correlation rec jet-gen jet",
294 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
295 SetProperties(fRGCorr50PhiH,"Rec Jet #phi","Gen Jet #phi");
299 void AliJetAnalysis::DefineShapH()
301 // leading reconstructed jet
303 fdEdrH = new TH2F("fdEdrH","dE/dr histo",20,0,1,40,0,200);
304 SetProperties(fdEdrH,"r","Rec Jet P_{T}");
305 fdEdrB = new TH2F("fdEdrB","dE/dr bkgdhisto",20,0,1,40,0,200);
306 SetProperties(fdEdrB,"r","Rec P_{T}");
307 fPtEneH2 = new TH2F("fPtEneH2","fPtEneH2",40,0,200,40,0,200);
308 SetProperties(fPtEneH2,"Rec Jet E","Rec Jet P_{T}");
309 fdEdrW = new TH1F("fdEdrW","weights for dE/dr",40,0,200);
310 SetProperties(fdEdrW,"Rec Jet P_{T}","weights");
312 fRShapSelH = new TH1F("RShapSel","Shape of generated jets (sel part)",20,0.,1.);
313 SetProperties(fRShapSelH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
314 fRShapRejH = new TH1F("RShapRej","Shape of generated jets (rej part)",20,0.,1.);
315 SetProperties(fRShapRejH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
316 fRShapAllH = new TH1F("RShapAll","Shape of generated jets (all part)",20,0.,1.);
317 SetProperties(fRShapAllH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
321 void AliJetAnalysis::DefineJtH()
323 // Define the histogram for J_T
325 fJtH = new TH2F("fJtH","J_{T} histogram",80,0.,4.,40,0.,200.);
326 SetProperties(fJtH,"J_{T}","Rec Jet P_{T}");
327 fJtB = new TH2F("fJtB","J_{T} bkgd histogram",80,0.,4.,40,0.,200.);
328 SetProperties(fJtB,"J_{T}","Rec P_{T}");
329 fJtW = new TH1F("fJtW","J_{T} weight",40,0,200);
330 SetProperties(fJtW,"J_{T}W","weight");
334 void AliJetAnalysis::DefinedNdxiH()
336 // Define the histogram for dN/dxi
338 fdNdxiH = new TH2F("fdNdxiH","dN/d#xi histo",200,0,10,40,0,200);
339 SetProperties(fdNdxiH,"xi","Rec Jet P_{T}");
340 fdNdxiB = new TH2F("fdNdxiB","dN/d#xi bkgd histo",200,0,10,40,0,200);
341 SetProperties(fdNdxiB,"xi","Rec P_{T}");
342 fdNdxiW = new TH1F("fdNdxiW","dN/d#xi histo",40,0,200);
343 SetProperties(fdNdxiW,"Rec Jet P_{T}","weights");
344 fPtEneH = new TH2F("fPtEneH","fPtEneH",40,0,200,40,0,200);
345 SetProperties(fPtEneH,"Rec Jet E","Rec Jet P_{T}");
349 void AliJetAnalysis::DefineFragH()
351 // leading reconstructed jet
353 fRFragSelH = new TH1F("RFragSel","Frag Fun of reconstructed jets (sel part)",20,0.,10.);
354 SetProperties(fRFragSelH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
355 fRFragRejH = new TH1F("RFragRej","Frag Fun of reconstructed jets (rej part)",20,0.,10.);
356 SetProperties(fRFragRejH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
357 fRFragAllH = new TH1F("RFragAll","Frag Fun of reconstructed jets (all part)",20,0.,10.);
358 SetProperties(fRFragAllH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
362 void AliJetAnalysis::DefineTrigH()
365 fGTriggerEneH = new TProfile("GTriggerEne","Generated energy (trigger bias)",
366 20, 0., 200., 0., 1., "S");
367 fGTriggerEneH->SetXTitle("E_{gen}");
368 fGTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
369 // reconstructed energy
370 fRTriggerEneH = new TProfile("RTriggerEne","Reconstructed energy (trigger bias)",
371 20, 0., 200., 0., 1., "S");
372 fRTriggerEneH->SetXTitle("E_{rec}");
373 fRTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
374 // generated energy vs generated/leading
375 fGPTriggerEneH = new TProfile("GPTriggerEne","Generated energy (trigger bias)",
376 20, 0., 200., 0., 1., "S");
377 fGPTriggerEneH->SetXTitle("E_{gen}");
378 fGPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
379 // leading particle energy
380 fPTriggerEneH = new TProfile("PTriggerEne","Leading particle energy (trigger bias)",
381 20, 0., 200., 0., 1., "S");
382 fPTriggerEneH->SetXTitle("E_{L}/0.2");
383 fPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
387 void AliJetAnalysis::SetProperties(TH1* h,const char* x, const char* y) const
389 //Set properties of histos (x and y title and error propagation)
395 ////////////////////////////////////////////////////////////////////////
396 // compute weights for dE/dr histogram
398 void AliJetAnalysis::SetdEdrWeight()
400 // Due to the limited acceptance, each bin in the dE/dr has a different
401 // weight given by the ratio of the area of a disk dR to the area
402 // within the eta acceptance. The weight depends on the eta position
403 // of the jet. Here a look up table for the weights is computed. It
404 // assumes |etaJet|<0.5 and |eta_lego|<0.9. It makes bins of 0.05
405 // units in eta. Note that this is fixed by the bin width chosen for
406 // the histogram --> this weights are tailored for the specific
407 // histogram definition used here!
409 // two dimensional table: first index, bin in eta of jet, second
410 // index bin of dE/dr histo
415 Float_t binSize = (xMax-xMin)/nBins;
417 Float_t da,ds,r1,r2,h1,h2,a1,a2,theta1,theta2,area1,area2;
419 for (Int_t i=0;i<(nBins/2);i++) {
420 // index of first histo bin needing a scale factor
422 for(Int_t j=0;j<nBins;j++) {
424 da = TMath::Pi()*(binSize*binSize)*(1.0+2.0*j);
426 if (j>=ji) { // compute missing area using segments of circle
431 a1=2.0*TMath::Sqrt(2.0*h1*r1-h1*h1);
432 a2=2.0*TMath::Sqrt(2.0*h2*r2-h2*h2);
433 theta1=2*TMath::ACos((r1-h1)/r1);
434 theta2=2*TMath::ACos((r2-h2)/r2);
435 area1=binSize*(r1*r1*theta1-a1*(r1-h1));
436 area2=binSize*(r2*r2*theta2-a2*(r2-h2));
437 ds = (da-(area2-area1))/da;
439 fWeightdEdr[i][j]=ds/da;
444 // get weight for dE/dr histogram
445 Float_t AliJetAnalysis::GetdEdrWeight(Float_t etaJet, Float_t r)
447 // Return the correponding weight for the dE/dr plot
451 Float_t binSize = (xMax-xMin)/nBins;
453 Float_t eta = TMath::Abs(etaJet);
454 if ((eta > 0.5) || (r > fdrdEdr)) return 0.0;
455 Int_t binJet = (Int_t) TMath::Floor(eta/binSize);
456 Int_t binR = (Int_t) TMath::Floor(r/binSize);
457 Float_t w = fWeightdEdr[binJet][binR];
462 ////////////////////////////////////////////////////////////////////////
465 void AliJetAnalysis::FillHistograms()
470 AliJetProductionDataPDC2004* runData = new AliJetProductionDataPDC2004();
475 for (Int_t iRun = fRunMin; iRun <= fRunMax; iRun++) {
478 sprintf(fn, "%s/%s.root", fDirectory, (runData->GetRunTitle(iRun)).Data());
479 jFile = new TFile(fn);
480 printf(" Analyzing run: %d %s\n", iRun,fn);
482 // Get reader header and events to be looped over
483 AliJetReaderHeader *jReaderH = (AliJetReaderHeader*)(jFile->Get(fReaderHeader));
484 if (fEventMin == -1) fEventMin = jReaderH->GetFirstEvent();
485 if (fEventMax == -1) {
486 fEventMax = jReaderH->GetLastEvent();
488 fEventMax = TMath::Min(fEventMax, jReaderH->GetLastEvent());
491 AliUA1JetHeader *jH = (AliUA1JetHeader *) (jFile->Get("AliUA1JetHeader"));
494 fWeight = runData->GetWeight(iRun)/ ( (Float_t) (fEventMax - fEventMin + 1));
497 for (Int_t i = fEventMin; i < fEventMax; i++) {
498 if (i%100 == 0) printf(" Analyzing event %d / %d \n",i,fEventMax);
500 // Get next tree with AliJet
502 sprintf(nameT, "TreeJ%d",i);
503 TTree *jetT =(TTree *)(jFile->Get(nameT));
504 if (fDoRecJ) jetT->SetBranchAddress("FoundJet", &fRecJ);
505 if (fDoGenJ) jetT->SetBranchAddress("GenJet", &fGenJ);
506 if (fDoPart) jetT->SetBranchAddress("LeadingPart", &fPart);
509 int nJets = fRecJ->GetNJets();
511 TArrayI inJet = fRecJ->GetInJet();
512 if(inJet.GetSize()>fminMult){ // removes events with low multiplicity
513 if (fDoKine) FillKineH();
514 if (fDoCorr) FillCorrH();
515 if (fDoCorr50) FillCorr50H();
516 if (fDoShap) FillShapH(jH->GetRadius());
517 if (fDoFrag) FillFragH();
518 if (fDoTrig) FillTrigH();
519 if (fDoJt) FillJtH();
520 if (fDodNdxi) FilldNdxiH();
522 delete jetT; // jet should be deleted before creating a new one
523 if(inJet.GetSize()>fminMult){ // removes events with low multiplicity
524 if (fDoBkgd && nJets>0) FillBkgd(i,iRun);
526 } // end loop over events in one file
527 if (jFile) jFile->Close();
529 } // end loop over files
532 void AliJetAnalysis::FillKineH()
535 if (fDoPart && fPart->LeadingFound()){
536 fPKineEneH->Fill(fPart->GetE(),fWeight);
537 fPKinePtH->Fill(fPart->GetPt(),fWeight);
538 fPKinePhiH->Fill(fPart->GetPhi(),fWeight);
539 fPKineEtaH->Fill(fPart->GetEta(),fWeight);
541 // leading generated jet
542 if (fDoGenJ && fGenJ->GetNJets()> 0){
543 fGKineEneH->Fill(fGenJ->GetE(0),fWeight);
544 fGKinePtH->Fill(fGenJ->GetPt(0),fWeight);
545 fGKinePhiH->Fill(fGenJ->GetPhi(0),fWeight);
546 fGKineEtaH->Fill(fGenJ->GetEta(0),fWeight);
548 // leading reconstructed jet
549 if (fDoRecJ && fRecJ->GetNJets()> 0) {
550 TArrayF p=fRecJ->GetPtFromSignal();
551 if (p[0]>fPercentage) {
552 fRKineEneH->Fill(fRecJ->GetE(0)/fEfactor,fWeight);
553 fRKinePtH->Fill(fRecJ->GetPt(0),fWeight);
554 fRKinePhiH->Fill(fRecJ->GetPhi(0),fWeight);
555 fRKineEtaH->Fill(fRecJ->GetEta(0),fWeight);
560 void AliJetAnalysis::FillCorrH()
562 // Fill correlation histograms
563 if (fDoPart && fPart->LeadingFound() && fDoGenJ && fGenJ->GetNJets()> 0)
564 Correlation(fPart->GetLeading(),fGenJ->GetJet(0),
565 fPGCorrEneH,fPGCorrPtH,fPGCorrEtaH,fPGCorrPhiH);
566 if (fDoPart && fPart->LeadingFound() && fDoRecJ && fRecJ->GetNJets()> 0) {
567 TArrayF p=fRecJ->GetPtFromSignal();
568 if (p[0]>fPercentage)
569 Correlation(fPart->GetLeading(),fRecJ->GetJet(0),
570 fPRCorrEneH,fPRCorrPtH,fPRCorrEtaH,fPRCorrPhiH);
572 if (fDoGenJ && fGenJ->GetNJets()> 0 && fDoRecJ && fRecJ->GetNJets()> 0) {
573 TArrayF p=fRecJ->GetPtFromSignal();
574 if (p[0]>fPercentage)
575 Correlation(fRecJ->GetJet(0), fGenJ->GetJet(0),
576 fRGCorrEneH,fRGCorrPtH,fRGCorrEtaH,fRGCorrPhiH);
580 void AliJetAnalysis::Correlation(TLorentzVector *lv1,TLorentzVector *lv2,
581 TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
583 // Correlation histograms
584 h1->Fill(lv1->E(),lv2->E(),fWeight);
585 h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
586 h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
588 p1 = ((lv1->Phi() < 0) ? (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
589 p2 = ((lv2->Phi() < 0) ? (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
590 h4->Fill(p1,p2,fWeight);
593 void AliJetAnalysis::FillCorr50H()
595 // Fill correlation histograms when one particle has > 50% of jet energy
596 if (fDoRecJ && fRecJ->GetNJets()> 0) {
597 TArrayF p = fRecJ->GetPtFromSignal();
598 if (p[0]>fPercentage) {
599 if (fDoPart && fPart->LeadingFound())
600 Correlation50(fRecJ, fPart->GetLeading(),fRecJ->GetJet(0),
601 fPRCorr50EneH,fPRCorr50PtH,fPRCorr50EtaH,fPRCorr50PhiH);
602 if (fDoGenJ && fGenJ->GetNJets()> 0)
603 Correlation50(fRecJ, fRecJ->GetJet(0), fGenJ->GetJet(0),
604 fRGCorr50EneH,fRGCorr50PtH,fRGCorr50EtaH,fRGCorr50PhiH);
609 void AliJetAnalysis::Correlation50(AliJet *j,TLorentzVector *lv1,TLorentzVector *lv2,
610 TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
612 // Correlation histograms when one particle has > 50% of jet energy
613 TArrayF ptin = j->GetPtIn();
614 TArrayF etain = j->GetEtaIn();
615 TArrayI inJet = j->GetInJet();
618 for(Int_t i=0;i<(inJet.GetSize());i++) {
620 Float_t t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
621 Float_t x = (ptin[i]*TMath::Sqrt(1.+1./(t1*t1)))/(j->GetE(0));
625 if (flag>1) cout << " Flag = " << flag << endl;
627 h1->Fill(lv1->E(),lv2->E(),fWeight);
628 h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
629 h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
631 p1 = ((lv1->Phi() < 0) ?
632 (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
633 p2 = ((lv2->Phi() < 0) ?
634 (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
635 h4->Fill(p1,p2,fWeight);
639 void AliJetAnalysis::FillJtH()
641 // Fill J_T histogram
642 if (fRecJ->GetNJets()> 0) {
643 fjv3X = 0.0; fjv3Y = 0.0; fjv3Z = 0.0;
644 TArrayF p=fRecJ->GetPtFromSignal();
645 if (p[0]>fPercentage) {
647 const TVector3 kJv3 = fRecJ->GetJet(0)->Vect();
648 fjv3X = kJv3.X(); fjv3Y = kJv3.Y(); fjv3Z = kJv3.Z();
650 TArrayI inJet = fRecJ->GetInJet();
651 TArrayF etain = fRecJ->GetEtaIn();
652 TArrayF ptin = fRecJ->GetPtIn();
653 TArrayF phiin = fRecJ->GetPhiIn();
654 Float_t deta, dphi,jt, dr;
655 for(Int_t i=0;i<inJet.GetSize();i++) {
656 deta = etain[i] - fRecJ->GetEta(0);
657 dphi = phiin[i] - fRecJ->GetPhi(0);
658 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
659 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
660 dr = TMath::Sqrt(deta * deta + dphi * dphi);
661 if ((dr<fdrJt) && (ptin[i] > fPartPtCut)) {
662 trk.SetPtEtaPhi(ptin[i],etain[i],phiin[i]);
663 jt = TMath::Sin(trk.Angle(kJv3))*trk.Mag();
664 fJtH->Fill(jt,fRecJ->GetPt(0),fWeight);
667 fJtW->Fill(fRecJ->GetPt(0),fWeight);
673 void AliJetAnalysis::FilldNdxiH()
675 // Fill dN/d#xi histograms
676 if (fRecJ->GetNJets()> 0) {
677 TArrayF p=fRecJ->GetPtFromSignal();
678 if (p[0]>fPercentage) {
680 TArrayI inJet = fRecJ->GetInJet();
681 TArrayF etain = fRecJ->GetEtaIn();
682 TArrayF ptin = fRecJ->GetPtIn();
683 TArrayF phiin = fRecJ->GetPhiIn();
684 Float_t xi,t1,ene,dr,deta,dphi;
685 for(Int_t i=0;i<inJet.GetSize();i++) {
686 deta = etain[i] - fRecJ->GetEta(0);
687 dphi = phiin[i] - fRecJ->GetPhi(0);
688 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
689 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
690 dr = TMath::Sqrt(deta * deta + dphi * dphi);
691 if ((dr<fdrdNdxi) && (ptin[i] > fPartPtCut)) {
692 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
693 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
694 xi = (Float_t) TMath::Log((fRecJ->GetE(0)/fEfactor)/ene);
695 fdNdxiH->Fill(xi,fRecJ->GetPt(0),fWeight);
698 fdNdxiW->Fill(fRecJ->GetPt(0),fWeight);
699 fPtEneH->Fill(fRecJ->GetE(0)/fEfactor,fRecJ->GetPt(0),fWeight);
705 void AliJetAnalysis::FillBkgd(Int_t eventN, Int_t runN)
707 // Background calculated with hijing events (no pythia)
708 if (fp0>fPercentage) {
709 fPtJ=0.,fEJ=0.,fEtaJ=0.,fPhiJ=0.;
710 // Background calculated with hijing events (no pythia)
711 AliJetProductionDataPDC2004* runDataB = new AliJetProductionDataPDC2004();
715 sprintf(fnB, "%s/%s.root", fBkgdDirectory, (runDataB->GetRunTitle(runN)).Data());
716 jFileB = new TFile(fnB);
719 sprintf(nameB, "TreeJ%d",eventN);
720 TTree *jetB =(TTree *)(jFileB->Get(nameB));
721 jetB->SetBranchAddress("FoundJet", &fRecB);
724 TArrayI inJetB = fRecB->GetInJet();
725 TArrayF etainB = fRecB->GetEtaIn();
726 TArrayF ptinB = fRecB->GetPtIn();
727 TArrayF phiinB = fRecB->GetPhiIn();
728 fPtJ = fRecJ->GetPt(0);
729 fEJ = fRecJ->GetE(0);
730 fEtaJ = fRecJ->GetEta(0);
731 fPhiJ = fRecJ->GetPhi(0);
732 Float_t t1,ene,xi,detaB,dphiB,drB,jt,wB;
735 jv3B.SetX(fjv3X); jv3B.SetY(fjv3Y); jv3B.SetZ(fjv3Z);
737 for(Int_t k=0;k<inJetB.GetSize();k++){
738 if(ptinB[k] > fPartPtCut){
739 detaB = etainB[k] - fEtaJ;
740 dphiB = phiinB[k] - fPhiJ;
741 if (dphiB < -TMath::Pi()) dphiB= -dphiB - 2.0 * TMath::Pi();
742 if (dphiB > TMath::Pi()) dphiB = 2.0 * TMath::Pi() - dphiB;
743 drB = TMath::Sqrt(detaB * detaB + dphiB * dphiB);
744 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etainB[k])));
745 ene = ptinB[k]*TMath::Sqrt(1.+1./(t1*t1));
746 trkB.SetPtEtaPhi(ptinB[k],etainB[k],phiinB[k]);
749 xi = (Float_t) TMath::Log(fEJ/ene);
750 fdNdxiB->Fill(xi,fPtJ,fWeight);
754 jt = TMath::Sin(trkB.Angle(jv3B))*(trkB.Mag());
755 fJtB->Fill(jt,fPtJ,fWeight);
759 wB = GetdEdrWeight(fEtaJ,drB)*fWeight*ene;
760 fdEdrB->Fill(drB,fPtJ,wB);
765 if (jFileB) jFileB->Close();
770 void AliJetAnalysis::FillShapH(Float_t r)
772 // Fill jet shape histograms
773 if (fDoRecJ && fRecJ->GetNJets()> 0) {
774 TArrayF p=fRecJ->GetPtFromSignal();
775 if (p[0]>fPercentage) {
776 Shape(fRecJ,fRShapSelH,fRShapRejH,fRShapAllH,fdEdrH,fPtEneH2,fdEdrW,r);
783 void AliJetAnalysis::Shape(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha,
784 TH2F *hdedr, TH2F *hptene, TH1F* wdedr, Float_t radius)
787 TArrayI inJet = j->GetInJet();
788 TArrayF etain = j->GetEtaIn();
789 TArrayF ptin = j->GetPtIn();
790 TArrayF phiin = j->GetPhiIn();
792 // first compute dE/dr histo
793 Float_t etaj = j->GetEta(0);
794 Float_t ene,w,deta,dphi,dr,t1;
795 for(Int_t i=0;i<inJet.GetSize();i++) {
796 deta = etain[i] - j->GetEta(0);
797 dphi = phiin[i] - j->GetPhi(0);
798 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
799 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
800 dr = TMath::Sqrt(deta * deta + dphi * dphi);
801 if ((dr<fdrdEdr) && (ptin[i] > fPartPtCut)) {
802 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
803 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
804 w = GetdEdrWeight(etaj,dr)*fWeight*ene;
805 hdedr->Fill(dr,j->GetPt(0),w);
808 hptene->Fill(fRecJ->GetE(0),fRecJ->GetPt(0),fWeight);
809 wdedr->Fill(j->GetPt(0),fWeight);
811 // now compute shape histos
812 Int_t nBins = ha->GetNbinsX();
813 Float_t xMin = ha->GetXaxis()->GetXmin();
814 Float_t xMax = ha->GetXaxis()->GetXmax();
815 Float_t binSize,halfBin;
816 binSize = (xMax-xMin)/nBins;
818 Float_t rptS[20], rptR[20], rptA[20];
819 for(Int_t i=0;i<nBins;i++) rptS[i]=rptR[i]=rptA[i]=0.0;
820 // fill bins in r for leading jet
821 for(Int_t i=0;i<inJet.GetSize();i++) {
822 if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
823 deta = etain[i] - j->GetEta(0);
824 dphi = phiin[i] - j->GetPhi(0);
825 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
826 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
827 dr = TMath::Sqrt(deta * deta + dphi * dphi);
828 Float_t rR = dr/radius;
829 Int_t bin = (Int_t) TMath::Floor(rR/binSize);
830 rptA[bin]+=ptin[i]/(j->GetPt(0));
831 if (inJet[i] == 1) rptS[bin]+=ptin[i]/(j->GetPt(0));
832 if (fPythia && inJet[i] == -1)
833 rptR[bin]+=ptin[i]/(j->GetPt(0));
837 // compute shape and fill histogram
838 Float_t ptS,ptR,ptA,r;
840 for (Int_t i=0;i<nBins;i++) {
844 r=(i+1)*binSize-halfBin;
845 hs->Fill(r,ptS*fWeight);
847 hr->Fill(r,ptR*fWeight);
848 ha->Fill(r,ptA*fWeight);
853 void AliJetAnalysis::FillFragH()
855 // Fill fragmentation histogram
856 if (fDoRecJ && fRecJ->GetNJets()> 0) {
857 TArrayF p=fRecJ->GetPtFromSignal();
858 if (p[0]>fPercentage) {
859 FragFun(fRecJ,fRFragSelH,fRFragRejH,fRFragAllH);
865 void AliJetAnalysis::FragFun(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha)
867 // Calculate de fragmentation function
868 TArrayI inJet = j->GetInJet();
869 TArrayF etain = j->GetEtaIn();
870 TArrayF ptin = j->GetPtIn();
874 for(Int_t i=0;i<(inJet.GetSize());i++) {
875 if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
876 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
877 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
878 xi = (Float_t) TMath::Log((j->GetE(0))/ene);
879 if (fPythia) ha->Fill(xi,fWeight);
880 if (inJet[i] == 1) hs->Fill(xi,fWeight);
881 if (fPythia && inJet[i] == -1) hr->Fill(xi,fWeight);
886 void AliJetAnalysis:: FillTrigH()
888 // Fill trigger bias histograms
889 if(fDoTrig && fGenJ->GetNJets()>0 && fRecJ->GetNJets()>0){
890 TArrayF p=fRecJ->GetPtFromSignal();
891 if (p[0]>fPercentage) {
892 float genEne = fGenJ->GetE(0);
893 float recEne = fRecJ->GetE(0);
894 float eneRatio = (recEne/fEfactor)/genEne;
896 fGTriggerEneH->Fill(genEne, eneRatio, fWeight);
897 fRTriggerEneH->Fill(recEne/fEfactor, eneRatio, fWeight);
899 if (fPart->LeadingFound()){
900 float leaEne = fPart->GetE();
901 float eneRatio2 = leaEne/genEne;
902 fGPTriggerEneH->Fill(genEne, eneRatio2, fWeight);
903 fPTriggerEneH->Fill(leaEne/0.2, eneRatio2, fWeight);
910 ////////////////////////////////////////////////////////////////////////
911 // Normalize histogrames
913 void AliJetAnalysis::NormHistograms()
915 // normalize shape histograms
917 if (fDoRecJ && fWShapR>0.) { // leading reconstructed jet
918 fRShapSelH->Scale(1.0/fWShapR);
919 fRShapRejH->Scale(1.0/fWShapR);
920 fRShapAllH->Scale(1.0/fWShapR);
924 // normalize fragmentation function histograms
926 if (fDoRecJ && fWFragR>0.) { // leading reconstructed jet
927 fRFragSelH->Scale(2.0/fWFragR);
928 fRFragRejH->Scale(2.0/fWFragR);
929 fRFragAllH->Scale(2.0/fWFragR);
934 ////////////////////////////////////////////////////////////////////////
936 void AliJetAnalysis::PlotHistograms()
938 // Plot histogramas (to be done...)
939 if (fDoKine) PlotKineH();
940 if (fDoCorr) PlotCorrH();
941 if (fDoCorr50) PlotCorr50H();
942 if (fDoShap) PlotShapH();
943 if (fDoFrag) PlotFragH();
944 if (fDoTrig) PlotTrigH();
947 void AliJetAnalysis::PlotKineH() const
952 void AliJetAnalysis::PlotCorrH() const
956 void AliJetAnalysis::PlotCorr50H() const
961 void AliJetAnalysis::PlotShapH() const
966 void AliJetAnalysis::PlotFragH() const
971 void AliJetAnalysis::PlotTrigH()
977 ////////////////////////////////////////////////////////////////////////
979 void AliJetAnalysis::SaveHistograms()
982 TFile *fOut = new TFile(fFile,"recreate");
984 if (fDoKine) SaveKineH();
985 if (fDoCorr) SaveCorrH();
986 if (fDoCorr50) SaveCorr50H();
987 if (fDoShap) SaveShapH();
988 if (fDoFrag) SaveFragH();
989 if (fDoTrig) SaveTrigH();
990 if (fDoJt) SaveJtH();
991 if (fDodNdxi) SavedNdxiH();
995 void AliJetAnalysis::SaveKineH()
997 // Save kinematic histograms
1001 fPKinePhiH->Write();
1002 fPKineEtaH->Write();
1006 fGKineEneH->Write();
1008 fGKinePhiH->Write();
1009 fGKineEtaH->Write();
1013 fRKineEneH->Write();
1015 fRKinePhiH->Write();
1016 fRKineEtaH->Write();
1020 void AliJetAnalysis::SaveCorrH()
1022 // Save correlation histograms
1023 if (fDoPart && fDoGenJ) {
1024 fPGCorrEneH->Write();
1025 fPGCorrPtH->Write();
1026 fPGCorrEtaH->Write();
1027 fPGCorrPhiH->Write();
1030 if (fDoPart && fDoRecJ) {
1031 fPRCorrEneH->Write();
1032 fPRCorrPtH->Write();
1033 fPRCorrEtaH->Write();
1034 fPRCorrPhiH->Write();
1037 if (fDoGenJ && fDoRecJ) {
1038 fRGCorrEneH->Write();
1039 fRGCorrPtH->Write();
1040 fRGCorrEtaH->Write();
1041 fRGCorrPhiH->Write();
1045 void AliJetAnalysis::SaveCorr50H()
1047 // Save correlation histograms (special case)
1048 if (fDoPart && fDoRecJ) {
1049 fPRCorr50EneH->Write();
1050 fPRCorr50PtH->Write();
1051 fPRCorr50EtaH->Write();
1052 fPRCorr50PhiH->Write();
1054 if (fDoGenJ && fDoRecJ) {
1055 fRGCorr50EneH->Write();
1056 fRGCorr50PtH->Write();
1057 fRGCorr50EtaH->Write();
1058 fRGCorr50PhiH->Write();
1062 void AliJetAnalysis::SaveShapH()
1064 // Save jet shape histograms
1066 fRShapSelH->Write();
1068 if(fDoBkgd) fdEdrB->Write();
1072 fRShapRejH->Write();
1073 fRShapAllH->Write();
1078 void AliJetAnalysis::SaveJtH()
1080 // Save J_T histograms
1083 if(fDoBkgd) fJtB->Write();
1088 void AliJetAnalysis::SavedNdxiH()
1090 // Save dN/d#xi histograms
1093 if(fDoBkgd) fdNdxiB->Write();
1099 void AliJetAnalysis::SaveFragH()
1101 // Save fragmentation histograms
1103 fRFragSelH->Write();
1105 fRFragRejH->Write();
1106 fRFragAllH->Write();
1111 void AliJetAnalysis::SaveTrigH()
1113 // Save trigger bias histograms
1115 fGTriggerEneH->Write();
1116 fRTriggerEneH->Write();
1117 fGPTriggerEneH->Write();
1118 fPTriggerEneH->Write();
1122 ////////////////////////////////////////////////////////////////////////
1123 // main Analysis function
1125 void AliJetAnalysis::Analyze()
1130 // leading generated jet
1131 // leading reconstructed jet
1133 // Correlations amd resolutions
1134 // a) correlations in energy, pt, phi, eta
1135 // b) resolutions in energy, pt, phi, eta, r
1136 // leading particle and leading generated jet
1137 // leading particle and leading reconstructed jet
1138 // leading generated jet and leading reconstructed jet
1140 // Fragmentation functions and Shapes
1141 // a) integrated over all pt
1143 // b.1) only for user selected particles in jet
1144 // b.2) only for user rejected particles in jet
1145 // b.3) for all particles in jet