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():
49 fFile("anaJets.root"),
146 // Default constructor
148 // initialize weight for dE/dr histo
152 AliJetAnalysis::~AliJetAnalysis()
158 ////////////////////////////////////////////////////////////////////////
159 // define histogrames
161 void AliJetAnalysis::DefineHistograms()
163 // Define the histograms to be filled
164 if (fDoKine) DefineKineH();
165 if (fDoCorr) DefineCorrH();
166 if (fDoCorr50) DefineCorr50H();
167 if (fDoShap) DefineShapH();
168 if (fDoFrag) DefineFragH();
169 if (fDoTrig) DefineTrigH();
170 if (fDoJt) DefineJtH();
171 if (fDodNdxi) DefinedNdxiH();
174 void AliJetAnalysis::DefineKineH()
178 fPKineEneH = new TH1F("PKineEne","Energy of leading particle",50,0,200);
179 SetProperties(fPKineEneH,"Energy (GeV)","Entries");
180 fPKinePtH = new TH1F("PKinePt","Pt of leading particle",50,0,200);
181 SetProperties(fPKinePtH,"P_{T} (GeV)","Entries");
182 fPKinePhiH = new TH1F("PKinePhiH","Azimuthal angle of leading particle",
183 90,0.,2.0*TMath::Pi());
184 SetProperties(fPKinePhiH,"#phi","Entries");
185 fPKineEtaH = new TH1F("PKineEtaH","Pseudorapidity of leading particle",
187 SetProperties(fPKineEtaH,"#eta","Entries");
189 // leading generated jet
191 fGKineEneH = new TH1F("GKineEne","Energy of generated jet",50,0,200);
192 SetProperties(fGKineEneH,"Energy (GeV)","Entries");
193 fGKinePtH = new TH1F("GKinePt","Pt of generated jet",50,0,200);
194 SetProperties(fGKinePtH,"P_{T} (GeV)","Entries");
195 fGKinePhiH = new TH1F("GKinePhiH","Azimuthal angle of generated jet",
196 90,0.,2.0*TMath::Pi());
197 SetProperties(fGKinePhiH,"#phi","Entries");
198 fGKineEtaH = new TH1F("GKineEtaH","Pseudorapidity of generated jet",
200 SetProperties(fGKineEtaH,"#eta","Entries");
202 // leading reconstructed jet
204 fRKineEneH = new TH1F("RKineEne","Energy of reconstructed jet",50,0,200);
205 SetProperties(fRKineEneH,"Energy (GeV)","Entries");
206 fRKinePtH = new TH1F("RKinePt","Pt of reconstructed jet",50,0,200);
207 SetProperties(fRKinePtH,"P_{T} (GeV)","Entries");
208 fRKinePhiH = new TH1F("RKinePhiH","Azimuthal angle of reconstructed jet",
209 90,0.,2.0*TMath::Pi());
210 SetProperties(fRKinePhiH,"#phi","Entries");
211 fRKineEtaH = new TH1F("RKineEtaH","Pseudorapidity of reconstructed jet",
213 SetProperties(fRKineEtaH,"#eta","Entries");
217 void AliJetAnalysis::DefineCorrH()
220 if (fDoPart && fDoGenJ) {
221 fPGCorrEneH = new TH2F("PGCorrEne","Energy correlation part-gen jet",
223 SetProperties(fPGCorrEneH,"Part Energy (GeV)","Gen Jet Energy (GeV)");
224 fPGCorrPtH = new TH2F("PGCorrPt","Pt correlation part-gen jet",
226 SetProperties(fPGCorrPtH,"Part P_{T} (GeV)","Gen Jet P_{T} (GeV)");
227 fPGCorrEtaH = new TH2F("PGCorrEta","Pseudorapidity correlation part-gen jet",
228 40,-1.0,1.0,40,-1.0,1.0);
229 SetProperties(fPGCorrEtaH,"Part #eta","Gen Jet #eta");
230 fPGCorrPhiH = new TH2F("PGCorrPhi","Azimuthal angle correlation part-gen jet",
231 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
232 SetProperties(fPGCorrPhiH,"Part #phi","Gen Jet #phi");
234 if (fDoPart && fDoRecJ) {
235 fPRCorrEneH = new TH2F("PRCorrEne","Energy correlation part-rec jet",
237 SetProperties(fPRCorrEneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
238 fPRCorrPtH = new TH2F("PRCorrPt","Pt correlation part-rec jet",
240 SetProperties(fPRCorrPtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
241 fPRCorrEtaH = new TH2F("PRCorrEta","Pseudorapidity correlation part-rec jet",
242 40,-1.0,1.0,40,-1.0,1.0);
243 SetProperties(fPRCorrEtaH,"part #eta","Rec Jet #eta");
244 fPRCorrPhiH = new TH2F("PRCorrPhi","Azimuthal angle correlation part-rec jet",
245 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
246 SetProperties(fPRCorrPhiH,"Part #phi","Rec Jet #phi");
248 if (fDoGenJ && fDoRecJ) {
249 fRGCorrEneH = new TH2F("RGCorrEne","Energy correlation rec jet-gen jet",
251 SetProperties(fRGCorrEneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
252 fRGCorrPtH = new TH2F("RGCorrPt","Pt correlation rec jet-gen jet",
254 SetProperties(fRGCorrPtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
255 fRGCorrEtaH = new TH2F("RGCorrEta","Pseudorapidity correlation rec jet-gen jet",
256 40,-1.0,1.0,40,-1.0,1.0);
257 SetProperties(fRGCorrEtaH,"Rec Jet #eta","Gen Jet #eta");
258 fRGCorrPhiH = new TH2F("RGCorrPhi","Azimuthal angle correlation rec jet-gen jet",
259 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
260 SetProperties(fRGCorrPhiH,"Rec Jet #phi","Gen Jet #phi");
264 void AliJetAnalysis::DefineCorr50H()
267 if (fDoPart && fDoRecJ) {
268 fPRCorr50EneH = new TH2F("PRCorr50Ene","Energy correlation part-rec jet",
270 SetProperties(fPRCorr50EneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
271 fPRCorr50PtH = new TH2F("PRCorr50Pt","Pt correlation part-rec jet",
273 SetProperties(fPRCorr50PtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
274 fPRCorr50EtaH = new TH2F("PRCorr50Eta","Pseudorapidity correlation part-rec jet",
275 40,-1.0,1.0,40,-1.0,1.0);
276 SetProperties(fPRCorr50EtaH,"part #eta","Rec Jet #eta");
277 fPRCorr50PhiH = new TH2F("PRCorr50Phi","Azimuthal angle correlation part-rec jet",
278 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
279 SetProperties(fPRCorr50PhiH,"Part #phi","Rec Jet #phi");
282 if (fDoGenJ && fDoRecJ) {
283 fRGCorr50EneH = new TH2F("RGCorr50Ene","Energy correlation rec jet-gen jet",
285 SetProperties(fRGCorr50EneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
286 fRGCorr50PtH = new TH2F("RGCorr50Pt","Pt correlation rec jet-gen jet",
288 SetProperties(fRGCorr50PtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
289 fRGCorr50EtaH = new TH2F("RGCorr50Eta","Pseudorapidity correlation rec jet-gen jet",
290 40,-1.0,1.0,40,-1.0,1.0);
291 SetProperties(fRGCorr50EtaH,"Rec Jet #eta","Gen Jet #eta");
292 fRGCorr50PhiH = new TH2F("RGCorr50Phi","Azimuthal angle correlation rec jet-gen jet",
293 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
294 SetProperties(fRGCorr50PhiH,"Rec Jet #phi","Gen Jet #phi");
298 void AliJetAnalysis::DefineShapH()
300 // leading reconstructed jet
302 fdEdrH = new TH2F("fdEdrH","dE/dr histo",20,0,1,40,0,200);
303 SetProperties(fdEdrH,"r","Rec Jet P_{T}");
304 fdEdrB = new TH2F("fdEdrB","dE/dr bkgdhisto",20,0,1,40,0,200);
305 SetProperties(fdEdrB,"r","Rec P_{T}");
306 fPtEneH2 = new TH2F("fPtEneH2","fPtEneH2",40,0,200,40,0,200);
307 SetProperties(fPtEneH2,"Rec Jet E","Rec Jet P_{T}");
308 fdEdrW = new TH1F("fdEdrW","weights for dE/dr",40,0,200);
309 SetProperties(fdEdrW,"Rec Jet P_{T}","weights");
311 fRShapSelH = new TH1F("RShapSel","Shape of generated jets (sel part)",20,0.,1.);
312 SetProperties(fRShapSelH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
313 fRShapRejH = new TH1F("RShapRej","Shape of generated jets (rej part)",20,0.,1.);
314 SetProperties(fRShapRejH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
315 fRShapAllH = new TH1F("RShapAll","Shape of generated jets (all part)",20,0.,1.);
316 SetProperties(fRShapAllH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
320 void AliJetAnalysis::DefineJtH()
322 // Define the histogram for J_T
324 fJtH = new TH2F("fJtH","J_{T} histogram",80,0.,4.,40,0.,200.);
325 SetProperties(fJtH,"J_{T}","Rec Jet P_{T}");
326 fJtB = new TH2F("fJtB","J_{T} bkgd histogram",80,0.,4.,40,0.,200.);
327 SetProperties(fJtB,"J_{T}","Rec P_{T}");
328 fJtW = new TH1F("fJtW","J_{T} weight",40,0,200);
329 SetProperties(fJtW,"J_{T}W","weight");
333 void AliJetAnalysis::DefinedNdxiH()
335 // Define the histogram for dN/dxi
337 fdNdxiH = new TH2F("fdNdxiH","dN/d#xi histo",200,0,10,40,0,200);
338 SetProperties(fdNdxiH,"xi","Rec Jet P_{T}");
339 fdNdxiB = new TH2F("fdNdxiB","dN/d#xi bkgd histo",200,0,10,40,0,200);
340 SetProperties(fdNdxiB,"xi","Rec P_{T}");
341 fdNdxiW = new TH1F("fdNdxiW","dN/d#xi histo",40,0,200);
342 SetProperties(fdNdxiW,"Rec Jet P_{T}","weights");
343 fPtEneH = new TH2F("fPtEneH","fPtEneH",40,0,200,40,0,200);
344 SetProperties(fPtEneH,"Rec Jet E","Rec Jet P_{T}");
348 void AliJetAnalysis::DefineFragH()
350 // leading reconstructed jet
352 fRFragSelH = new TH1F("RFragSel","Frag Fun of reconstructed jets (sel part)",20,0.,10.);
353 SetProperties(fRFragSelH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
354 fRFragRejH = new TH1F("RFragRej","Frag Fun of reconstructed jets (rej part)",20,0.,10.);
355 SetProperties(fRFragRejH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
356 fRFragAllH = new TH1F("RFragAll","Frag Fun of reconstructed jets (all part)",20,0.,10.);
357 SetProperties(fRFragAllH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
361 void AliJetAnalysis::DefineTrigH()
364 fGTriggerEneH = new TProfile("GTriggerEne","Generated energy (trigger bias)",
365 20, 0., 200., 0., 1., "S");
366 fGTriggerEneH->SetXTitle("E_{gen}");
367 fGTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
368 // reconstructed energy
369 fRTriggerEneH = new TProfile("RTriggerEne","Reconstructed energy (trigger bias)",
370 20, 0., 200., 0., 1., "S");
371 fRTriggerEneH->SetXTitle("E_{rec}");
372 fRTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
373 // generated energy vs generated/leading
374 fGPTriggerEneH = new TProfile("GPTriggerEne","Generated energy (trigger bias)",
375 20, 0., 200., 0., 1., "S");
376 fGPTriggerEneH->SetXTitle("E_{gen}");
377 fGPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
378 // leading particle energy
379 fPTriggerEneH = new TProfile("PTriggerEne","Leading particle energy (trigger bias)",
380 20, 0., 200., 0., 1., "S");
381 fPTriggerEneH->SetXTitle("E_{L}/0.2");
382 fPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
386 void AliJetAnalysis::SetProperties(TH1* h,const char* x, const char* y) const
388 //Set properties of histos (x and y title and error propagation)
394 ////////////////////////////////////////////////////////////////////////
395 // compute weights for dE/dr histogram
397 void AliJetAnalysis::SetdEdrWeight()
399 // Due to the limited acceptance, each bin in the dE/dr has a different
400 // weight given by the ratio of the area of a disk dR to the area
401 // within the eta acceptance. The weight depends on the eta position
402 // of the jet. Here a look up table for the weights is computed. It
403 // assumes |etaJet|<0.5 and |eta_lego|<0.9. It makes bins of 0.05
404 // units in eta. Note that this is fixed by the bin width chosen for
405 // the histogram --> this weights are tailored for the specific
406 // histogram definition used here!
408 // two dimensional table: first index, bin in eta of jet, second
409 // index bin of dE/dr histo
414 Float_t binSize = (xMax-xMin)/nBins;
416 Float_t da,ds,r1,r2,h1,h2,a1,a2,theta1,theta2,area1,area2;
418 for (Int_t i=0;i<(nBins/2);i++) {
419 // index of first histo bin needing a scale factor
421 for(Int_t j=0;j<nBins;j++) {
423 da = TMath::Pi()*(binSize*binSize)*(1.0+2.0*j);
425 if (j>=ji) { // compute missing area using segments of circle
430 a1=2.0*TMath::Sqrt(2.0*h1*r1-h1*h1);
431 a2=2.0*TMath::Sqrt(2.0*h2*r2-h2*h2);
432 theta1=2*TMath::ACos((r1-h1)/r1);
433 theta2=2*TMath::ACos((r2-h2)/r2);
434 area1=binSize*(r1*r1*theta1-a1*(r1-h1));
435 area2=binSize*(r2*r2*theta2-a2*(r2-h2));
436 ds = (da-(area2-area1))/da;
438 fWeightdEdr[i][j]=ds/da;
443 // get weight for dE/dr histogram
444 Float_t AliJetAnalysis::GetdEdrWeight(Float_t etaJet, Float_t r)
446 // Return the correponding weight for the dE/dr plot
450 Float_t binSize = (xMax-xMin)/nBins;
452 Float_t eta = TMath::Abs(etaJet);
453 if ((eta > 0.5) || (r > fdrdEdr)) return 0.0;
454 Int_t binJet = (Int_t) TMath::Floor(eta/binSize);
455 Int_t binR = (Int_t) TMath::Floor(r/binSize);
456 Float_t w = fWeightdEdr[binJet][binR];
461 ////////////////////////////////////////////////////////////////////////
464 void AliJetAnalysis::FillHistograms()
469 AliJetProductionDataPDC2004* runData = new AliJetProductionDataPDC2004();
474 for (Int_t iRun = fRunMin; iRun <= fRunMax; iRun++) {
477 sprintf(fn, "%s/%s.root", fDirectory, (runData->GetRunTitle(iRun)).Data());
478 jFile = new TFile(fn);
479 printf(" Analyzing run: %d %s\n", iRun,fn);
481 // Get reader header and events to be looped over
482 AliJetReaderHeader *jReaderH = (AliJetReaderHeader*)(jFile->Get(fReaderHeader));
483 if (fEventMin == -1) fEventMin = jReaderH->GetFirstEvent();
484 if (fEventMax == -1) {
485 fEventMax = jReaderH->GetLastEvent();
487 fEventMax = TMath::Min(fEventMax, jReaderH->GetLastEvent());
490 AliUA1JetHeader *jH = (AliUA1JetHeader *) (jFile->Get("AliUA1JetHeader"));
493 fWeight = runData->GetWeight(iRun)/ ( (Float_t) (fEventMax - fEventMin + 1));
496 for (Int_t i = fEventMin; i < fEventMax; i++) {
497 if (i%100 == 0) printf(" Analyzing event %d / %d \n",i,fEventMax);
499 // Get next tree with AliJet
501 sprintf(nameT, "TreeJ%d",i);
502 TTree *jetT =(TTree *)(jFile->Get(nameT));
503 if (fDoRecJ) jetT->SetBranchAddress("FoundJet", &fRecJ);
504 if (fDoGenJ) jetT->SetBranchAddress("GenJet", &fGenJ);
505 if (fDoPart) jetT->SetBranchAddress("LeadingPart", &fPart);
508 int nJets = fRecJ->GetNJets();
510 TArrayI inJet = fRecJ->GetInJet();
511 if(inJet.GetSize()>fminMult){ // removes events with low multiplicity
512 if (fDoKine) FillKineH();
513 if (fDoCorr) FillCorrH();
514 if (fDoCorr50) FillCorr50H();
515 if (fDoShap) FillShapH(jH->GetRadius());
516 if (fDoFrag) FillFragH();
517 if (fDoTrig) FillTrigH();
518 if (fDoJt) FillJtH();
519 if (fDodNdxi) FilldNdxiH();
521 delete jetT; // jet should be deleted before creating a new one
522 if(inJet.GetSize()>fminMult){ // removes events with low multiplicity
523 if (fDoBkgd && nJets>0) FillBkgd(i,iRun);
525 } // end loop over events in one file
526 if (jFile) jFile->Close();
528 } // end loop over files
531 void AliJetAnalysis::FillKineH()
534 if (fDoPart && fPart->LeadingFound()){
535 fPKineEneH->Fill(fPart->GetE(),fWeight);
536 fPKinePtH->Fill(fPart->GetPt(),fWeight);
537 fPKinePhiH->Fill(fPart->GetPhi(),fWeight);
538 fPKineEtaH->Fill(fPart->GetEta(),fWeight);
540 // leading generated jet
541 if (fDoGenJ && fGenJ->GetNJets()> 0){
542 fGKineEneH->Fill(fGenJ->GetE(0),fWeight);
543 fGKinePtH->Fill(fGenJ->GetPt(0),fWeight);
544 fGKinePhiH->Fill(fGenJ->GetPhi(0),fWeight);
545 fGKineEtaH->Fill(fGenJ->GetEta(0),fWeight);
547 // leading reconstructed jet
548 if (fDoRecJ && fRecJ->GetNJets()> 0) {
549 TArrayF p=fRecJ->GetPtFromSignal();
550 if (p[0]>fPercentage) {
551 fRKineEneH->Fill(fRecJ->GetE(0)/fEfactor,fWeight);
552 fRKinePtH->Fill(fRecJ->GetPt(0),fWeight);
553 fRKinePhiH->Fill(fRecJ->GetPhi(0),fWeight);
554 fRKineEtaH->Fill(fRecJ->GetEta(0),fWeight);
559 void AliJetAnalysis::FillCorrH()
561 // Fill correlation histograms
562 if (fDoPart && fPart->LeadingFound() && fDoGenJ && fGenJ->GetNJets()> 0)
563 Correlation(fPart->GetLeading(),fGenJ->GetJet(0),
564 fPGCorrEneH,fPGCorrPtH,fPGCorrEtaH,fPGCorrPhiH);
565 if (fDoPart && fPart->LeadingFound() && fDoRecJ && fRecJ->GetNJets()> 0) {
566 TArrayF p=fRecJ->GetPtFromSignal();
567 if (p[0]>fPercentage)
568 Correlation(fPart->GetLeading(),fRecJ->GetJet(0),
569 fPRCorrEneH,fPRCorrPtH,fPRCorrEtaH,fPRCorrPhiH);
571 if (fDoGenJ && fGenJ->GetNJets()> 0 && fDoRecJ && fRecJ->GetNJets()> 0) {
572 TArrayF p=fRecJ->GetPtFromSignal();
573 if (p[0]>fPercentage)
574 Correlation(fRecJ->GetJet(0), fGenJ->GetJet(0),
575 fRGCorrEneH,fRGCorrPtH,fRGCorrEtaH,fRGCorrPhiH);
579 void AliJetAnalysis::Correlation(TLorentzVector *lv1,TLorentzVector *lv2,
580 TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
582 // Correlation histograms
583 h1->Fill(lv1->E(),lv2->E(),fWeight);
584 h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
585 h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
587 p1 = ((lv1->Phi() < 0) ? (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
588 p2 = ((lv2->Phi() < 0) ? (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
589 h4->Fill(p1,p2,fWeight);
592 void AliJetAnalysis::FillCorr50H()
594 // Fill correlation histograms when one particle has > 50% of jet energy
595 if (fDoRecJ && fRecJ->GetNJets()> 0) {
596 TArrayF p = fRecJ->GetPtFromSignal();
597 if (p[0]>fPercentage) {
598 if (fDoPart && fPart->LeadingFound())
599 Correlation50(fRecJ, fPart->GetLeading(),fRecJ->GetJet(0),
600 fPRCorr50EneH,fPRCorr50PtH,fPRCorr50EtaH,fPRCorr50PhiH);
601 if (fDoGenJ && fGenJ->GetNJets()> 0)
602 Correlation50(fRecJ, fRecJ->GetJet(0), fGenJ->GetJet(0),
603 fRGCorr50EneH,fRGCorr50PtH,fRGCorr50EtaH,fRGCorr50PhiH);
608 void AliJetAnalysis::Correlation50(AliJet *j,TLorentzVector *lv1,TLorentzVector *lv2,
609 TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
611 // Correlation histograms when one particle has > 50% of jet energy
612 TArrayF ptin = j->GetPtIn();
613 TArrayF etain = j->GetEtaIn();
614 TArrayI inJet = j->GetInJet();
617 for(Int_t i=0;i<(inJet.GetSize());i++) {
619 Float_t t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
620 Float_t x = (ptin[i]*TMath::Sqrt(1.+1./(t1*t1)))/(j->GetE(0));
624 if (flag>1) cout << " Flag = " << flag << endl;
626 h1->Fill(lv1->E(),lv2->E(),fWeight);
627 h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
628 h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
630 p1 = ((lv1->Phi() < 0) ?
631 (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
632 p2 = ((lv2->Phi() < 0) ?
633 (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
634 h4->Fill(p1,p2,fWeight);
638 void AliJetAnalysis::FillJtH()
640 // Fill J_T histogram
641 if (fRecJ->GetNJets()> 0) {
642 fjv3X = 0.0; fjv3Y = 0.0; fjv3Z = 0.0;
643 TArrayF p=fRecJ->GetPtFromSignal();
644 if (p[0]>fPercentage) {
646 const TVector3 kJv3 = fRecJ->GetJet(0)->Vect();
647 fjv3X = kJv3.X(); fjv3Y = kJv3.Y(); fjv3Z = kJv3.Z();
649 TArrayI inJet = fRecJ->GetInJet();
650 TArrayF etain = fRecJ->GetEtaIn();
651 TArrayF ptin = fRecJ->GetPtIn();
652 TArrayF phiin = fRecJ->GetPhiIn();
653 Float_t deta, dphi,jt, dr;
654 for(Int_t i=0;i<inJet.GetSize();i++) {
655 deta = etain[i] - fRecJ->GetEta(0);
656 dphi = phiin[i] - fRecJ->GetPhi(0);
657 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
658 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
659 dr = TMath::Sqrt(deta * deta + dphi * dphi);
660 if ((dr<fdrJt) && (ptin[i] > fPartPtCut)) {
661 trk.SetPtEtaPhi(ptin[i],etain[i],phiin[i]);
662 jt = TMath::Sin(trk.Angle(kJv3))*trk.Mag();
663 fJtH->Fill(jt,fRecJ->GetPt(0),fWeight);
666 fJtW->Fill(fRecJ->GetPt(0),fWeight);
672 void AliJetAnalysis::FilldNdxiH()
674 // Fill dN/d#xi histograms
675 if (fRecJ->GetNJets()> 0) {
676 TArrayF p=fRecJ->GetPtFromSignal();
677 if (p[0]>fPercentage) {
679 TArrayI inJet = fRecJ->GetInJet();
680 TArrayF etain = fRecJ->GetEtaIn();
681 TArrayF ptin = fRecJ->GetPtIn();
682 TArrayF phiin = fRecJ->GetPhiIn();
683 Float_t xi,t1,ene,dr,deta,dphi;
684 for(Int_t i=0;i<inJet.GetSize();i++) {
685 deta = etain[i] - fRecJ->GetEta(0);
686 dphi = phiin[i] - fRecJ->GetPhi(0);
687 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
688 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
689 dr = TMath::Sqrt(deta * deta + dphi * dphi);
690 if ((dr<fdrdNdxi) && (ptin[i] > fPartPtCut)) {
691 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
692 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
693 xi = (Float_t) TMath::Log((fRecJ->GetE(0)/fEfactor)/ene);
694 fdNdxiH->Fill(xi,fRecJ->GetPt(0),fWeight);
697 fdNdxiW->Fill(fRecJ->GetPt(0),fWeight);
698 fPtEneH->Fill(fRecJ->GetE(0)/fEfactor,fRecJ->GetPt(0),fWeight);
704 void AliJetAnalysis::FillBkgd(Int_t eventN, Int_t runN)
706 // Background calculated with hijing events (no pythia)
707 if (fp0>fPercentage) {
708 fPtJ=0.,fEJ=0.,fEtaJ=0.,fPhiJ=0.;
709 // Background calculated with hijing events (no pythia)
710 AliJetProductionDataPDC2004* runDataB = new AliJetProductionDataPDC2004();
714 sprintf(fnB, "%s/%s.root", fBkgdDirectory, (runDataB->GetRunTitle(runN)).Data());
715 jFileB = new TFile(fnB);
718 sprintf(nameB, "TreeJ%d",eventN);
719 TTree *jetB =(TTree *)(jFileB->Get(nameB));
720 jetB->SetBranchAddress("FoundJet", &fRecB);
723 TArrayI inJetB = fRecB->GetInJet();
724 TArrayF etainB = fRecB->GetEtaIn();
725 TArrayF ptinB = fRecB->GetPtIn();
726 TArrayF phiinB = fRecB->GetPhiIn();
727 fPtJ = fRecJ->GetPt(0);
728 fEJ = fRecJ->GetE(0);
729 fEtaJ = fRecJ->GetEta(0);
730 fPhiJ = fRecJ->GetPhi(0);
731 Float_t t1,ene,xi,detaB,dphiB,drB,jt,wB;
734 jv3B.SetX(fjv3X); jv3B.SetY(fjv3Y); jv3B.SetZ(fjv3Z);
736 for(Int_t k=0;k<inJetB.GetSize();k++){
737 if(ptinB[k] > fPartPtCut){
738 detaB = etainB[k] - fEtaJ;
739 dphiB = phiinB[k] - fPhiJ;
740 if (dphiB < -TMath::Pi()) dphiB= -dphiB - 2.0 * TMath::Pi();
741 if (dphiB > TMath::Pi()) dphiB = 2.0 * TMath::Pi() - dphiB;
742 drB = TMath::Sqrt(detaB * detaB + dphiB * dphiB);
743 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etainB[k])));
744 ene = ptinB[k]*TMath::Sqrt(1.+1./(t1*t1));
745 trkB.SetPtEtaPhi(ptinB[k],etainB[k],phiinB[k]);
748 xi = (Float_t) TMath::Log(fEJ/ene);
749 fdNdxiB->Fill(xi,fPtJ,fWeight);
753 jt = TMath::Sin(trkB.Angle(jv3B))*(trkB.Mag());
754 fJtB->Fill(jt,fPtJ,fWeight);
758 wB = GetdEdrWeight(fEtaJ,drB)*fWeight*ene;
759 fdEdrB->Fill(drB,fPtJ,wB);
764 if (jFileB) jFileB->Close();
769 void AliJetAnalysis::FillShapH(Float_t r)
771 // Fill jet shape histograms
772 if (fDoRecJ && fRecJ->GetNJets()> 0) {
773 TArrayF p=fRecJ->GetPtFromSignal();
774 if (p[0]>fPercentage) {
775 Shape(fRecJ,fRShapSelH,fRShapRejH,fRShapAllH,fdEdrH,fPtEneH2,fdEdrW,r);
782 void AliJetAnalysis::Shape(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha,
783 TH2F *hdedr, TH2F *hptene, TH1F* wdedr, Float_t radius)
786 TArrayI inJet = j->GetInJet();
787 TArrayF etain = j->GetEtaIn();
788 TArrayF ptin = j->GetPtIn();
789 TArrayF phiin = j->GetPhiIn();
791 // first compute dE/dr histo
792 Float_t etaj = j->GetEta(0);
793 Float_t ene,w,deta,dphi,dr,t1;
794 for(Int_t i=0;i<inJet.GetSize();i++) {
795 deta = etain[i] - j->GetEta(0);
796 dphi = phiin[i] - j->GetPhi(0);
797 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
798 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
799 dr = TMath::Sqrt(deta * deta + dphi * dphi);
800 if ((dr<fdrdEdr) && (ptin[i] > fPartPtCut)) {
801 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
802 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
803 w = GetdEdrWeight(etaj,dr)*fWeight*ene;
804 hdedr->Fill(dr,j->GetPt(0),w);
807 hptene->Fill(fRecJ->GetE(0),fRecJ->GetPt(0),fWeight);
808 wdedr->Fill(j->GetPt(0),fWeight);
810 // now compute shape histos
811 Int_t nBins = ha->GetNbinsX();
812 Float_t xMin = ha->GetXaxis()->GetXmin();
813 Float_t xMax = ha->GetXaxis()->GetXmax();
814 Float_t binSize,halfBin;
815 binSize = (xMax-xMin)/nBins;
817 Float_t rptS[20], rptR[20], rptA[20];
818 for(Int_t i=0;i<nBins;i++) rptS[i]=rptR[i]=rptA[i]=0.0;
819 // fill bins in r for leading jet
820 for(Int_t i=0;i<inJet.GetSize();i++) {
821 if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
822 deta = etain[i] - j->GetEta(0);
823 dphi = phiin[i] - j->GetPhi(0);
824 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
825 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
826 dr = TMath::Sqrt(deta * deta + dphi * dphi);
827 Float_t rR = dr/radius;
828 Int_t bin = (Int_t) TMath::Floor(rR/binSize);
829 rptA[bin]+=ptin[i]/(j->GetPt(0));
830 if (inJet[i] == 1) rptS[bin]+=ptin[i]/(j->GetPt(0));
831 if (fPythia && inJet[i] == -1)
832 rptR[bin]+=ptin[i]/(j->GetPt(0));
836 // compute shape and fill histogram
837 Float_t ptS,ptR,ptA,r;
839 for (Int_t i=0;i<nBins;i++) {
843 r=(i+1)*binSize-halfBin;
844 hs->Fill(r,ptS*fWeight);
846 hr->Fill(r,ptR*fWeight);
847 ha->Fill(r,ptA*fWeight);
852 void AliJetAnalysis::FillFragH()
854 // Fill fragmentation histogram
855 if (fDoRecJ && fRecJ->GetNJets()> 0) {
856 TArrayF p=fRecJ->GetPtFromSignal();
857 if (p[0]>fPercentage) {
858 FragFun(fRecJ,fRFragSelH,fRFragRejH,fRFragAllH);
864 void AliJetAnalysis::FragFun(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha)
866 // Calculate de fragmentation function
867 TArrayI inJet = j->GetInJet();
868 TArrayF etain = j->GetEtaIn();
869 TArrayF ptin = j->GetPtIn();
873 for(Int_t i=0;i<(inJet.GetSize());i++) {
874 if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
875 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
876 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
877 xi = (Float_t) TMath::Log((j->GetE(0))/ene);
878 if (fPythia) ha->Fill(xi,fWeight);
879 if (inJet[i] == 1) hs->Fill(xi,fWeight);
880 if (fPythia && inJet[i] == -1) hr->Fill(xi,fWeight);
885 void AliJetAnalysis:: FillTrigH()
887 // Fill trigger bias histograms
888 if(fDoTrig && fGenJ->GetNJets()>0 && fRecJ->GetNJets()>0){
889 TArrayF p=fRecJ->GetPtFromSignal();
890 if (p[0]>fPercentage) {
891 float genEne = fGenJ->GetE(0);
892 float recEne = fRecJ->GetE(0);
893 float eneRatio = (recEne/fEfactor)/genEne;
895 fGTriggerEneH->Fill(genEne, eneRatio, fWeight);
896 fRTriggerEneH->Fill(recEne/fEfactor, eneRatio, fWeight);
898 if (fPart->LeadingFound()){
899 float leaEne = fPart->GetE();
900 float eneRatio2 = leaEne/genEne;
901 fGPTriggerEneH->Fill(genEne, eneRatio2, fWeight);
902 fPTriggerEneH->Fill(leaEne/0.2, eneRatio2, fWeight);
909 ////////////////////////////////////////////////////////////////////////
910 // Normalize histogrames
912 void AliJetAnalysis::NormHistograms()
914 // normalize shape histograms
916 if (fDoRecJ && fWShapR>0.) { // leading reconstructed jet
917 fRShapSelH->Scale(1.0/fWShapR);
918 fRShapRejH->Scale(1.0/fWShapR);
919 fRShapAllH->Scale(1.0/fWShapR);
923 // normalize fragmentation function histograms
925 if (fDoRecJ && fWFragR>0.) { // leading reconstructed jet
926 fRFragSelH->Scale(2.0/fWFragR);
927 fRFragRejH->Scale(2.0/fWFragR);
928 fRFragAllH->Scale(2.0/fWFragR);
933 ////////////////////////////////////////////////////////////////////////
935 void AliJetAnalysis::PlotHistograms()
937 // Plot histogramas (to be done...)
938 if (fDoKine) PlotKineH();
939 if (fDoCorr) PlotCorrH();
940 if (fDoCorr50) PlotCorr50H();
941 if (fDoShap) PlotShapH();
942 if (fDoFrag) PlotFragH();
943 if (fDoTrig) PlotTrigH();
946 void AliJetAnalysis::PlotKineH() const
951 void AliJetAnalysis::PlotCorrH() const
955 void AliJetAnalysis::PlotCorr50H() const
960 void AliJetAnalysis::PlotShapH() const
965 void AliJetAnalysis::PlotFragH() const
970 void AliJetAnalysis::PlotTrigH()
976 ////////////////////////////////////////////////////////////////////////
978 void AliJetAnalysis::SaveHistograms()
981 TFile *fOut = new TFile(fFile,"recreate");
983 if (fDoKine) SaveKineH();
984 if (fDoCorr) SaveCorrH();
985 if (fDoCorr50) SaveCorr50H();
986 if (fDoShap) SaveShapH();
987 if (fDoFrag) SaveFragH();
988 if (fDoTrig) SaveTrigH();
989 if (fDoJt) SaveJtH();
990 if (fDodNdxi) SavedNdxiH();
994 void AliJetAnalysis::SaveKineH()
996 // Save kinematic histograms
1000 fPKinePhiH->Write();
1001 fPKineEtaH->Write();
1005 fGKineEneH->Write();
1007 fGKinePhiH->Write();
1008 fGKineEtaH->Write();
1012 fRKineEneH->Write();
1014 fRKinePhiH->Write();
1015 fRKineEtaH->Write();
1019 void AliJetAnalysis::SaveCorrH()
1021 // Save correlation histograms
1022 if (fDoPart && fDoGenJ) {
1023 fPGCorrEneH->Write();
1024 fPGCorrPtH->Write();
1025 fPGCorrEtaH->Write();
1026 fPGCorrPhiH->Write();
1029 if (fDoPart && fDoRecJ) {
1030 fPRCorrEneH->Write();
1031 fPRCorrPtH->Write();
1032 fPRCorrEtaH->Write();
1033 fPRCorrPhiH->Write();
1036 if (fDoGenJ && fDoRecJ) {
1037 fRGCorrEneH->Write();
1038 fRGCorrPtH->Write();
1039 fRGCorrEtaH->Write();
1040 fRGCorrPhiH->Write();
1044 void AliJetAnalysis::SaveCorr50H()
1046 // Save correlation histograms (special case)
1047 if (fDoPart && fDoRecJ) {
1048 fPRCorr50EneH->Write();
1049 fPRCorr50PtH->Write();
1050 fPRCorr50EtaH->Write();
1051 fPRCorr50PhiH->Write();
1053 if (fDoGenJ && fDoRecJ) {
1054 fRGCorr50EneH->Write();
1055 fRGCorr50PtH->Write();
1056 fRGCorr50EtaH->Write();
1057 fRGCorr50PhiH->Write();
1061 void AliJetAnalysis::SaveShapH()
1063 // Save jet shape histograms
1065 fRShapSelH->Write();
1067 if(fDoBkgd) fdEdrB->Write();
1071 fRShapRejH->Write();
1072 fRShapAllH->Write();
1077 void AliJetAnalysis::SaveJtH()
1079 // Save J_T histograms
1082 if(fDoBkgd) fJtB->Write();
1087 void AliJetAnalysis::SavedNdxiH()
1089 // Save dN/d#xi histograms
1092 if(fDoBkgd) fdNdxiB->Write();
1098 void AliJetAnalysis::SaveFragH()
1100 // Save fragmentation histograms
1102 fRFragSelH->Write();
1104 fRFragRejH->Write();
1105 fRFragAllH->Write();
1110 void AliJetAnalysis::SaveTrigH()
1112 // Save trigger bias histograms
1114 fGTriggerEneH->Write();
1115 fRTriggerEneH->Write();
1116 fGPTriggerEneH->Write();
1117 fPTriggerEneH->Write();
1121 ////////////////////////////////////////////////////////////////////////
1122 // main Analysis function
1124 void AliJetAnalysis::Analyze()
1129 // leading generated jet
1130 // leading reconstructed jet
1132 // Correlations amd resolutions
1133 // a) correlations in energy, pt, phi, eta
1134 // b) resolutions in energy, pt, phi, eta, r
1135 // leading particle and leading generated jet
1136 // leading particle and leading reconstructed jet
1137 // leading generated jet and leading reconstructed jet
1139 // Fragmentation functions and Shapes
1140 // a) integrated over all pt
1142 // b.1) only for user selected particles in jet
1143 // b.2) only for user rejected particles in jet
1144 // b.3) for all particles in jet