PWG4/UserTasks/PHOS_embedding/AliPHOSEmbedding.cxx

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