JETAN/AliCdfJetFinder.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   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  **************************************************************************/
  15
  16 //---------------------------------------------------------------------
  17 //
  18 //
  19 //---------------------------------------------------------------------
  20 #include <vector>
  21 #include <iostream>
  22 #include <Riostream.h>
  23 #include <TFile.h>
  24 #include <TCanvas.h>
  25 #include <TROOT.h>
  26 #include <TClonesArray.h>
  27 #include <TLorentzVector.h>
  28 #include <TH1F.h>
  29 #include <TH2F.h>
  30 #include <TBits.h>
  31 #include <TArrayF.h>
  32 #include "AliCdfJetFinder.h"
  33 #include "AliCdfJetHeader.h"
  34 #include "AliJetReader.h"
  35 #include "AliJetReaderHeader.h"
  36 #include "AliJet.h"
  37 #include "AliAODJet.h"
  38 #include "AliAODEvent.h"
  39 #include "TProfile.h"
  40
  41 using namespace std ;
  42
  43 struct varContainer // container for Particle variables
  44   { // variables of container struct
  45   Double_t  pt; Double_t eta; Double_t phi;
  46   Int_t njet; // if jets are stored in varContainer nr_jet is multiplicity of jet
  47   } ;
  48
  49 ClassImp ( AliCdfJetFinder )
  50
  51 //______________________________________________________________________________
  52 AliCdfJetFinder::AliCdfJetFinder():
  53    AliJetFinder(),
  54    fHistos(0)
  55   {  /* Constructor */  }
  56
  57 //______________________________________________________________________________
  58 AliCdfJetFinder::~AliCdfJetFinder()
  59
  60   {
  61   // destructor
  62   cout << "Calling the destructor ... " << endl ;
  63   Reset();
  64   cout << "Destructor called!" << endl;
  65   }
  66
  67 //______________________________________________________________________________
  68 void AliCdfJetFinder::CreateOutputObjects(TList *histos)
  69 {
  70   // Create the list of histograms. Only the list is owned.
  71   fHistos = histos;
  72
  73   TH1F *h1 = new TH1F ("histo1", "no jets/event", 100, 0,20);
  74   h1->SetStats(kTRUE);
  75   h1->GetXaxis()->SetTitle("N_{jets}");
  76   h1->GetYaxis()->SetTitle("#frac{dN}{dN_{jets}}");
  77   h1->GetXaxis()->SetTitleColor(1);
  78   h1->SetMarkerStyle(kFullCircle);
  79   fHistos->Add(h1);
  80
  81   TH1F *h2 = new TH1F ("histo2", "no part/leading jet", 40, 0,40);
  82   h2->SetStats(kTRUE);
  83   h2->GetXaxis()->SetTitle("N_{leading}");
  84   h2->GetYaxis()->SetTitle("#frac{dN}{dN_{leading}}");
  85   h2->GetXaxis()->SetTitleColor(1);
  86   h2->SetMarkerStyle(kFullCircle);
  87   fHistos->Add(h2);
  88
  89   TProfile * h3 = new TProfile ("histo3","ProfileX of (pt,npart) of leading jets", 25, 0. ,25. , 0.,25. ) ;
  90   h3->SetStats(kTRUE);
  91   h3->GetXaxis()->SetTitle("P_{t}");
  92   h3->GetYaxis()->SetTitle("N_{leading}");
  93   h3->GetXaxis()->SetTitleColor(1);
  94   h3->SetMarkerStyle(kFullCircle);
  95   fHistos->Add(h3);
  96
  97   TH1F *h4 = new TH1F ("histo4", "Charge momentum distribution for leading jet", 500, 0 , 5);
  98   fHistos->Add(h4);
  99
 100   TProfile *h5 = new TProfile ("histo5", "ProfileX of N_{charge} vs dphi leading jet", 100 , 0. , 200. , 0 , 100 );
 101   fHistos->Add(h5);
 102
 103   TH1F *h6 = new TH1F ("histo6", " \"Transverse\" Pt Distribution ", 70, 0 , 14);
 104   fHistos->Add(h6);
 105
 106
 107 //cout << "CreateOutputObjetcs done!" << endl ;
 108 }
 109
 110 //______________________________________________________________________________
 111 void AliCdfJetFinder::FindJets()
 112   {
 113   //cout << "Begining of FindJets ..." << endl ;
 114   //1) Fill 1 array with momentum and initialise another array for indexes
 115   //2) Run algorithm
 116   //   2.1) Sort momentum array
 117   //   2.2) loop over arrays of TLorentzVectors
 118   //   2.3) flag as a possible jet
 119   //3) fill AliAODJet list
 120   Bool_t debug = kFALSE;
 121
 122   TRefArray *refs = 0;
 123   Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
 124   if (fromAod) { refs = fReader->GetReferences(); }
 125
 126   AliCdfJetHeader *header = (AliCdfJetHeader*)fHeader;
 127   Double_t radius = header->GetRadius(); // get Radius from header
 128 //  cout << "Radius is : " << radius << endl ;
 129
 130   TClonesArray * vectArray = fReader->GetMomentumArray() ;
 131   if ( vectArray == 0 ) { cout << "Could not get the momentum array" << endl; return; }
 132
 133   Int_t nPart = vectArray->GetEntries()  ; // n particles in this event;
 134   if ( nPart == 0 ) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; }
 135
 136   Int_t nJets = 0; // n jets in this event
 137
 138   fJets->SetNinput ( nPart ) ; // number of input objects
 139
 140   varContainer **vectParticle = new varContainer* [nPart]; // container for Particles
 141   varContainer **vectJet      = new varContainer* [nPart]; // container for Jets
 142
 143   Double_t *ptArray = new Double_t [nPart] ; // momentum array
 144   Int_t   *idxArray = new Int_t    [nPart] ; // index array of sorted pts
 145
 146   // initialisation of momentum and index arrays
 147   //  cout << "Filling idxArray && momArray for sorting " << endl;
 148   for ( Int_t i = 0 ; i < nPart ; i++ )
 149     {// SORTING STEP :: ptArray with data from TClonesArray of TLorentzVector
 150     TLorentzVector * lv = (TLorentzVector*) vectArray->At(i);
 151     // INITIALISATION of local arrays for temporary storage
 152     varContainer *aParticle = new varContainer;
 153     aParticle->pt   = lv->Pt();
 154     aParticle->eta  = lv->Eta();
 155     aParticle->phi  = lv->Phi();
 156     aParticle->njet = -1;
 157     vectParticle[i] = aParticle;
 158
 159     // initializing arrays
 160     idxArray [i] = 0 ;
 161      ptArray [i] = aParticle->pt ;
 162     }
 163
 164   TMath::Sort ( nPart, ptArray, idxArray ) ; // get a sorted array of indexes with size = size of TClonesArray
 165
 166   Double_t  pt_jet = 0. ,
 167            eta_jet = 0. ,
 168            phi_jet = 0. ,
 169            sum_eta = 0. ,
 170            sum_phi = 0. ;  // jet variables
 171
 172   Int_t idxPtSort = 0 ,  // index of array of sorted pt indexes
 173          npartJet = 0 ;  // number of particles in curent jet
 174
 175   TBits lkup_table ( nPart ) ;  // bit container of size npart
 176
 177   while ( lkup_table.CountBits(0) != (UInt_t)nPart )
 178     { // loop over particles in event until all flags are set
 179     UInt_t first_non_flagged = lkup_table.FirstNullBit(0) ; // set the index to the first NON flagged bit ; less conditions
 180
 181     npartJet = 0 ; // reseting number of particles in jet
 182     pt_jet   = 0.;
 183     eta_jet  = 0.;
 184     phi_jet  = 0.;
 185     sum_eta  = 0.;
 186     sum_phi  = 0.;
 187
 188     for ( UInt_t i_part = first_non_flagged ; i_part < (UInt_t)nPart ; i_part++ )
 189       {// iteration over particles in event
 190       // reseting variables
 191
 192       // the loop is done over sorted array of pt
 193       idxPtSort = idxArray[i_part] ;  // index entry of TLorentzVector from vect_arr
 194
 195       if ( lkup_table.TestBitNumber(i_part) == 1 ) { continue; } // if 4vector is already flagged skip it
 196
 197       //taking info from vectParticle ;
 198       Double_t  pt_tmp = vectParticle[idxPtSort]->pt ;
 199       Double_t eta_tmp = vectParticle[idxPtSort]->eta ;
 200       Double_t phi_tmp = vectParticle[idxPtSort]->phi ;
 201       if (debug) printf("   particle %d: pt=%g\n", i_part, pt_tmp);
 202
 203       // all angles are expressed in rad
 204
 205       if ( TMath::Abs(eta_tmp) > 0.9 )
 206          {
 207          lkup_table.SetBitNumber ( i_part, 1 ) ; // mark particle as used to be skipped
 208          continue ;
 209          }
 210
 211
 212       if ( i_part == first_non_flagged )
 213         {// this is first particle in event; leading particle
 214         // initialise the jet variables with leading particle numbers
 215
 216         npartJet++; // incrementing counter of particles in jet (nparJet == 1) (leading particle)
 217         pt_jet = pt_tmp ;
 218         eta_jet = eta_tmp ;
 219         phi_jet = phi_tmp ;
 220         if (debug) printf("  first part in jet: npartjet=%d  pt_jet=%g\n", npartJet, pt_jet);
 221
 222         lkup_table.SetBitNumber ( i_part, 1 ) ; // flag the 4vector
 223
 224         vectParticle[idxPtSort]->njet = nJets ; // associate particle with current jet
 225
 226         continue ; // skip to next particle
 227         }
 228
 229       // condition to be in jet
 230       Double_t deta = eta_jet - eta_tmp ;
 231       Double_t dphi = phi_jet - phi_tmp ;
 232
 233       if ( dphi < -TMath::Pi() ) { dphi = -dphi - 2.0 * TMath::Pi() ; }
 234       if ( dphi >  TMath::Pi() ) { dphi = -dphi + 2.0 * TMath::Pi() ; }
 235
 236       Double_t r_computed = TMath::Sqrt ( deta * deta + dphi * dphi ) ;
 237
 238       Bool_t in_jet = ( r_computed <= radius ) ? 1 : 0 ; // if r_computed is within jet_r in_jet == 1 else 0
 239
 240       if ( in_jet )
 241         { // calculus of jet variables
 242         npartJet++;  // incrementing counter of particles in jet
 243
 244         pt_jet += pt_tmp ;
 245         sum_eta += (pt_tmp * eta_tmp) ;
 246         sum_phi += (pt_tmp * phi_tmp) ;
 247         eta_jet = sum_eta / pt_jet ;
 248         phi_jet = sum_phi / pt_jet ;
 249
 250         lkup_table.SetBitNumber ( i_part, 1 ) ;  // flag the 4vector
 251         vectParticle[idxPtSort]->njet = nJets ; // setting in particle list the associated jet
 252         if (debug) printf("   particle added to jet: npartjet=%d  pt_jet=%g\n", npartJet, pt_jet);
 253         continue ; // skip to next particle
 254         }
 255
 256       }
 257       // end of iteration over event; one jet definition
 258
 259     varContainer *aJet = new varContainer;
 260     aJet->pt = pt_jet;
 261     aJet->eta = eta_jet;
 262     aJet->phi = phi_jet;
 263     aJet->njet = npartJet;
 264     vectJet[nJets++] = aJet;   // store the jet (from jet 0 to nJets-1 )
 265     if (debug) printf("=== current jet: npartjet=%d  ptjet=%g\n", npartJet,pt_jet);
 266
 267        // writing aod information
 268     Double_t px = 0., py = 0., pz = 0., en = 0.; // convert to 4-vector
 269     px = pt_jet * TMath::Cos ( phi_jet ) ;
 270     py = pt_jet * TMath::Sin ( phi_jet ) ;
 271     pz = pt_jet / TMath::Tan ( 2.0 * TMath::ATan ( TMath::Exp ( -eta_jet ) ) ) ;
 272     en = TMath::Sqrt ( px * px + py * py + pz * pz );
 273
 274     if (npartJet<2) continue;         // do not add jets with less than 2 particles
 275     if (debug) cout << "Jet 4vect : " << "px = " << px << " ; py = " << py << " ; pz = " << pz << " ; E = " << en << endl;
 276
 277     AliAODJet jet (px, py, pz, en);
 278 //    cout << "Printing jet " << endl;
 279     if (debug) jet.Print("");
 280
 281 //    cout << "Adding jet ... " ;
 282     AddJet(jet);
 283 //    cout << "added \n" << endl;
 284
 285 //     if (fromAod)
 286 //       {
 287 //       for (Int_t parts = 0; parts < nPart; parts++ )
 288 //          { if (idx_jetT[parts] == nr_jet) {jet.AddTrack(refs->At(parts));} }
 289 //       }
 290
 291
 292     }
 293     // end of while loop over particles ; ends when all particles were flagged as used and all jets defined
 294
 295   /////////////////////////////////
 296   ////   END OF EVENT PARSING   ///
 297   /////////////////////////////////
 298
 299   Int_t   *jets_pt_idx = 0;     // sorted array of jets pt
 300   Double_t    *jets_pt = 0;     // array of jets pts
 301   Int_t leading_jet_index = -1 ;   // index of leading jet from vectJet
 302   Int_t part_leadjet = 0 ;         // number of particles in leading jet
 303   Double_t   pt_leadjet = 0. ; // pt  of leading jet
 304   Double_t  eta_leadjet = 0. ; // eta of leading jet
 305   Double_t  phi_leadjet = 0. ; // phi of leading jet
 306   Int_t * idx_leadjet_part = 0;
 307
 308   if (nJets > 0)
 309     { // if there is at least one jet in event
 310     jets_pt_idx = new Int_t    [nJets] ;
 311     jets_pt     = new Double_t [nJets] ;
 312
 313     // filing the idx_ptjets array
 314     if (debug) printf("List of unsorted jets:\n");
 315     for( Int_t i = 0 ; i < nJets ; i++ )
 316       {
 317       jets_pt_idx [i] = 0 ;
 318       jets_pt [i] = vectJet[i]->pt ;
 319       if (debug) cout << "   jet found: " << i << " npartjet=" << vectJet[i]->njet << "  jets_pt[i]= " << jets_pt [i] << endl;
 320       }
 321     TMath::Sort ( nJets, jets_pt , jets_pt_idx ) ; // sorting pt of jets
 322
 323     // selection of leading jet with nr of particles > 1
 324     for( Int_t i = 0 ; i < nJets ; i++ )
 325       {
 326       if ( vectJet[ jets_pt_idx[i] ]->njet > 1 )
 327         {
 328         leading_jet_index = jets_pt_idx[i] ;
 329         part_leadjet = vectJet[ leading_jet_index ]->njet ; // number of particles in leading jet
 330           pt_leadjet = vectJet[ leading_jet_index ]->pt   ; // pt  of leading jet
 331          eta_leadjet = vectJet[ leading_jet_index ]->eta  ; // eta of leading jet
 332          phi_leadjet = vectJet[ leading_jet_index ]->phi  ; // phi of leading jet
 333
 334         cout << "Leading jet: npart = " << part_leadjet
 335              << " ; pt  = " << pt_leadjet
 336              << " ; phi = " << phi_leadjet
 337              << " ; eta = " << eta_leadjet
 338              << endl ;
 339         break ;
 340         }
 341       }
 342       // end of selection of leading jet
 343
 344    // Filling an array with indexes of leading jet particles
 345    idx_leadjet_part = new Int_t [part_leadjet] ;
 346    Int_t counter = 0;
 347    if (debug) printf("   Searching particles with jet index %d\n", leading_jet_index);
 348    for( Int_t i = 0 ; i <nPart ; i++ )
 349      {
 350      if ( vectParticle[i]->njet == leading_jet_index )
 351        {
 352        idx_leadjet_part[counter++] = i ;
 353        if (debug) cout << "   " << counter-1 << ": index=" << i << "  pt=" <<  vectParticle[i]->pt << endl;
 354        }
 355      }
 356    if ( (counter-1) > part_leadjet ) { cout << " Counter > part_leadjet !!!!" << endl;}
 357
 358    }
 359     // end of loop over DETECTED jets
 360
 361
 362
 363 // Calculus of part distribution in leading jet
 364 Double_t z = 0. ;
 365 Double_t *z_part_ljet = new Double_t [ part_leadjet ] ; // array of z of particles in leading jet
 366 for( Int_t j = 0 ; j < part_leadjet ; j++ )
 367   {
 368   Double_t z_j = vectParticle[idx_leadjet_part[j]]->pt ;
 369   z =  z_j / pt_leadjet ;
 370   z_part_ljet [j] = z ;
 371   cout << "idx_leadjet_part[j] = " << idx_leadjet_part[j]
 372       << " p of particle = " << z_j
 373       << " pt lead jet = " << pt_leadjet
 374       << " Z = " << z << endl;
 375   }
 376
 377
 378 /*
 379 Double_t dphi_partLJ = 0. ;
 380 Double_t dphiArray [nPart];
 381 for( Int_t part = 0 ; part < nPart ; part++ )
 382   {
 383   dphi_partLJ = phi_leadjet - phi_partT [part] ;
 384
 385   if ( dphi_partLJ < -TMath::Pi() ) { dphi_partLJ = -dphi_partLJ - 2.0 * TMath::Pi() ; }
 386   if ( dphi_partLJ >  TMath::Pi() ) { dphi_partLJ = -dphi_partLJ + 2.0 * TMath::Pi() ; }
 387
 388   dphi_partLJ = TMath::Abs(dphi_partLJ) ;
 389
 390   dphiArray [part] = dphi_partLJ ;
 391   }
 392 */
 393   // Fill nr_jet histogram
 394   if (fHistos)
 395     {
 396 //    printf("FILLING histograms\n");
 397     TH1F *h_jets = (TH1F*)fHistos->FindObject("histo1");
 398     if (h_jets) h_jets->Fill(nJets);
 399
 400     // There may be no leading jet with more that 2 particles !
 401     TH1F *h_leadpart = (TH1F*)fHistos->FindObject("histo2");
 402     if (part_leadjet && h_leadpart) h_leadpart->Fill(part_leadjet);
 403
 404     TProfile * h_prof = (TProfile*)fHistos->FindObject("histo3");
 405     if (part_leadjet && h_prof) h_prof->Fill(pt_leadjet,part_leadjet);
 406
 407     TH1F *h_MD = (TH1F*)fHistos->FindObject("histo4");
 408       for( Int_t k = 0  ; k < part_leadjet ; k++)
 409         { h_MD->Fill( z_part_ljet[k] ); }
 410
 411
 412 //     TProfile *h_phi = (TProfile*)fHistos->At(4);
 413 //     for( Int_t k = 0  ; k < nPart ; k++)
 414 //       { h_phi->Fill( TMath::RadToDeg() * dphiArray [k] , nPart ) ; }
 415
 416
 417 //    TH1F *h_tpd = (TH1F*)fHistos->At(5);
 418 //     for( Int_t k = 0  ; k < nPart ; k++)
 419 //       {
 420 //       if ( (phi_partT[k] > (TMath::Pi()/3.)) && (phi_partT[k] < (2. * TMath::Pi()/3.)) )
 421 //         { h_tpd->Fill( pt_partT [k] ) ; }
 422 //       }
 423
 424
 425     }
 426
 427
 428   // CLEANING SECTION
 429   for (Int_t i=0; i<nPart; i++) delete vectParticle[i];
 430   delete [] vectParticle;
 431   for (Int_t i=0; i<nJets; i++) delete vectJet[i];
 432   delete [] vectJet;
 433   delete [] ptArray;
 434   delete [] idxArray;
 435
 436   if (z_part_ljet) delete [] z_part_ljet;
 437   if (idx_leadjet_part) delete [] idx_leadjet_part;
 438   if (jets_pt_idx) delete [] jets_pt_idx ;
 439   if (jets_pt) delete [] jets_pt ;
 440
 441   vectArray->Delete() ; // TClonesArray of lorentz vectors
 442   }
 443 // end of FindJets
 444
 445
 446 //______________________________________________________________________________
 447 void AliCdfJetFinder::FinishRun()
 448 {}
 449