PWGLF/SPECTRA/MultEvShape/AliAnalysisCentralCutEvtMC.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 //  MC event level cuts for azimuthal isotropic
  18 //  expansion in highly central collisions analysis
  19 //  author: Cristian ANDREI
  20 //          acristian@niham.nipne.ro
  21 // ----------------------------------------------------
  22
  23
  24 #include <TParticle.h>
  25
  26 #include "AliMCEvent.h"
  27 #include "AliStack.h"
  28
  29 #include "AliAnalysisCentralCutEvtMC.h"
  30
  31
  32 class TObject;
  33
  34 //____________________________________________________________________
  35 ClassImp(AliAnalysisCentralCutEvtMC)
  36
  37 //____________________________________________________________________
  38 AliAnalysisCentralCutEvtMC::AliAnalysisCentralCutEvtMC(const Char_t* name, const Char_t* title)
  39     :AliAnalysisCuts(name,title)
  40     ,fReqMult(kFALSE)
  41     ,fReqDir(kFALSE)
  42     ,fReqDirUnit(kFALSE)
  43     ,fMultMin(0)
  44     ,fMultMax(0)
  45     ,fDirMin(0)
  46     ,fDirMax(0)
  47     ,fDirUMin(0)
  48     ,fDirUMax(0)
  49
  50 {
  51     //constructor
  52
  53 }
  54
  55 AliAnalysisCentralCutEvtMC::~AliAnalysisCentralCutEvtMC(){
  56 //destructor
  57
  58 }
  59
  60 //___________________________________________________________________________
  61 Bool_t AliAnalysisCentralCutEvtMC::IsSelected(TObject *obj){
  62 // Checks whether the event passes the cuts
  63
  64     if (!obj){
  65         printf("Pointer to MC Event is null! \n");
  66         return kFALSE;
  67     }
  68
  69     AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*> (obj);
  70                 if(!mcEvent){
  71                 printf("AliAnalysisCentralCutEvtMC -> Can't get MCEvt!\n");
  72                 return kFALSE;
  73     }
  74
  75
  76     if(fReqMult){
  77                 Int_t mult = CalcMult(mcEvent);
  78                 if((mult<fMultMin)||(mult>fMultMax)){
  79                         return kFALSE;
  80                 }
  81     }
  82
  83     if(fReqDir){
  84                 Double_t dir = CalcDir(mcEvent);
  85                 if((dir<fDirMin)||(dir>fDirMax)){
  86                         return kFALSE;
  87                 }
  88     }
  89
  90     if(fReqDirUnit){
  91                 Double_t dirU = CalcDirUnit(mcEvent);
  92                 if((dirU<fDirUMin)||(dirU>fDirUMax)){
  93                         return kFALSE;
  94                 }
  95     }
  96
  97     return kTRUE;
  98 }
  99
 100 //___________________________________________________________________________
 101 Int_t AliAnalysisCentralCutEvtMC::CalcMult(AliMCEvent* const mcEvent) {
 102 // Computes the multiplicity of the event
 103
 104     AliStack *stack=mcEvent->Stack();
 105     if(!stack) return -1;
 106
 107     Int_t nPrim = stack->GetNprimary();
 108
 109     Int_t charged = 0;
 110     Double_t eta;
 111
 112     for(Int_t i = 0; i < nPrim; i++){//track loop -> compute multiplicity
 113
 114                 TParticle *particle=stack->Particle(i);
 115                 if(!particle){
 116                         continue;
 117                 }
 118
 119                 eta = particle->Eta();
 120
 121                 if((eta > 0.5)||(eta < -0.5)){         //------------mid-rapidity
 122                         continue;
 123                 }
 124
 125                 TParticlePDG *particlePDG = particle->GetPDG(0);
 126                 if(!particlePDG){
 127                         continue;
 128                 }
 129
 130                 if(TMath::Abs(particlePDG->Charge()) < 3){
 131                         continue;
 132                 }
 133
 134                 charged++;
 135
 136     }//end track loop
 137
 138     return charged;
 139
 140 }
 141
 142 //____________________________________________________________________________
 143 Double_t AliAnalysisCentralCutEvtMC::CalcDir(AliMCEvent* const mcEvent) {
 144 // computes the directivity of the event
 145
 146     AliStack *stack=mcEvent->Stack();
 147     if(!stack) return -1;
 148
 149     Int_t nPrim = stack->GetNprimary();
 150
 151     Int_t goodtrack = 0;
 152
 153     Double_t eta;
 154     Double_t pt;
 155
 156     Double_t dir;
 157     Double_t px,py;
 158     Double_t sumaPt = 0;
 159     Double_t sumaPx = 0;
 160     Double_t sumaPy = 0;
 161
 162     for(Int_t i = 0; i < nPrim; i++){//track loop -> compute directivity
 163
 164                 TParticle *particle=stack->Particle(i);
 165                 if(!particle){
 166                         continue;
 167                 }
 168
 169                 eta = particle->Eta();
 170
 171                 if((eta > 1.9)||(eta < 0.0)){   // half of the SPD coverage -> directivity
 172                         continue;
 173                 }
 174
 175                 TParticlePDG *particlePDG = particle->GetPDG(0);
 176                 if(!particlePDG){
 177                         continue;
 178                 }
 179
 180                 if(TMath::Abs(particlePDG->Charge()) < 3){
 181                         continue;
 182                 }
 183
 184                 px = particle->Px();
 185                 py = particle->Py();
 186
 187                 pt = particle->Pt();
 188
 189                 sumaPx = sumaPx + px;
 190                 sumaPy = sumaPy + py;
 191
 192                 sumaPt = sumaPt + pt;
 193
 194                 goodtrack++;
 195
 196     }//end track loop
 197
 198
 199     if(sumaPt < 0.0000001){
 200                 return -1;
 201     }
 202
 203     dir = (sqrt(pow(sumaPx,2)+pow(sumaPy,2)))/sumaPt;
 204
 205     return dir;
 206 }
 207
 208 //__________________________________________________________________
 209 Double_t AliAnalysisCentralCutEvtMC::CalcDirUnit(AliMCEvent* const mcEvent) {
 210 // computes the directivity using only the SPD unity vectors
 211
 212     AliStack *stack=mcEvent->Stack();
 213     if(!stack){
 214                 return -1;
 215         }
 216
 217     Int_t nPrim = stack->GetNprimary();
 218
 219     Int_t goodtrack = 0;
 220
 221     Double_t eta;
 222     Double_t pt;
 223
 224     Double_t dirU;
 225     Double_t px,py,pxU,pyU;
 226     Double_t sumaPxU = 0;
 227     Double_t sumaPyU = 0;
 228
 229
 230     for(Int_t i = 0; i < nPrim; i++){//track loop -> compute directivity
 231
 232                 TParticle *particle=stack->Particle(i);
 233                 if(!particle){
 234                         continue;
 235                 }
 236
 237                 eta = particle->Eta();
 238
 239                 if((eta > 1.9)||(eta < 0.0)){ // half of the SPD coverage -> directivity
 240                         continue;
 241                 }
 242
 243                 TParticlePDG *particlePDG = particle->GetPDG(0);
 244                 if(!particlePDG){
 245                         continue;
 246                 }
 247
 248                 if(TMath::Abs(particlePDG->Charge()) < 3){
 249                         continue;
 250                 }
 251
 252                 px = particle->Px();
 253                 py = particle->Py();
 254                 pt = particle->Pt();
 255
 256                 pxU = px/pt;
 257                 pyU = py/pt;
 258
 259                 sumaPxU = sumaPxU + pxU;
 260                 sumaPyU = sumaPyU + pyU;
 261
 262                 goodtrack++;
 263
 264     }//end track loop
 265
 266     if(goodtrack == 0){
 267                 return -1;
 268         }
 269
 270     dirU = (sqrt(pow(sumaPxU,2)+pow(sumaPyU,2)))/goodtrack;
 271
 272
 273     return dirU;
 274 }