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
 112   for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {
 113       TParticle* trackmc = stack->Particle(iMCTracks);
 114       if (!trackmc) continue;
 115       Double_t etamc = trackmc ->Eta();
 116       Double_t pxmc  = trackmc->Px();
 117       Double_t pymc  = trackmc->Py();
 118       Double_t ptmc  = trackmc->Pt();
 119       Int_t pdgCode  = TMath::Abs(trackmc->GetPdgCode());
 120       if (TMath::Abs(etamc) > etacutoff) continue;
 121       if(ptmc < ptcutoff) continue;
 122       Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
 123       if(isprimary==0) continue;
 124       TParticlePDG* pdgPart =trackmc ->GetPDG();
 125
 126       if(chom==1){
 127           if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16)continue;
 128       } else {
 129           if (pdgPart->Charge() == 0) continue;
 130       }
 131
 132       ptT[j] = ptmc;
 133       pxT[j] = pxmc;
 134       pyT[j] = pymc;
 135       ptsuma += ptmc;
 136       pxsuma+=pxmc;
 137       pysuma+=pymc;
 138       j++;
 139   }
 140
 141   Double_t numerador = 0;
 142   Double_t numerador2 = 0;
 143   Double_t phimax = -1;
 144   Double_t pFull = -1;
 145   Double_t pMax = 0;
 146   Double_t phi = 0;
 147   Double_t thrust = 80;
 148   Double_t thrustminor = 80;
 149   Double_t nx = 0;
 150   Double_t ny = 0;
 151   Double_t phiparam = 0;
 152   //Getting thrust
 153   for(Int_t i = 0; i < 360; ++i){
 154       numerador = 0;
 155       phiparam  = 0;
 156       nx = 0;
 157       ny = 0;
 158       phiparam=((TMath::Pi()) * i) / 180; // parametrization of the angle
 159       nx = TMath::Cos(phiparam);            // x component of an unitary vector n
 160       ny = TMath::Sin(phiparam);            // y component of an unitary vector n
 161       for(Int_t i1 = 0; i1 < nmctracks; ++i1){
 162           numerador += TMath::Abs(nx * pxT[i1] + ny * pyT[i1]);//product between momentum proyection in XY plane and the unitari vector.
 163       }
 164       pFull=numerador / ptsuma;
 165       if(pFull > pMax)//maximization of pFull
 166       {
 167           pMax = pFull;
 168           phi = phiparam;
 169       }
 170   }
 171
 172   phimax=(phi * 180) / TMath::Pi();//angular parameter of the unitary vector which maximiza thrust
 173   //if n vector and beam axis form a plane, then we can calculate a second unitary vector perpendicular to that plane
 174   Double_t nx1 = TMath::Cos(phi);
 175   Double_t ny1 = TMath::Sin(phi);
 176   for(Int_t i2 =0; i2 < nmctracks; ++i2){
 177       numerador2 += TMath::Abs(pxT[i2] * ny1 - nx1 * pyT[i2]);//cross product: P_{i} X n, P_{i}=(px_{i},py_{i})
 178   }
 179   thrust = 1 - pMax;//this is the value of thrust
 180   thrustminor = numerador2 / ptsuma;//this is the value of thrust minor
 181   Double_t recoil = TMath::Abs(TMath::Sqrt(pxsuma * pxsuma + pysuma * pysuma)) / (ptsuma);//factor sentsitive to radiation outside from acceptance
 182
 183   evsh->AddAt(thrust, 0);
 184   evsh->AddAt(thrustminor, 1);
 185   evsh->AddAt(recoil, 2);
 186
 187
 188   delete evsh;
 189   delete ptT;
 190   delete pxT;
 191   delete pyT;
 192
 193   return evsh;
 194 }
 195
 196
 197 Double_t AliEventShape::GetCircularityMC(AliMCEvent* mcEvent, Int_t  nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
 198 {
 199   /*
 200     This function returns the circularity value of the event
 201
 202     The calculus uses  primary  particles. The input parameters:
 203     1. nstudymin, is the minumum number of particles which you want that participate in the calculus.(default:nstudymin =3)
 204     2. ptcutoff, is the cut in pt applied to participants to calculate the variables.(default: ptcutoff=1)
 205     3. etacutoff, is the cut in acceptance applied to participants to calculate the variables.(default: etacutoff=1)
 206     4. if chom=kTRUE, then the calculus includes neutral particles (rejecting photons and neutrinos).
 207        if chom=kFALSE, then the calculus includes only charged particles (rejecting photons and neutrinos).
 208        1>=circularity>=0
 209   */
 210
 211
 212   AliStack* stack = 0;
 213
 214   stack = mcEvent->Stack();
 215
 216   TMatrixDSym s(2);
 217
 218   Double_t s00 =  0;
 219   Double_t s01 =  0;
 220   Double_t s10 =  0;
 221   Double_t s11 =  0;
 222   Double_t ptot = 0;
 223   Double_t circularity = -2;
 224   Int_t  nmctracks = 0;
 225   Int_t nPrim  = stack->GetNprimary();
 226
 227   for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {
 228     TParticle* trackmc = stack->Particle(iMCTracks);
 229     if (!trackmc) continue;
 230     Double_t etamc = trackmc ->Eta();
 231     Double_t ptmc  = trackmc->Pt();
 232     Double_t pxmc  = trackmc->Px();
 233     Double_t pymc  = trackmc->Py();
 234     Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
 235     if (TMath::Abs(etamc) > etacutoff) continue;
 236     if (ptmc < ptcutoff) continue;
 237     Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
 238     if (isprimary == 0) continue;
 239     TParticlePDG* pdgPart = trackmc ->GetPDG();
 240     if(chom == kTRUE){
 241       // skip photons and neutrinos
 242       if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16) continue;
 243     }
 244     else{
 245       if (pdgPart->Charge() == 0)continue;
 246     }
 247
 248     ptot = ptot + (ptmc * ptmc);
 249     s00 = s00 + (pxmc * pxmc);
 250     s01 = s01 + (pxmc * pymc);
 251     s10 = s10 + (pymc * pxmc);
 252     s11 = s11 + (pymc * pymc);
 253     nmctracks++;
 254   } //track loop
 255
 256
 257   if (nmctracks < nstudymin) {
 258       Printf("Too few particles, abort");
 259       return -2;
 260   }
 261
 262
 263
 264   if(ptot != 0){
 265     s(0,0) = s00 / ptot;
 266     s(0,1) = s01 / ptot;
 267     s(1,0) = s10 / ptot;
 268     s(1,1) = s11 / ptot;
 269     const TMatrixDSymEigen eigen(s);
 270     const TVectorD eigenVal=eigen.GetEigenValues();
 271     circularity = 2 * (1 - eigenVal(0));
 272   }
 273   return circularity;
 274 }
 275
 276
 277
 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