JETAN/AliEventShape.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 /* $Id$ */
  17 //---------------------------------------------------------------------
  18 // Event shape utility class
  19 // Circularity, Thrust, ...
  20 // Authors: Antonio Ortiz Velasquez <Antonio.Ortiz.Velasquez@cern.ch>
  21 //
  22 //---------------------------------------------------------------------
  23
  24
  25 #include "AliEventShape.h"
  26
  27 #include "AliStack.h"
  28 #include "AliLog.h"
  29 #include "AliMCEventHandler.h"
  30 #include "AliMCEvent.h"
  31
  32 #include <TMatrixDSym.h>
  33 #include <TVectorD.h>
  34 #include <TMatrixDSymEigen.h>
  35 #include <TParticle.h>
  36 #include <TParticlePDG.h>
  37
  38
  39
  40
  41 //____________________________________________________________________
  42 ClassImp(AliEventShape)
  43
  44 //___________________________________________________________________
  45 TArrayD * AliEventShape::GetThrustParamMC(AliMCEvent* mcEvent, Int_t  nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
  46 {
  47   /*
  48     This function returns an array of values of thrust. To get these values you have to do:
  49     TArrayD* eventshapes = 0;
  50     eventshapes = AliShapeEvent::GetThrustParamMC(mcEvent, 3, 1, 1, kFALSE);
  51     Double_t thrust=eventshapes->GetAt(0);
  52     Double_t thrustmin=eventshapes->GetAt(1);
  53     Double_t recoil=eventshapes->GetAt(2);
  54     The calculus uses  primary  particles. The input parameters:
  55     1. nstudymin, is the minumum number of particles which you want that participate in the calculus.(default:nstudymin =3)
  56     2. ptcutoff, is the cut in pt applied to participants to calculate the variables.(default: ptcutoff=1)
  57     3. etacutoff, is the cut in acceptance applied to participants to calculate the variables.(default: etacutoff=1)
  58     4. if chom=kTRUE, then the calculus includes neutral particles (rejecting photons and neutrinos).
  59        if chom=kFALSE, then the calculus includes only charged particles (rejecting photons and neutrinos).
  60        Returned values: thrust->0: 2-jet event, thrust->0.5: isotropic event
  61        Recoil is a term which is sensitive to radiation outside from acceptance, 1>=recoil>=0,
  62        thrustmin, is a measure of the radiation which is perpendicular to the plane formed by beam axis and thrust axis, 2/TMath::Pi()>thrustmin>0. In the limit of 2 back-to-back jets thrusmin->0, while in the case of a uniformly distributed event thrustmin->2/TMath::Pi();
  63   */
  64
  65   AliStack* stack = 0;
  66
  67   stack = mcEvent->Stack();
  68   Double_t * ptT = 0;
  69   Double_t * pxT = 0;
  70   Double_t * pyT = 0;
  71   Double_t ptsuma = 0;
  72   Double_t pxsuma = 0;
  73   Double_t pysuma = 0;
  74
  75   TArrayD* evsh = new TArrayD(3);
  76   Int_t nPrim  = stack->GetNprimary();
  77   Int_t  nmctracks = 0;
  78   for (Int_t iMCTracks = 0; iMCTracks < nPrim; iMCTracks++) {
  79       TParticle* trackmc = stack->Particle(iMCTracks);
  80       if (!trackmc) continue;
  81       Double_t etamc =trackmc ->Eta();
  82       Double_t ptmc=trackmc->Pt();
  83       Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
  84       if (TMath::Abs(etamc) > etacutoff) continue; //only particles in |eta|<=etacutoff
  85       if(ptmc < ptcutoff) continue;  // PT cut
  86       Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);    // Check if particle is charged, and primary
  87       if(isprimary == 0) continue;  // only primary particles
  88       TParticlePDG* pdgPart =trackmc ->GetPDG();
  89       if(chom == 1){//include neutral particles
  90           // skip photons and neutrinos
  91           if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16) continue;
  92           nmctracks++;
  93       }
  94       else{ //only charged particles
  95           if (pdgPart->Charge() == 0)continue;
  96           nmctracks++;
  97       }
  98   }
  99   // Minimum number of particles used in the analysis
 100   if(nmctracks < nstudymin){
 101       evsh->AddAt(-2,0);
 102       evsh->AddAt(-2,1);
 103       evsh->AddAt(-2,2);
 104       return evsh;
 105   }
 106
 107   Int_t j=0;
 108   pxT = new Double_t[nmctracks];
 109   pyT = new Double_t[nmctracks];
 110   ptT = new Double_t[nmctracks];
 111   for (Int_t i = 0; i < nmctracks; i++)
 112     {
 113       pxT[i] = 0;
 114       pyT[i] = 0;
 115       ptT[i] = 0;
 116     }
 117   for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {
 118       TParticle* trackmc = stack->Particle(iMCTracks);
 119       if (!trackmc) continue;
 120       Double_t etamc = trackmc ->Eta();
 121       Double_t pxmc  = trackmc->Px();
 122       Double_t pymc  = trackmc->Py();
 123       Double_t ptmc  = trackmc->Pt();
 124       Int_t pdgCode  = TMath::Abs(trackmc->GetPdgCode());
 125       if (TMath::Abs(etamc) > etacutoff) continue;
 126       if(ptmc < ptcutoff) continue;
 127       Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
 128       if(isprimary==0) continue;
 129       TParticlePDG* pdgPart =trackmc ->GetPDG();
 130
 131       if(chom==1){
 132           if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16)continue;
 133       } else {
 134           if (pdgPart->Charge() == 0) continue;
 135       }
 136
 137       ptT[j] = ptmc;
 138       pxT[j] = pxmc;
 139       pyT[j] = pymc;
 140       ptsuma += ptmc;
 141       pxsuma+=pxmc;
 142       pysuma+=pymc;
 143       j++;
 144   }
 145
 146   Double_t numerador = 0;
 147   Double_t numerador2 = 0;
 148   Double_t phimax = -1;
 149   Double_t pFull = -1;
 150   Double_t pMax = 0;
 151   Double_t phi = 0;
 152   Double_t thrust = 80;
 153   Double_t thrustminor = 80;
 154   Double_t nx = 0;
 155   Double_t ny = 0;
 156   Double_t phiparam = 0;
 157   //Getting thrust
 158   for(Int_t i = 0; i < 360; ++i){
 159       numerador = 0;
 160       phiparam  = 0;
 161       nx = 0;
 162       ny = 0;
 163       phiparam=((TMath::Pi()) * i) / 180; // parametrization of the angle
 164       nx = TMath::Cos(phiparam);            // x component of an unitary vector n
 165       ny = TMath::Sin(phiparam);            // y component of an unitary vector n
 166       for(Int_t i1 = 0; i1 < nmctracks; ++i1){
 167           numerador += TMath::Abs(nx * pxT[i1] + ny * pyT[i1]);//product between momentum proyection in XY plane and the unitari vector.
 168       }
 169       pFull=numerador / ptsuma;
 170       if(pFull > pMax)//maximization of pFull
 171       {
 172           pMax = pFull;
 173           phi = phiparam;
 174       }
 175   }
 176
 177   phimax=(phi * 180) / TMath::Pi();//angular parameter of the unitary vector which maximiza thrust
 178   //if n vector and beam axis form a plane, then we can calculate a second unitary vector perpendicular to that plane
 179   Double_t nx1 = TMath::Cos(phi);
 180   Double_t ny1 = TMath::Sin(phi);
 181   for(Int_t i2 =0; i2 < nmctracks; ++i2){
 182       numerador2 += TMath::Abs(pxT[i2] * ny1 - nx1 * pyT[i2]);//cross product: P_{i} X n, P_{i}=(px_{i},py_{i})
 183   }
 184   thrust = 1 - pMax;//this is the value of thrust
 185   thrustminor = numerador2 / ptsuma;//this is the value of thrust minor
 186   Double_t recoil = TMath::Abs(TMath::Sqrt(pxsuma * pxsuma + pysuma * pysuma)) / (ptsuma);//factor sentsitive to radiation outside from acceptance
 187
 188   evsh->AddAt(thrust, 0);
 189   evsh->AddAt(thrustminor, 1);
 190   evsh->AddAt(recoil, 2);
 191
 192
 193   delete [] ptT;
 194   delete [] pxT;
 195   delete [] pyT;
 196
 197   return evsh;
 198 }
 199
 200
 201 Double_t AliEventShape::GetCircularityMC(AliMCEvent* mcEvent, Int_t  nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
 202 {
 203   /*
 204     This function returns the circularity value of the event
 205
 206     The calculus uses  primary  particles. The input parameters:
 207     1. nstudymin, is the minumum number of particles which you want that participate in the calculus.(default:nstudymin =3)
 208     2. ptcutoff, is the cut in pt applied to participants to calculate the variables.(default: ptcutoff=1)
 209     3. etacutoff, is the cut in acceptance applied to participants to calculate the variables.(default: etacutoff=1)
 210     4. if chom=kTRUE, then the calculus includes neutral particles (rejecting photons and neutrinos).
 211        if chom=kFALSE, then the calculus includes only charged particles (rejecting photons and neutrinos).
 212        1>=circularity>=0
 213   */
 214
 215
 216   AliStack* stack = 0;
 217
 218   stack = mcEvent->Stack();
 219
 220   TMatrixDSym s(2);
 221
 222   Double_t s00 =  0;
 223   Double_t s01 =  0;
 224   Double_t s10 =  0;
 225   Double_t s11 =  0;
 226   Double_t ptot = 0;
 227   Double_t circularity = -2;
 228   Int_t  nmctracks = 0;
 229   Int_t nPrim  = stack->GetNprimary();
 230
 231   for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {
 232     TParticle* trackmc = stack->Particle(iMCTracks);
 233     if (!trackmc) continue;
 234     Double_t etamc = trackmc ->Eta();
 235     Double_t ptmc  = trackmc->Pt();
 236     Double_t pxmc  = trackmc->Px();
 237     Double_t pymc  = trackmc->Py();
 238     Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
 239     if (TMath::Abs(etamc) > etacutoff) continue;
 240     if (ptmc < ptcutoff) continue;
 241     Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
 242     if (isprimary == 0) continue;
 243     TParticlePDG* pdgPart = trackmc ->GetPDG();
 244     if(chom == kTRUE){
 245       // skip photons and neutrinos
 246       if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16) continue;
 247     }
 248     else{
 249       if (pdgPart->Charge() == 0)continue;
 250     }
 251
 252     ptot = ptot + (ptmc * ptmc);
 253     s00 = s00 + (pxmc * pxmc);
 254     s01 = s01 + (pxmc * pymc);
 255     s10 = s10 + (pymc * pxmc);
 256     s11 = s11 + (pymc * pymc);
 257     nmctracks++;
 258   } //track loop
 259
 260
 261   if (nmctracks < nstudymin) {
 262       Printf("Too few particles, stopping");
 263       return -2;
 264   }
 265
 266
 267
 268   if(ptot != 0){
 269     s(0,0) = s00 / ptot;
 270     s(0,1) = s01 / ptot;
 271     s(1,0) = s10 / ptot;
 272     s(1,1) = s11 / ptot;
 273     const TMatrixDSymEigen eigen(s);
 274     const TVectorD eigenVal=eigen.GetEigenValues();
 275     circularity = 2 * (1 - eigenVal(0));
 276   }
 277   return circularity;
 278 }
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311