PWGUD/DIFFRACTIVE/example/AliCDMesonUtilsStripped.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 // AliCDMesonUtilsStripped
  17 //
  18 //  Author:
  19 //  Felix Reidt <Felix.Reidt@cern.ch>
  20
  21 #include <TH1.h>
  22 #include <TH2.h>
  23 #include <TGeoMatrix.h>
  24 #include <THnSparse.h>
  25 #include <TLorentzVector.h>
  26 #include <TRandom3.h>
  27 #include <TParticle.h>
  28
  29 #include "AliCDBManager.h"
  30 #include "AliCDBEntry.h"
  31 #include "AliCDBStorage.h"
  32 #include "AliESDEvent.h"
  33 #include "AliPIDResponse.h"
  34 #include "AliVTrack.h"
  35 #include "AliVParticle.h"
  36 #include "AliESDtrack.h"
  37 #include "AliESDFMD.h"
  38 #include "AliGeomManager.h"
  39 #include "AliITSAlignMille2Module.h"
  40 #include "AliITSsegmentationSPD.h"
  41 #include "AliMultiplicity.h"
  42 #include "AliPIDResponse.h"
  43 #include "AliSPDUtils.h"
  44 #include "AliTriggerAnalysis.h"
  45
  46 #include "AliAODTracklets.h"
  47 #include "AliAODEvent.h"
  48
  49 #include "AliCDMesonBaseStripped.h"
  50 #include "AliCDMesonTracks.h"
  51
  52 #include "AliCDMesonUtilsStripped.h"
  53
  54
  55 //==============================================================================
  56 //------------------------------------------------------------------------------
  57 Bool_t AliCDMesonUtilsStripped::CutEvent(const AliESDEvent *ESDEvent)
  58 {
  59         //
  60         // CutEvent
  61         //
  62
  63         AliTriggerAnalysis triggerAnalysis;
  64
  65         // collision vertex cut
  66         // A cut in XY is implicitly done during the reconstruction by constraining
  67         // the vertex to the beam diamond.
  68
  69         // Primary vertex
  70         Bool_t kpr0 = kTRUE;
  71         const AliESDVertex *vertex = ESDEvent->GetPrimaryVertexTracks();
  72         if(vertex->GetNContributors()<1) {
  73                 // SPD vertex
  74                 vertex = ESDEvent->GetPrimaryVertexSPD();
  75                 if(vertex->GetNContributors()<1) {
  76                         // NO GOOD VERTEX, SKIP EVENT
  77                         kpr0 = kFALSE;
  78                 }
  79         }
  80         const Bool_t kpriv = kpr0 && (fabs(ESDEvent->GetPrimaryVertex()->GetZ())<10.);
  81         // 10 is the common value, unit: cm
  82         if(!kpriv)
  83                 return kFALSE;
  84
  85         return kTRUE;
  86 }
  87
  88 //------------------------------------------------------------------------------
  89 Bool_t AliCDMesonUtilsStripped::CutEvent(const AliAODEvent *AODEvent)
  90 {
  91         //
  92         // Cut Event for AOD Events, to be combined with the ESD Track Cut
  93         //
  94
  95         // TODO: no idea about fast or yet, to be thought of
  96
  97         // Primary vertex
  98         Bool_t kpr0 = kTRUE;
  99         const AliAODVertex *vertex = AODEvent->GetPrimaryVertex();
 100         if(vertex->GetNContributors()<1) {
 101                 // SPD vertex
 102                 vertex = AODEvent->GetPrimaryVertexSPD();
 103                 if(vertex->GetNContributors()<1) {
 104                         // NO GOOD VERTEX, SKIP EVENT
 105                         kpr0 = kFALSE;
 106                 }
 107         }
 108         const Bool_t kpriv = kpr0 && (fabs(AODEvent->GetPrimaryVertex()->GetZ())<10.);
 109         // 10 is the common value, unit: cm
 110
 111         if(!kpriv)
 112                 return kFALSE;
 113
 114         return kTRUE;
 115 }
 116
 117
 118 //------------------------------------------------------------------------------
 119 Int_t AliCDMesonUtilsStripped::GetGapConfig(const AliESDEvent *ESDEvent)
 120 {
 121         //
 122         // GetGapConfigAndTracks
 123         //
 124         // retrieves the gap configuration of a track and returns it as
 125         // an bit vector
 126         // kBaseLine ensures, that this event is valid
 127         // + is equivalent to | in this case
 128         return AliCDMesonBaseStripped::kBitBaseLine
 129                 + GetV0(ESDEvent) + GetFMD(ESDEvent) + GetSPD(ESDEvent) + GetTPC(ESDEvent);
 130 }
 131
 132
 133 //==============================================================================
 134 //------------------------------------------------------------------------------
 135 Int_t AliCDMesonUtilsStripped::GetV0(const AliESDEvent * ESDEvent)
 136 {
 137         //
 138         //GetV0
 139         //
 140
 141         AliTriggerAnalysis triggerAnalysis;
 142         const Bool_t khw = kFALSE;
 143         const Bool_t v0A =
 144                 (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kASide, khw) ==
 145                  AliTriggerAnalysis::kV0BB);
 146         const Bool_t v0C =
 147                 (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kCSide, khw) ==
 148                  AliTriggerAnalysis::kV0BB);
 149
 150         return v0A * AliCDMesonBaseStripped::kBitV0A
 151                 + v0C * AliCDMesonBaseStripped::kBitV0C;
 152 }
 153
 154
 155 //------------------------------------------------------------------------------
 156 Int_t AliCDMesonUtilsStripped::GetFMD(const AliESDEvent *ESDEvent)
 157 {
 158         //
 159         // GetFMD
 160         //
 161
 162         AliTriggerAnalysis triggerAnalysis;
 163         triggerAnalysis.SetFMDThreshold(0.3, 0.5); // parameters got from FMD
 164         const Bool_t fmdA =
 165                 triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kASide);
 166         const Bool_t fmdC =
 167                 triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kCSide);
 168
 169         return fmdA * AliCDMesonBaseStripped::kBitFMDA
 170                 + fmdC * AliCDMesonBaseStripped::kBitFMDC;
 171 }
 172
 173
 174 //------------------------------------------------------------------------------
 175 Int_t AliCDMesonUtilsStripped::GetSPD(const AliESDEvent *ESDEvent)
 176 {
 177         //
 178         // GetSPD
 179         //
 180
 181         Int_t nfoctr[10];
 182         GetNFO(ESDEvent, "]0.9[", nfoctr);
 183         // get multiplicity from fastOR and fill corresponding hit maps
 184
 185
 186         const Int_t ipA = nfoctr[kIPA]; // inner layer A side
 187         const Int_t ipC = nfoctr[kIPC]; // inner layer C side
 188         const Int_t opA = nfoctr[kOPA]; // outer layer A side
 189         const Int_t opC = nfoctr[kOPC]; // outer layer C side
 190
 191         const Bool_t spdA = ipA + opA; // A side hit?
 192         const Bool_t spdC = ipC + opC; // C side hit?
 193
 194         return spdA * AliCDMesonBaseStripped::kBitSPDA
 195                 + spdC * AliCDMesonBaseStripped::kBitSPDC;
 196 }
 197
 198
 199 //------------------------------------------------------------------------------
 200 Int_t AliCDMesonUtilsStripped::GetTPC(const AliESDEvent * ESDEvent)
 201 {
 202         //
 203         //GetTPC
 204         //
 205
 206         const Double_t etacut = 0.9;
 207         Int_t nA = 0;
 208         Int_t nC = 0;
 209         for(Int_t itrack = 0; itrack < ESDEvent->GetNumberOfTracks(); itrack++){
 210                 const AliESDtrack* esdtrack = ESDEvent->GetTrack(itrack);
 211                 if( esdtrack->Eta() > etacut ){
 212                         nA ++;
 213                 }
 214                 else if( esdtrack->Eta() < -etacut ){
 215                         nC ++;
 216                 }
 217         }
 218
 219         const Bool_t tpcA = nA;
 220         const Bool_t tpcC = nC;
 221
 222         return tpcA * AliCDMesonBaseStripped::kBitTPCA
 223                 + tpcC * AliCDMesonBaseStripped::kBitTPCC;
 224 }
 225
 226 //==============================================================================
 227 //------------------------------------------------------------------------------
 228 void AliCDMesonUtilsStripped::SPDLoadGeom(const Int_t run)
 229 {
 230         // method to get the gGeomanager
 231         // it is called at the CreatedOutputObject stage
 232         // to comply with the CAF environment
 233
 234         AliCDBManager *man = AliCDBManager::Instance();
 235         // WARNING THE OCDB PATH SHOULD BE ADJUSTED TO THE RUNNING CONDITIONS
 236
 237
 238         TString cdbpath;
 239         if (man->IsDefaultStorageSet()) {
 240                 const AliCDBStorage *dsto = man->GetDefaultStorage();
 241                 cdbpath = TString(dsto->GetBaseFolder());
 242         }
 243         else {
 244                 man->SetDefaultStorage("raw://");
 245                 cdbpath = "raw://";
 246         }
 247
 248         man->SetSpecificStorage("ITS/Align/Data",cdbpath);
 249         man->SetSpecificStorage("GRP/Geometry/Data",cdbpath);
 250         man->SetRun(run);
 251
 252         AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
 253         if (!obj) {
 254                 printf("freidtlog failed loading geometry object\n");
 255                 return;
 256         }
 257         AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
 258         AliGeomManager::ApplyAlignObjsFromCDB("ITS");
 259 }
 260
 261 //------------------------------------------------------------------------------
 262 Bool_t AliCDMesonUtilsStripped::SPDLoc2Glo(const Int_t id, const Double_t *loc,
 263                                            Double_t *glo)
 264 {
 265         //
 266         //SPDLoc2Glo, do not touch
 267         //
 268
 269         static TGeoHMatrix mat;
 270         Int_t vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
 271         if (vid<0) {
 272                 printf("freidtlog Did not find module with such ID %d\n",id);
 273                 return kFALSE;
 274         }
 275         AliITSAlignMille2Module::SensVolMatrix(vid,&mat);
 276         mat.LocalToMaster(loc,glo);
 277         return kTRUE;
 278 }
 279
 280
 281 //------------------------------------------------------------------------------
 282 Int_t AliCDMesonUtilsStripped::CheckChipEta(const Int_t chipKey,
 283                                             const TString scut,
 284                                             const Double_t vtxPos[])
 285 {
 286         //
 287         //CheckChipEta
 288         //
 289
 290         // retrieves the position in eta for a given chip and applies the cut
 291         // results:
 292         // 0 <= out of range
 293         // -1 <= negative pseudo-rapidity position, in range (C-Side)
 294         // 1 <= positive pseudo-rapidity position, in range (A-Side)
 295         //
 296         // scut: "[0.9" or "]0.9", only 3 digits for the value!!
 297
 298
 299         const Bool_t kincl = (scut[0] == '[');
 300         const TString cutval = scut(1,3);
 301         const Double_t etacut = fabs(cutval.Atof());
 302
 303         //no eta cut, save time
 304         if(kincl && etacut>=2)
 305                 return kTRUE;
 306
 307         Int_t etaside = 1;
 308         //------------------------------- NOT TO TOUCH ------------------------>>
 309         UInt_t module=999, offchip=999;
 310         AliSPDUtils::GetOfflineFromOfflineChipKey(chipKey,module,offchip);
 311         UInt_t hs = AliSPDUtils::GetOnlineHSFromOffline(module);
 312         if(hs<2) offchip = 4 - offchip; // inversion  in the inner layer...
 313
 314         const Int_t col[]={
 315                 hs<2? 0 : 31,
 316                 hs<2? 31 : 0,
 317                 hs<2? 31 : 0,
 318                 hs<2? 0 : 31};
 319         const Int_t aa[]={0, 0, 255, 255};
 320         const AliITSsegmentationSPD seg;
 321
 322         for(Int_t ic=0; ic<4; ic++){
 323                 Float_t localchip[3]={0.,0.,0.};
 324                 seg.DetToLocal(aa[ic],col[ic]+32*offchip,localchip[0],localchip[2]);
 325                 // local coordinate of the chip center
 326                 //printf("local coordinates %d %d: %f %f \n",chipKey, ic, localchip[0],localchip[2]);
 327                 const Double_t local[3] = {localchip[0],localchip[1],localchip[2]};
 328                 Double_t glochip[3]={0.,0.,0.};
 329                 if(!SPDLoc2Glo(module,local,glochip)){
 330                         return kFALSE;
 331                 }
 332
 333                 //-------------------------------------------------------------------<<
 334
 335                 const TVector3 pos(glochip[0]-vtxPos[0], glochip[1]-vtxPos[1],
 336                                    glochip[2]-vtxPos[2]);
 337                 //pos.Print();
 338
 339                 if( kincl && fabs(pos.Eta()) > etacut)
 340                         return kFALSE;
 341
 342                 if(!kincl){
 343                         if(fabs(pos.Eta()) < etacut)
 344                                 return kFALSE;
 345                         else if(pos.Eta()<0)
 346                                 etaside = -1;
 347                         else
 348                                 etaside = 1;
 349                 }
 350         }
 351
 352         return etaside;
 353 }
 354
 355
 356 //------------------------------------------------------------------------------
 357 void AliCDMesonUtilsStripped::GetNFO(const AliESDEvent *ESDEvent,
 358                                      const TString etacut, Int_t ctr[])
 359 {
 360         //
 361         // GetNFO
 362         //
 363         // analyzes the SPD fastOR for a given eta range and returns
 364         // an array with the number of hits in:
 365
 366         Int_t ninner=0; // inner layer
 367         Int_t nouter=0; // outer layer
 368         Int_t ipA = 0; // inner layer A side
 369         Int_t ipC = 0; // inner layer C side
 370         Int_t opA = 0; // outer layer A side
 371         Int_t opC = 0; // outer layer C side
 372
 373         const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
 374
 375         // position of the primary vertex
 376         Double_t tmp[3] = { 0., 0., 0. };
 377         ESDEvent->GetPrimaryVertex()->GetXYZ(tmp);
 378         Double_t vtxPos[3] = { tmp[0], tmp[1], tmp[2] };
 379
 380
 381         for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
 382                 if(mult->TestFastOrFiredChips(iChipKey)){
 383                         // here you check if the FastOr bit is 1 or 0
 384                         const Int_t iseta = CheckChipEta(iChipKey, etacut, vtxPos);
 385                         if(iseta==0)
 386                                 continue;
 387
 388                         if(iChipKey<400) {
 389                                 ninner++;  // here you count the FastOr bits in the inner layer
 390                                 if(iseta>0)
 391                                         ipA ++;
 392                                 else
 393                                         ipC ++;
 394                         }
 395                         else {
 396                                 nouter++;  // here you count the FastOr bits in the outer layer
 397                                 if(iseta>0)
 398                                         opA ++;
 399                                 else
 400                                         opC ++;
 401                         }
 402                 }
 403         }
 404
 405         ctr[kInnerPixel]= ninner;
 406         ctr[kOuterPixel]= nouter;
 407         ctr[kIPA]= ipA;
 408         ctr[kIPC]= ipC;
 409         ctr[kOPA]= opA;
 410         ctr[kOPC]= opC;
 411
 412         return;
 413 }