-/**************************************************************************\r
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- * *\r
- * Author: The ALICE Off-line Project. *\r
- * Contributors are mentioned in the code where appropriate. *\r
- * *\r
- * Permission to use, copy, modify and distribute this software and its *\r
- * documentation strictly for non-commercial purposes is hereby granted *\r
- * without fee, provided that the above copyright notice appears in all *\r
- * copies and that both the copyright notice and this permission notice *\r
- * appear in the supporting documentation. The authors make no claims *\r
- * about the suitability of this software for any purpose. It is *\r
- * provided "as is" without express or implied warranty. *\r
- **************************************************************************/\r
-\r
-//---------------------------------------------------------------------\r
-//\r
-//\r
-//---------------------------------------------------------------------\r
-#include <vector>\r
-#include <iostream>\r
-#include <Riostream.h>\r
-#include <TFile.h>\r
-#include <TCanvas.h>\r
-#include <TROOT.h>\r
-#include <TClonesArray.h>\r
-#include <TLorentzVector.h>\r
-#include <TH1F.h>\r
-#include <TH2F.h>\r
-#include <TBits.h>\r
-#include <TArrayF.h>\r
-#include "AliCdfJetFinder.h"\r
-#include "AliCdfJetHeader.h"\r
-#include "AliJetReader.h"\r
-#include "AliJetReaderHeader.h"\r
-#include "AliJet.h"\r
-#include "AliAODJet.h"\r
-#include "AliAODEvent.h"\r
-#include "TProfile.h"\r
-\r
-using namespace std ;\r
-\r
-struct varContainer // container for Particle variables\r
- { // variables of container struct\r
- Double_t pt; Double_t eta; Double_t phi;\r
- Int_t njet; // if jets are stored in varContainer nr_jet is multiplicity of jet\r
- } ;\r
-\r
-ClassImp ( AliCdfJetFinder )\r
-\r
-//______________________________________________________________________________\r
-AliCdfJetFinder::AliCdfJetFinder():\r
- AliJetFinder(),\r
- fHistos(0)\r
- { /* Constructor */ }\r
-\r
-//______________________________________________________________________________\r
-AliCdfJetFinder::~AliCdfJetFinder()\r
-\r
- {\r
- // destructor\r
- cout << "Calling the destructor ... " << endl ;\r
- Reset();\r
- cout << "Destructor called!" << endl;\r
- }\r
-\r
-//______________________________________________________________________________\r
-void AliCdfJetFinder::CreateOutputObjects(TList *histos)\r
-{\r
- // Create the list of histograms. Only the list is owned.\r
- fHistos = histos;\r
-\r
- TH1F *h1 = new TH1F ("histo1", "no jets/event", 100, 0,20);\r
- h1->SetStats(kTRUE);\r
- h1->GetXaxis()->SetTitle("N_{jets}");\r
- h1->GetYaxis()->SetTitle("#frac{dN}{dN_{jets}}");\r
- h1->GetXaxis()->SetTitleColor(1);\r
- h1->SetMarkerStyle(kFullCircle);\r
- fHistos->Add(h1);\r
-\r
- TH1F *h2 = new TH1F ("histo2", "no part/leading jet", 40, 0,40);\r
- h2->SetStats(kTRUE);\r
- h2->GetXaxis()->SetTitle("N_{leading}");\r
- h2->GetYaxis()->SetTitle("#frac{dN}{dN_{leading}}");\r
- h2->GetXaxis()->SetTitleColor(1);\r
- h2->SetMarkerStyle(kFullCircle);\r
- fHistos->Add(h2);\r
-\r
- TProfile * h3 = new TProfile ("histo3","ProfileX of (pt,npart) of leading jets", 25, 0. ,25. , 0.,25. ) ;\r
- h3->SetStats(kTRUE);\r
- h3->GetXaxis()->SetTitle("P_{t}");\r
- h3->GetYaxis()->SetTitle("N_{leading}");\r
- h3->GetXaxis()->SetTitleColor(1);\r
- h3->SetMarkerStyle(kFullCircle);\r
- fHistos->Add(h3);\r
-\r
- TH1F *h4 = new TH1F ("histo4", "Charge momentum distribution for leading jet", 500, 0 , 5);\r
- fHistos->Add(h4);\r
-\r
- TProfile *h5 = new TProfile ("histo5", "ProfileX of N_{charge} vs dphi leading jet", 100 , 0. , 200. , 0 , 100 );\r
- fHistos->Add(h5);\r
-\r
- TH1F *h6 = new TH1F ("histo6", " \"Transverse\" Pt Distribution ", 70, 0 , 14);\r
- fHistos->Add(h6);\r
-\r
-\r
-//cout << "CreateOutputObjetcs done!" << endl ;\r
-}\r
-\r
-//______________________________________________________________________________\r
-void AliCdfJetFinder::FindJets()\r
- {\r
- //cout << "Begining of FindJets ..." << endl ;\r
- //1) Fill 1 array with momentum and initialise another array for indexes\r
- //2) Run algorithm\r
- // 2.1) Sort momentum array\r
- // 2.2) loop over arrays of TLorentzVectors\r
- // 2.3) flag as a possible jet\r
- //3) fill AliAODJet list\r
- Bool_t debug = kFALSE;\r
-\r
- TRefArray *refs = 0;\r
- Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");\r
- if (fromAod) { refs = fReader->GetReferences(); }\r
-\r
- AliCdfJetHeader *header = (AliCdfJetHeader*)fHeader;\r
- Double_t radius = header->GetRadius(); // get Radius from header\r
-// cout << "Radius is : " << radius << endl ;\r
-\r
- TClonesArray * vectArray = fReader->GetMomentumArray() ;\r
- if ( vectArray == 0 ) { cout << "Could not get the momentum array" << endl; return; }\r
-\r
- Int_t nPart = vectArray->GetEntries() ; // n particles in this event;\r
- if ( nPart == 0 ) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; }\r
-\r
- Int_t nJets = 0; // n jets in this event\r
-\r
- fJets->SetNinput ( nPart ) ; // number of input objects\r
-\r
- varContainer **vectParticle = new varContainer* [nPart]; // container for Particles\r
- varContainer **vectJet = new varContainer* [nPart]; // container for Jets\r
-\r
- Double_t *ptArray = new Double_t [nPart] ; // momentum array\r
- Int_t *idxArray = new Int_t [nPart] ; // index array of sorted pts\r
-\r
- // initialisation of momentum and index arrays\r
- // cout << "Filling idxArray && momArray for sorting " << endl;\r
- for ( Int_t i = 0 ; i < nPart ; i++ )\r
- {// SORTING STEP :: ptArray with data from TClonesArray of TLorentzVector\r
- TLorentzVector * lv = (TLorentzVector*) vectArray->At(i);\r
- // INITIALISATION of local arrays for temporary storage\r
- varContainer *aParticle = new varContainer;\r
- aParticle->pt = lv->Pt();\r
- aParticle->eta = lv->Eta();\r
- aParticle->phi = lv->Phi();\r
- aParticle->njet = -1;\r
- vectParticle[i] = aParticle;\r
-\r
- // initializing arrays\r
- idxArray [i] = 0 ;\r
- ptArray [i] = aParticle->pt ;\r
- }\r
-\r
- TMath::Sort ( nPart, ptArray, idxArray ) ; // get a sorted array of indexes with size = size of TClonesArray\r
-\r
- Double_t pt_jet = 0. ,\r
- eta_jet = 0. ,\r
- phi_jet = 0. ,\r
- sum_eta = 0. ,\r
- sum_phi = 0. ; // jet variables\r
-\r
- Int_t idxPtSort = 0 , // index of array of sorted pt indexes\r
- npartJet = 0 ; // number of particles in curent jet\r
-\r
- TBits lkup_table ( nPart ) ; // bit container of size npart\r
-\r
- while ( lkup_table.CountBits(0) != (UInt_t)nPart )\r
- { // loop over particles in event until all flags are set\r
- UInt_t first_non_flagged = lkup_table.FirstNullBit(0) ; // set the index to the first NON flagged bit ; less conditions\r
-\r
- npartJet = 0 ; // reseting number of particles in jet\r
- pt_jet = 0.;\r
- eta_jet = 0.;\r
- phi_jet = 0.;\r
- sum_eta = 0.;\r
- sum_phi = 0.;\r
-\r
- for ( UInt_t i_part = first_non_flagged ; i_part < (UInt_t)nPart ; i_part++ )\r
- {// iteration over particles in event\r
- // reseting variables\r
-\r
- // the loop is done over sorted array of pt\r
- idxPtSort = idxArray[i_part] ; // index entry of TLorentzVector from vect_arr\r
-\r
- if ( lkup_table.TestBitNumber(i_part) == 1 ) { continue; } // if 4vector is already flagged skip it\r
-\r
- //taking info from vectParticle ;\r
- Double_t pt_tmp = vectParticle[idxPtSort]->pt ;\r
- Double_t eta_tmp = vectParticle[idxPtSort]->eta ;\r
- Double_t phi_tmp = vectParticle[idxPtSort]->phi ;\r
- if (debug) printf(" particle %d: pt=%g\n", i_part, pt_tmp);\r
-\r
- // all angles are expressed in rad\r
-\r
- if ( TMath::Abs(eta_tmp) > 0.9 )\r
- {\r
- lkup_table.SetBitNumber ( i_part, 1 ) ; // mark particle as used to be skipped\r
- continue ;\r
- }\r
-\r
-\r
- if ( i_part == first_non_flagged )\r
- {// this is first particle in event; leading particle\r
- // initialise the jet variables with leading particle numbers\r
-\r
- npartJet++; // incrementing counter of particles in jet (nparJet == 1) (leading particle)\r
- pt_jet = pt_tmp ;\r
- eta_jet = eta_tmp ;\r
- phi_jet = phi_tmp ;\r
- if (debug) printf(" first part in jet: npartjet=%d pt_jet=%g\n", npartJet, pt_jet);\r
-\r
- lkup_table.SetBitNumber ( i_part, 1 ) ; // flag the 4vector\r
-\r
- vectParticle[idxPtSort]->njet = nJets ; // associate particle with current jet\r
-\r
- continue ; // skip to next particle\r
- }\r
-\r
- // condition to be in jet\r
- Double_t deta = eta_jet - eta_tmp ;\r
- Double_t dphi = phi_jet - phi_tmp ;\r
-\r
- if ( dphi < -TMath::Pi() ) { dphi = -dphi - 2.0 * TMath::Pi() ; }\r
- if ( dphi > TMath::Pi() ) { dphi = -dphi + 2.0 * TMath::Pi() ; }\r
-\r
- Double_t r_computed = TMath::Sqrt ( deta * deta + dphi * dphi ) ;\r
-\r
- Bool_t in_jet = ( r_computed <= radius ) ? 1 : 0 ; // if r_computed is within jet_r in_jet == 1 else 0\r
-\r
- if ( in_jet )\r
- { // calculus of jet variables\r
- npartJet++; // incrementing counter of particles in jet\r
-\r
- pt_jet += pt_tmp ;\r
- sum_eta += (pt_tmp * eta_tmp) ;\r
- sum_phi += (pt_tmp * phi_tmp) ;\r
- eta_jet = sum_eta / pt_jet ;\r
- phi_jet = sum_phi / pt_jet ;\r
-\r
- lkup_table.SetBitNumber ( i_part, 1 ) ; // flag the 4vector\r
- vectParticle[idxPtSort]->njet = nJets ; // setting in particle list the associated jet\r
- if (debug) printf(" particle added to jet: npartjet=%d pt_jet=%g\n", npartJet, pt_jet);\r
- continue ; // skip to next particle\r
- }\r
-\r
- }\r
- // end of iteration over event; one jet definition\r
-\r
- varContainer *aJet = new varContainer;\r
- aJet->pt = pt_jet;\r
- aJet->eta = eta_jet;\r
- aJet->phi = phi_jet;\r
- aJet->njet = npartJet;\r
- vectJet[nJets++] = aJet; // store the jet (from jet 0 to nJets-1 )\r
- if (debug) printf("=== current jet: npartjet=%d ptjet=%g\n", npartJet,pt_jet);\r
-\r
- // writing aod information\r
- Double_t px = 0., py = 0., pz = 0., en = 0.; // convert to 4-vector\r
- px = pt_jet * TMath::Cos ( phi_jet ) ;\r
- py = pt_jet * TMath::Sin ( phi_jet ) ;\r
- pz = pt_jet / TMath::Tan ( 2.0 * TMath::ATan ( TMath::Exp ( -eta_jet ) ) ) ;\r
- en = TMath::Sqrt ( px * px + py * py + pz * pz );\r
-\r
- if (npartJet<2) continue; // do not add jets with less than 2 particles\r
- if (debug) cout << "Jet 4vect : " << "px = " << px << " ; py = " << py << " ; pz = " << pz << " ; E = " << en << endl;\r
-\r
- AliAODJet jet (px, py, pz, en);\r
-// cout << "Printing jet " << endl;\r
- if (debug) jet.Print("");\r
-\r
- if (fromAod) {\r
- for (Int_t parts = 0; parts < nPart; parts++ ) { \r
- if (vectParticle[parts]->njet == nJets-1) \r
- jet.AddTrack(refs->At(parts));\r
- }\r
- }\r
-\r
-// cout << "Adding jet ... " ;\r
- AddJet(jet);\r
-// cout << "added \n" << endl;\r
-\r
- }\r
- // end of while loop over particles ; ends when all particles were flagged as used and all jets defined\r
-\r
- /////////////////////////////////\r
- //// END OF EVENT PARSING ///\r
- /////////////////////////////////\r
-\r
- Int_t *jets_pt_idx = 0; // sorted array of jets pt\r
- Double_t *jets_pt = 0; // array of jets pts\r
- Int_t leading_jet_index = -1 ; // index of leading jet from vectJet\r
- Int_t part_leadjet = 0 ; // number of particles in leading jet\r
- Double_t pt_leadjet = 0. ; // pt of leading jet\r
- Double_t eta_leadjet = 0. ; // eta of leading jet\r
- Double_t phi_leadjet = 0. ; // phi of leading jet\r
- Int_t * idx_leadjet_part = 0;\r
-\r
- if (nJets > 0)\r
- { // if there is at least one jet in event\r
- jets_pt_idx = new Int_t [nJets] ;\r
- jets_pt = new Double_t [nJets] ;\r
-\r
- // filing the idx_ptjets array\r
- if (debug) printf("List of unsorted jets:\n");\r
- for( Int_t i = 0 ; i < nJets ; i++ )\r
- {\r
- jets_pt_idx [i] = 0 ;\r
- jets_pt [i] = vectJet[i]->pt ;\r
- if (debug) cout << " jet found: " << i << " npartjet=" << vectJet[i]->njet << " jets_pt[i]= " << jets_pt [i] << endl;\r
- }\r
- TMath::Sort ( nJets, jets_pt , jets_pt_idx ) ; // sorting pt of jets\r
-\r
- // selection of leading jet with nr of particles > 1\r
- for( Int_t i = 0 ; i < nJets ; i++ )\r
- {\r
- if ( vectJet[ jets_pt_idx[i] ]->njet > 1 )\r
- {\r
- leading_jet_index = jets_pt_idx[i] ;\r
- part_leadjet = vectJet[ leading_jet_index ]->njet ; // number of particles in leading jet\r
- pt_leadjet = vectJet[ leading_jet_index ]->pt ; // pt of leading jet\r
- eta_leadjet = vectJet[ leading_jet_index ]->eta ; // eta of leading jet\r
- phi_leadjet = vectJet[ leading_jet_index ]->phi ; // phi of leading jet\r
-\r
- cout << "Leading jet: npart = " << part_leadjet\r
- << " ; pt = " << pt_leadjet\r
- << " ; phi = " << phi_leadjet\r
- << " ; eta = " << eta_leadjet\r
- << endl ;\r
- break ;\r
- }\r
- }\r
- // end of selection of leading jet\r
-\r
- // Filling an array with indexes of leading jet particles\r
- idx_leadjet_part = new Int_t [part_leadjet] ;\r
- Int_t counter = 0;\r
- if (debug) printf(" Searching particles with jet index %d\n", leading_jet_index);\r
- for( Int_t i = 0 ; i <nPart ; i++ )\r
- {\r
- if ( vectParticle[i]->njet == leading_jet_index )\r
- {\r
- idx_leadjet_part[counter++] = i ;\r
- if (debug) cout << " " << counter-1 << ": index=" << i << " pt=" << vectParticle[i]->pt << endl;\r
- }\r
- }\r
- if ( (counter-1) > part_leadjet ) { cout << " Counter > part_leadjet !!!!" << endl;}\r
-\r
- }\r
- // end of loop over DETECTED jets\r
-\r
-\r
-\r
-// Calculus of part distribution in leading jet\r
-Double_t z = 0. ;\r
-Double_t *z_part_ljet = new Double_t [ part_leadjet ] ; // array of z of particles in leading jet\r
-for( Int_t j = 0 ; j < part_leadjet ; j++ )\r
- {\r
- Double_t z_j = vectParticle[idx_leadjet_part[j]]->pt ;\r
- z = z_j / pt_leadjet ;\r
- z_part_ljet [j] = z ;\r
- cout << "idx_leadjet_part[j] = " << idx_leadjet_part[j]\r
- << " p of particle = " << z_j\r
- << " pt lead jet = " << pt_leadjet\r
- << " Z = " << z << endl;\r
- }\r
-\r
-\r
-/*\r
-Double_t dphi_partLJ = 0. ;\r
-Double_t dphiArray [nPart];\r
-for( Int_t part = 0 ; part < nPart ; part++ )\r
- {\r
- dphi_partLJ = phi_leadjet - phi_partT [part] ;\r
-\r
- if ( dphi_partLJ < -TMath::Pi() ) { dphi_partLJ = -dphi_partLJ - 2.0 * TMath::Pi() ; }\r
- if ( dphi_partLJ > TMath::Pi() ) { dphi_partLJ = -dphi_partLJ + 2.0 * TMath::Pi() ; }\r
-\r
- dphi_partLJ = TMath::Abs(dphi_partLJ) ;\r
-\r
- dphiArray [part] = dphi_partLJ ;\r
- }\r
-*/\r
- // Fill nr_jet histogram\r
- if (fHistos)\r
- {\r
-// printf("FILLING histograms\n");\r
- TH1F *h_jets = (TH1F*)fHistos->FindObject("histo1");\r
- if (h_jets) h_jets->Fill(nJets);\r
-\r
- // There may be no leading jet with more that 2 particles !\r
- TH1F *h_leadpart = (TH1F*)fHistos->FindObject("histo2");\r
- if (part_leadjet && h_leadpart) h_leadpart->Fill(part_leadjet);\r
-\r
- TProfile * h_prof = (TProfile*)fHistos->FindObject("histo3");\r
- if (part_leadjet && h_prof) h_prof->Fill(pt_leadjet,part_leadjet);\r
-\r
- TH1F *h_MD = (TH1F*)fHistos->FindObject("histo4");\r
- for( Int_t k = 0 ; k < part_leadjet ; k++)\r
- { h_MD->Fill( z_part_ljet[k] ); }\r
-\r
-\r
-// TProfile *h_phi = (TProfile*)fHistos->At(4);\r
-// for( Int_t k = 0 ; k < nPart ; k++)\r
-// { h_phi->Fill( TMath::RadToDeg() * dphiArray [k] , nPart ) ; }\r
-\r
-\r
-// TH1F *h_tpd = (TH1F*)fHistos->At(5);\r
-// for( Int_t k = 0 ; k < nPart ; k++)\r
-// {\r
-// if ( (phi_partT[k] > (TMath::Pi()/3.)) && (phi_partT[k] < (2. * TMath::Pi()/3.)) )\r
-// { h_tpd->Fill( pt_partT [k] ) ; }\r
-// }\r
-\r
-\r
- }\r
-\r
-\r
- // CLEANING SECTION\r
- for (Int_t i=0; i<nPart; i++) delete vectParticle[i];\r
- delete [] vectParticle;\r
- for (Int_t i=0; i<nJets; i++) delete vectJet[i];\r
- delete [] vectJet;\r
- delete [] ptArray;\r
- delete [] idxArray;\r
-\r
- if (z_part_ljet) delete [] z_part_ljet;\r
- if (idx_leadjet_part) delete [] idx_leadjet_part;\r
- if (jets_pt_idx) delete [] jets_pt_idx ;\r
- if (jets_pt) delete [] jets_pt ;\r
-\r
- vectArray->Delete() ; // TClonesArray of lorentz vectors\r
- }\r
-// end of FindJets\r
-\r
-\r
-//______________________________________________________________________________\r
-void AliCdfJetFinder::FinishRun()\r
-{}\r
-\r
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//---------------------------------------------------------------------
+// Jet Finder based on CDF algortihm
+// Charged jet evolution and the underlying event in proton-antiproton collisions at 1.8 TeV
+// Physical Review D, vol. 65, Issue 9, id. 092002
+// http://www.phys.ufl.edu/~rfield/cdf/chgjet/chgjet_intro.html
+// Authors : Adrian.Sevcenco@cern.ch (adriansev@spacescience.ro )
+// Daniel.Felea@cern.ch (dfelea@spacescience.ro)
+// Ciprian.Mihai.Mitu@cern.ch (mcm@spacescience.ro)
+//---------------------------------------------------------------------
+
+/*
+Changelog
+
+
+
+*/
+
+#include <Riostream.h>
+#include <TROOT.h>
+#include <TMath.h>
+#include <TBits.h>
+#include <TFile.h>
+#include <TCanvas.h>
+#include <TClonesArray.h>
+#include <TLorentzVector.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TProfile.h>
+#include <TArrayF.h>
+#include <TVector2.h>
+
+#include "AliJetReader.h"
+#include "AliJetReaderHeader.h"
+#include "AliJet.h"
+#include "AliAODJet.h"
+#include "AliAODEvent.h"
+#include "AliJetFinder.h"
+
+#include "AliCdfJetFinder.h"
+#include "AliCdfJetHeader.h"
+
+ClassImp(AliCdfJetFinder)
+
+//______________________________________________________________________________
+AliCdfJetFinder::AliCdfJetFinder():
+ AliJetFinder(),
+ fHistos(0),
+ fDebug(0),
+ fFromAod(0),
+ fAODwrite(0),
+ fAODtracksWrite(0),
+ fAnalyseJets(0),
+ fRefArr (NULL),
+ fNJets(-9999),
+ fNPart(-9999),
+ fRadius(0.7),
+ fMinJetParticles(1),
+ fJetPtCut(0.),
+ fVectParticle(NULL),
+ fVectJet(NULL),
+ fPtArray(NULL),
+ fIdxArray(NULL)
+ { /* Constructor */ }
+
+//______________________________________________________________________________
+AliCdfJetFinder::~AliCdfJetFinder()
+ {
+ // destructor
+ }
+
+//______________________________________________________________________________
+void AliCdfJetFinder::CreateOutputObjects(TList * const histos)
+{
+ // Create the list of histograms. Only the list is owned.
+ fHistos = histos;
+
+// gStyle->SetOptStat(11111111);
+
+ TH1F *h1 = new TH1F ("histo1", "Pt distribution of jets", 200, 0,200); // 1GeV/bin
+ h1->SetStats(kTRUE);
+ h1->GetXaxis()->SetTitle("P_{T} of jets");
+ h1->GetYaxis()->SetTitle("Number of jets");
+ h1->GetXaxis()->SetTitleColor(1);
+ h1->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h1);
+
+ TH1F *h2 = new TH1F ("histo2", "Eta distribution of jets", 240, -1.2,1.2); // 1 unit of rapidity / 100 bin
+ h2->SetStats(kTRUE);
+ h2->GetXaxis()->SetTitle("Eta of jets");
+ h2->GetYaxis()->SetTitle("Number of jets");
+ h2->GetXaxis()->SetTitleColor(1);
+ h2->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h2);
+
+ TH1F *h3 = new TH1F ("histo3", "Phi distribution of jets", 400, -4,4);
+ h3->SetStats(kTRUE);
+ h3->GetXaxis()->SetTitle("Phi of jets");
+ h3->GetYaxis()->SetTitle("Number of jets");
+ h3->GetXaxis()->SetTitleColor(1);
+ h3->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h3);
+
+ TH1F *h4 = new TH1F ("histo4", "Multiplicity of jets", 40, 0,40); // 1 unit of multiplicity /bin
+ h4->SetStats(kTRUE);
+ h4->GetXaxis()->SetTitle("Particles in jets");
+ h4->GetYaxis()->SetTitle("Number of jets");
+ h4->GetXaxis()->SetTitleColor(1);
+ h4->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h4);
+
+ TH1F *h5 = new TH1F ("histo5", "Distribution of jets in events", 100, 0,100);
+ h5->SetStats(kTRUE);
+ h5->GetXaxis()->SetTitle("Number of jets");
+ h5->GetYaxis()->SetTitle("Number of events");
+ h5->GetXaxis()->SetTitleColor(1);
+ h5->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h5);
+
+ TH1F *h6 = new TH1F ("histo6", "Jet1 Charged Multiplicity Distribution", 30, 0,30);
+ h6->SetStats(kTRUE);
+ h6->GetXaxis()->SetTitle("N_{chg}");
+ h6->GetYaxis()->SetTitle("Number of jets");
+ h6->GetXaxis()->SetTitleColor(1);
+ h6->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h6);
+
+ TProfile * h7 = new TProfile ("histo7","N_{chg}(jet1) vs P_{T}(charged jet1)", 200, 0. ,200. , 0.,200. ) ;
+ h7->SetStats(kTRUE);
+ h7->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h7->GetYaxis()->SetTitle("<N_{chg}(jet1)> in 1 GeV/c bin");
+ h7->GetXaxis()->SetTitleColor(1);
+ h7->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h7);
+
+ TH1F *h8 = new TH1F ("histo8", "Charge momentum distribution for leading jet", 120, 0 , 1.2);
+ h8->SetStats(kTRUE);
+ h8->GetXaxis()->SetTitle("Jets");
+ h8->GetYaxis()->SetTitle("Particle distribution");
+ h8->GetXaxis()->SetTitleColor(1);
+ h8->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h8);
+
+ TProfile *h9 = new TProfile ("histo9", "N_{chg} vs the Azimuthal Angle from Charged Jet1", 50 , 0. , 180. , 0 , 20 );
+ h9->SetStats(kTRUE);
+ h9->GetXaxis()->SetTitle("#Delta#phi (degrees)");
+ h9->GetYaxis()->SetTitle("<N_{chg}> in 3.6 degree bin");
+ h9->GetXaxis()->SetTitleColor(1);
+ h9->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h9);
+
+ TProfile *h10 = new TProfile ("histo10", "P_{T} sum vs the Azimuthal Angle from Charged Jet1", 50 , 0. , 180. , 0 , 100 );
+ h10->SetStats(kTRUE);
+ h10->GetXaxis()->SetTitle("#Delta#phi (degrees)");
+ h10->GetYaxis()->SetTitle("<P_{T} sum> in 3.6 degree bin");
+ h10->GetXaxis()->SetTitleColor(1);
+ h10->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h10);
+
+ TH1F *h11 = new TH1F ("histo11", " \"Transverse\" Pt Distribution ", 70, 0 , 14);
+ h11->SetStats(kTRUE);
+ h11->GetXaxis()->SetTitle("P_{T} (GeV/c)");
+ h11->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
+ h11->GetXaxis()->SetTitleColor(1);
+ h11->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h11);
+
+ TH1F *h20 = new TH1F ("histo20", "Distribution of R in leading jet", 400, 0.,4.);
+ h20->SetStats(kTRUE);
+ h20->GetXaxis()->SetTitle("R [formula]");
+ h20->GetYaxis()->SetTitle("dN/dR");
+ h20->GetXaxis()->SetTitleColor(1);
+ h20->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h20);
+
+ TProfile * h21 = new TProfile ("histo21","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 30. ) ;
+ h21->SetStats(kTRUE);
+ h21->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h21->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
+ h21->GetXaxis()->SetTitleColor(1);
+ h21->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h21);
+
+ TProfile * h22 = new TProfile ("histo22","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
+ h22->SetStats(kTRUE);
+ h22->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h22->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
+ h22->GetXaxis()->SetTitleColor(1);
+ h22->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h22);
+
+ TProfile * h21Toward = new TProfile ("histo21_toward","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
+ h21Toward->SetStats(kTRUE);
+ h21Toward->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h21Toward->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
+ h21Toward->GetXaxis()->SetTitleColor(1);
+ h21Toward->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h21Toward);
+
+ TProfile * h21Transverse = new TProfile ("histo21_transverse","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
+ h21Transverse->SetStats(kTRUE);
+ h21Transverse->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h21Transverse->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
+ h21Transverse->GetXaxis()->SetTitleColor(1);
+ h21Transverse->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h21Transverse);
+
+ TProfile * h21Away = new TProfile ("histo21_away","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
+ h21Away->SetStats(kTRUE);
+ h21Away->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h21Away->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
+ h21Away->GetXaxis()->SetTitleColor(1);
+ h21Away->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h21Away);
+
+ TProfile * h22Toward = new TProfile ("histo22_toward","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
+ h22Toward->SetStats(kTRUE);
+ h22Toward->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h22Toward->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
+ h22Toward->GetXaxis()->SetTitleColor(1);
+ h22Toward->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h22Toward);
+
+ TProfile * h22Transverse = new TProfile ("histo22_transverse","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
+ h22Transverse->SetStats(kTRUE);
+ h22Transverse->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h22Transverse->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
+ h22Transverse->GetXaxis()->SetTitleColor(1);
+ h22Transverse->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h22Transverse);
+
+ TProfile * h22Away = new TProfile ("histo22_away","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
+ h22Away->SetStats(kTRUE);
+ h22Away->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h22Away->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
+ h22Away->GetXaxis()->SetTitleColor(1);
+ h22Away->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h22Away);
+
+ TH1F *h23Toward = new TH1F ("histo23_toward","'Toward' Pt Distribution of charged particles", 200, 0., 14.);
+ h23Toward->SetStats(kTRUE);
+ h23Toward->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
+ h23Toward->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
+ h23Toward->GetXaxis()->SetTitleColor(1);
+ h23Toward->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h23Toward);
+
+ TH1F *h23Transverse = new TH1F ("histo23_transverse","'Transverse' Pt Distribution of charged particles", 200, 0., 14.);
+ h23Transverse->SetStats(kTRUE);
+ h23Transverse->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
+ h23Transverse->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
+ h23Transverse->GetXaxis()->SetTitleColor(1);
+ h23Transverse->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h23Transverse);
+
+ TH1F *h23Away = new TH1F ("histo23_away","'Away' Pt Distribution of charged particles", 200, 0., 14.);
+ h23Away->SetStats(kTRUE);
+ h23Away->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
+ h23Away->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
+ h23Away->GetXaxis()->SetTitleColor(1);
+ h23Away->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h23Away);
+
+ TProfile * h24 = new TProfile ("histo24","Jet1 Size vs P_{T}(charged jet1)", 200, 0., 50. , 0., 0.5) ;
+ h24->SetStats(kTRUE);
+ h24->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h24->GetYaxis()->SetTitle("<R(chgjet1)> in 1 GeV/c bin");
+ h24->GetXaxis()->SetTitleColor(1);
+ h24->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h24);
+
+ TProfile * h25 = new TProfile ("histo25","Jet1 Size vs P_{T}(charged jet1)", 200, 0., 50. , 0., 0.5) ;
+ h25->SetStats(kTRUE);
+ h25->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
+ h25->GetYaxis()->SetTitle("<R(chgjet1)> in 1 GeV/c bin");
+ h25->GetXaxis()->SetTitleColor(1);
+ h25->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h25);
+
+ TProfile *h26 = new TProfile ("histo26", "N_{chg} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0., 0.8);
+ h26->SetStats(kTRUE);
+ h26->GetXaxis()->SetTitle("Distance R");
+ h26->GetYaxis()->SetTitle("<N_{chg}> in 0.02 bin");
+ h26->GetXaxis()->SetTitleColor(1);
+ h26->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h26);
+
+ TProfile *h27 = new TProfile ("histo27", "N_{chg} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0., 0.8);
+ h27->SetStats(kTRUE);
+ h27->GetXaxis()->SetTitle("Distance R");
+ h27->GetYaxis()->SetTitle("<N_{chg}> in 0.02 bin");
+ h27->GetXaxis()->SetTitleColor(1);
+ h27->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h27);
+
+ TProfile *h28 = new TProfile ("histo28", "PT_{sum} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0.01, 10.);
+ h28->SetStats(kTRUE);
+ h28->GetXaxis()->SetTitle("Distance R");
+ h28->GetYaxis()->SetTitle("<PT_{sum} (GeV/c)> in 0.02 bin");
+ h28->GetXaxis()->SetTitleColor(1);
+ h28->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h28);
+
+ TProfile *h29 = new TProfile ("histo29", "PT_{sum} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0.01, 10.);
+ h29->SetStats(kTRUE);
+ h29->GetXaxis()->SetTitle("Distance R");
+ h29->GetYaxis()->SetTitle("<PT_{sum} (GeV/c)> in 0.02 bin");
+ h29->GetXaxis()->SetTitleColor(1);
+ h29->SetMarkerStyle(kFullCircle);
+ fHistos->Add(h29);
+
+}
+
+//______________________________________________________________________________
+void AliCdfJetFinder::FindJets()
+{
+// Jet Algorithm:
+// * Order all charged particles according to their PT.
+// * Start with the highest PT particle and include in the "jet" all particles within the "radius" R = 0.7
+// (considering each particle in the order of decreasing PT and recalculating the centroid of the jet after
+// each new particle is added to the jet).
+// * Go to the next highest PT particle (not already included in a jet) and include in the "jet" all particles
+// (not already included in a jet) within the radius R =0.7.
+// * Continue until all particles are in a "jet".
+
+AliCdfJetHeader *header = (AliCdfJetHeader*)fHeader;
+
+ if (header)
+ {
+ fDebug = header->IsDebugCDF();
+ fAODwrite = header->IsAODwrite() ; // write jets to AOD
+ fAODtracksWrite = header->IsAODtracksWrite() ; // write jet tracks to AOD
+ fRadius = header->GetRadius(); // get Radius from jet finder header
+ fMinJetParticles = header->GetMinPartJet (); // get minimum multiplicity of an jet
+ fJetPtCut = header->GetJetPtCut (); // get minimum of jet pt
+ }
+ else
+ { cout << "Header not found" << endl; return; }
+
+if (fAODwrite)
+ {
+ fFromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
+ if (fFromAod) { fRefArr = fReader->GetReferences(); }
+ }
+
+InitData();
+
+if (!fNPart) { if (fDebug) {cout << "entries = 0 ; Event empty !!!" << endl ;} return; } // if event empty then exit
+
+FindCones();
+
+ComputeConesWeight();
+
+if (fAODwrite) { cout << "Writing AOD" << endl ; WriteJets(); }
+
+if (fAnalyseJets) AnalizeJets();
+
+Clean();
+
+}
+
+//______________________________________________________________________________
+void AliCdfJetFinder::InitData()
+{
+// initialisation of variables and data members
+
+ TClonesArray * vectArray = fReader->GetMomentumArray() ;
+ if ( vectArray == 0 ) { cout << "Could not get the momentum array" << endl; return; }
+
+ fNPart = vectArray->GetEntries() ; // n particles in this event;
+
+ if ( fNPart == 0 ) { return ; } // if event empty then exit
+
+ fJets->SetNinput ( fNPart ) ; // number of input objects
+
+ fVectParticle = new varContainer* [fNPart]; // container for Particles
+
+ fPtArray = new Double_t [fNPart] ; // momentum array
+ fIdxArray = new Int_t [fNPart] ; // index array of sorted pts
+
+ // initialisation of momentum and index arrays
+ for ( Int_t i = 0 ; i < fNPart ; i++ )
+ {// SORTING STEP :: fPtArray with data from TClonesArray of TLorentzVector
+ TLorentzVector * lv = (TLorentzVector*) vectArray->At(i);
+
+ // INITIALISATION of local arrays for temporary storage
+ varContainer *aParticle = new varContainer;
+ aParticle->pt = lv->Pt();
+ aParticle->eta = lv->Eta();
+ aParticle->phi = TVector2::Phi_mpi_pi ( lv->Phi() ); // normalize to -pi,pi
+ aParticle->njet = -999;
+
+ fVectParticle[i] = aParticle; // vector of Particles
+
+ // initializing arrays
+ fIdxArray [i] = -999 ;
+ fPtArray [i] = aParticle->pt ;
+ }
+
+ TMath::Sort ( fNPart, fPtArray, fIdxArray ) ; // get a sorted array of indexes with TClonesArray.Size()
+
+}
+
+
+//______________________________________________________________________________
+void AliCdfJetFinder::FindCones()
+{
+// parsing of particles in event and estlabish jets (label them with jet index)
+
+ Double_t ptSeed = 0. , etaSeed = 0. , phiSeed = 0. ; // leading particle params
+ Double_t pttmp = 0. , etatmp = 0. , phitmp = 0. ; // temporary variables to be used in various calculations
+ Double_t deta = 0. , dphi = 0. , dcomputed = 0. ;
+ Bool_t injet = 0 ;
+
+ fNJets = -1 ; // n jets in this event
+ Int_t idxPtSort = -1 ; // index of array of sorted pt indexes
+
+ if (fDebug) { cout << "\n\n\n\n\n\n------------------\nBegin Event Analysis\n------------------\n\n" << endl ;}
+
+ cout << "fNPart = " << fNPart << endl;
+
+ TBits lkupTable ( fNPart ) ; // bit container ; 1-to-1 corespondence with fIdxArray
+
+ while ( lkupTable.CountBits() != (UInt_t)fNPart )
+ { // loop over particles in event until all flags are set
+ UInt_t firstnonflagged = lkupTable.FirstNullBit() ; // set the index to the first NON flagged bit ; less conditions
+
+ cout << "\n\nfirst_non_flagged : " << firstnonflagged << endl;
+
+ ++fNJets; // incrementing the jet counter
+ if (fDebug) { printf("JET %d \n", fNJets); }
+
+ ptSeed = 0. ; etaSeed = 0. ; phiSeed = 0. ; // reseting leading particle params
+
+ for ( UInt_t ipart = firstnonflagged ; ipart < (UInt_t)fNPart ; ipart++ )
+ {// iteration over particles in event
+ // the loop is done over sorted array of pt
+ idxPtSort = fIdxArray[ipart] ; // index of particle ! fIdxArray is an index list pt sorted
+
+ if ( lkupTable.TestBitNumber(ipart) ) { continue; } // if 4vector is already flagged skip it
+
+ //init computed and used vars
+ pttmp = 0. ; etatmp = 0. ; phitmp = 0. ;
+ deta = 0. ; dphi = 0. ; dcomputed = 0. ; injet = 0 ;
+
+ //taking info from fVectParticle ;
+ pttmp = fVectParticle[idxPtSort]->pt ;
+ etatmp = fVectParticle[idxPtSort]->eta ;
+ phitmp = fVectParticle[idxPtSort]->phi ;
+
+ if ( ipart == firstnonflagged )
+ {// this is first particle in event; leading particle
+ // begin the search around this particle in a fRadius
+
+ // CENTRE OF THE JET
+ ptSeed = pttmp ; etaSeed = etatmp ; phiSeed = phitmp ; // seeding the jet with first particle idxPtSort
+
+ lkupTable.SetBitNumber ( ipart ) ; // flag the index of particle in lkup_table
+ fVectParticle[idxPtSort]->njet = fNJets ; // associate particle with current jet number
+
+ if (fDebug) { printf("\nLeading particle :: particle index = %d ; at sorted index %d ; in jet %d \n", idxPtSort, ipart, fNJets); }
+ if (fDebug) { printf("pt= %g ; eta= %g ; phi = %g \n", pttmp, etatmp, phitmp) ; }
+ if (fDebug) { lkupTable.Print() ;}
+
+ continue ; // skip to next particle
+ }
+
+ // condition to be in jet
+ deta = etatmp - etaSeed ;
+ dphi = TVector2::Phi_mpi_pi (phitmp - phiSeed) ; // computing dphi and normalizing to (0,2pi) interval in one step
+
+ dcomputed = TMath::Hypot(deta, dphi) ; // Distance(fRadius) to (eta,phi) seed
+
+ injet = ( ( fRadius - dcomputed ) >= 0.000000001 ) ? 1 : 0 ; // if r_computed is within jet_r in_jet == 1 else 0
+
+ if ( injet )
+ { // calculus of jet variables
+ lkupTable.SetBitNumber ( ipart ) ; // flag the index of particle in lkup_table
+ fVectParticle[idxPtSort]->njet = fNJets ; // setting in particle list the associated jet
+
+ if (fDebug) { printf("\njet particle :: particle index = %d ; at sorted index %d ; in jet %d ; found at radius %g ; \n", idxPtSort, ipart, fNJets, dcomputed); }
+ if (fDebug) { printf("pt= %g ; eta= %g ; phi = %g \n", pttmp, etatmp, phitmp) ; }
+ if (fDebug) { lkupTable.Print() ;}
+
+ continue ; // skip to next particle
+ }
+
+ }
+ // end of iteration over event; one jet definition of content ; jet parameters to be computed later
+ }
+}
+
+
+//______________________________________________________________________________
+void AliCdfJetFinder::ComputeConesWeight()
+{
+// computing of jets Pt, Eta and Phi (centre of weight in (eta,phi) plane)
+// rescan the vector of particles by identify them by asociate jet number for computing of weight centre
+
+// JET CONTAINER
+fVectJet = new varContainer* [fNJets]; // container for Jets
+
+Double_t ptJet, ptJet2 , etaJet , phiJet ; Int_t npartJet ;
+Double_t pttmp = 0. , etatmp = 0. , phitmp = 0. ; // temporary variables to be used in various calculations
+Int_t idxPtSort = -999 ; // index of array of sorted pt indexes
+
+for( Int_t jet = 0 ; jet < fNJets ; jet++ )
+ {
+ if (fDebug) { printf("\n\n--- Computing weight of Jet %d \n", jet ); }
+ npartJet = 0 ; ptJet = 0. ; etaJet = 0. ; phiJet = 0. ; // reset variables for a new computation
+
+ for ( Int_t ipart = 0 ; ipart < fNPart ; ipart++ )
+ {// iteration over particles in event
+ // the loop is done over sorted array of pt
+ idxPtSort = fIdxArray[ipart] ; // index of particle ! fIdxArray is an index list pt sorted
+
+ if ( fVectParticle[idxPtSort]->njet == jet )
+ {
+ ++npartJet; // incrementing the counter of jet particles
+
+ //taking info from fVectParticle ;
+ pttmp = fVectParticle[idxPtSort]->pt ;
+ etatmp = fVectParticle[idxPtSort]->eta ;
+ phitmp = TVector2::Phi_mpi_pi (fVectParticle[idxPtSort]->phi) ;
+
+// jet_new_angular_coordinate = jet_old_angular_coordinate * jet_old_pt / jet_new_pt +
+// part[i]_angular_coordinate * part[i]_pt/jet_new_pt
+
+ ptJet2 = ptJet + pttmp ;
+
+ etaJet = etaJet * ptJet / ptJet2 + etatmp * pttmp / ptJet2 ;
+ phiJet = phiJet * ptJet / ptJet2 + phitmp * pttmp / ptJet2 ;
+
+ ptJet = ptJet2 ;
+
+ }
+ // add a particle and recalculation of centroid
+ }
+ // end of 1 jet computation
+
+ varContainer *aJet = new varContainer; // Jet container
+ aJet->pt = ptJet; aJet->eta = etaJet; aJet->phi = phiJet; aJet->njet = npartJet; // setting jet vars in container
+ fVectJet[jet] = aJet; // store the number of the jet(fNJets) and increment afterwards
+
+ if (fDebug) { printf ("=== current jet %d : npartjet= %d ; pt_jet= %g ; eta_jet = %g ; phi_jet = %g \n\n\n",
+ jet, npartJet, ptJet, etaJet, phiJet ) ; }
+
+ }
+ //end loop over jets
+
+}
+
+
+//______________________________________________________________________________
+void AliCdfJetFinder::WriteJets()
+{
+// Writing AOD jets and AOD tracks
+
+for( Int_t jetnr = 0 ; jetnr < fNJets ; jetnr++ )
+ {
+ Double_t pt = 0., eta = 0., phi = 0., // jet variables
+ px = 0., py = 0., pz = 0., en = 0.; // convert to 4-vector
+ pt = fVectJet[ jetnr ]->pt ; // pt of jet
+ eta = fVectJet[ jetnr ]->eta ; // eta of jet
+ phi = fVectJet[ jetnr ]->phi ; // phi of jet
+
+ px = pt * TMath::Cos ( phi ) ;
+ py = pt * TMath::Sin ( phi ) ;
+ pz = pt / TMath::Tan ( 2.0 * TMath::ATan ( TMath::Exp ( -eta ) ) ) ;
+ en = TMath::Sqrt ( px * px + py * py + pz * pz );
+
+ AliAODJet jet (px, py, pz, en);
+
+
+ if (fDebug) jet.Print("");
+
+ if (fFromAod && fAODtracksWrite)
+ {
+ for ( Int_t jetTrack = 0; jetTrack < fNPart; jetTrack++ )
+ { if ( fVectParticle[jetTrack]->njet == jetnr ) { jet.AddTrack(fRefArr->At(jetTrack)) ; } }
+ }
+ // tracks REFs written in AOD
+ AddJet(jet);
+ }
+//jets vector parsed and written to AOD
+}
+
+
+//______________________________________________________________________________
+void AliCdfJetFinder::AnalizeJets()
+{
+// analyzing of jets and filling of histograms
+
+ const Double_t kPI = TMath::Pi();
+
+ //persistent pointer to histo20
+ TH1F *hR = (TH1F*)fHistos->FindObject("histo20");
+
+ Int_t *jetsptidx = 0; // sorted array of jets pt
+ Double_t *jetspt = 0; // array of jets pts
+ Int_t leadingjetindex = -1 ; // index of leading jet from fVectJet
+ Int_t partleadjet = 0 ; // number of particles in leading jet
+ Double_t ptleadjet = 0. ; // pt of leading jet
+ Double_t etaleadjet = 0. ; // eta of leading jet
+ Double_t phileadjet = 0. ; // phi of leading jet
+
+ jetsptidx = new Int_t [fNJets] ;
+ jetspt = new Double_t [fNJets] ;
+
+//________________________________________________________________________________________
+// Jet sorting and finding the leading jet that coresponds to cuts in pt and multiplicity
+//________________________________________________________________________________________
+
+ // filing the idx_ptjets array
+ if (fDebug) printf("List of unsorted jets:\n");
+ for( Int_t i = 0 ; i < fNJets ; i++ )
+ {
+ jetsptidx [i] = 0 ;
+ jetspt [i] = fVectJet[i]->pt ;
+ if (fDebug) { cout << " jet found: " << i << " npartjet=" << fVectJet[i]->njet << " ; jets_pt = " << jetspt[i] << endl; }
+ }
+
+ TMath::Sort ( fNJets, jetspt , jetsptidx ) ; // sorting pt of jets
+
+ // selection of leading jet
+ // looping over jets searching for __first__ one that coresponds to cuts
+ for( Int_t i = 0 ; i < fNJets ; i++ )
+ {
+ if ( ( fVectJet[ jetsptidx[i] ]->njet >= fMinJetParticles ) && ( fVectJet[ jetsptidx[i] ]->pt >= fJetPtCut ) )
+ {
+ leadingjetindex = jetsptidx[i] ;
+ partleadjet = fVectJet[ leadingjetindex ]->njet ; // number of particles in leading jet
+ ptleadjet = fVectJet[ leadingjetindex ]->pt ; // pt of leading jet
+ etaleadjet = fVectJet[ leadingjetindex ]->eta ; // eta of leading jet
+ phileadjet = fVectJet[ leadingjetindex ]->phi ; // phi of leading jet
+
+ if (fDebug)
+ { printf("Leading jet %d : npart= %d ; pt= %g ; eta = %g ; phi = %g \n", leadingjetindex, partleadjet, ptleadjet, etaleadjet, phileadjet ); }
+
+ break ;
+ }
+ }
+ // end of selection of leading jet
+
+
+
+//////////////////////////////////////////////////
+//// Computing of values used in histograms
+//////////////////////////////////////////////////
+
+//___________________________________________________________________________
+// pt_sum of all particles in event
+//___________________________________________________________________________
+cout << "Computing sum of pt in event" << endl ;
+Double_t ptsumevent = 0.;
+for ( Int_t i = 0 ; i< fNPart ; i++ ) { ptsumevent += fVectParticle[i]->pt ; }
+printf ("Sum of all Pt in event : pt_sum_event = %g", ptsumevent) ;
+
+//___________________________________________________________________________
+// Filling an array with indexes of leading jet particles
+//___________________________________________________________________________
+Int_t * idxpartLJ = new Int_t [partleadjet] ;
+Int_t counterpartleadjet = 0;
+
+cout << "Filling an array with indexes of leading jet particles" << endl;
+
+for( Int_t i = 0 ; i < fNPart ; i++ )
+ {
+ if ( fVectParticle[i]->njet == leadingjetindex )
+ { idxpartLJ[counterpartleadjet++] = i ; }
+ }
+
+if ( (counterpartleadjet-1) > partleadjet ) { cout << " Counter_part_lead_jet > part_leadjet !!!!" << endl;}
+
+
+//___________________________________________________________________________
+// Calculus of part distribution in leading jet
+//___________________________________________________________________________
+Double_t z = 0. ;
+Double_t *zpartljet = new Double_t [ partleadjet ] ; // array of z of particles in leading jet
+
+cout << "Entering loop of calculus of part distribution in leading jet" << endl ;
+
+for( Int_t j = 0 ; j < partleadjet ; j++ )
+ {
+ Double_t zj = fVectParticle[idxpartLJ[j]]->pt ;
+ z = zj / ptleadjet ;
+ zpartljet [j] = z ;
+ cout << "idx_leadjet_part[j] = " << idxpartLJ[j]
+ << " p of particle = " << zj
+ << " pt lead jet = " << ptleadjet
+ << " Z = " << z << endl;
+ }
+
+
+//___________________________________________________________________________
+// array of delta phi's between phi of particles and leading jet phi
+//___________________________________________________________________________
+cout << "array of delta phi's between phi of particles and leading jet phi" << endl;
+Double_t dphipartLJ = 0. ;
+Double_t *dphipartljet = new Double_t [fNPart];
+for( Int_t part = 0 ; part < fNPart ; part++ )
+ {
+ dphipartLJ = fVectParticle[part]->phi - phileadjet ;
+ dphipartLJ = TVector2::Phi_mpi_pi (dphipartLJ) ; // restrict the delta phi to (-pi,pi) interval
+ dphipartljet [part] = dphipartLJ ;
+ printf("part= %d ; dphi_partLJ = %g \n", part, dphipartLJ );
+ }
+
+
+//______________________________________________________________________________
+// Pt distribution for all particles
+//______________________________________________________________________________
+TH1F * hpt = (TH1F*)fHistos->FindObject("histo11");
+if ( hpt ) { for ( Int_t i = 0 ; i < fNPart ; i++ ) { hpt->Fill( fVectParticle[i]->pt ); } }
+
+//___________________________________________________________________________
+// Recomputing of radius of particles in leading jet
+//___________________________________________________________________________
+if (fDebug) { printf(" Searching particles with jet index %d\n", leadingjetindex); }
+
+Double_t ddeta = 0. , ddphi = 0. , rpart = 0. ;
+
+for( Int_t j = 0 ; j < partleadjet ; j++ )
+ {
+ ddeta = etaleadjet - fVectParticle[idxpartLJ[j]]->eta;
+
+ Double_t phitmp = fVectParticle[idxpartLJ[j]]->phi ;
+
+ ddphi = TVector2::Phi_mpi_pi ( phileadjet - phitmp ) ; // restrict the delta phi to (-pi,pi) interval
+
+ rpart = TMath::Hypot (ddeta, ddphi) ;
+
+ printf ("Particle %d with Re-Computed radius = %f ", idxpartLJ[j], rpart) ;
+ if ( (rpart - fRadius) >= 0.00000001 )
+ { printf (" bigger than selected radius of %f\n", fRadius ); }
+ else
+ { printf ("\n") ; }
+
+ if (hR) hR->Fill(rpart);
+
+ }
+
+
+
+//_______________________________________________________________________
+// Computing of radius that contain 80% of Leading Jet ( PT and multiplicity )
+//_______________________________________________________________________
+Double_t corepartleadjet = 0.8 * partleadjet ;
+Double_t coreptleadjet = 0.8 * ptleadjet ;
+Int_t countercorepart = 0 ;
+Double_t countercorept = 0. ;
+Int_t sortedindex = -1 ;
+
+TProfile * hprof24 = (TProfile*)fHistos->FindObject("histo24");
+TProfile * hprof25 = (TProfile*)fHistos->FindObject("histo25");
+
+TProfile * hprof26 = (TProfile*)fHistos->FindObject("histo26");
+TProfile * hprof27 = (TProfile*)fHistos->FindObject("histo27");
+TProfile * hprof28 = (TProfile*)fHistos->FindObject("histo28");
+TProfile * hprof29 = (TProfile*)fHistos->FindObject("histo29");
+
+
+if ((hprof24) && (hprof25) && (hprof26) && (hprof27) && (hprof28) && (hprof29) )
+{
+for( Int_t part = 0 ; part < fNPart ; part++ )
+ {
+ Double_t pttmp = 0. ; Double_t etatmp = 0. ; Double_t phitmp = 0. ; // temporary variables
+ Double_t dpart = 0. ;
+ sortedindex = fIdxArray[part] ;
+
+ if ( fVectParticle [ sortedindex ]->njet == leadingjetindex )
+ {
+ pttmp = fVectParticle[sortedindex]->pt ;
+ etatmp = fVectParticle[sortedindex]->eta ;
+ phitmp = fVectParticle[sortedindex]->phi ;
+
+ ++countercorepart ;
+ countercorept += pttmp ;
+
+ dpart = TMath::Hypot ( etaleadjet - etatmp, TVector2::Phi_mpi_pi (phileadjet - phitmp) ) ;
+
+ if ( countercorepart <= corepartleadjet ) { hprof24->Fill(ptleadjet, dpart); }
+ if ( countercorept <= coreptleadjet ) { hprof25->Fill(ptleadjet, dpart); }
+
+ if (ptleadjet > 5.) { hprof26->Fill(dpart, countercorepart); hprof28->Fill(dpart, countercorept); }
+ if (ptleadjet > 30.) { hprof27->Fill(dpart, countercorepart); hprof29->Fill(dpart, countercorept); }
+
+ }
+ }
+}
+
+ TH1F *hjetpt = (TH1F*)fHistos->FindObject("histo1");
+ TH1F *hjeteta = (TH1F*)fHistos->FindObject("histo2");
+ TH1F *hjetphi = (TH1F*)fHistos->FindObject("histo3");
+ TH1F *hjetnjet = (TH1F*)fHistos->FindObject("histo4");
+
+ for( Int_t jet = 0 ; jet < fNJets ; jet++ )
+ {
+ if (hjetpt) hjetpt ->Fill ( fVectJet[jet]->pt ) ;
+ if (hjeteta) hjeteta ->Fill ( fVectJet[jet]->eta ) ;
+ if (hjetphi) hjetphi ->Fill ( fVectJet[jet]->phi ) ;
+ if (hjetnjet) hjetnjet ->Fill ( fVectJet[jet]->njet ) ;
+ }
+
+ TH1F *hjets = (TH1F*)fHistos->FindObject("histo5");
+ if (hjets) hjets->Fill(fNJets);
+
+ TH1F *hleadpart = (TH1F*)fHistos->FindObject("histo6");
+ if (hleadpart) hleadpart->Fill(partleadjet);
+
+ TProfile * hprof = (TProfile*)fHistos->FindObject("histo7");
+ if (hprof) hprof->Fill(ptleadjet,partleadjet);
+
+ TH1F *hMD = (TH1F*)fHistos->FindObject("histo8");
+ for( Int_t k = 0 ; k < partleadjet ; k++)
+ { hMD->Fill( zpartljet[k] ); }
+
+ TProfile * hphi = (TProfile*)fHistos->FindObject("histo9");
+ for( Int_t k = 0 ; k < partleadjet ; k++)
+ { hphi->Fill( TMath::RadToDeg() * dphipartljet [k] , fNPart ) ; }
+
+ TProfile * htpd = (TProfile*)fHistos->FindObject("histo10");
+ for( Int_t k = 0 ; k < fNPart ; k++)
+ { htpd->Fill( TMath::RadToDeg() * dphipartljet [k] , ptsumevent ) ; }
+
+
+ TProfile * hprof1 = (TProfile*)fHistos->FindObject("histo21");
+ if (hprof1) hprof1->Fill(ptleadjet, fNPart);
+
+ TProfile * hprof2 = (TProfile*)fHistos->FindObject("histo22");
+ if (hprof2) hprof2->Fill(ptleadjet, ptsumevent);
+
+ TProfile * hprof1toward = (TProfile*)fHistos->FindObject("histo21_toward");
+ TProfile * hprof1transverse = (TProfile*)fHistos->FindObject("histo21_transverse");
+ TProfile * hprof1away = (TProfile*)fHistos->FindObject("histo21_away");
+ TProfile * hprof2toward = (TProfile*)fHistos->FindObject("histo22_toward");
+ TProfile * hprof2transverse = (TProfile*)fHistos->FindObject("histo22_transverse");
+ TProfile * hprof2away = (TProfile*)fHistos->FindObject("histo22_away");
+ TH1F * hpttoward = (TH1F*)fHistos->FindObject("histo23_toward");
+ TH1F * hpttransverse = (TH1F*)fHistos->FindObject("histo23_transverse");
+ TH1F * hptaway = (TH1F*)fHistos->FindObject("histo23_away");
+
+ if ( (hprof1toward) && (hprof1transverse) && (hprof1away) && (hprof2toward) && (hprof2transverse) && (hprof2away) )
+ {
+ for( Int_t part = 0 ; part < fNPart ; part++)
+ {
+ Double_t ptpart = fVectParticle[part]->pt ; // pt of particle
+ if ( ( dphipartljet[part] >=0.) && ( dphipartljet[part] < kPI/3. ) )
+ {
+ hprof1toward->Fill( ptleadjet, fNPart );
+ hprof2toward->Fill( ptleadjet, ptsumevent);
+ hpttoward->Fill( ptpart );
+ }
+ else
+ if ( ( dphipartljet[part] >= (kPI/3.)) && ( dphipartljet[part] < (2.*kPI/3.)) )
+ {
+ hprof1transverse->Fill( ptleadjet, fNPart );
+ hprof2transverse->Fill( ptleadjet, ptsumevent);
+ hpttransverse->Fill( ptpart );
+ }
+ else
+ if ( ( dphipartljet[part] >= ( 2.*kPI/3.)) && ( dphipartljet[part] < kPI ) )
+ {
+ hprof1away->Fill( ptleadjet, fNPart );
+ hprof2away->Fill( ptleadjet, ptsumevent);
+ hptaway->Fill( ptpart );
+ }
+ }
+ }
+
+}
+
+
+//______________________________________________________________________________
+void AliCdfJetFinder::Clean()
+{
+// CLEANING SECTION
+ delete [] fVectParticle;
+ delete [] fVectJet;
+ delete [] fPtArray;
+ delete [] fIdxArray;
+
+ Reset();
+}
+
+//______________________________________________________________________________
+void AliCdfJetFinder::FinishRun()
+{
+// do i need this?
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff