ANALYSIS/TenderSupplies/AliPHOSTenderSupply.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
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //  PHOS tender, recalibrate PHOS clusters                                   //
  20 //  and do track matching                                                    //
  21 //  Author : Dmitri Peressounko (RRC KI)                                     //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25
  26 #include <AliLog.h>
  27 #include <AliESDEvent.h>
  28 #include <AliAnalysisManager.h>
  29 #include <AliTender.h>
  30 #include <AliCDBManager.h>
  31 #include "AliMagF.h"
  32 #include "TGeoGlobalMagField.h"
  33
  34 #include "AliESDCaloCluster.h"
  35 #include "AliPHOSTenderSupply.h"
  36 #include "AliPHOSCalibData.h"
  37 #include "AliPHOSGeometry.h"
  38 #include "AliPHOSEsdCluster.h"
  39
  40 ClassImp(AliPHOSTenderSupply)
  41
  42 AliPHOSTenderSupply::AliPHOSTenderSupply() :
  43   AliTenderSupply()
  44   ,fOCDBpass("local://OCDB")
  45   ,fNonlinearityVersion("Default")
  46   ,fPHOSGeo(0x0)
  47   ,fPHOSCalibData(0x0)
  48 {
  49         //
  50         // default ctor
  51         //
  52    for(Int_t i=0;i<10;i++)fNonlinearityParams[i]=0. ;
  53 }
  54
  55 //_____________________________________________________
  56 AliPHOSTenderSupply::AliPHOSTenderSupply(const char *name, const AliTender *tender) :
  57   AliTenderSupply(name,tender)
  58   ,fOCDBpass("alien:///alice/cern.ch/user/p/prsnko/PHOSrecalibrations/")
  59   ,fNonlinearityVersion("Default")
  60   ,fPHOSGeo(0x0)
  61   ,fPHOSCalibData(0x0)
  62 {
  63         //
  64         // named ctor
  65         //
  66    for(Int_t i=0;i<10;i++)fNonlinearityParams[i]=0. ;
  67 }
  68
  69 //_____________________________________________________
  70 AliPHOSTenderSupply::~AliPHOSTenderSupply()
  71 {
  72   //Destructor
  73   if(fPHOSCalibData)delete fPHOSCalibData;
  74 }
  75
  76 //_____________________________________________________
  77 void AliPHOSTenderSupply::Init()
  78 {
  79   //
  80   // Initialise PHOS tender
  81   //
  82
  83
  84
  85 }
  86
  87 //_____________________________________________________
  88 void AliPHOSTenderSupply::ProcessEvent()
  89 {
  90   //Choose PHOS clusters and recalibrate them
  91   //that it recalculate energy, position and distance
  92   //to closest track extrapolation
  93
  94   AliESDEvent *event=fTender->GetEvent();
  95   if (!event) return;
  96
  97   // Init goemetry
  98   if(!fPHOSGeo){
  99     fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP") ;
 100     for(Int_t mod=0; mod<5; mod++) {
 101       if(!event->GetPHOSMatrix(mod)) continue;
 102       fPHOSGeo->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
 103     }
 104   }
 105
 106
 107   if(!fPHOSCalibData || fTender->RunChanged()){
 108     AliCDBManager * man = AliCDBManager::Instance();
 109     man->SetRun(event->GetRunNumber()) ;
 110     //    man->SetDefaultStorage("local://OCDB");
 111     man->SetSpecificStorage("PHOS/Calib/EmcGainPedestals",fOCDBpass);
 112     if(fPHOSCalibData) delete fPHOSCalibData;
 113     fPHOSCalibData = new AliPHOSCalibData();
 114   }
 115
 116
 117   const AliESDVertex *esdVertex = event->GetPrimaryVertex();
 118   AliESDCaloCells * cells = event->GetPHOSCells() ;
 119   TVector3 vertex(esdVertex->GetX(),esdVertex->GetY(),esdVertex->GetZ());
 120   if(vertex.Mag()>99.) //vertex not defined?
 121     vertex.SetXYZ(0.,0.,0.) ;
 122
 123   //For re-calibration
 124   const Double_t logWeight=4.5 ;
 125
 126   Int_t multClust = event->GetNumberOfCaloClusters();
 127   for (Int_t i=0; i<multClust; i++) {
 128     AliESDCaloCluster *clu = event->GetCaloCluster(i);
 129     if ( !clu->IsPHOS()) continue;
 130
 131     //Apply re-Calibreation
 132     AliPHOSEsdCluster cluPHOS1(*clu);
 133     cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies
 134     cluPHOS1.EvalAll(logWeight,vertex);         // recalculate the cluster parameters
 135     cluPHOS1.SetE(CorrectNonlinearity(cluPHOS1.E()));// Users's nonlinearity
 136
 137     Float_t  xyz[3];
 138     cluPHOS1.GetPosition(xyz);
 139     clu->SetPosition(xyz);                       //rec.point position in MARS
 140     clu->SetE(cluPHOS1.E());                           //total or core particle energy
 141     clu->SetDispersion(cluPHOS1.GetDispersion());  //cluster dispersion
 142     //    ec->SetPID(rp->GetPID()) ;            //array of particle identification
 143     clu->SetM02(cluPHOS1.GetM02()) ;               //second moment M2x
 144     clu->SetM20(cluPHOS1.GetM20()) ;               //second moment M2z
 145     Double_t r=0.,dx=0.,dz=0. ;
 146     TVector3 locPos;
 147     TVector3 globPos(xyz) ;
 148     Int_t relId[4] ;
 149     fPHOSGeo->GlobalPos2RelId(globPos,relId) ;
 150     Int_t mod  = relId[0] ;
 151     fPHOSGeo->Global2Local(locPos,globPos,mod) ;
 152
 153     FindTrackMatching(mod,&locPos,r,dx,dz) ;
 154     clu->SetEmcCpvDistance(r);
 155     clu->SetTrackDistance(dx,dz);
 156     //    clu->SetChi2(-1);                     //not yet implemented
 157     clu->SetTOF(cluPHOS1.GetTOF());
 158
 159   }
 160
 161
 162 }
 163 //___________________________________________________________________________________________________
 164 void AliPHOSTenderSupply::FindTrackMatching(Int_t mod,TVector3 *locpos,Double_t &r,Double_t &dx, Double_t &dz){
 165   //Find track with closest extrapolation to cluster
 166
 167   AliESDEvent *event= fTender->GetEvent();
 168   Double_t  magF = event->GetMagneticField();
 169   Double_t magSign = 1.0;
 170   if(magF<0)magSign = -1.0;
 171
 172   if (!TGeoGlobalMagField::Instance()->GetField()) {
 173     AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG);
 174     TGeoGlobalMagField::Instance()->SetField(field);
 175   }
 176
 177   // *** Start the matching
 178   Int_t nt=event->GetNumberOfTracks();
 179   //Calculate actual distance to PHOS module
 180   TVector3 globaPos ;
 181   fPHOSGeo->Local2Global(mod, 0.,0., globaPos) ;
 182   const Double_t rPHOS = globaPos.Pt() ; //Distance to center of  PHOS module
 183   const Double_t kYmax = 72.+10. ; //Size of the module (with some reserve) in phi direction
 184   const Double_t kZmax = 64.+10. ; //Size of the module (with some reserve) in z direction
 185   const Double_t kAlpha0=330./180.*TMath::Pi() ; //First PHOS module angular direction
 186   const Double_t kAlpha= 20./180.*TMath::Pi() ; //PHOS module angular size
 187   Double_t minDistance = 1.e6;
 188
 189
 190   Double_t gposTrack[3] ;
 191
 192   Double_t bz = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField();
 193   bz = TMath::Sign(0.5*kAlmost0Field,bz) + bz;
 194
 195   Double_t b[3];
 196   for (Int_t i=0; i<nt; i++) {
 197     AliESDtrack *esdTrack=event->GetTrack(i);
 198
 199     // Skip the tracks having "wrong" status (has to be checked/tuned)
 200     ULong_t status = esdTrack->GetStatus();
 201     if ((status & AliESDtrack::kTPCout)   == 0) continue;
 202     //     if ((status & AliESDtrack::kTRDout)   == 0) continue;
 203     //     if ((status & AliESDtrack::kTRDrefit) == 1) continue;
 204
 205     //Continue extrapolation from TPC outer surface
 206     const AliExternalTrackParam *outerParam=esdTrack->GetOuterParam();
 207     if (!outerParam) continue;
 208     AliExternalTrackParam t(*outerParam);
 209
 210     t.GetBxByBz(b) ;
 211     //Direction to the current PHOS module
 212     Double_t phiMod=kAlpha0-kAlpha*mod ;
 213     if(!t.Rotate(phiMod))
 214       continue ;
 215
 216     Double_t y;                       // Some tracks do not reach the PHOS
 217     if (!t.GetYAt(rPHOS,bz,y)) continue; //    because of the bending
 218
 219     Double_t z;
 220     if(!t.GetZAt(rPHOS,bz,z))
 221       continue ;
 222     if (TMath::Abs(z) > kZmax)
 223       continue; // Some tracks miss the PHOS in Z
 224     if(TMath::Abs(y) < kYmax){
 225       t.PropagateToBxByBz(rPHOS,b);        // Propagate to the matching module
 226       //t.CorrectForMaterial(...); // Correct for the TOF material, if needed
 227       t.GetXYZ(gposTrack) ;
 228       TVector3 globalPositionTr(gposTrack) ;
 229       TVector3 localPositionTr ;
 230       fPHOSGeo->Global2Local(localPositionTr,globalPositionTr,mod) ;
 231       Double_t ddx = locpos->X()-localPositionTr.X();
 232       Double_t ddz = locpos->Z()-localPositionTr.Z();
 233       Double_t d2 = ddx*ddx + ddz*ddz;
 234       if(d2 < minDistance) {
 235         dx = ddx ;
 236         dz = ddz ;
 237         minDistance=d2 ;
 238       }
 239     }
 240   } //Scanned all tracks
 241   r=TMath::Sqrt(minDistance) ;
 242
 243 }
 244 //____________________________________________________________
 245 Float_t AliPHOSTenderSupply::CorrectNonlinearity(Float_t en){
 246
 247   //For backward compatibility, if no RecoParameters found
 248   if(fNonlinearityVersion=="Default"){
 249     return 0.0241+1.0504*en+0.000249*en*en ;
 250   }
 251
 252   if(fNonlinearityVersion=="NoCorrection"){
 253     return en ;
 254   }
 255   if(fNonlinearityVersion=="Gustavo2005"){
 256     return fNonlinearityParams[0]+fNonlinearityParams[1]*en + fNonlinearityParams[2]*en*en ;
 257   }
 258   if(fNonlinearityVersion=="Henrik2010"){
 259     return en*(fNonlinearityParams[0]+fNonlinearityParams[1]*TMath::Exp(-en*fNonlinearityParams[2]))*(1.+fNonlinearityParams[3]*TMath::Exp(-en*fNonlinearityParams[4]))*(1.+fNonlinearityParams[6]/(en*en+fNonlinearityParams[5])) ;
 260   }
 261
 262   return en ;
 263 }
 264
 265
 266
 267