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 "AliJetKineReaderHeader.h"
42 #include "AliJetESDReaderHeader.h"
43 #include "AliUA1JetHeader.h"
44 #include "AliLeading.h"
47 AliJetAnalysis::AliJetAnalysis()
52 fFile = "anajets.root";
111 // correlation histograms
125 // correlation histograms when one particle
126 // has more than 50% of the energy of the jet
142 // fragmentation function histos
148 // trigger bias histos
154 // dE/dr histo and its dr
161 // jt histo and its dr
167 // dN/dxi histo and its dr
174 // initialize weight for dE/dr histo
178 AliJetAnalysis::~AliJetAnalysis()
184 ////////////////////////////////////////////////////////////////////////
185 // define histogrames
187 void AliJetAnalysis::DefineHistograms()
189 // Define the histograms to be filled
190 if (fDoKine) DefineKineH();
191 if (fDoCorr) DefineCorrH();
192 if (fDoCorr50) DefineCorr50H();
193 if (fDoShap) DefineShapH();
194 if (fDoFrag) DefineFragH();
195 if (fDoTrig) DefineTrigH();
196 if (fDoJt) DefineJtH();
197 if (fDodNdxi) DefinedNdxiH();
200 void AliJetAnalysis::DefineKineH()
204 fPKineEneH = new TH1F("PKineEne","Energy of leading particle",50,0,200);
205 SetProperties(fPKineEneH,"Energy (GeV)","Entries");
206 fPKinePtH = new TH1F("PKinePt","Pt of leading particle",50,0,200);
207 SetProperties(fPKinePtH,"P_{T} (GeV)","Entries");
208 fPKinePhiH = new TH1F("PKinePhiH","Azimuthal angle of leading particle",
209 90,0.,2.0*TMath::Pi());
210 SetProperties(fPKinePhiH,"#phi","Entries");
211 fPKineEtaH = new TH1F("PKineEtaH","Pseudorapidity of leading particle",
213 SetProperties(fPKineEtaH,"#eta","Entries");
215 // leading generated jet
217 fGKineEneH = new TH1F("GKineEne","Energy of generated jet",50,0,200);
218 SetProperties(fGKineEneH,"Energy (GeV)","Entries");
219 fGKinePtH = new TH1F("GKinePt","Pt of generated jet",50,0,200);
220 SetProperties(fGKinePtH,"P_{T} (GeV)","Entries");
221 fGKinePhiH = new TH1F("GKinePhiH","Azimuthal angle of generated jet",
222 90,0.,2.0*TMath::Pi());
223 SetProperties(fGKinePhiH,"#phi","Entries");
224 fGKineEtaH = new TH1F("GKineEtaH","Pseudorapidity of generated jet",
226 SetProperties(fGKineEtaH,"#eta","Entries");
228 // leading reconstructed jet
230 fRKineEneH = new TH1F("RKineEne","Energy of reconstructed jet",50,0,200);
231 SetProperties(fRKineEneH,"Energy (GeV)","Entries");
232 fRKinePtH = new TH1F("RKinePt","Pt of reconstructed jet",50,0,200);
233 SetProperties(fRKinePtH,"P_{T} (GeV)","Entries");
234 fRKinePhiH = new TH1F("RKinePhiH","Azimuthal angle of reconstructed jet",
235 90,0.,2.0*TMath::Pi());
236 SetProperties(fRKinePhiH,"#phi","Entries");
237 fRKineEtaH = new TH1F("RKineEtaH","Pseudorapidity of reconstructed jet",
239 SetProperties(fRKineEtaH,"#eta","Entries");
243 void AliJetAnalysis::DefineCorrH()
246 if (fDoPart && fDoGenJ) {
247 fPGCorrEneH = new TH2F("PGCorrEne","Energy correlation part-gen jet",
249 SetProperties(fPGCorrEneH,"Part Energy (GeV)","Gen Jet Energy (GeV)");
250 fPGCorrPtH = new TH2F("PGCorrPt","Pt correlation part-gen jet",
252 SetProperties(fPGCorrPtH,"Part P_{T} (GeV)","Gen Jet P_{T} (GeV)");
253 fPGCorrEtaH = new TH2F("PGCorrEta","Pseudorapidity correlation part-gen jet",
254 40,-1.0,1.0,40,-1.0,1.0);
255 SetProperties(fPGCorrEtaH,"Part #eta","Gen Jet #eta");
256 fPGCorrPhiH = new TH2F("PGCorrPhi","Azimuthal angle correlation part-gen jet",
257 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
258 SetProperties(fPGCorrPhiH,"Part #phi","Gen Jet #phi");
260 if (fDoPart && fDoRecJ) {
261 fPRCorrEneH = new TH2F("PRCorrEne","Energy correlation part-rec jet",
263 SetProperties(fPRCorrEneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
264 fPRCorrPtH = new TH2F("PRCorrPt","Pt correlation part-rec jet",
266 SetProperties(fPRCorrPtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
267 fPRCorrEtaH = new TH2F("PRCorrEta","Pseudorapidity correlation part-rec jet",
268 40,-1.0,1.0,40,-1.0,1.0);
269 SetProperties(fPRCorrEtaH,"part #eta","Rec Jet #eta");
270 fPRCorrPhiH = new TH2F("PRCorrPhi","Azimuthal angle correlation part-rec jet",
271 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
272 SetProperties(fPRCorrPhiH,"Part #phi","Rec Jet #phi");
274 if (fDoGenJ && fDoRecJ) {
275 fRGCorrEneH = new TH2F("RGCorrEne","Energy correlation rec jet-gen jet",
277 SetProperties(fRGCorrEneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
278 fRGCorrPtH = new TH2F("RGCorrPt","Pt correlation rec jet-gen jet",
280 SetProperties(fRGCorrPtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
281 fRGCorrEtaH = new TH2F("RGCorrEta","Pseudorapidity correlation rec jet-gen jet",
282 40,-1.0,1.0,40,-1.0,1.0);
283 SetProperties(fRGCorrEtaH,"Rec Jet #eta","Gen Jet #eta");
284 fRGCorrPhiH = new TH2F("RGCorrPhi","Azimuthal angle correlation rec jet-gen jet",
285 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
286 SetProperties(fRGCorrPhiH,"Rec Jet #phi","Gen Jet #phi");
290 void AliJetAnalysis::DefineCorr50H()
293 if (fDoPart && fDoRecJ) {
294 fPRCorr50EneH = new TH2F("PRCorr50Ene","Energy correlation part-rec jet",
296 SetProperties(fPRCorr50EneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
297 fPRCorr50PtH = new TH2F("PRCorr50Pt","Pt correlation part-rec jet",
299 SetProperties(fPRCorr50PtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
300 fPRCorr50EtaH = new TH2F("PRCorr50Eta","Pseudorapidity correlation part-rec jet",
301 40,-1.0,1.0,40,-1.0,1.0);
302 SetProperties(fPRCorr50EtaH,"part #eta","Rec Jet #eta");
303 fPRCorr50PhiH = new TH2F("PRCorr50Phi","Azimuthal angle correlation part-rec jet",
304 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
305 SetProperties(fPRCorr50PhiH,"Part #phi","Rec Jet #phi");
308 if (fDoGenJ && fDoRecJ) {
309 fRGCorr50EneH = new TH2F("RGCorr50Ene","Energy correlation rec jet-gen jet",
311 SetProperties(fRGCorr50EneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
312 fRGCorr50PtH = new TH2F("RGCorr50Pt","Pt correlation rec jet-gen jet",
314 SetProperties(fRGCorr50PtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
315 fRGCorr50EtaH = new TH2F("RGCorr50Eta","Pseudorapidity correlation rec jet-gen jet",
316 40,-1.0,1.0,40,-1.0,1.0);
317 SetProperties(fRGCorr50EtaH,"Rec Jet #eta","Gen Jet #eta");
318 fRGCorr50PhiH = new TH2F("RGCorr50Phi","Azimuthal angle correlation rec jet-gen jet",
319 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
320 SetProperties(fRGCorr50PhiH,"Rec Jet #phi","Gen Jet #phi");
324 void AliJetAnalysis::DefineShapH()
326 // leading reconstructed jet
328 fdEdrH = new TH2F("fdEdrH","dE/dr histo",20,0,1,40,0,200);
329 SetProperties(fdEdrH,"r","Rec Jet P_{T}");
330 fdEdrB = new TH2F("fdEdrB","dE/dr bkgdhisto",20,0,1,40,0,200);
331 SetProperties(fdEdrB,"r","Rec P_{T}");
332 fPtEneH2 = new TH2F("fPtEneH2","fPtEneH2",40,0,200,40,0,200);
333 SetProperties(fPtEneH2,"Rec Jet E","Rec Jet P_{T}");
334 fdEdrW = new TH1F("fdEdrW","weights for dE/dr",40,0,200);
335 SetProperties(fdEdrW,"Rec Jet P_{T}","weights");
337 fRShapSelH = new TH1F("RShapSel","Shape of generated jets (sel part)",20,0.,1.);
338 SetProperties(fRShapSelH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
339 fRShapRejH = new TH1F("RShapRej","Shape of generated jets (rej part)",20,0.,1.);
340 SetProperties(fRShapRejH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
341 fRShapAllH = new TH1F("RShapAll","Shape of generated jets (all part)",20,0.,1.);
342 SetProperties(fRShapAllH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
346 void AliJetAnalysis::DefineJtH()
348 // Define the histogram for J_T
350 fJtH = new TH2F("fJtH","J_{T} histogram",80,0.,4.,40,0.,200.);
351 SetProperties(fJtH,"J_{T}","Rec Jet P_{T}");
352 fJtB = new TH2F("fJtB","J_{T} bkgd histogram",80,0.,4.,40,0.,200.);
353 SetProperties(fJtB,"J_{T}","Rec P_{T}");
354 fJtW = new TH1F("fJtW","J_{T} weight",40,0,200);
355 SetProperties(fJtW,"J_{T}W","weight");
359 void AliJetAnalysis::DefinedNdxiH()
361 // Define the histogram for dN/dxi
363 fdNdxiH = new TH2F("fdNdxiH","dN/d#xi histo",200,0,10,40,0,200);
364 SetProperties(fdNdxiH,"xi","Rec Jet P_{T}");
365 fdNdxiB = new TH2F("fdNdxiB","dN/d#xi bkgd histo",200,0,10,40,0,200);
366 SetProperties(fdNdxiB,"xi","Rec P_{T}");
367 fdNdxiW = new TH1F("fdNdxiW","dN/d#xi histo",40,0,200);
368 SetProperties(fdNdxiW,"Rec Jet P_{T}","weights");
369 fPtEneH = new TH2F("fPtEneH","fPtEneH",40,0,200,40,0,200);
370 SetProperties(fPtEneH,"Rec Jet E","Rec Jet P_{T}");
374 void AliJetAnalysis::DefineFragH()
376 // leading reconstructed jet
378 fRFragSelH = new TH1F("RFragSel","Frag Fun of reconstructed jets (sel part)",20,0.,10.);
379 SetProperties(fRFragSelH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
380 fRFragRejH = new TH1F("RFragRej","Frag Fun of reconstructed jets (rej part)",20,0.,10.);
381 SetProperties(fRFragRejH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
382 fRFragAllH = new TH1F("RFragAll","Frag Fun of reconstructed jets (all part)",20,0.,10.);
383 SetProperties(fRFragAllH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
387 void AliJetAnalysis::DefineTrigH()
390 fGTriggerEneH = new TProfile("GTriggerEne","Generated energy (trigger bias)",
391 20, 0., 200., 0., 1., "S");
392 fGTriggerEneH->SetXTitle("E_{gen}");
393 fGTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
394 // reconstructed energy
395 fRTriggerEneH = new TProfile("RTriggerEne","Reconstructed energy (trigger bias)",
396 20, 0., 200., 0., 1., "S");
397 fRTriggerEneH->SetXTitle("E_{rec}");
398 fRTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
399 // generated energy vs generated/leading
400 fGPTriggerEneH = new TProfile("GPTriggerEne","Generated energy (trigger bias)",
401 20, 0., 200., 0., 1., "S");
402 fGPTriggerEneH->SetXTitle("E_{gen}");
403 fGPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
404 // leading particle energy
405 fPTriggerEneH = new TProfile("PTriggerEne","Leading particle energy (trigger bias)",
406 20, 0., 200., 0., 1., "S");
407 fPTriggerEneH->SetXTitle("E_{L}/0.2");
408 fPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
412 void AliJetAnalysis::SetProperties(TH1* h,const char* x, const char* y) const
414 //Set properties of histos (x and y title and error propagation)
420 ////////////////////////////////////////////////////////////////////////
421 // compute weights for dE/dr histogram
423 void AliJetAnalysis::SetdEdrWeight()
425 // Due to the limited acceptance, each bin in the dE/dr has a different
426 // weight given by the ratio of the area of a disk dR to the area
427 // within the eta acceptance. The weight depends on the eta position
428 // of the jet. Here a look up table for the weights is computed. It
429 // assumes |etaJet|<0.5 and |eta_lego|<0.9. It makes bins of 0.05
430 // units in eta. Note that this is fixed by the bin width chosen for
431 // the histogram --> this weights are tailored for the specific
432 // histogram definition used here!
434 // two dimensional table: first index, bin in eta of jet, second
435 // index bin of dE/dr histo
440 Float_t binSize = (xMax-xMin)/nBins;
442 Float_t da,ds,r1,r2,h1,h2,a1,a2,theta1,theta2,area1,area2;
444 for (Int_t i=0;i<(nBins/2);i++) {
445 // index of first histo bin needing a scale factor
447 for(Int_t j=0;j<nBins;j++) {
449 da = TMath::Pi()*(binSize*binSize)*(1.0+2.0*j);
451 if (j>=ji) { // compute missing area using segments of circle
456 a1=2.0*TMath::Sqrt(2.0*h1*r1-h1*h1);
457 a2=2.0*TMath::Sqrt(2.0*h2*r2-h2*h2);
458 theta1=2*TMath::ACos((r1-h1)/r1);
459 theta2=2*TMath::ACos((r2-h2)/r2);
460 area1=binSize*(r1*r1*theta1-a1*(r1-h1));
461 area2=binSize*(r2*r2*theta2-a2*(r2-h2));
462 ds = (da-(area2-area1))/da;
464 fWeightdEdr[i][j]=ds/da;
469 // get weight for dE/dr histogram
470 Float_t AliJetAnalysis::GetdEdrWeight(Float_t etaJet, Float_t r)
472 // Return the correponding weight for the dE/dr plot
476 Float_t binSize = (xMax-xMin)/nBins;
478 Float_t eta = TMath::Abs(etaJet);
479 if ((eta > 0.5) || (r > fdrdEdr)) return 0.0;
480 Int_t binJet = (Int_t) TMath::Floor(eta/binSize);
481 Int_t binR = (Int_t) TMath::Floor(r/binSize);
482 Float_t w = fWeightdEdr[binJet][binR];
487 ////////////////////////////////////////////////////////////////////////
490 void AliJetAnalysis::FillHistograms()
495 AliJetProductionDataPDC2004* runData = new AliJetProductionDataPDC2004();
500 for (Int_t iRun = fRunMin; iRun <= fRunMax; iRun++) {
503 sprintf(fn, "%s/%s.root", fDirectory, (runData->GetRunTitle(iRun)).Data());
504 jFile = new TFile(fn);
505 printf(" Analyzing run: %d %s\n", iRun,fn);
507 // Get reader header and events to be looped over
508 AliJetReaderHeader *jReaderH = (AliJetReaderHeader*)(jFile->Get(fReaderHeader));
509 if (fEventMin == -1) fEventMin = jReaderH->GetFirstEvent();
510 if (fEventMax == -1) {
511 fEventMax = jReaderH->GetLastEvent();
513 fEventMax = TMath::Min(fEventMax, jReaderH->GetLastEvent());
516 AliUA1JetHeader *jH = (AliUA1JetHeader *) (jFile->Get("AliUA1JetHeader"));
519 fWeight = runData->GetWeight(iRun)/ ( (Float_t) (fEventMax - fEventMin + 1));
522 for (Int_t i = fEventMin; i < fEventMax; i++) {
523 if (i%100 == 0) printf(" Analyzing event %d / %d \n",i,fEventMax);
525 // Get next tree with AliJet
527 sprintf(nameT, "TreeJ%d",i);
528 TTree *jetT =(TTree *)(jFile->Get(nameT));
529 if (fDoRecJ) jetT->SetBranchAddress("FoundJet", &fRecJ);
530 if (fDoGenJ) jetT->SetBranchAddress("GenJet", &fGenJ);
531 if (fDoPart) jetT->SetBranchAddress("LeadingPart", &fPart);
534 int nJets = fRecJ->GetNJets();
536 TArrayI inJet = fRecJ->GetInJet();
537 if(inJet.GetSize()>fminMult){ // removes events with low multiplicity
538 if (fDoKine) FillKineH();
539 if (fDoCorr) FillCorrH();
540 if (fDoCorr50) FillCorr50H();
541 if (fDoShap) FillShapH(jH->GetRadius());
542 if (fDoFrag) FillFragH();
543 if (fDoTrig) FillTrigH();
544 if (fDoJt) FillJtH();
545 if (fDodNdxi) FilldNdxiH();
547 delete jetT; // jet should be deleted before creating a new one
548 if(inJet.GetSize()>fminMult){ // removes events with low multiplicity
549 if (fDoBkgd && nJets>0) FillBkgd(i,iRun);
551 } // end loop over events in one file
552 if (jFile) jFile->Close();
554 } // end loop over files
557 void AliJetAnalysis::FillKineH()
560 if (fDoPart && fPart->LeadingFound()){
561 fPKineEneH->Fill(fPart->GetE(),fWeight);
562 fPKinePtH->Fill(fPart->GetPt(),fWeight);
563 fPKinePhiH->Fill(fPart->GetPhi(),fWeight);
564 fPKineEtaH->Fill(fPart->GetEta(),fWeight);
566 // leading generated jet
567 if (fDoGenJ && fGenJ->GetNJets()> 0){
568 fGKineEneH->Fill(fGenJ->GetE(0),fWeight);
569 fGKinePtH->Fill(fGenJ->GetPt(0),fWeight);
570 fGKinePhiH->Fill(fGenJ->GetPhi(0),fWeight);
571 fGKineEtaH->Fill(fGenJ->GetEta(0),fWeight);
573 // leading reconstructed jet
574 if (fDoRecJ && fRecJ->GetNJets()> 0) {
575 TArrayF p=fRecJ->GetPtFromSignal();
576 if (p[0]>fPercentage) {
577 fRKineEneH->Fill(fRecJ->GetE(0)/fEfactor,fWeight);
578 fRKinePtH->Fill(fRecJ->GetPt(0),fWeight);
579 fRKinePhiH->Fill(fRecJ->GetPhi(0),fWeight);
580 fRKineEtaH->Fill(fRecJ->GetEta(0),fWeight);
585 void AliJetAnalysis::FillCorrH()
587 // Fill correlation histograms
588 if (fDoPart && fPart->LeadingFound() && fDoGenJ && fGenJ->GetNJets()> 0)
589 Correlation(fPart->GetLeading(),fGenJ->GetJet(0),
590 fPGCorrEneH,fPGCorrPtH,fPGCorrEtaH,fPGCorrPhiH);
591 if (fDoPart && fPart->LeadingFound() && fDoRecJ && fRecJ->GetNJets()> 0) {
592 TArrayF p=fRecJ->GetPtFromSignal();
593 if (p[0]>fPercentage)
594 Correlation(fPart->GetLeading(),fRecJ->GetJet(0),
595 fPRCorrEneH,fPRCorrPtH,fPRCorrEtaH,fPRCorrPhiH);
597 if (fDoGenJ && fGenJ->GetNJets()> 0 && fDoRecJ && fRecJ->GetNJets()> 0) {
598 TArrayF p=fRecJ->GetPtFromSignal();
599 if (p[0]>fPercentage)
600 Correlation(fRecJ->GetJet(0), fGenJ->GetJet(0),
601 fRGCorrEneH,fRGCorrPtH,fRGCorrEtaH,fRGCorrPhiH);
605 void AliJetAnalysis::Correlation(TLorentzVector *lv1,TLorentzVector *lv2,
606 TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
608 // Correlation histograms
609 h1->Fill(lv1->E(),lv2->E(),fWeight);
610 h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
611 h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
613 p1 = ((lv1->Phi() < 0) ? (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
614 p2 = ((lv2->Phi() < 0) ? (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
615 h4->Fill(p1,p2,fWeight);
618 void AliJetAnalysis::FillCorr50H()
620 // Fill correlation histograms when one particle has > 50% of jet energy
621 if (fDoRecJ && fRecJ->GetNJets()> 0) {
622 TArrayF p = fRecJ->GetPtFromSignal();
623 if (p[0]>fPercentage) {
624 if (fDoPart && fPart->LeadingFound())
625 Correlation50(fRecJ, fPart->GetLeading(),fRecJ->GetJet(0),
626 fPRCorr50EneH,fPRCorr50PtH,fPRCorr50EtaH,fPRCorr50PhiH);
627 if (fDoGenJ && fGenJ->GetNJets()> 0)
628 Correlation50(fRecJ, fRecJ->GetJet(0), fGenJ->GetJet(0),
629 fRGCorr50EneH,fRGCorr50PtH,fRGCorr50EtaH,fRGCorr50PhiH);
634 void AliJetAnalysis::Correlation50(AliJet *j,TLorentzVector *lv1,TLorentzVector *lv2,
635 TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
637 // Correlation histograms when one particle has > 50% of jet energy
638 TArrayF ptin = j->GetPtIn();
639 TArrayF etain = j->GetEtaIn();
640 TArrayI inJet = j->GetInJet();
643 for(Int_t i=0;i<(inJet.GetSize());i++) {
645 Float_t t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
646 Float_t x = (ptin[i]*TMath::Sqrt(1.+1./(t1*t1)))/(j->GetE(0));
650 if (flag>1) cout << " Flag = " << flag << endl;
652 h1->Fill(lv1->E(),lv2->E(),fWeight);
653 h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
654 h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
656 p1 = ((lv1->Phi() < 0) ?
657 (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
658 p2 = ((lv2->Phi() < 0) ?
659 (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
660 h4->Fill(p1,p2,fWeight);
664 void AliJetAnalysis::FillJtH()
666 // Fill J_T histogram
667 if (fRecJ->GetNJets()> 0) {
668 fjv3X = 0.0; fjv3Y = 0.0; fjv3Z = 0.0;
669 TArrayF p=fRecJ->GetPtFromSignal();
670 if (p[0]>fPercentage) {
672 const TVector3 kJv3 = fRecJ->GetJet(0)->Vect();
673 fjv3X = kJv3.X(); fjv3Y = kJv3.Y(); fjv3Z = kJv3.Z();
675 TArrayI inJet = fRecJ->GetInJet();
676 TArrayF etain = fRecJ->GetEtaIn();
677 TArrayF ptin = fRecJ->GetPtIn();
678 TArrayF phiin = fRecJ->GetPhiIn();
679 Float_t deta, dphi,jt, dr;
680 for(Int_t i=0;i<inJet.GetSize();i++) {
681 deta = etain[i] - fRecJ->GetEta(0);
682 dphi = phiin[i] - fRecJ->GetPhi(0);
683 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
684 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
685 dr = TMath::Sqrt(deta * deta + dphi * dphi);
686 if ((dr<fdrJt) && (ptin[i] > fPartPtCut)) {
687 trk.SetPtEtaPhi(ptin[i],etain[i],phiin[i]);
688 jt = TMath::Sin(trk.Angle(kJv3))*trk.Mag();
689 fJtH->Fill(jt,fRecJ->GetPt(0),fWeight);
692 fJtW->Fill(fRecJ->GetPt(0),fWeight);
698 void AliJetAnalysis::FilldNdxiH()
700 // Fill dN/d#xi histograms
701 if (fRecJ->GetNJets()> 0) {
702 TArrayF p=fRecJ->GetPtFromSignal();
703 if (p[0]>fPercentage) {
705 TArrayI inJet = fRecJ->GetInJet();
706 TArrayF etain = fRecJ->GetEtaIn();
707 TArrayF ptin = fRecJ->GetPtIn();
708 TArrayF phiin = fRecJ->GetPhiIn();
709 Float_t xi,t1,ene,dr,deta,dphi;
710 for(Int_t i=0;i<inJet.GetSize();i++) {
711 deta = etain[i] - fRecJ->GetEta(0);
712 dphi = phiin[i] - fRecJ->GetPhi(0);
713 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
714 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
715 dr = TMath::Sqrt(deta * deta + dphi * dphi);
716 if ((dr<fdrdNdxi) && (ptin[i] > fPartPtCut)) {
717 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
718 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
719 xi = (Float_t) TMath::Log((fRecJ->GetE(0)/fEfactor)/ene);
720 fdNdxiH->Fill(xi,fRecJ->GetPt(0),fWeight);
723 fdNdxiW->Fill(fRecJ->GetPt(0),fWeight);
724 fPtEneH->Fill(fRecJ->GetE(0)/fEfactor,fRecJ->GetPt(0),fWeight);
730 void AliJetAnalysis::FillBkgd(Int_t eventN, Int_t runN)
732 // Background calculated with hijing events (no pythia)
733 if (fp0>fPercentage) {
734 fPtJ=0.,fEJ=0.,fEtaJ=0.,fPhiJ=0.;
735 // Background calculated with hijing events (no pythia)
736 AliJetProductionDataPDC2004* runDataB = new AliJetProductionDataPDC2004();
740 sprintf(fnB, "%s/%s.root", fBkgdDirectory, (runDataB->GetRunTitle(runN)).Data());
741 jFileB = new TFile(fnB);
744 sprintf(nameB, "TreeJ%d",eventN);
745 TTree *jetB =(TTree *)(jFileB->Get(nameB));
746 jetB->SetBranchAddress("FoundJet", &fRecB);
749 TArrayI inJetB = fRecB->GetInJet();
750 TArrayF etainB = fRecB->GetEtaIn();
751 TArrayF ptinB = fRecB->GetPtIn();
752 TArrayF phiinB = fRecB->GetPhiIn();
753 fPtJ = fRecJ->GetPt(0);
754 fEJ = fRecJ->GetE(0);
755 fEtaJ = fRecJ->GetEta(0);
756 fPhiJ = fRecJ->GetPhi(0);
757 Float_t t1,ene,xi,detaB,dphiB,drB,jt,wB;
760 jv3B.SetX(fjv3X); jv3B.SetY(fjv3Y); jv3B.SetZ(fjv3Z);
762 for(Int_t k=0;k<inJetB.GetSize();k++){
763 if(ptinB[k] > fPartPtCut){
764 detaB = etainB[k] - fEtaJ;
765 dphiB = phiinB[k] - fPhiJ;
766 if (dphiB < -TMath::Pi()) dphiB= -dphiB - 2.0 * TMath::Pi();
767 if (dphiB > TMath::Pi()) dphiB = 2.0 * TMath::Pi() - dphiB;
768 drB = TMath::Sqrt(detaB * detaB + dphiB * dphiB);
769 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etainB[k])));
770 ene = ptinB[k]*TMath::Sqrt(1.+1./(t1*t1));
771 trkB.SetPtEtaPhi(ptinB[k],etainB[k],phiinB[k]);
774 xi = (Float_t) TMath::Log(fEJ/ene);
775 fdNdxiB->Fill(xi,fPtJ,fWeight);
779 jt = TMath::Sin(trkB.Angle(jv3B))*(trkB.Mag());
780 fJtB->Fill(jt,fPtJ,fWeight);
784 wB = GetdEdrWeight(fEtaJ,drB)*fWeight*ene;
785 fdEdrB->Fill(drB,fPtJ,wB);
790 if (jFileB) jFileB->Close();
795 void AliJetAnalysis::FillShapH(Float_t r)
797 // Fill jet shape histograms
798 if (fDoRecJ && fRecJ->GetNJets()> 0) {
799 TArrayF p=fRecJ->GetPtFromSignal();
800 if (p[0]>fPercentage) {
801 Shape(fRecJ,fRShapSelH,fRShapRejH,fRShapAllH,fdEdrH,fPtEneH2,fdEdrW,r);
808 void AliJetAnalysis::Shape(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha,
809 TH2F *hdedr, TH2F *hptene, TH1F* wdedr, Float_t radius)
812 TArrayI inJet = j->GetInJet();
813 TArrayF etain = j->GetEtaIn();
814 TArrayF ptin = j->GetPtIn();
815 TArrayF phiin = j->GetPhiIn();
817 // first compute dE/dr histo
818 Float_t etaj = j->GetEta(0);
819 Float_t ene,w,deta,dphi,dr,t1;
820 for(Int_t i=0;i<inJet.GetSize();i++) {
821 deta = etain[i] - j->GetEta(0);
822 dphi = phiin[i] - j->GetPhi(0);
823 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
824 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
825 dr = TMath::Sqrt(deta * deta + dphi * dphi);
826 if ((dr<fdrdEdr) && (ptin[i] > fPartPtCut)) {
827 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
828 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
829 w = GetdEdrWeight(etaj,dr)*fWeight*ene;
830 hdedr->Fill(dr,j->GetPt(0),w);
833 hptene->Fill(fRecJ->GetE(0),fRecJ->GetPt(0),fWeight);
834 wdedr->Fill(j->GetPt(0),fWeight);
836 // now compute shape histos
837 Int_t nBins = ha->GetNbinsX();
838 Float_t xMin = ha->GetXaxis()->GetXmin();
839 Float_t xMax = ha->GetXaxis()->GetXmax();
840 Float_t binSize,halfBin;
841 binSize = (xMax-xMin)/nBins;
843 Float_t rptS[20], rptR[20], rptA[20];
844 for(Int_t i=0;i<nBins;i++) rptS[i]=rptR[i]=rptA[i]=0.0;
845 // fill bins in r for leading jet
846 for(Int_t i=0;i<inJet.GetSize();i++) {
847 if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
848 deta = etain[i] - j->GetEta(0);
849 dphi = phiin[i] - j->GetPhi(0);
850 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
851 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
852 dr = TMath::Sqrt(deta * deta + dphi * dphi);
853 Float_t rR = dr/radius;
854 Int_t bin = (Int_t) TMath::Floor(rR/binSize);
855 rptA[bin]+=ptin[i]/(j->GetPt(0));
856 if (inJet[i] == 1) rptS[bin]+=ptin[i]/(j->GetPt(0));
857 if (fPythia && inJet[i] == -1)
858 rptR[bin]+=ptin[i]/(j->GetPt(0));
862 // compute shape and fill histogram
863 Float_t ptS,ptR,ptA,r;
865 for (Int_t i=0;i<nBins;i++) {
869 r=(i+1)*binSize-halfBin;
870 hs->Fill(r,ptS*fWeight);
872 hr->Fill(r,ptR*fWeight);
873 ha->Fill(r,ptA*fWeight);
878 void AliJetAnalysis::FillFragH()
880 // Fill fragmentation histogram
881 if (fDoRecJ && fRecJ->GetNJets()> 0) {
882 TArrayF p=fRecJ->GetPtFromSignal();
883 if (p[0]>fPercentage) {
884 FragFun(fRecJ,fRFragSelH,fRFragRejH,fRFragAllH);
890 void AliJetAnalysis::FragFun(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha)
892 // Calculate de fragmentation function
893 TArrayI inJet = j->GetInJet();
894 TArrayF etain = j->GetEtaIn();
895 TArrayF ptin = j->GetPtIn();
899 for(Int_t i=0;i<(inJet.GetSize());i++) {
900 if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
901 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
902 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
903 xi = (Float_t) TMath::Log((j->GetE(0))/ene);
904 if (fPythia) ha->Fill(xi,fWeight);
905 if (inJet[i] == 1) hs->Fill(xi,fWeight);
906 if (fPythia && inJet[i] == -1) hr->Fill(xi,fWeight);
911 void AliJetAnalysis:: FillTrigH()
913 // Fill trigger bias histograms
914 if(fDoTrig && fGenJ->GetNJets()>0 && fRecJ->GetNJets()>0){
915 TArrayF p=fRecJ->GetPtFromSignal();
916 if (p[0]>fPercentage) {
917 float genEne = fGenJ->GetE(0);
918 float recEne = fRecJ->GetE(0);
919 float eneRatio = (recEne/fEfactor)/genEne;
921 fGTriggerEneH->Fill(genEne, eneRatio, fWeight);
922 fRTriggerEneH->Fill(recEne/fEfactor, eneRatio, fWeight);
924 if (fPart->LeadingFound()){
925 float leaEne = fPart->GetE();
926 float eneRatio2 = leaEne/genEne;
927 fGPTriggerEneH->Fill(genEne, eneRatio2, fWeight);
928 fPTriggerEneH->Fill(leaEne/0.2, eneRatio2, fWeight);
935 ////////////////////////////////////////////////////////////////////////
936 // Normalize histogrames
938 void AliJetAnalysis::NormHistograms()
940 // normalize shape histograms
942 if (fDoRecJ && fWShapR>0.) { // leading reconstructed jet
943 fRShapSelH->Scale(1.0/fWShapR);
944 fRShapRejH->Scale(1.0/fWShapR);
945 fRShapAllH->Scale(1.0/fWShapR);
949 // normalize fragmentation function histograms
951 if (fDoRecJ && fWFragR>0.) { // leading reconstructed jet
952 fRFragSelH->Scale(2.0/fWFragR);
953 fRFragRejH->Scale(2.0/fWFragR);
954 fRFragAllH->Scale(2.0/fWFragR);
959 ////////////////////////////////////////////////////////////////////////
961 void AliJetAnalysis::PlotHistograms()
963 // Plot histogramas (to be done...)
964 if (fDoKine) PlotKineH();
965 if (fDoCorr) PlotCorrH();
966 if (fDoCorr50) PlotCorr50H();
967 if (fDoShap) PlotShapH();
968 if (fDoFrag) PlotFragH();
969 if (fDoTrig) PlotTrigH();
972 void AliJetAnalysis::PlotKineH() const
980 void AliJetAnalysis::PlotCorrH() const
983 if (fDoPart && fDoGenJ) ;
984 if (fDoPart && fDoRecJ) ;
985 if (fDoGenJ && fDoRecJ) ;
987 void AliJetAnalysis::PlotCorr50H() const
990 if (fDoPart && fDoGenJ) ;
991 if (fDoPart && fDoRecJ) ;
992 if (fDoGenJ && fDoRecJ) ;
995 void AliJetAnalysis::PlotShapH() const
1002 void AliJetAnalysis::PlotFragH() const
1009 void AliJetAnalysis::PlotTrigH()
1015 ////////////////////////////////////////////////////////////////////////
1017 void AliJetAnalysis::SaveHistograms()
1020 TFile *fOut = new TFile(fFile,"recreate");
1022 if (fDoKine) SaveKineH();
1023 if (fDoCorr) SaveCorrH();
1024 if (fDoCorr50) SaveCorr50H();
1025 if (fDoShap) SaveShapH();
1026 if (fDoFrag) SaveFragH();
1027 if (fDoTrig) SaveTrigH();
1028 if (fDoJt) SaveJtH();
1029 if (fDodNdxi) SavedNdxiH();
1033 void AliJetAnalysis::SaveKineH()
1035 // Save kinematic histograms
1037 fPKineEneH->Write();
1039 fPKinePhiH->Write();
1040 fPKineEtaH->Write();
1044 fGKineEneH->Write();
1046 fGKinePhiH->Write();
1047 fGKineEtaH->Write();
1051 fRKineEneH->Write();
1053 fRKinePhiH->Write();
1054 fRKineEtaH->Write();
1058 void AliJetAnalysis::SaveCorrH()
1060 // Save correlation histograms
1061 if (fDoPart && fDoGenJ) {
1062 fPGCorrEneH->Write();
1063 fPGCorrPtH->Write();
1064 fPGCorrEtaH->Write();
1065 fPGCorrPhiH->Write();
1068 if (fDoPart && fDoRecJ) {
1069 fPRCorrEneH->Write();
1070 fPRCorrPtH->Write();
1071 fPRCorrEtaH->Write();
1072 fPRCorrPhiH->Write();
1075 if (fDoGenJ && fDoRecJ) {
1076 fRGCorrEneH->Write();
1077 fRGCorrPtH->Write();
1078 fRGCorrEtaH->Write();
1079 fRGCorrPhiH->Write();
1083 void AliJetAnalysis::SaveCorr50H()
1085 // Save correlation histograms (special case)
1086 if (fDoPart && fDoRecJ) {
1087 fPRCorr50EneH->Write();
1088 fPRCorr50PtH->Write();
1089 fPRCorr50EtaH->Write();
1090 fPRCorr50PhiH->Write();
1092 if (fDoGenJ && fDoRecJ) {
1093 fRGCorr50EneH->Write();
1094 fRGCorr50PtH->Write();
1095 fRGCorr50EtaH->Write();
1096 fRGCorr50PhiH->Write();
1100 void AliJetAnalysis::SaveShapH()
1102 // Save jet shape histograms
1104 fRShapSelH->Write();
1106 if(fDoBkgd) fdEdrB->Write();
1110 fRShapRejH->Write();
1111 fRShapAllH->Write();
1116 void AliJetAnalysis::SaveJtH()
1118 // Save J_T histograms
1121 if(fDoBkgd) fJtB->Write();
1126 void AliJetAnalysis::SavedNdxiH()
1128 // Save dN/d#xi histograms
1131 if(fDoBkgd) fdNdxiB->Write();
1137 void AliJetAnalysis::SaveFragH()
1139 // Save fragmentation histograms
1141 fRFragSelH->Write();
1143 fRFragRejH->Write();
1144 fRFragAllH->Write();
1149 void AliJetAnalysis::SaveTrigH()
1151 // Save trigger bias histograms
1153 fGTriggerEneH->Write();
1154 fRTriggerEneH->Write();
1155 fGPTriggerEneH->Write();
1156 fPTriggerEneH->Write();
1160 ////////////////////////////////////////////////////////////////////////
1161 // main Analysis function
1163 void AliJetAnalysis::Analyze()
1168 // leading generated jet
1169 // leading reconstructed jet
1171 // Correlations amd resolutions
1172 // a) correlations in energy, pt, phi, eta
1173 // b) resolutions in energy, pt, phi, eta, r
1174 // leading particle and leading generated jet
1175 // leading particle and leading reconstructed jet
1176 // leading generated jet and leading reconstructed jet
1178 // Fragmentation functions and Shapes
1179 // a) integrated over all pt
1181 // b.1) only for user selected particles in jet
1182 // b.2) only for user rejected particles in jet
1183 // b.3) for all particles in jet