PWG2/SPECTRA/AliAnalysisCentralCutEvtESD.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 //  ESD 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 #include "TMath.h"
  24 #include <TObjArray.h>
  25 #include <TList.h>
  26
  27 #include "AliESDEvent.h"
  28 #include "AliESDtrack.h"
  29 #include "AliESDtrackCuts.h"
  30  #include "AliMultiplicity.h"
  31
  32 #include "AliAnalysisCentralCutESD.h"
  33 #include "AliAnalysisCentralCutEvtESD.h"
  34
  35 class TObject;
  36
  37
  38 //____________________________________________________________________
  39 ClassImp(AliAnalysisCentralCutEvtESD)
  40
  41 //____________________________________________________________________
  42 AliAnalysisCentralCutEvtESD::AliAnalysisCentralCutEvtESD(const Char_t* name, const Char_t* title)
  43     :AliAnalysisCuts(name,title)
  44     ,fReqMult(kFALSE)
  45     ,fReqDir(kFALSE)
  46     ,fReqSPDMult(kFALSE)
  47     ,fReqSPDDir(kFALSE)
  48     ,fMultMin(0)
  49     ,fMultMax(0)
  50     ,fDirMin(0)
  51     ,fDirMax(0)
  52     ,fSPDMultMin(0)
  53     ,fSPDMultMax(0)
  54     ,fSPDDirMin(0)
  55     ,fSPDDirMax(0)
  56
  57 {
  58 //constructor
  59     for(Int_t i=0; i<10; i++){
  60                 fCutsList[i] = 0;
  61     }
  62
  63         InitCuts();
  64         if (!fCutsList) {
  65                 Printf("ERROR: fCutsList not available");
  66                 return;
  67     }
  68
  69 }
  70
  71 AliAnalysisCentralCutEvtESD::~AliAnalysisCentralCutEvtESD() {
  72 //Destructor
  73         if(fCutsList){
  74                 delete [] fCutsList;
  75         }
  76
  77 }
  78
  79
  80
  81 Bool_t AliAnalysisCentralCutEvtESD::IsSelected(TObject *obj){
  82 // check whether the event passes the cuts
  83     AliESDEvent *esdEvent = dynamic_cast<AliESDEvent *>(obj);
  84     if(!esdEvent){
  85                 printf("AliAnalysisCentralCutEvtESD:IsSelected ->Can't get ESD event!\n");
  86                 return kFALSE;
  87     }
  88
  89         if(fReqMult){
  90                 Int_t mult = CalcMult(esdEvent);
  91                 if((mult<fMultMin)||(mult>fMultMax)){
  92                         return kFALSE;
  93                 }
  94     }
  95
  96     if(fReqDir){
  97                 Double_t dir = CalcDir(esdEvent);
  98                 if((dir<fDirMin)||(dir>fDirMax)){
  99                         return kFALSE;
 100                 }
 101     }
 102
 103     if(fReqSPDMult){
 104                 Double_t spdMult = CalcSPDMult(esdEvent);
 105                 if((spdMult<fSPDMultMin)||(spdMult>fSPDMultMax)){
 106                         return kFALSE;
 107                 }
 108     }
 109
 110     if(fReqSPDDir){
 111                 Double_t spdDir = CalcSPDDir(esdEvent);
 112                 if((spdDir<fSPDDirMin)||(spdDir>fSPDDirMax)){
 113                         return kFALSE;
 114                 }
 115     }
 116
 117     return kTRUE;
 118 }
 119
 120 //_________________________________________________________________________
 121 void AliAnalysisCentralCutEvtESD::InitCuts(){
 122 //Initialize internal cuts
 123
 124 ////////////////  FULL ALICE ///////////////////
 125 //------------General ESD Cuts------------------
 126     AliESDtrackCuts *esdCutsGen = new AliESDtrackCuts("AliESDtrackCuts", "CutEvtESDInternal");
 127     esdCutsGen->SetMinNClustersTPC(50);
 128     esdCutsGen->SetMaxChi2PerClusterTPC(2.2);
 129     esdCutsGen->SetMaxCovDiagonalElements(0.5,0.5,0.5,0.5,0.5);
 130     esdCutsGen->SetRequireTPCRefit(kTRUE);
 131     esdCutsGen->SetAcceptKinkDaughters(kFALSE);
 132     esdCutsGen->SetMaxNsigmaToVertex(2.0);
 133     esdCutsGen->SetRequireSigmaToVertex(kTRUE);
 134
 135     AliAnalysisCentralCutESD *esdCutsGen1 = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
 136     esdCutsGen1->SetReqIsCharged();
 137
 138 //-------------Specific ESD Cuts------------------
 139     AliESDtrackCuts *esdCutsMult = new AliESDtrackCuts("AliESDCutsMult", "midrapidity");
 140     esdCutsMult->SetEtaRange(-0.5,0.5);
 141     AliESDtrackCuts *esdCutsDir = new AliESDtrackCuts("AliESDCutsDir", "SPD coverage");
 142     esdCutsDir->SetEtaRange(0.0,1.9);
 143
 144 ///////////////////  SPD ONLY  ////////////////////
 145     AliESDtrackCuts *esdCutsSPD = new AliESDtrackCuts("AliESDtrackCuts", "CutEvtESDInternal");
 146     esdCutsSPD->SetAcceptKinkDaughters(kFALSE);
 147     esdCutsSPD->SetMaxNsigmaToVertex(2.0);
 148     esdCutsSPD->SetRequireSigmaToVertex(kTRUE);
 149
 150 //--------------set the cuts ----------------------
 151     TObjArray* esdListMult = new TObjArray();
 152     esdListMult->AddLast(esdCutsGen);
 153     esdListMult->AddLast(esdCutsGen1);
 154     esdListMult->AddLast(esdCutsMult);
 155
 156     TObjArray* esdListDir = new TObjArray();
 157     esdListDir->AddLast(esdCutsGen);
 158     esdListDir->AddLast(esdCutsGen1);
 159     esdListDir->AddLast(esdCutsDir);
 160
 161     TObjArray* esdListSPD = new TObjArray();
 162     esdListSPD->AddLast(esdCutsSPD);
 163
 164     fCutsList[0]=esdListDir;
 165     fCutsList[1]=esdListMult;
 166     fCutsList[2]=esdListSPD;
 167
 168 }
 169
 170
 171 //__________________________________________________________________________
 172 Bool_t AliAnalysisCentralCutEvtESD::CheckIntCuts(Int_t no, TObject *obj) const{
 173 // Check if the particle passes the internal cuts
 174
 175     if(no > 9){
 176                 printf("\n AliAnalysisCentralCutEvtESD::CheckIntCuts -> Cut number is not ok! \n");
 177                 return kFALSE;
 178     }
 179
 180     if(!fCutsList[no]){
 181                 printf("AliAnalysisCentralCutEvtESD::CheckIntCuts -> cuts list problem! \n");
 182                 return kFALSE;
 183     }
 184
 185
 186     TObjArrayIter iter(fCutsList[no]);
 187     AliAnalysisCuts *cut = 0;
 188
 189
 190     while((cut = (AliAnalysisCuts*)iter.Next())){
 191
 192                 if(!cut->IsSelected(obj)) return kFALSE;
 193     }
 194
 195     return kTRUE;
 196 }
 197
 198
 199 //__________________________________________________________________________
 200 Double_t AliAnalysisCentralCutEvtESD::CalcDir(AliESDEvent* const esdEv) {
 201
 202 //Compute the directivity - FULL ALICE
 203
 204     Double_t dir;
 205     Double_t px,py;
 206     Double_t sumaPt = 0;
 207     Double_t sumaPx = 0;
 208     Double_t sumaPy = 0;
 209
 210     Double_t pt;
 211
 212     if (!esdEv){
 213                 printf("NULL tree\n");
 214                 return -1;
 215     }
 216
 217     Int_t totTracks=esdEv->GetNumberOfTracks();
 218
 219     for(Int_t itrack = 0; itrack < totTracks; itrack++){//first loop->compute event directivity
 220
 221                 AliESDtrack* track = esdEv->GetTrack(itrack);
 222
 223                 if (!track) {
 224                                 Printf("ERROR: Could not receive track %d", itrack);
 225                         continue;
 226                 }
 227
 228                 if(!CheckIntCuts(0, track)) continue;
 229
 230                 px = track->Px();
 231                 py = track->Py();
 232                 pt = track->Pt();
 233
 234                 sumaPx = sumaPx + px;
 235                 sumaPy = sumaPy + py;
 236
 237                 sumaPt = sumaPt + pt;
 238
 239     }//end track loop
 240
 241
 242         if(sumaPt < 0.0000001){
 243                 return -1;
 244     }
 245
 246     dir = (sqrt(pow(sumaPx,2)+pow(sumaPy,2)))/sumaPt;
 247
 248     return dir;
 249 }
 250
 251 //__________________________________________________________________________
 252 Int_t AliAnalysisCentralCutEvtESD::CalcMult(AliESDEvent* const esdEv) {
 253
 254 //Compute multiplicity - FULL ALICE
 255
 256     Int_t charged = 0;
 257
 258     if (!esdEv){
 259                 printf("NULL tree\n");
 260                 return -1;
 261     }
 262
 263     Int_t totTracks=esdEv->GetNumberOfTracks();
 264
 265     for(Int_t iTrack = 0; iTrack < totTracks; iTrack++){//second track loop -> compute event multiplicity
 266
 267                 AliESDtrack* track = esdEv->GetTrack(iTrack);
 268
 269                 if (!track) {
 270                         Printf("ERROR: Could not receive track %d", iTrack);
 271                         continue;
 272                 }
 273
 274                 if(!CheckIntCuts(1, track)) continue;
 275
 276
 277                 charged++; //multiplicity
 278
 279     }//end second track loop
 280
 281
 282     return charged;
 283 }
 284
 285
 286 //__________________________________________________________________________
 287 Double_t AliAnalysisCentralCutEvtESD::CalcSPDDir(AliESDEvent* const esdEv) {
 288
 289 //Compute directivity - SPD ONLY
 290
 291     Double_t dirU;
 292     Double_t pxU,pyU;
 293     Double_t sumaPxU = 0;
 294     Double_t sumaPyU = 0;
 295
 296     Double_t goodtrack = 0;
 297     Double_t spdEta = 0.;
 298
 299     const AliMultiplicity *spdMult=esdEv->GetMultiplicity();
 300     if(!spdMult){
 301         printf("Unable to get multiplicity! \n");
 302         return -1;
 303     }
 304
 305     Int_t spdTracks = spdMult->GetNumberOfTracklets();  //SPD multiplicity
 306
 307     for(Int_t iTrack = 0; iTrack< spdTracks; iTrack++){  //SPD track loop -> Directivity
 308
 309                 AliESDtrack* track = esdEv->GetTrack(iTrack);
 310
 311                 if (!track) {
 312                                 Printf("ERROR: Could not receive track %d", iTrack);
 313                                 continue;
 314                 }
 315
 316                 if(!CheckIntCuts(2, track)) continue;
 317
 318                 spdEta = spdMult->GetEta(iTrack);
 319
 320                 if((spdEta<0.0)||(spdEta>1.9)) continue;
 321
 322                 Double_t phi = spdMult->GetPhi(iTrack);
 323
 324                 pxU = TMath::Cos(phi);
 325                 pyU = TMath::Sin(phi);
 326
 327                 sumaPxU = sumaPxU + pxU;
 328                 sumaPyU = sumaPyU + pyU;
 329
 330                 goodtrack++;
 331
 332     }//end SPD track loop
 333
 334     if(goodtrack < 1.) return -1;
 335
 336     dirU = (sqrt(pow(sumaPxU,2)+pow(sumaPyU,2)))/goodtrack;
 337
 338     return dirU;
 339 }
 340
 341 //__________________________________________________________________________
 342
 343 Int_t AliAnalysisCentralCutEvtESD::CalcSPDMult(AliESDEvent* const esdEv) {
 344
 345         //Compute multiplicity - SPD ONLY
 346
 347     const AliMultiplicity *spdMult=esdEv->GetMultiplicity();
 348     if(!spdMult){
 349         printf("Unable to get multiplicity! \n");
 350         return -1;
 351     }
 352
 353     Int_t spdTracks = spdMult->GetNumberOfTracklets();  //SPD multiplicity
 354
 355     Double_t spdEta;
 356     Int_t charged = 0;
 357
 358     for(Int_t iTrack = 0; iTrack< spdTracks; iTrack++){  //second SPD track loop -> Multiplicity
 359
 360                 AliESDtrack* track = esdEv->GetTrack(iTrack);
 361
 362                 if (!track) {
 363                                 Printf("ERROR: Could not receive track %d", iTrack);
 364                                 continue;
 365                 }
 366
 367                 if(!CheckIntCuts(2, track)) continue;
 368
 369                 spdEta = spdMult->GetEta(iTrack);
 370
 371                 if((spdEta<-0.5)||(spdEta>0.5)) continue;
 372
 373                 charged++;
 374
 375     }//end second SPD track loop
 376
 377     return charged;
 378 }