PWGGA/PHOSTasks/PHOS_embedding/AliPHOSEmbedding.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 // Analysis task to perform embedding of the simulated AOD
  17 // into real data (ESD)
  18 // Output conntainer contain AOD with
  19 //   0. Data before embedding
  20 //   1. Data+with embedded signal
  21 //   2. MC information from signal
  22 // Data for embedding are obtained using standard Manager
  23 // Signal is read from TChain with aodTree
  24 //
  25 // Authors: Dmitri Peressounko
  26 // Date   : 28.05.2011
  27
  28 #include "TChain.h"
  29 #include "TFile.h"
  30 #include "TROOT.h"
  31 #include "TClonesArray.h"
  32 #include "TH2.h"
  33 #include "TGeoManager.h"
  34 #include "TGeoGlobalMagField.h"
  35
  36 #include "AliGeomManager.h"
  37 #include "AliGRPObject.h"
  38 #include "AliCDBPath.h"
  39 #include "AliCDBEntry.h"
  40 #include "AliESDEvent.h"
  41 #include "AliAODEvent.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliVEventHandler.h"
  44 #include "AliInputEventHandler.h"
  45 #include "AliMultiplicity.h"
  46 #include "AliAODMCParticle.h"
  47 #include "AliAODMCHeader.h"
  48 #include "AliCDBManager.h"
  49 #include "AliAODHandler.h"
  50 #include "AliMagF.h"
  51
  52 #include "AliPHOSReconstructor.h"
  53 #include "AliPHOSClusterizerv1.h"
  54 #include "AliPHOSDigit.h"
  55 #include "AliPHOSGeometry.h"
  56
  57 #include "AliPHOSEmbedding.h"
  58
  59 ClassImp(AliPHOSEmbedding)
  60
  61
  62 //________________________________________________________________________
  63 AliPHOSEmbedding::AliPHOSEmbedding(const char *name)
  64   : AliAnalysisTaskSE(name),
  65     fAODChain(0x0),
  66     fDigitsTree(0x0),
  67     fClustersTree(0x0),
  68     fTreeOut(0x0),
  69     fDigitsArr(0x0),
  70     fEmbeddedClusters(0x0),
  71     fEmbeddedCells(0x0),
  72     fCellsPHOS(0x0),
  73     fClusterizer(0x0),
  74     fPHOSReconstructor(0x0),
  75     fNSignal(0),
  76     fNCaloClustersOld(0),
  77     fInitialized(0)
  78 {
  79   // Constructor
  80   for(int i=0;i<5;i++)fOldPHOSCalibration[i]=0x0 ;
  81
  82 }
  83
  84 //________________________________________________________________________
  85 void AliPHOSEmbedding::UserCreateOutputObjects()
  86 {
  87  //prepare output containers
  88
  89   //Create branch for output MC particles
  90   TClonesArray *mcparticles = new TClonesArray("AliAODMCParticle",0);
  91   mcparticles->SetName(AliAODMCParticle::StdBranchName());
  92   AddAODBranch("TClonesArray", &mcparticles);
  93
  94   //Branch with clusters after embedding
  95   fEmbeddedClusters = new TClonesArray("AliAODCaloCluster",0);
  96   fEmbeddedClusters->SetName("EmbeddedCaloClusters");
  97   AddAODBranch("TClonesArray", &fEmbeddedClusters);
  98
  99
 100   //Branch with cells after embedding
 101   fEmbeddedCells = new AliAODCaloCells();
 102   fEmbeddedCells->SetName("EmbeddedPHOScells");
 103   AddAODBranch("AliAODCaloCells", &fEmbeddedCells);
 104
 105
 106   fDigitsArr = new TClonesArray("AliPHOSDigit",3*56*64) ;
 107
 108   AliAODHandler* handler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 109   if (handler) {
 110       fTreeOut = handler->GetTree();
 111   }
 112   else {
 113        AliWarning("No output AOD Event Handler connected.") ;
 114   }
 115
 116 }
 117
 118 //________________________________________________________________________
 119 void AliPHOSEmbedding::UserExec(Option_t *) {
 120   // Main loop, called for each event
 121   // Perform embedding
 122
 123   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 124   if (!event) {
 125     Printf("ERROR: Could not retrieve event");
 126     PostData(0, fTreeOut);
 127     return;
 128   }
 129
 130   AliCentrality *centrality = event->GetCentrality();
 131   Float_t acentrality=centrality->GetCentralityPercentileUnchecked("V0M");
 132   if( acentrality <= 0. || acentrality>80.){
 133     return;
 134   }
 135
 136   //Read next AOD event
 137   //If necesary method checks if AOD event is good to embed
 138   //e.g. there are PHOS clusters etc.
 139   AliAODEvent * signal = GetNextSignalEvent() ;
 140   if (!signal) {
 141     Printf("ERROR: Could not retrieve signal");
 142     PostData(0, fTreeOut);
 143     return;
 144   }
 145
 146   Init() ;
 147
 148   //Remember PHOS digits before any modification
 149   CopyRecalibrateDigits() ;
 150
 151   //First re-reconstruct existing digits to ensure identical reconstruction before and
 152   //after embeding
 153   AliAODEvent * nullSignal=0x0 ;
 154   MakeEmbedding(event,nullSignal) ;
 155
 156   //Convert ESD + MC information from signal to AOD
 157   ConvertESDtoAOD(event);
 158
 159   //Now perform real embedding
 160   //Embed signal clusters into event
 161   MakeEmbedding(event,signal) ;
 162
 163   //Convert ESD to AOD
 164   ConvertEmbeddedClusters(event) ;
 165
 166   //Fill MC information
 167   ConvertMCParticles(signal) ;
 168
 169   delete signal ;
 170
 171   PostData(0, fTreeOut);
 172 }
 173 //______________________________________________________________________________
 174 void AliPHOSEmbedding::CopyRecalibrateDigits(){
 175   //Recalibrate digits if there is better calibration ("new")
 176   //exists than one used in reconstruction ("ESD")
 177
 178   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 179   if(!event){
 180     AliError("No ESD event!");
 181     return ;
 182   }
 183   if(fCellsPHOS)
 184     delete fCellsPHOS ;
 185   fCellsPHOS = new AliESDCaloCells() ;
 186   fCellsPHOS->CreateContainer(event->GetPHOSCells()->GetNumberOfCells());
 187
 188   //Apply recalibration if necessary to account for
 189   //difference between calibration used in reconstruction and
 190   //final calibration
 191 //  if(fOldPHOSCalibration[0]){ //there is a difference in calibrations
 192     for (Short_t icell = 0; icell < event->GetPHOSCells()->GetNumberOfCells(); icell++) {
 193       Short_t id=0;
 194       Double_t time=0., amp=0. ;
 195       Int_t mclabel;
 196       Double_t efrac =0. ;
 197 //      if (fCellsPHOS->GetCell(icell, id, amp,time,mclabel,efrac) != kTRUE)
 198       if (event->GetPHOSCells()->GetCell(icell, id, amp,time,mclabel,efrac) != kTRUE)
 199         break;
 200       Int_t relId[4] ;
 201       AliPHOSGeometry::GetInstance()->AbsToRelNumbering(id,relId);
 202       Float_t calibESD=fOldPHOSCalibration[relId[0]-1]->GetBinContent(relId[2],relId[3]) ;
 203       Float_t calibNew=fPHOSReconstructor->Calibrate(1.,id) ;
 204       if(calibESD>0.)
 205         amp=amp*calibNew/calibESD ;
 206       fCellsPHOS->SetCell(icell, id, amp, time,mclabel,efrac);
 207     }
 208 //  }
 209 }
 210 //______________________________________________________________________________
 211 void AliPHOSEmbedding::Init(){
 212   if(fInitialized)
 213     return ;
 214
 215   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 216   if(!event){
 217     AliError("Can not obtain InputEvent!") ;
 218     return ;
 219   }
 220   Int_t runNum = event->GetRunNumber();
 221   AliCDBManager::Instance()->SetRun(runNum);
 222
 223   fPHOSReconstructor = new AliPHOSReconstructor() ;
 224   AliCDBPath path("PHOS","Calib","RecoParam");
 225   AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
 226   if(!entry){
 227     AliError(Form("Can not get OCDB entry %s",path.GetPath().Data())) ;
 228     return ;
 229   }
 230
 231   TObjArray* recoParamArray = (TObjArray*)entry->GetObject();
 232   AliPHOSRecoParam* recoParam = (AliPHOSRecoParam*)recoParamArray->At(2);
 233   fPHOSReconstructor->SetRecoParam(recoParam) ;
 234
 235   InitMF() ;
 236   InitGeometry() ;
 237
 238   fInitialized=kTRUE ;
 239 }
 240
 241 //________________________________________________________________________
 242 AliAODEvent* AliPHOSEmbedding::GetNextSignalEvent(){
 243   //Read signal AOD event from the chain
 244
 245   if(fAODChain==0){
 246     AliError("No chain to read signal events") ;
 247     return 0 ;
 248   }
 249
 250   AliAODEvent* aod = new AliAODEvent;
 251   aod->ReadFromTree(fAODChain);
 252
 253   if(fNSignal>=fAODChain->GetEntries()){
 254     delete aod ;
 255     return 0 ;
 256   }
 257   fAODChain->GetEvent(fNSignal) ;
 258   fNSignal++ ;
 259
 260   //check if event is read and there is at least 2 PHOS clusters
 261   //Otherwise fill histogram and read next event
 262   //TODO
 263
 264   return aod ;
 265
 266 }
 267
 268 //________________________________________________________________________
 269 void AliPHOSEmbedding::MakeEmbedding(AliESDEvent *event,AliAODEvent * signal){
 270   //Perform embedding of the signal to the event
 271
 272   gROOT->cd(); //make sure that the digits and RecPoints Trees are memory resident
 273   if(fDigitsTree){
 274     delete fDigitsTree ;
 275   }
 276   fDigitsTree = new TTree("digitstree","digitstree");
 277   fDigitsTree->Branch("PHOS","TClonesArray", &fDigitsArr, 32000);
 278
 279   if(fClustersTree){
 280     delete fClustersTree ;
 281   }
 282   fClustersTree = new TTree("clustertree","clustertree");
 283
 284   //Remember number of Clusters before we added new ones...
 285   fNCaloClustersOld=event->GetNumberOfCaloClusters() ;
 286
 287
 288 /*
 289   printf("---------------Before---------\n") ;
 290   Int_t n=event->GetNumberOfCaloClusters() ;
 291   for (Int_t iClust=0; iClust<event->GetNumberOfCaloClusters(); ++iClust) {
 292     AliESDCaloCluster * cluster = event->GetCaloCluster(iClust);
 293      if(!cluster->IsPHOS())
 294       continue;
 295    Float_t pos[3] = { 0.};
 296     cluster->GetPosition(pos);
 297     printf("Cluster(%d): E=%f, x=%f, z=%f, CPV=%f \n",iClust,cluster->E(),pos[0],pos[2],cluster->GetEmcCpvDistance()) ;
 298     UShort_t * index    = cluster->GetCellsAbsId() ;
 299     for(Int_t ic=0; ic < cluster->GetNCells(); ic++ )
 300     printf("Dig(%d)=%d, ",ic,index[ic]) ;
 301   printf("\n") ;
 302
 303
 304   }
 305   printf("---------------END Before---------\n") ;
 306 */
 307
 308
 309   //create digits
 310   MakeDigits(signal);
 311
 312 /*
 313 for(Int_t i=0;i<fDigitsArr->GetEntriesFast();i++){
 314    AliPHOSDigit * d = (AliPHOSDigit *) fDigitsArr->At(i) ;
 315    printf("  Digit(%d) = %d, E=%f \n",i,d->GetId(),d->GetEnergy()) ;
 316 }
 317 */
 318
 319   //clusterize and make tracking
 320   fPHOSReconstructor->Reconstruct(fDigitsTree,fClustersTree) ;
 321
 322   //Note that the cirrent ESDEvent will be modified!
 323   fPHOSReconstructor->FillESD(fDigitsTree, fClustersTree, event) ;
 324
 325
 326 /*
 327    printf("---------------After---------\n") ;
 328   for (Int_t iClust=n; iClust<event->GetNumberOfCaloClusters(); ++iClust) {
 329     AliESDCaloCluster * cluster = event->GetCaloCluster(iClust);
 330     if(!cluster->IsPHOS())
 331       continue;
 332     Float_t pos[3] = { 0.};
 333     cluster->GetPosition(pos);
 334     printf("Cluster(%d): E=%f, x=%f, z=%f, CPV=%f, Label=%d \n",iClust-n,cluster->E(),pos[0],pos[2],cluster->GetEmcCpvDistance(),cluster->GetLabel()) ;
 335     UShort_t * index    = cluster->GetCellsAbsId() ;
 336     for(Int_t ic=0; ic < cluster->GetNCells(); ic++ )
 337       printf("Dig(%d)=%d, ",ic,index[ic]) ;
 338   printf("\n") ;
 339   }
 340   printf("---------------END After---------\n") ;
 341 */
 342
 343
 344 }
 345 //______________________________________________________________________________
 346 void AliPHOSEmbedding::ConvertESDtoAOD(AliESDEvent* esd)
 347 {
 348   // ESD Filter analysis task executed for each event
 349
 350   if(!esd)return;
 351
 352   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
 353   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
 354
 355   ConvertHeader(*esd);
 356
 357   // Fetch Stack for debuggging if available
 358
 359   // Update the header
 360   Int_t nTracks    = esd->GetNumberOfTracks();
 361   Int_t nV0s      = esd->GetNumberOfV0s();
 362   Int_t nCascades = esd->GetNumberOfCascades();
 363   Int_t nKinks    = esd->GetNumberOfKinks();
 364   Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;
 365   Int_t nPileSPDVertices=1+esd->GetNumberOfPileupVerticesSPD(); // also SPD main vertex
 366   Int_t nPileTrkVertices=esd->GetNumberOfPileupVerticesTracks();
 367   nVertices+=nPileSPDVertices;
 368   nVertices+=nPileTrkVertices;
 369   Int_t nJets     = 0;
 370   Int_t nCaloClus = esd->GetNumberOfCaloClusters();
 371   Int_t nFmdClus  = 0;
 372   Int_t nPmdClus  = esd->GetNumberOfPmdTracks();
 373
 374
 375   AODEvent()->ResetStd(nTracks, nVertices, nV0s, nCascades, nJets, nCaloClus, nFmdClus, nPmdClus);
 376
 377   ConvertPrimaryVertices(*esd);
 378
 379   ConvertCaloClusters(*esd);
 380
 381   ConvertEMCALCells(*esd);
 382
 383   ConvertPHOSCells(*esd);
 384
 385 }
 386 //______________________________________________________________________________
 387 void AliPHOSEmbedding::ConvertPHOSCells(const AliESDEvent& esd)
 388 {
 389   // fill PHOS cell info
 390   if (esd.GetPHOSCells()) { // protection against missing ESD information
 391     AliESDCaloCells &esdPHcells = *(esd.GetPHOSCells());
 392     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
 393
 394     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
 395     aodPHcells.CreateContainer(nPHcell);
 396     aodPHcells.SetType(AliAODCaloCells::kPHOSCell);
 397     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
 398       aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell),esdPHcells.GetTime(iCell),-1,0.);
 399     }
 400     aodPHcells.Sort();
 401   }
 402 }
 403
 404 //______________________________________________________________________________
 405 void AliPHOSEmbedding::ConvertEMCALCells(const AliESDEvent& esd)
 406 {
 407   // fill EMCAL cell info
 408   if (esd.GetEMCALCells()) { // protection against missing ESD information
 409     AliESDCaloCells &esdEMcells = *(esd.GetEMCALCells());
 410     Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
 411
 412     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 413     aodEMcells.CreateContainer(nEMcell);
 414     aodEMcells.SetType(AliAODCaloCells::kEMCALCell);
 415     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
 416       aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell),esdEMcells.GetTime(iCell));
 417     }
 418     aodEMcells.Sort();
 419   }
 420 }
 421 //______________________________________________________________________________
 422 void AliPHOSEmbedding::ConvertCaloClusters(const AliESDEvent& esd)
 423 {
 424
 425   // Access to the AOD container of clusters
 426   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
 427   Int_t jClusters(0);
 428
 429   for (Int_t iClust=fNCaloClustersOld; iClust<esd.GetNumberOfCaloClusters(); ++iClust) {
 430
 431     AliESDCaloCluster * cluster = esd.GetCaloCluster(iClust);
 432
 433     Int_t  id        = cluster->GetID();
 434     Int_t  nLabel    = cluster->GetNLabels();
 435     Int_t *labels    = cluster->GetLabels();
 436
 437     Float_t energy = cluster->E();
 438     Float_t posF[3] = { 0.};
 439     cluster->GetPosition(posF);
 440
 441     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,
 442                                                                                       nLabel,
 443                                                                                       labels,
 444                                                                                       energy,
 445                                                                                       posF,
 446                                                                                       NULL,
 447                                                                                       cluster->GetType(),0);
 448
 449     Double_t dx=cluster->GetTrackDx() ;
 450     Double_t dz=cluster->GetTrackDz() ;
 451     Float_t cpv=99. ; //No track matched by default
 452     TArrayI * itracks = cluster->GetTracksMatched() ;
 453     if(itracks->GetSize()>0){
 454       Int_t iTr = itracks->At(0);
 455       if(iTr>=0 && iTr<esd.GetNumberOfTracks()){
 456         AliESDtrack *track = esd.GetTrack(iTr) ;
 457         Double_t pt = track->Pt() ;
 458         Short_t charge = track->Charge() ;
 459         cpv=TestCPV(dx, dz, pt,charge) ;
 460       }
 461     }
 462
 463     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
 464                                 cluster->GetDispersion(),
 465                                 cluster->GetM20(), cluster->GetM02(),
 466                                 cpv,
 467                                 cluster->GetNExMax(),cluster->GetTOF()) ;
 468
 469     caloCluster->SetPIDFromESD(cluster->GetPID());
 470     caloCluster->SetNCells(cluster->GetNCells());
 471     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
 472     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
 473
 474     /*
 475     TArrayI* matchedT =         cluster->GetTracksMatched();
 476     if (fNumberOfTracks>0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {
 477       for (Int_t im = 0; im < matchedT->GetSize(); im++) {
 478         Int_t iESDtrack = matchedT->At(im);;
 479         if (fAODTrackRefs->At(iESDtrack) != 0) {
 480           caloCluster->AddTrackMatched((AliAODTrack*)fAODTrackRefs->At(iESDtrack));
 481         }
 482       }
 483     }
 484     */
 485
 486   }
 487   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
 488 }
 489 //______________________________________________________________________________
 490 void AliPHOSEmbedding::ConvertEmbeddedClusters(const AliESDEvent* esd)
 491 {
 492    //Copy PHOS clusters and cells after embedding
 493
 494   // Access to the AOD container of clusters
 495   Int_t jClusters(0);
 496   fEmbeddedClusters->Clear();
 497   fEmbeddedClusters->Expand(esd->GetNumberOfCaloClusters()-fNCaloClustersOld) ;
 498   for (Int_t iClust=fNCaloClustersOld; iClust<esd->GetNumberOfCaloClusters(); ++iClust) {
 499
 500     AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
 501
 502     Int_t  id        = cluster->GetID();
 503     Int_t  nLabel    = cluster->GetNLabels();
 504     Int_t *labels    = cluster->GetLabels();
 505
 506     Float_t energy = cluster->E();
 507     Float_t posF[3] = { 0.};
 508     cluster->GetPosition(posF);
 509
 510     AliAODCaloCluster *caloCluster = new((*fEmbeddedClusters)[jClusters++]) AliAODCaloCluster(id,
 511                                                                                       nLabel,
 512                                                                                       labels,
 513                                                                                       energy,
 514                                                                                       posF,
 515                                                                                       NULL,
 516                                                                                       cluster->GetType(),0);
 517
 518
 519     Double_t dx=cluster->GetTrackDx() ;
 520     Double_t dz=cluster->GetTrackDz() ;
 521     Float_t cpv=99. ; //No track matched by default
 522     TArrayI * itracks = cluster->GetTracksMatched() ;
 523     if(itracks->GetSize()>0){
 524       Int_t iTr = itracks->At(0);
 525       if(iTr>=0 && iTr<esd->GetNumberOfTracks()){
 526         AliESDtrack *track = esd->GetTrack(iTr) ;
 527         Double_t pt = track->Pt() ;
 528         Short_t charge = track->Charge() ;
 529         cpv=TestCPV(dx, dz, pt,charge) ;
 530       }
 531     }
 532
 533     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
 534                                 cluster->GetDispersion(),
 535                                 cluster->GetM20(), cluster->GetM02(),
 536                                 cpv,
 537                                 cluster->GetNExMax(),cluster->GetTOF()) ;
 538
 539     caloCluster->SetPIDFromESD(cluster->GetPID());
 540     caloCluster->SetNCells(cluster->GetNCells());
 541     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
 542     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
 543   }
 544
 545   for(Int_t i=0;i<fEmbeddedClusters->GetEntriesFast();i++){
 546     AliAODCaloCluster *caloCluster =(AliAODCaloCluster *)fEmbeddedClusters->At(i) ;
 547     caloCluster->GetID() ;
 548
 549   }
 550
 551   //Now cells
 552   if (esd->GetPHOSCells()) { // protection against missing ESD information
 553     AliESDCaloCells &esdPHcells = *(esd->GetPHOSCells());
 554     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
 555     fEmbeddedCells->CreateContainer(nPHcell);
 556     fEmbeddedCells->SetType(AliAODCaloCells::kPHOSCell);
 557     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
 558       fEmbeddedCells->SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell),esdPHcells.GetTime(iCell));
 559     }
 560     fEmbeddedCells->Sort();
 561   }
 562
 563
 564 }
 565 //______________________________________________________________________________
 566 void AliPHOSEmbedding::ConvertMCParticles(const AliAODEvent* aod)
 567 {
 568   //Copy MC branches to new AOD
 569
 570   TClonesArray *mcArray = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
 571   TClonesArray * aodMcParticles = (TClonesArray*)AODEvent()->FindListObject(AliAODMCParticle::StdBranchName());
 572   for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
 573
 574     AliAODMCParticle* aodpart =  (AliAODMCParticle*) mcArray->At(i);
 575 //    AliAODMCParticle * mcp =new ((*aodMcParticles)[i]) AliAODMCParticle(*aodpart);
 576     new ((*aodMcParticles)[i]) AliAODMCParticle(*aodpart);
 577
 578   }
 579
 580   /*
 581   // MC header
 582   AliAODMCHeader * mcHeader = (AliAODMCHeader *)aod.FindListObject(AliAODMCHeader::StdBranchName());
 583   AliAODMCHeader * aodMcHeader = new AliAODMCHeader(*mcHeader);
 584   aodMcHeader->SetName(AliAODMCHeader::StdBranchName());
 585   AddAODBranch("AliAODMCHeader",&aodMcHeader);
 586   */
 587
 588 }
 589 //______________________________________________________________________________
 590 void AliPHOSEmbedding::ConvertHeader(AliESDEvent & esd){
 591
 592   AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
 593   if(!header) AliFatal("Not a standard AOD");
 594
 595   header->SetRunNumber(esd.GetRunNumber());
 596   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 597
 598   TTree* tree = fInputHandler->GetTree();
 599   if (tree) {
 600     TFile* file = tree->GetCurrentFile();
 601     if (file) header->SetESDFileName(file->GetName());
 602   }
 603
 604   header->SetBunchCrossNumber(esd.GetBunchCrossNumber());
 605   header->SetOrbitNumber(esd.GetOrbitNumber());
 606   header->SetPeriodNumber(esd.GetPeriodNumber());
 607   header->SetEventType(esd.GetEventType());
 608
 609   header->SetEventNumberESDFile(esd.GetHeader()->GetEventNumberInFile());
 610   if(const_cast<AliESDEvent&>(esd).GetCentrality()){
 611     header->SetCentrality(const_cast<AliESDEvent&>(esd).GetCentrality());
 612   }
 613   else{
 614     header->SetCentrality(0);
 615   }
 616   if(const_cast<AliESDEvent&>(esd).GetEventplane()){
 617     header->SetEventplane(const_cast<AliESDEvent&>(esd).GetEventplane());
 618   }
 619   else{
 620     header->SetEventplane(0);
 621   }
 622
 623   // Trigger
 624   header->SetFiredTriggerClasses(esd.GetFiredTriggerClasses());
 625   header->SetTriggerMask(esd.GetTriggerMask());
 626   header->SetTriggerCluster(esd.GetTriggerCluster());
 627   header->SetL0TriggerInputs(esd.GetHeader()->GetL0TriggerInputs());
 628   header->SetL1TriggerInputs(esd.GetHeader()->GetL1TriggerInputs());
 629   header->SetL2TriggerInputs(esd.GetHeader()->GetL2TriggerInputs());
 630
 631   header->SetMagneticField(esd.GetMagneticField());
 632   header->SetMuonMagFieldScale(esd.GetCurrentDip()/6000.);
 633   header->SetZDCN1Energy(esd.GetZDCN1Energy());
 634   header->SetZDCP1Energy(esd.GetZDCP1Energy());
 635   header->SetZDCN2Energy(esd.GetZDCN2Energy());
 636   header->SetZDCP2Energy(esd.GetZDCP2Energy());
 637   header->SetZDCEMEnergy(esd.GetZDCEMEnergy(0),esd.GetZDCEMEnergy(1));
 638
 639   // ITS Cluster Multiplicty
 640   const AliMultiplicity *mult = esd.GetMultiplicity();
 641   for (Int_t ilay = 0; ilay < 6; ilay++) header->SetITSClusters(ilay, mult->GetNumberOfITSClusters(ilay));
 642
 643   // TPC only Reference Multiplicty
 644   //  Int_t refMult  = fTPCaloneTrackCuts ? (Short_t)fTPCaloneTrackCuts->GetReferenceMultiplicity(&esd, kTRUE) : -1;
 645   header->SetTPConlyRefMultiplicity(-1);
 646
 647   //
 648   Float_t diamxy[2]={static_cast<Float_t>(esd.GetDiamondX()),static_cast<Float_t>(esd.GetDiamondY())};
 649   Float_t diamcov[3];
 650   esd.GetDiamondCovXY(diamcov);
 651   header->SetDiamond(diamxy,diamcov);
 652   header->SetDiamondZ(esd.GetDiamondZ(),esd.GetSigma2DiamondZ());
 653
 654
 655   //PHOS matrixes
 656   char path[255] ;
 657   TGeoHMatrix * m ;
 658   for(Int_t mod=0; mod<5; mod++){
 659     snprintf(path,255,"/ALIC_1/PHOS_%d",mod+1) ; //In Geometry modules numbered 1,2,.,5
 660     if (gGeoManager->cd(path)){
 661       m = gGeoManager->GetCurrentMatrix() ;
 662       header->SetPHOSMatrix(new TGeoHMatrix(*m),mod) ;
 663     }
 664     else{
 665       header->SetPHOSMatrix(NULL,mod) ;
 666     }
 667   }
 668
 669   //EventPlane
 670   AliEventplane *eventplane = esd.GetEventplane();
 671   Double_t epQ = eventplane->GetEventplane("Q");
 672
 673   //Standard Eventplane setter is too complicated...
 674   //use durty hack
 675   header->SetZDCN1Energy(epQ) ;
 676
 677   AliCentrality *centrality = esd.GetCentrality();
 678   Double_t acentrality=centrality->GetCentralityPercentileUnchecked("V0M");
 679   header->SetZDCN2Energy(acentrality) ;
 680
 681
 682 }//______________________________________________________________________________
 683 void AliPHOSEmbedding::ConvertPrimaryVertices(AliESDEvent const&esd){
 684   // Access to the AOD container of vertices
 685
 686   // Access to the AOD container of vertices
 687   Int_t fNumberOfVertices = 0;
 688
 689   Double_t pos[3] = { 0. };
 690   Double_t covVtx[6] = { 0. };
 691
 692   // Add primary vertex. The primary tracks will be defined
 693   // after the loops on the composite objects (V0, cascades, kinks)
 694   const AliESDVertex *vtx = esd.GetPrimaryVertex();
 695
 696   vtx->GetXYZ(pos); // position
 697   vtx->GetCovMatrix(covVtx); //covariance matrix
 698
 699   TClonesArray * vertexes=AODEvent()->GetVertices();
 700
 701   AliAODVertex* primaryVertex = new((*vertexes)[fNumberOfVertices++])
 702   AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
 703   primaryVertex->SetName(vtx->GetName());
 704   primaryVertex->SetTitle(vtx->GetTitle());
 705
 706   TString vtitle = vtx->GetTitle();
 707   if (!vtitle.Contains("VertexerTracks"))
 708     primaryVertex->SetNContributors(vtx->GetNContributors());
 709
 710   if (fDebug > 0) primaryVertex->Print();
 711
 712   // Add SPD "main" vertex
 713   const AliESDVertex *vtxS = esd.GetPrimaryVertexSPD();
 714   vtxS->GetXYZ(pos); // position
 715   vtxS->GetCovMatrix(covVtx); //covariance matrix
 716   AliAODVertex * mVSPD = new((*vertexes)[fNumberOfVertices++])
 717   AliAODVertex(pos, covVtx, vtxS->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainSPD);
 718   mVSPD->SetName(vtxS->GetName());
 719   mVSPD->SetTitle(vtxS->GetTitle());
 720   mVSPD->SetNContributors(vtxS->GetNContributors());
 721
 722   // Add SPD pileup vertices
 723   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesSPD(); ++iV)
 724   {
 725     const AliESDVertex *vtxP = esd.GetPileupVertexSPD(iV);
 726     vtxP->GetXYZ(pos); // position
 727     vtxP->GetCovMatrix(covVtx); //covariance matrix
 728     AliAODVertex * pVSPD =  new((*vertexes)[fNumberOfVertices++])
 729     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupSPD);
 730     pVSPD->SetName(vtxP->GetName());
 731     pVSPD->SetTitle(vtxP->GetTitle());
 732     pVSPD->SetNContributors(vtxP->GetNContributors());
 733   }
 734
 735   // Add TRK pileup vertices
 736   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesTracks(); ++iV)
 737   {
 738     const AliESDVertex *vtxP = esd.GetPileupVertexTracks(iV);
 739     vtxP->GetXYZ(pos); // position
 740     vtxP->GetCovMatrix(covVtx); //covariance matrix
 741     AliAODVertex * pVTRK = new((*vertexes)[fNumberOfVertices++])
 742     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupTracks);
 743     pVTRK->SetName(vtxP->GetName());
 744     pVTRK->SetTitle(vtxP->GetTitle());
 745     pVTRK->SetNContributors(vtxP->GetNContributors());
 746   }
 747
 748
 749 }
 750 //__________________________________________________________________________________
 751 void AliPHOSEmbedding::MakeDigits(AliAODEvent * signal){
 752
 753   //-------------------------------------------------------------------------------------
 754   //Transform CaloCells into Digits which can be used for standard reconstruction
 755   //Add signal digits to the event
 756   //-------------------------------------------------------------------------------------
 757
 758   fDigitsArr->Clear() ;
 759   fDigitsTree->Branch("PHOS","TClonesArray", &fDigitsArr, 32000);
 760
 761   //First copy data digits
 762   Int_t ndigit=0 ;
 763   for (Short_t icell = 0; icell < fCellsPHOS->GetNumberOfCells(); icell++) {
 764     Short_t id=0;
 765     Double_t time=0., amp=0. ;
 766     Int_t mclabel;
 767     Double_t efrac =0. ;
 768     if (fCellsPHOS->GetCell(icell, id, amp, time,mclabel, efrac) != kTRUE)
 769       break;
 770
 771     new((*fDigitsArr)[ndigit]) AliPHOSDigit(-1,id,float(amp),float(time),ndigit);
 772     ndigit++;
 773   }
 774
 775
 776   //Add Digits from Signal
 777   TClonesArray sdigits("AliPHOSDigit",0) ;
 778   Int_t isdigit=0 ;
 779   if(signal){
 780     AliAODCaloCells* cellsS = signal->GetPHOSCells();
 781     Int_t cellLabels[1000]={0} ;       //1000 should be enough for simulated
 782 //    Int_t cellSecondLabels[1000]={0} ; //low-statistics event.
 783     for(Int_t i=0;i<cellsS->GetNumberOfCells();i++){
 784       cellLabels[i]=-1 ;
 785 //      cellSecondLabels[i]=-1;
 786     }
 787     //------------------------------------------------------------------------------------
 788     //Ancestry information
 789     //Celect digits contributing to signal clusters and add primary information
 790     //(it is not stored in CaloCells)
 791     //------------------------------------------------------------------------------------
 792     sdigits.Expand(cellsS->GetNumberOfCells());
 793     for(Int_t i=0; i<signal->GetNumberOfCaloClusters(); i++) {
 794       //cluster from embedded signal
 795       AliVCluster *clus = signal->GetCaloCluster(i);
 796
 797       if(!clus->IsPHOS())
 798         continue;
 799
 800       Int_t label = clus->GetLabel();
 801 //      Int_t label2 = -1 ;
 802 //      if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
 803
 804       UShort_t * index    = clus->GetCellsAbsId() ;
 805       for(Int_t ic=0; ic < clus->GetNCells(); ic++ ){
 806         for (Int_t icell = 0; icell < cellsS->GetNumberOfCells(); icell++){
 807            Short_t cellNumber;
 808            Double_t cellAmplitude=0., cellTime=0. ;
 809            Int_t mclabel;
 810            Double_t efrac =0. ;
 811            cellsS->GetCell(icell, cellNumber, cellAmplitude, cellTime,mclabel,efrac) ;
 812            if(cellNumber==index[ic]){
 813               cellLabels[icell]=label;
 814 //              cellSecondLabels[icell]=label2;
 815               break ;
 816           }
 817         }
 818       }
 819     }
 820
 821     for (Int_t icell = 0; icell < cellsS->GetNumberOfCells(); icell++) {
 822       Short_t cellNumber;
 823       Double_t cellAmplitude=0., cellTime=0. ;
 824       Int_t mclabel;
 825       Double_t efrac =0. ;
 826       if (cellsS->GetCell(icell, cellNumber, cellAmplitude, cellTime,mclabel,efrac) != kTRUE)
 827         break;
 828       //Add only digits related to the cluster, no noisy digits...
 829       if(cellLabels[icell]==-1){
 830         continue ;
 831       }
 832       new(sdigits[isdigit]) AliPHOSDigit(cellLabels[icell],cellNumber,float(cellAmplitude),float(cellTime),isdigit);
 833       isdigit++;
 834     }
 835   }
 836
 837
 838   //Merge digits
 839   Int_t icurrent = 0 ; //index of the last used digit in underlying event
 840   fDigitsArr->Expand(ndigit+isdigit) ;
 841   for(Int_t i=0; i<isdigit;i++){
 842     AliPHOSDigit * sdigit=static_cast<AliPHOSDigit*>(sdigits.At(i)) ;
 843     Bool_t added=kFALSE ;
 844     for(Int_t id=icurrent;id<ndigit;id++){
 845       AliPHOSDigit * digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(id)) ;
 846       if(sdigit->GetId() ==  digit->GetId() ){
 847         *digit += *sdigit ;  //add energies
 848         icurrent=id+1 ;
 849         added=kTRUE ;
 850         break ; //no more digits with same ID in the list
 851       }
 852       if(sdigit->GetId() < digit->GetId() ){
 853         icurrent=id ;
 854         break ; //no more digits with same ID in the list
 855       }
 856     }
 857     if(!added){
 858       new((*fDigitsArr)[ndigit]) AliPHOSDigit(*sdigit) ;
 859       ndigit++ ;
 860     }
 861   }
 862
 863   //Change Amp back from Energy to ADC counts
 864   //Note that Reconstructor uses best ("new") calibration
 865   for(Int_t i=0; i<ndigit;i++){
 866      AliPHOSDigit * digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
 867      Float_t calib =fPHOSReconstructor->Calibrate(1.,digit->GetId()) ;
 868      if(calib>0.)
 869        digit->SetEnergy(digit->GetEnergy()/calib) ;
 870   }
 871
 872   fDigitsArr->Sort() ;
 873   for (Int_t i = 0 ; i < ndigit ; i++) {
 874     AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(i) ) ;
 875     digit->SetIndexInList(i) ;
 876   }
 877   fDigitsTree->Fill() ;
 878
 879 }
 880 //__________________________________________________________________________________
 881 void AliPHOSEmbedding::SetOldCalibration(TH2F **calibH){
 882   //Copy histograms with calibration coeff
 883   gROOT->cd() ;
 884   char name[55] ;
 885   for(Int_t i=0;i<5;i++){
 886     fOldPHOSCalibration[i] = new TH2F(*((TH2F*)calibH[i])) ;
 887     snprintf(name,55,"OldPHOSCalibrationMod%d",i) ;
 888     fOldPHOSCalibration[i]->SetName(name) ;
 889   }
 890 }
 891 //____________________________________________________________________________
 892 void AliPHOSEmbedding::InitMF(){
 893
 894   printf("............Init MF \n") ;
 895   //------------------------------------
 896   // Initialization of the Mag.Fiels from GRP entry
 897   // Copied from AliReconstruction::InitGRP()
 898   //------------------------------------
 899   AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
 900   AliGRPObject * aGRPData=0 ;
 901   if (entry) {
 902
 903     TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
 904
 905     if (m) {
 906        AliInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
 907        m->Print();
 908        aGRPData = new AliGRPObject();
 909        aGRPData->ReadValuesFromMap(m);
 910     }
 911
 912     else {
 913        AliInfo("Found an AliGRPObject in GRP/GRP/Data, reading it");
 914        aGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
 915        entry->SetOwner(0);
 916     }
 917   }
 918
 919   if (!aGRPData) {
 920      AliError("No GRP entry found in OCDB!");
 921      return ;
 922   }
 923   //*** Dealing with the magnetic field map
 924
 925   TString lhcState = aGRPData->GetLHCState();
 926   if (lhcState==AliGRPObject::GetInvalidString()) {
 927     AliError("GRP/GRP/Data entry:  missing value for the LHC state ! Using UNKNOWN");
 928     lhcState = "UNKNOWN";
 929   }
 930
 931   TString beamType = aGRPData->GetBeamType();
 932   if (beamType==AliGRPObject::GetInvalidString()) {
 933     AliError("GRP/GRP/Data entry:  missing value for the beam type ! Using UNKNOWN");
 934     beamType = "UNKNOWN";
 935   }
 936
 937   Float_t beamEnergy = aGRPData->GetBeamEnergy();
 938   if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
 939     AliError("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
 940     beamEnergy = 0;
 941   }
 942
 943   TString runType = aGRPData->GetRunType();
 944   if (runType==AliGRPObject::GetInvalidString()) {
 945     AliError("GRP/GRP/Data entry:  missing value for the run type ! Using UNKNOWN");
 946     runType = "UNKNOWN";
 947   }
 948
 949   if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
 950     if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
 951       AliInfo("PHOSEmbedding: MF is locked - GRP information will be ignored !");
 952       AliInfo("Running with the externally locked B field !");
 953     }
 954     else {
 955       AliInfo("Destroying existing B field instance!");
 956       delete TGeoGlobalMagField::Instance();
 957     }
 958   }
 959   if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
 960     // Construct the field map out of the information retrieved from GRP.
 961     Bool_t ok = kTRUE;
 962     // L3
 963     Float_t l3Current = aGRPData->GetL3Current((AliGRPObject::Stats)0);
 964     if (l3Current == AliGRPObject::GetInvalidFloat()) {
 965       AliError("GRP/GRP/Data entry:  missing value for the L3 current !");
 966       ok = kFALSE;
 967     }
 968
 969     Char_t l3Polarity = aGRPData->GetL3Polarity();
 970     if (l3Polarity == AliGRPObject::GetInvalidChar()) {
 971       AliError("GRP/GRP/Data entry:  missing value for the L3 polarity !");
 972       ok = kFALSE;
 973     }
 974
 975     // Dipole
 976     Float_t diCurrent = aGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
 977     if (diCurrent == AliGRPObject::GetInvalidFloat()) {
 978       AliError("GRP/GRP/Data entry:  missing value for the dipole current !");
 979       ok = kFALSE;
 980     }
 981
 982     Char_t diPolarity = aGRPData->GetDipolePolarity();
 983     if (diPolarity == AliGRPObject::GetInvalidChar()) {
 984       AliError("GRP/GRP/Data entry:  missing value for the dipole polarity !");
 985       ok = kFALSE;
 986     }
 987
 988     // read special bits for the polarity convention and map type
 989     Int_t  polConvention = aGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
 990     Bool_t uniformB = aGRPData->IsUniformBMap();
 991
 992     if (ok) {
 993       AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1),
 994                                              TMath::Abs(diCurrent) * (diPolarity ? -1:1),
 995                                              polConvention,uniformB,beamEnergy, beamType.Data());
 996       if (fld) {
 997         TGeoGlobalMagField::Instance()->SetField( fld );
 998         TGeoGlobalMagField::Instance()->Lock();
 999         AliInfo("Running with the B field constructed out of GRP !");
1000       }
1001       else AliFatal("Failed to create a B field map !");
1002     }
1003     else AliFatal("B field is neither set nor constructed from GRP ! Exitig...");
1004   }
1005
1006 }
1007 //____________________________________________________________________________
1008 void AliPHOSEmbedding::InitGeometry(){
1009
1010   // Import ideal TGeo geometry and apply misalignment
1011   if (!gGeoManager) {
1012     AliGeomManager::LoadGeometry("geometry.root");
1013     if (!gGeoManager) {
1014       AliFatal("Can not load geometry");
1015     }
1016     if(!AliGeomManager::CheckSymNamesLUT("PHOS")) {
1017       AliFatal("CheckSymNamesLUT");
1018     }
1019   }
1020
1021
1022   TString detStr = "PHOS";
1023   TString loadAlObjsListOfDets = "PHOS";
1024
1025   if(AliGeomManager::GetNalignable("GRP") != 0)
1026       loadAlObjsListOfDets.Prepend("GRP "); //add alignment objects for non-sensitive modules
1027 //  AliGeomManager::ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
1028
1029   AliCDBManager::Instance()->UnloadFromCache("*/Align/*");
1030   AliCDBManager::Instance()->UnloadFromCache("GRP/Geometry/Data");
1031
1032 }
1033 //____________________________________________________________________________
1034 Float_t AliPHOSEmbedding::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1035   //Parameterization of LHC10h period
1036   //_true if neutral_
1037
1038   Double_t meanX=0;
1039   Double_t meanZ=0.;
1040   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1041               6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
1042   Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
1043
1044   AliESDEvent *event = static_cast<AliESDEvent*>(InputEvent());
1045   Double_t mf = event->GetMagneticField(); //Positive for ++ and negative for --
1046
1047   if(mf<0.){ //field --
1048     meanZ = -0.468318 ;
1049     if(charge>0)
1050       meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
1051     else
1052       meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
1053   }
1054   else{ //Field ++
1055     meanZ= -0.468318;
1056     if(charge>0)
1057       meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
1058     else
1059       meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
1060   }
1061
1062   Double_t rz=(dz-meanZ)/sz ;
1063   Double_t rx=(dx-meanX)/sx ;
1064   return TMath::Sqrt(rx*rx+rz*rz) ;
1065 }
1066