Event embedding tasks for PHOS are added (D.Peressounko)
authorkharlov <kharlov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Sep 2011 07:58:31 +0000 (07:58 +0000)
committerkharlov <kharlov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Sep 2011 07:58:31 +0000 (07:58 +0000)
19 files changed:
PWG4/CMakelibPWG4UserTasks.pkg
PWG4/PWG4UserTasksLinkDef.h
PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0DiffEfficiency.cxx [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0DiffEfficiency.h [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0Efficiency.cxx [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0Efficiency.h [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/AliPHOSDigitDecalibrate.cxx [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/AliPHOSDigitDecalibrate.h [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/AliPHOSEmbedding.cxx [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/AliPHOSEmbedding.h [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/Analyze.C [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/AnalyzeDiff.C [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/Config.C [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/CreateAOD.C [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/Embedding.C [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/Readme [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/rec.C [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/sim.C [new file with mode: 0644]
PWG4/UserTasks/PHOS_embedding/simrun.C [new file with mode: 0644]

index 85d43bea4ccfabe0db9a4be7bf12886bb65f4eaf..fd9679df5f6434a1dc99f23eed1bf5e6baf75f2b 100644 (file)
 set ( SRCS 
  UserTasks/PHOS_pp_pi0/AliCaloPhoton.cxx
  UserTasks/PHOS_pp_pi0/AliAnalysisTaskPi0.cxx
+ UserTasks/PHOS_embedding/AliPHOSEmbedding.cxx
+ UserTasks/PHOS_embedding/AliAnalysisTaskPi0Efficiency.cxx
+ UserTasks/PHOS_embedding/AliAnalysisTaskPi0DiffEfficiency.cxx
+ UserTasks/PHOS_embedding/AliPHOSDigitDecalibrate.cxx
  UserTasks/CaloCellQA/AliCaloCellsQA.cxx
  UserTasks/CaloCellQA/AliAnalysisTaskCaloCellsQA.cxx
  UserTasks/EmcalTasks/AliEmcalPhysicsSelection.cxx
index 7bcd3c5ad82017856ed6461d258eb50624823742..6067c9482dd80b7d75853fd57b3fd7460d18048b 100644 (file)
@@ -8,6 +8,12 @@
 #pragma link C++ class AliCaloPhoton+;
 #pragma link C++ class AliAnalysisTaskPi0+;
 
+// PHOS_embedding
+#pragma link C++ class AliPHOSEmbedding+;
+#pragma link C++ class AliAnalysisTaskPi0Efficiency+;
+#pragma link C++ class AliAnalysisTaskPi0DiffEfficiency+;
+#pragma link C++ class AliPHOSDigitDecalibrate+;
+
 // CaloCellQA
 #pragma link C++ class AliCaloCellsQA+;
 #pragma link C++ class AliAnalysisTaskCaloCellsQA+;
diff --git a/PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0DiffEfficiency.cxx b/PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0DiffEfficiency.cxx
new file mode 100644 (file)
index 0000000..049527e
--- /dev/null
@@ -0,0 +1,913 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TObjArray.h"
+#include "TF1.h"
+#include "TFile.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH2I.h"
+#include "TH3F.h"
+#include "TParticle.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TRandom.h"
+
+#include "AliAODMCParticle.h"
+#include "AliAnalysisManager.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisTaskPi0DiffEfficiency.h"
+#include "AliCaloPhoton.h"
+#include "AliPHOSGeometry.h"
+#include "AliPHOSEsdCluster.h"
+#include "AliPHOSCalibData.h"
+#include "AliAODEvent.h"
+#include "AliAODCaloCluster.h"
+#include "AliAODVertex.h"
+#include "AliESDtrackCuts.h"
+#include "AliLog.h"
+#include "AliPID.h"
+#include "AliCDBManager.h"
+#include "AliCentrality.h" 
+
+// Analysis task to fill histograms with PHOS ESD clusters and cells
+// Authors: Yuri Kharlov
+// Date   : 28.05.2009
+
+ClassImp(AliAnalysisTaskPi0DiffEfficiency)
+
+//________________________________________________________________________
+AliAnalysisTaskPi0DiffEfficiency::AliAnalysisTaskPi0DiffEfficiency(const char *name) 
+: AliAnalysisTaskSE(name),
+  fStack(0),
+  fOutputContainer(0),
+  fPHOSEvent1(0),
+  fPHOSEvent2(0),
+  fPHOSCalibData(0),
+  fNonLinCorr(0),
+  fRPfull(0),
+  fRPA(0),
+  fRPC(0),
+  fRPFar(0),
+  fRPAFar(0),
+  fRPCFar(0),
+  fCentrality(0),
+  fCenBin(0),
+  fPHOSGeo(0),
+  fEventCounter(0)
+{
+  // Constructor
+  for(Int_t i=0;i<1;i++){
+    for(Int_t j=0;j<10;j++)
+      for(Int_t k=0;k<11;k++)
+       fPHOSEvents[i][j][k]=0 ;
+  }
+  
+  // Output slots #0 write into a TH1 container
+  DefineOutput(1,TList::Class());
+
+  // Set bad channel map
+  char key[55] ;
+  for(Int_t i=0; i<6; i++){
+    sprintf(key,"PHOS_BadMap_mod%d",i) ;
+    fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
+  }
+  // Initialize the PHOS geometry
+  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
+
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPi0DiffEfficiency::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+
+  // ESD histograms
+  if(fOutputContainer != NULL){
+    delete fOutputContainer;
+  }
+  fOutputContainer = new TList();
+  fOutputContainer->SetOwner(kTRUE);
+
+  //Event selection
+  fOutputContainer->Add(new TH1F("hSelEvents","Event celection", 10,0.,10.)) ;
+
+  //vertex distribution
+  fOutputContainer->Add(new TH1F("hZvertex","Z vertex position", 50,-25.,25.)) ;
+
+  //Centrality
+  fOutputContainer->Add(new TH1F("hCentrality","Event centrality", 100,0.,100.)) ;
+
+  //QA histograms                      
+  fOutputContainer->Add(new TH2F("hCluM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
+  fOutputContainer->Add(new TH2F("hCluM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
+  fOutputContainer->Add(new TH2F("hCluM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
+
+  Int_t nM       = 500;
+  Double_t mMin  = 0.0;
+  Double_t mMax  = 1.0;
+  Int_t nPt      = 200;
+  Double_t ptMin = 0;
+  Double_t ptMax = 20;
+
+  char key[55] ;
+  for(Int_t cent=0; cent<6; cent++){
+    //Single photon
+    snprintf(key,55,"hPhotAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    ((TH1F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hPhotCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    ((TH1F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hPhotDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    ((TH1F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hPhotBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    ((TH1F*)fOutputContainer->Last())->Sumw2() ;
+
+    snprintf(key,55,"hOldMassPtAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hOldMassPtCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hOldMassPtDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hOldMassPtBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    
+    snprintf(key,55,"hNewMassPtAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hNewMassPtCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hNewMassPtDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hNewMassPtBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    
+    snprintf(key,55,"hMassPtAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hMassPtCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hMassPtDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hMassPtBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    ((TH2F*)fOutputContainer->Last())->Sumw2() ;
+    
+    //Mixed
+    snprintf(key,55,"hMiMassPtAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMiMassPtCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMiMassPtDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMiMassPtBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+
+    snprintf(key,55,"hMCMassPtAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMCMassPtCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMCMassPtDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMCMassPtBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+
+    //Single photon
+    snprintf(key,55,"hMCPhotAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    ((TH1F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hMCPhotCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    ((TH1F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hMCPhotDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    ((TH1F*)fOutputContainer->Last())->Sumw2() ;
+    snprintf(key,55,"hMCPhotBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    ((TH1F*)fOutputContainer->Last())->Sumw2() ;
+
+    
+    
+  }
+
+
+  //MC
+  for(Int_t cent=0; cent<6; cent++){
+    snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
+    snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
+    snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F("hMC_rap_eta","Rapidity eta",200,-1.,1.)) ;
+    snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
+    snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
+    snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
+    snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
+    snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
+    snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
+    snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
+    snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
+    snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
+  }
+  
+  PostData(1, fOutputContainer);
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPi0DiffEfficiency::UserExec(Option_t *) 
+{
+  // Main loop, called for each event
+  // Analyze ESD/AOD
+
+  FillHistogram("hSelEvents",0.5) ;  
+  
+  AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
+  if (!event) {
+    Printf("ERROR: Could not retrieve event");
+    PostData(1, fOutputContainer);
+    return;
+  }
+  
+  FillHistogram("hSelEvents",1.5) ;
+  AliAODHeader *header = event->GetHeader() ;
+  
+  // Checks if we have a primary vertex
+  // Get primary vertices form ESD
+  const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
+
+ // don't rely on ESD vertex, assume (0,0,0)
+  Double_t vtx0[3] ={0.,0.,0.};
+  Double_t vtx5[3] ={0.,0.,0.};
+  
+  vtx5[0] = esdVertex5->GetX();
+  vtx5[1] = esdVertex5->GetY();
+  vtx5[2] = esdVertex5->GetZ();
+  
+  
+  FillHistogram("hZvertex",esdVertex5->GetZ());
+  if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
+    PostData(1, fOutputContainer);
+    return;
+  }
+  FillHistogram("hSelEvents",2.5) ;
+
+  //Vtx class z-bin
+  //  Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
+  //  if(zvtx<0)zvtx=0 ;
+  //  if(zvtx>9)zvtx=9 ;
+  Int_t zvtx=0 ;
+
+//  fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile, 
+//                                                          //a float from 0 to 100 (or to the trigger efficiency)
+   fCentrality=header->GetZDCN2Energy() ;
+
+  if( fCentrality < 0. || fCentrality>80.){
+    PostData(1, fOutputContainer);
+    return;
+  }
+  FillHistogram("hSelEvents",3.5) ;
+  Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
+  fCenBin=0 ;
+  while(fCenBin<6 && fCentrality > bins[fCenBin+1])
+    fCenBin++ ; 
+
+  //reaction plain
+  fRPfull= header->GetZDCN1Energy() ;
+  if(fRPfull==999){ //reaction plain was not defined
+    PostData(1, fOutputContainer);
+    return;
+  } 
+
+  FillHistogram("hSelEvents",4.5) ;
+  //All event selections done
+  FillHistogram("hCentrality",fCentrality) ;
+  //Reaction plain is defined in the range (-pi/2;pi/2)
+  //We have 10 bins
+  Int_t irp=Int_t(10.*(fRPfull+TMath::PiOver2())/TMath::Pi());
+  if(irp>9)irp=9 ;
+
+  if(!fPHOSEvents[zvtx][fCenBin][irp]) 
+    fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
+  TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
+
+  // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
+  if(fEventCounter == 0) {
+    for(Int_t mod=0; mod<5; mod++) {
+      const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
+      fPHOSGeo->SetMisalMatrix(m,mod) ;
+      Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
+    }
+    fEventCounter++ ;
+  }
+
+  ProcessMC() ;
+
+  if(fPHOSEvent1){
+    fPHOSEvent1->Clear() ;
+    fPHOSEvent2->Clear() ;
+  }
+  else{
+    fPHOSEvent1 = new TClonesArray("AliCaloPhoton",200) ;
+    fPHOSEvent2 = new TClonesArray("AliCaloPhoton",200) ;
+  }
+
+  TClonesArray * clustersEmb = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
+  TClonesArray * clustersOld = event->GetCaloClusters() ;
+  TVector3 vertex(vtx0);
+  char key[55] ;
+  //Before Embedding
+  Int_t multClustOld = clustersOld->GetEntriesFast();
+  Int_t multClustEmb = clustersEmb->GetEntriesFast();
+  Int_t inPHOSold=0 ;
+  Int_t inPHOSemb=0 ;
+  for (Int_t i=0; i<multClustOld; i++) {
+    AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersOld->At(i);
+    if ( !clu->IsPHOS() || clu->E()<0.3) continue;
+
+    Bool_t survive=kFALSE ;
+    for(Int_t ii=0;(ii<multClustEmb)&&(!survive);ii++){
+       AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersEmb->At(ii);
+       survive=IsSameCluster(clu,clu2);
+    }
+
+
+    Float_t  position[3];
+    clu->GetPosition(position);
+    TVector3 global(position) ;
+    Int_t relId[4] ;
+    fPHOSGeo->GlobalPos2RelId(global,relId) ;
+    Int_t mod  = relId[0] ;
+    Int_t cellX = relId[2];
+    Int_t cellZ = relId[3] ;
+    if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) 
+      continue ;
+    if(clu->GetNCells()<3)
+      continue ;
+
+    sprintf(key,"hCluM%d",mod) ;
+    FillHistogram(key,cellX,cellZ,1.);
+
+    TLorentzVector pv1 ;
+    clu->GetMomentum(pv1 ,vtx0);
+    
+    if(inPHOSold>=fPHOSEvent1->GetSize()){
+      fPHOSEvent1->Expand(inPHOSold+50) ;
+    }
+    AliCaloPhoton * ph = new((*fPHOSEvent1)[inPHOSold]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
+    ph->SetModule(mod) ;
+    ph->SetMomV2(&pv1) ;
+    ph->SetNCells(clu->GetNCells());
+    ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
+    ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
+    if(!survive) //this cluster found in list after embedding, skipping it
+     ph->SetTagged(1) ;
+
+    inPHOSold++ ;
+  }
+
+  for (Int_t i=0; i<multClustEmb; i++) {
+    AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersEmb->At(i);
+    if ( !clu->IsPHOS() || clu->E()<0.3) continue;
+
+    Bool_t survive=kFALSE ;
+    for(Int_t ii=0;(ii<multClustOld)&&(!survive);ii++){
+       AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersOld->At(ii);
+       survive=IsSameCluster(clu,clu2);
+    }
+    
+    Float_t  position[3];
+    clu->GetPosition(position);
+    TVector3 global(position) ;
+    Int_t relId[4] ;
+    fPHOSGeo->GlobalPos2RelId(global,relId) ;
+    Int_t mod  = relId[0] ;
+    Int_t cellX = relId[2];
+    Int_t cellZ = relId[3] ;
+    if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) 
+      continue ;
+    if(clu->GetNCells()<3)
+      continue ;
+
+    sprintf(key,"hCluM%d",mod) ;
+    FillHistogram(key,cellX,cellZ,1.);
+
+    TLorentzVector pv1 ;
+    clu->GetMomentum(pv1 ,vtx0);
+    
+    if(inPHOSemb>=fPHOSEvent2->GetSize()){
+      fPHOSEvent2->Expand(inPHOSemb+50) ;
+    }
+    AliCaloPhoton * ph = new((*fPHOSEvent2)[inPHOSemb]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
+    ph->SetModule(mod) ;
+    ph->SetMomV2(&pv1) ;
+    ph->SetNCells(clu->GetNCells());
+    ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
+    ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
+    if(!survive) //this cluster found in list after embedding, skipping it
+     ph->SetTagged(1) ;
+
+    inPHOSemb++ ;
+  }
+  
+  //Single photon
+  for (Int_t i1=0; i1<inPHOSold; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
+    if(!ph1->IsTagged())
+      continue ;
+    snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
+    FillHistogram(key,ph1->Pt(),-1.) ;
+    if(ph1->IsCPVOK() ){
+      snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
+      FillHistogram(key,ph1->Pt(),-1.) ;
+    }
+    if(ph1->IsDispOK()){
+      snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
+      FillHistogram(key,ph1->Pt(),-1.) ;
+      if(ph1->IsCPVOK()){
+       snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
+        FillHistogram(key,ph1->Pt(),-1.) ;
+      }
+    } // end of loop i2
+  } // end of loop i1 
+
+  for (Int_t i1=0; i1<inPHOSemb; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
+    if(!ph1->IsTagged())
+      continue ;
+    snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
+    FillHistogram(key,ph1->Pt(),1.) ;
+    if(ph1->IsCPVOK() ){
+      snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
+      FillHistogram(key,ph1->Pt(),1.) ;
+    }
+    if(ph1->IsDispOK()){
+      snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
+      FillHistogram(key,ph1->Pt(),1.) ;
+      if(ph1->IsCPVOK()){
+       snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
+        FillHistogram(key,ph1->Pt(),1.) ;
+      }
+    } // end of loop i2
+  } // end of loop i1 
+
+
+
+  // Fill Real disribution:
+  // Disappeared clusters enter with negative contribution
+  // In addition fill control histogam with Real before embedding
+  for (Int_t i1=0; i1<inPHOSold-1; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
+    for (Int_t i2=i1+1; i2<inPHOSold; i2++) {
+      AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent1->At(i2) ;
+
+      TLorentzVector p12  = *ph1  + *ph2;
+      TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());      
+      //Fill Controll histogram: Real before embedding
+      snprintf(key,55,"hOldMassPtAll_cen%d",fCenBin) ;
+      FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
+      if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+       snprintf(key,55,"hOldMassPtCPV_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
+      }
+      if(ph1->IsDispOK() && ph2->IsDispOK()){
+       snprintf(key,55,"hOldMassPtDisp_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
+
+       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+         snprintf(key,55,"hOldMassPtBoth_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
+        }
+      }
+      
+      //Now fill mail histograms with negative contributions
+      if(!(ph1->IsTagged() || ph2->IsTagged()) )
+        continue ;
+      snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
+      FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
+      if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+       snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
+      }
+      if(ph1->IsDispOK() && ph2->IsDispOK()){
+       snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
+
+       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+         snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
+        }
+      }
+    } // end of loop i2
+  } // end of loop i1 
+
+
+  // Further fill Real disribution
+  // now with positive contribution from new clusters
+  // ass well fill controll histogram
+  for (Int_t i1=0; i1<inPHOSemb-1; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
+    for (Int_t i2=i1+1; i2<inPHOSemb; i2++) {
+      AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent2->At(i2) ;
+
+      TLorentzVector p12  = *ph1  + *ph2;
+      TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());      
+
+      // Controll histogram: Real after embedding
+      snprintf(key,55,"hNewMassPtAll_cen%d",fCenBin) ;
+      FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
+      if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+       snprintf(key,55,"hNewMassPtCPV_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt(),1) ;
+      }
+      if(ph1->IsDispOK() && ph2->IsDispOK()){
+       snprintf(key,55,"hNewMassPtDisp_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt(),1) ;
+
+       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+         snprintf(key,55,"hNewMassPtBoth_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
+        }
+      }
+     
+      //Now fill main histogamm
+      //new clusters with positive contribution
+      if(!(ph1->IsTagged() || ph2->IsTagged()) )
+        continue ;
+      snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
+      FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
+      if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+       snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt(),1) ;
+      }
+      if(ph1->IsDispOK() && ph2->IsDispOK()){
+       snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt(),1) ;
+
+       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+         snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
+        }
+      }
+    } // end of loop i2
+  } // end of loop i1 
+
+
+  //now mixed, does not really matter old or new list of clusters
+  for (Int_t i1=0; i1<inPHOSemb; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
+    for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
+      TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
+      for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
+       AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
+       TLorentzVector p12  = *ph1  + *ph2;
+       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
+       
+       snprintf(key,55,"hMiMassPtAll_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
+       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+         snprintf(key,55,"hMiMassPtCPV_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
+       }
+       if(ph1->IsDispOK() && ph2->IsDispOK()){
+         snprintf(key,55,"hMiMassPtDisp_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
+         
+         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+           snprintf(key,55,"hMiMassPtBoth_cen%d",fCenBin) ;
+           FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
+         }
+       }
+      } // end of loop i2
+    }
+  } // end of loop i1
+  
+  
+  //Now we either add current events to stack or remove
+  //If no photons in current event - no need to add it to mixed
+  if(fPHOSEvent2->GetEntriesFast()>0){
+    prevPHOS->AddFirst(fPHOSEvent2) ;
+    fPHOSEvent2=0;
+    delete fPHOSEvent1;
+    fPHOSEvent1=0;
+    if(prevPHOS->GetSize()>100){//Remove redundant events
+      TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
+      prevPHOS->RemoveLast() ;
+      delete tmp ;
+    }
+  }
+  // Post output data.
+  PostData(1, fOutputContainer);
+  fEventCounter++;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPi0DiffEfficiency::Terminate(Option_t *)
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+  
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskPi0DiffEfficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
+{
+  //Check if this channel belogs to the good ones
+
+  if(strcmp(det,"PHOS")==0){
+    if(mod>5 || mod<1){
+      AliError(Form("No bad map for PHOS module %d ",mod)) ;
+      return kTRUE ;
+    }
+    if(!fPHOSBadMap[mod]){
+      AliError(Form("No Bad map for PHOS module %d",mod)) ;
+      return kTRUE ;
+    }
+    if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
+      return kFALSE ;
+    else
+      return kTRUE ;
+  }
+  else{
+    AliError(Form("Can not find bad channels for detector %s ",det)) ;
+  }
+  return kTRUE ;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x)const{
+  //FillHistogram
+  TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
+  if(tmpI){
+    tmpI->Fill(x) ;
+    return ;
+  }
+  TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
+  if(tmpF){
+    tmpF->Fill(x) ;
+    return ;
+  }
+  TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
+  if(tmpD){
+    tmpD->Fill(x) ;
+    return ;
+  }
+  AliInfo(Form("can not find histogram <%s> ",key)) ;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
+  //FillHistogram
+  TObject * tmp = fOutputContainer->FindObject(key) ;
+  if(!tmp){
+    AliInfo(Form("can not find histogram <%s> ",key)) ;
+    return ;
+  }
+  if(tmp->IsA() == TClass::GetClass("TH1F")){
+    ((TH1F*)tmp)->Fill(x,y) ;
+    return ;
+  }
+  if(tmp->IsA() == TClass::GetClass("TH2F")){
+    ((TH2F*)tmp)->Fill(x,y) ;
+    return ;
+  }
+  AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
+  //Fills 1D histograms with key
+  TObject * tmp = fOutputContainer->FindObject(key) ;
+  if(!tmp){
+    AliInfo(Form("can not find histogram <%s> ",key)) ;
+    return ;
+  }
+  if(tmp->IsA() == TClass::GetClass("TH2F")){
+    ((TH2F*)tmp)->Fill(x,y,z) ;
+    return ;
+  }
+  if(tmp->IsA() == TClass::GetClass("TH3F")){
+    ((TH3F*)tmp)->Fill(x,y,z) ;
+    return ;
+  }
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskPi0DiffEfficiency::TestLambda(Double_t l1,Double_t l2){
+  Double_t l1Mean=1.22 ;
+  Double_t l2Mean=2.0 ;
+  Double_t l1Sigma=0.42 ;
+  Double_t l2Sigma=0.71 ;
+  Double_t c=-0.59 ;
+  Double_t R2=(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma+(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma-c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
+  return (R2<9.) ;
+  
+}
+//___________________________________________________________________________
+void AliAnalysisTaskPi0DiffEfficiency::ProcessMC(){
+  //fill histograms for efficiensy etc. calculation
+  const Double_t rcut = 1. ; //cut for primary particles
+  //---------First pi0/eta-----------------------------
+  char partName[10] ;
+  char hkey[55] ;
+
+  AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
+  TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
+  for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
+     AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(i);
+    if(particle->GetPdgCode() == 111)
+      sprintf(partName,"pi0") ;
+    else
+      if(particle->GetPdgCode() == 221)
+        sprintf(partName,"eta") ;
+      else
+        if(particle->GetPdgCode() == 22)
+           sprintf(partName,"gamma") ;
+       else
+           continue ;
+
+    //Primary particle
+    Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
+    if(r >rcut)
+      continue ;
+
+    Double_t pt = particle->Pt() ;
+    //Total number of pi0 with creation radius <1 cm
+    sprintf(hkey,"hMC_all_%s_cen%d",partName,fCenBin) ;
+    FillHistogram(hkey,pt) ;
+    if(TMath::Abs(particle->Y())<0.12){
+      sprintf(hkey,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
+      FillHistogram(hkey,pt) ;
+    }
+
+    sprintf(hkey,"hMC_rap_%s_cen%d",partName,fCenBin) ;
+    FillHistogram(hkey,particle->Y()) ;
+    
+    Double_t phi=particle->Phi() ;
+    while(phi<0.)phi+=TMath::TwoPi() ;
+    while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
+    sprintf(hkey,"hMC_phi_%s_cen%d",partName,fCenBin) ;
+    FillHistogram(hkey,phi) ;
+   
+
+    //Check if one of photons converted
+    if(particle->GetNDaughters()!=2)
+     continue ; //Do not account Dalitz decays
+
+/*
+    TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
+    TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
+    //Number of pi0s decayed into acceptance
+    Int_t mod1,mod2 ;
+    Double_t x=0.,z=0. ;
+    Bool_t hitPHOS1 = fPHOSGeo->ImpactOnEmc(gamma1, mod1, z,x) ;
+    Bool_t hitPHOS2 = fPHOSGeo->ImpactOnEmc(gamma2, mod2, z,x) ;
+
+    Bool_t goodPair=kFALSE ;
+    if( hitPHOS1 && hitPHOS2){
+      sprintf(hkey,"hMC_PHOSacc_%s",partName) ;
+      FillHistogram(hkey,pt) ;
+      goodPair=kTRUE ;
+    }
+
+*/
+  }
+  //Now calculate "Real" distribution of clusters with primary
+  TClonesArray cluPrim("AliCaloPhoton",200) ; //clusters with primary
+  Int_t multClust = event->GetNumberOfCaloClusters();
+  Int_t inPHOS=0 ;
+  Double_t vtx0[3] = {0,0,0}; 
+  for (Int_t i=0; i<multClust; i++) {
+    AliVCluster *clu = event->GetCaloCluster(i);
+    if ( !clu->IsPHOS() || clu->E()<0.3) continue;
+    if(clu->GetLabel()<0) continue ;
+
+    Float_t  position[3];
+    clu->GetPosition(position);
+    TVector3 global(position) ;
+    Int_t relId[4] ;
+    fPHOSGeo->GlobalPos2RelId(global,relId) ;
+    Int_t mod  = relId[0] ;
+    Int_t cellX = relId[2];
+    Int_t cellZ = relId[3] ;
+    if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) 
+      continue ;
+    if(clu->GetNCells()<3)
+      continue ;
+
+    TLorentzVector pv1 ;
+    clu->GetMomentum(pv1 ,vtx0);
+    
+    if(inPHOS>=cluPrim.GetSize()){
+      cluPrim.Expand(inPHOS+50) ;
+    }
+    AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
+    //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
+    ph->SetModule(mod) ;
+    ph->SetMomV2(&pv1) ;
+    ph->SetNCells(clu->GetNCells());
+    ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
+    ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
+
+    inPHOS++ ;
+
+  }
+  
+  //Single photon
+  char key[55] ;
+  for (Int_t i1=0; i1<inPHOS; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
+    snprintf(key,55,"hMCPhotAll_cen%d",fCenBin) ;
+    FillHistogram(key,ph1->Pt()) ;
+    if(ph1->IsCPVOK() ){
+      snprintf(key,55,"hMCPhotCPV_cen%d",fCenBin) ;
+      FillHistogram(key,ph1->Pt()) ;
+      
+    }
+    if(ph1->IsDispOK()){
+      snprintf(key,55,"hMCPhotDisp_cen%d",fCenBin) ;
+      FillHistogram(key,ph1->Pt()) ;
+      if(ph1->IsCPVOK()){
+       snprintf(key,55,"hMCPhotBoth_cen%d",fCenBin) ;
+       FillHistogram(key,ph1->Pt()) ;
+      }
+    } // end of loop i2
+  } // end of loop i1 
+
+  // Fill Real disribution
+  for (Int_t i1=0; i1<inPHOS-1; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
+    for (Int_t i2=i1+1; i2<inPHOS; i2++) {
+      AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
+      TLorentzVector p12  = *ph1  + *ph2;
+      TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());      
+  
+      snprintf(key,55,"hMCMassPtAll_cen%d",fCenBin) ;
+      FillHistogram(key,p12.M() ,p12.Pt()) ;
+      if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+       snprintf(key,55,"hMCMassPtCPV_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt()) ;
+      }
+      if(ph1->IsDispOK() && ph2->IsDispOK()){
+       snprintf(key,55,"hMCMassPtDisp_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt()) ;
+
+       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+         snprintf(key,55,"hMCMassPtBoth_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt()) ;
+        }
+      }
+    } // end of loop i2
+  } // end of loop i1
+
+  
+}
+//___________________________________________________________________________
+Bool_t AliAnalysisTaskPi0DiffEfficiency::IsSameCluster(AliAODCaloCluster * c1, AliAODCaloCluster * c2)const{
+ //Compare clusters before and after embedding
+ //clusters are the same if 
+ // - Energy changed less than 0.1%  (numerical accuracy in reconstruction)
+ // - lists of digits are the same
+  
+ if(c1->GetNCells() != c2->GetNCells())
+   return kFALSE ;
+ if(TMath::Abs(c1->E()-c2->E())>0.001*c1->E())
+   return kFALSE ;
+
+ UShort_t *list1 = c1->GetCellsAbsId() ; 
+ UShort_t *list2 = c2->GetCellsAbsId() ; 
+ for(Int_t i=0; i<c1->GetNCells(); i++){
+  if(list1[i] != list2[i])
+    return kFALSE ;
+ }
+ return kTRUE ; 
+  
+}
+
+
+
diff --git a/PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0DiffEfficiency.h b/PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0DiffEfficiency.h
new file mode 100644 (file)
index 0000000..63c0047
--- /dev/null
@@ -0,0 +1,79 @@
+#ifndef ALIANALYSISTASKPI0DIFFEFFICIENCY_H\r
+#define ALIANALYSISTASKPI0DIFFEFFICIENCY_H\r
+\r
+// example of an analysis task creating a p_t spectrum\r
+// Authors: Panos Cristakoglou, Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing\r
+\r
+class TObjArray;\r
+class TH1F;\r
+class TH2I;\r
+class TH2F;\r
+class TH3F;\r
+class TF1 ;\r
+class AliStack ;\r
+class AliESDtrackCuts;\r
+class AliPHOSGeometry;\r
+class AliAODEvent ;\r
+class AliPHOSCalibData;\r
+class AliAODTrack ;\r
+class AliAODCaloCluster ;\r
+\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class AliAnalysisTaskPi0DiffEfficiency : public AliAnalysisTaskSE {\r
+public:\r
+  AliAnalysisTaskPi0DiffEfficiency(const char *name = "AliAnalysisTaskPi0DiffEfficiency");\r
+  virtual ~AliAnalysisTaskPi0DiffEfficiency() {}\r
+  \r
+  virtual void   UserCreateOutputObjects();\r
+  virtual void   UserExec(Option_t *option);\r
+  virtual void   Terminate(Option_t *);\r
+  void SetPHOSBadMap(Int_t mod,TH2I * h)\r
+  {\r
+    if(fPHOSBadMap[mod]) delete fPHOSBadMap[mod] ;\r
+    fPHOSBadMap[mod]=new TH2I(*h) ;\r
+    printf("Set %s \n",fPHOSBadMap[mod]->GetName());\r
+  }\r
+  \r
+private:\r
+  AliAnalysisTaskPi0DiffEfficiency(const AliAnalysisTaskPi0DiffEfficiency&); // not implemented\r
+  AliAnalysisTaskPi0DiffEfficiency& operator=(const AliAnalysisTaskPi0DiffEfficiency&); // not implemented\r
+  Bool_t IsSameCluster(AliAODCaloCluster * c1,AliAODCaloCluster * c2)const ;\r
+  Bool_t IsGoodChannel(const char * det, Int_t mod,Int_t ix, Int_t iz); //Use addisional bad map for PHOS\r
+  void FillHistogram(const char * key,Double_t x) const ; //Fill 1D histogram witn name key\r
+  void FillHistogram(const char * key,Double_t x, Double_t y) const ; //Fill 2D histogram witn name key\r
+  void FillHistogram(const char * key,Double_t x, Double_t y, Double_t z) const ; //Fill 3D histogram witn name key\r
+  Bool_t TestLambda(Double_t l1,Double_t l2) ;  //Evaluate Dispersion cuts for photons\r
+  void ProcessMC() ;\r
\r
+private:\r
+  AliStack * fStack ;              //stack of MC tracks\r
+  TList * fOutputContainer;        //final histogram container\r
+  TList * fPHOSEvents[1][10][11] ; //Containers for events with PHOS photons\r
+  TClonesArray * fPHOSEvent1 ;      //PHOS photons in current event\r
+  TClonesArray * fPHOSEvent2 ;      //PHOS photons in current event\r
+  AliPHOSCalibData *fPHOSCalibData; // PHOS calibration object\r
+  TF1 *fNonLinCorr;          // Non-linearity correction\r
\r
+  //Reaction plain for v2\r
+  Float_t fRPfull ; //!Reaction plain calculated with full TPC \r
+  Float_t fRPA ;    //!Reaction plain calculated with A-side TPC: eta>0.15 \r
+  Float_t fRPC ;    //!Reaction plain calculated with C-side TPC: eta<-0.15\r
+  Float_t fRPFar ;  //!Reaction plain calculated with TPC: eta>0.6 \r
+  Float_t fRPAFar ; //!Reaction plain calculated with A-side TPC: eta>0.6 \r
+  Float_t fRPCFar ; //!Reaction plain calculated with C-side TPC: eta<-0.6\r
+\r
+  Float_t fCentrality ; //!Centrality of the currecnt event\r
+\r
+  Int_t fCenBin ;       //! Current centrality bin\r
+\r
+  TH2I *fPHOSBadMap[6] ;    //Container for PHOS bad channels map\r
+\r
+  AliPHOSGeometry  *fPHOSGeo;  //! PHOS geometry\r
+  Int_t fEventCounter;         // number of analyzed events\r
+\r
+  ClassDef(AliAnalysisTaskPi0DiffEfficiency, 1); // PHOS analysis task\r
+};\r
+\r
+#endif\r
diff --git a/PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0Efficiency.cxx b/PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0Efficiency.cxx
new file mode 100644 (file)
index 0000000..3236b3c
--- /dev/null
@@ -0,0 +1,691 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TObjArray.h"
+#include "TF1.h"
+#include "TFile.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH2I.h"
+#include "TH3F.h"
+#include "TParticle.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TRandom.h"
+
+#include "AliAODMCParticle.h"
+#include "AliAnalysisManager.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisTaskPi0Efficiency.h"
+#include "AliCaloPhoton.h"
+#include "AliPHOSGeometry.h"
+#include "AliPHOSEsdCluster.h"
+#include "AliPHOSCalibData.h"
+#include "AliAODEvent.h"
+#include "AliAODCaloCluster.h"
+#include "AliAODVertex.h"
+#include "AliESDtrackCuts.h"
+#include "AliLog.h"
+#include "AliPID.h"
+#include "AliCDBManager.h"
+#include "AliCentrality.h" 
+
+// Analysis task to fill histograms with PHOS ESD clusters and cells
+// Authors: Yuri Kharlov
+// Date   : 28.05.2009
+
+ClassImp(AliAnalysisTaskPi0Efficiency)
+
+//________________________________________________________________________
+AliAnalysisTaskPi0Efficiency::AliAnalysisTaskPi0Efficiency(const char *name) 
+: AliAnalysisTaskSE(name),
+  fStack(0),
+  fOutputContainer(0),
+  fPHOSEvent(0),
+  fPHOSCalibData(0),
+  fNonLinCorr(0),
+  fRPfull(0),
+  fRPA(0),
+  fRPC(0),
+  fRPFar(0),
+  fRPAFar(0),
+  fRPCFar(0),
+  fCentrality(0),
+  fCenBin(0),
+  fPHOSGeo(0),
+  fEventCounter(0)
+{
+  // Constructor
+  for(Int_t i=0;i<1;i++){
+    for(Int_t j=0;j<10;j++)
+      for(Int_t k=0;k<11;k++)
+       fPHOSEvents[i][j][k]=0 ;
+  }
+  
+  // Output slots #0 write into a TH1 container
+  DefineOutput(1,TList::Class());
+
+  // Set bad channel map
+  char key[55] ;
+  for(Int_t i=0; i<6; i++){
+    sprintf(key,"PHOS_BadMap_mod%d",i) ;
+    fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
+  }
+  // Initialize the PHOS geometry
+  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
+
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPi0Efficiency::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+
+  // ESD histograms
+  if(fOutputContainer != NULL){
+    delete fOutputContainer;
+  }
+  fOutputContainer = new TList();
+  fOutputContainer->SetOwner(kTRUE);
+
+  //Event selection
+  fOutputContainer->Add(new TH1F("hSelEvents","Event celection", 10,0.,10.)) ;
+
+  //vertex distribution
+  fOutputContainer->Add(new TH1F("hZvertex","Z vertex position", 50,-25.,25.)) ;
+
+  //Centrality
+  fOutputContainer->Add(new TH1F("hCentrality","Event centrality", 100,0.,100.)) ;
+
+  //QA histograms                      
+  fOutputContainer->Add(new TH2F("hCluM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
+  fOutputContainer->Add(new TH2F("hCluM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
+  fOutputContainer->Add(new TH2F("hCluM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
+
+  Int_t nM       = 500;
+  Double_t mMin  = 0.0;
+  Double_t mMax  = 1.0;
+  Int_t nPt      = 200;
+  Double_t ptMin = 0;
+  Double_t ptMax = 20;
+
+  char key[55] ;
+  for(Int_t cent=0; cent<6; cent++){
+    //Single photon
+    snprintf(key,55,"hPhotAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    snprintf(key,55,"hPhotCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    snprintf(key,55,"hPhotDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    snprintf(key,55,"hPhotBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    
+    snprintf(key,55,"hMassPtAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMassPtCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMassPtDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMassPtBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    
+    //Mixed
+    snprintf(key,55,"hMiMassPtAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMiMassPtCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMiMassPtDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMiMassPtBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+
+    snprintf(key,55,"hMCMassPtAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMCMassPtCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMCMassPtDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMCMassPtBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
+
+    //Single photon
+    snprintf(key,55,"hMCPhotAll_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMCPhotCPV_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMCPhotDisp_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+    snprintf(key,55,"hMCPhotBoth_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
+
+    
+    
+  }
+
+
+  //MC
+  for(Int_t cent=0; cent<6; cent++){
+    snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
+    snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
+    snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F("hMC_rap_eta","Rapidity eta",200,-1.,1.)) ;
+    snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
+    snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
+    snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
+    snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
+    snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
+    snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
+    snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
+    snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
+    snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
+    fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
+  }
+  
+  PostData(1, fOutputContainer);
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPi0Efficiency::UserExec(Option_t *) 
+{
+  // Main loop, called for each event
+  // Analyze ESD/AOD
+
+  FillHistogram("hSelEvents",0.5) ;  
+  
+  AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
+  if (!event) {
+    Printf("ERROR: Could not retrieve event");
+    PostData(1, fOutputContainer);
+    return;
+  }
+  
+  FillHistogram("hSelEvents",1.5) ;
+  AliAODHeader *header = event->GetHeader() ;
+  
+  // Checks if we have a primary vertex
+  // Get primary vertices form ESD
+  const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
+
+ // don't rely on ESD vertex, assume (0,0,0)
+  Double_t vtx0[3] ={0.,0.,0.};
+  Double_t vtx5[3] ={0.,0.,0.};
+  
+  vtx5[0] = esdVertex5->GetX();
+  vtx5[1] = esdVertex5->GetY();
+  vtx5[2] = esdVertex5->GetZ();
+  
+  
+  FillHistogram("hZvertex",esdVertex5->GetZ());
+  if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
+    PostData(1, fOutputContainer);
+    return;
+  }
+  FillHistogram("hSelEvents",2.5) ;
+
+  //Vtx class z-bin
+  //  Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
+  //  if(zvtx<0)zvtx=0 ;
+  //  if(zvtx>9)zvtx=9 ;
+  Int_t zvtx=0 ;
+
+//  fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile, 
+//                                                          //a float from 0 to 100 (or to the trigger efficiency)
+   fCentrality=header->GetZDCN2Energy() ;
+
+  if( fCentrality < 0. || fCentrality>80.){
+    PostData(1, fOutputContainer);
+    return;
+  }
+  FillHistogram("hSelEvents",3.5) ;
+  Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
+  fCenBin=0 ;
+  while(fCenBin<6 && fCentrality > bins[fCenBin+1])
+    fCenBin++ ; 
+
+  //reaction plain
+  fRPfull= header->GetZDCN1Energy() ;
+  if(fRPfull==999){ //reaction plain was not defined
+    PostData(1, fOutputContainer);
+    return;
+  } 
+
+  FillHistogram("hSelEvents",4.5) ;
+  //All event selections done
+  FillHistogram("hCentrality",fCentrality) ;
+  //Reaction plain is defined in the range (-pi/2;pi/2)
+  //We have 10 bins
+  Int_t irp=Int_t(10.*(fRPfull+TMath::PiOver2())/TMath::Pi());
+  if(irp>9)irp=9 ;
+
+  if(!fPHOSEvents[zvtx][fCenBin][irp]) 
+    fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
+  TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
+
+  // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
+  if(fEventCounter == 0) {
+    for(Int_t mod=0; mod<5; mod++) {
+      const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
+      fPHOSGeo->SetMisalMatrix(m,mod) ;
+      Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
+    }
+    fEventCounter++ ;
+  }
+
+  ProcessMC() ;
+
+  if(fPHOSEvent)
+    fPHOSEvent->Clear() ;
+  else
+    fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
+
+
+  char key[55] ;
+  Int_t inPHOS=0 ;
+  TVector3 vertex(vtx0);
+  TClonesArray * clusters = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
+  Int_t multClust = clusters->GetEntriesFast();
+  for (Int_t i=0; i<multClust; i++) {
+    AliVCluster *clu = (AliVCluster*) clusters->At(i);
+    if ( !clu->IsPHOS() || clu->E()<0.3) continue;
+
+    Float_t  position[3];
+    clu->GetPosition(position);
+    TVector3 global(position) ;
+    Int_t relId[4] ;
+    fPHOSGeo->GlobalPos2RelId(global,relId) ;
+    Int_t mod  = relId[0] ;
+    Int_t cellX = relId[2];
+    Int_t cellZ = relId[3] ;
+    if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) 
+      continue ;
+    if(clu->GetNCells()<3)
+      continue ;
+
+    sprintf(key,"hCluM%d",mod) ;
+    FillHistogram(key,cellX,cellZ,1.);
+
+    TLorentzVector pv1 ;
+    clu->GetMomentum(pv1 ,vtx0);
+    
+    if(inPHOS>=fPHOSEvent->GetSize()){
+      fPHOSEvent->Expand(inPHOS+50) ;
+    }
+    new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
+    AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
+    ph->SetModule(mod) ;
+    ph->SetMomV2(&pv1) ;
+    ph->SetNCells(clu->GetNCells());
+    ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
+    ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
+
+    inPHOS++ ;
+  }
+  //Single photon
+  for (Int_t i1=0; i1<inPHOS; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
+    snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
+    FillHistogram(key,ph1->Pt()) ;
+    if(ph1->IsCPVOK() ){
+      snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
+      FillHistogram(key,ph1->Pt()) ;
+    }
+    if(ph1->IsDispOK()){
+      snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
+      FillHistogram(key,ph1->Pt()) ;
+      if(ph1->IsCPVOK()){
+       snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
+       FillHistogram(key,ph1->Pt()) ;
+      }
+    } // end of loop i2
+  } // end of loop i1 
+
+  // Fill Real disribution
+  for (Int_t i1=0; i1<inPHOS-1; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
+    for (Int_t i2=i1+1; i2<inPHOS; i2++) {
+      AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
+      TLorentzVector p12  = *ph1  + *ph2;
+      TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());      
+
+      snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
+      FillHistogram(key,p12.M() ,p12.Pt()) ;
+      if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+       snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt()) ;
+      }
+      if(ph1->IsDispOK() && ph2->IsDispOK()){
+       snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt()) ;
+
+       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+         snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt()) ;
+        }
+      }
+    } // end of loop i2
+  } // end of loop i1 
+
+  //now mixed
+  for (Int_t i1=0; i1<inPHOS; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
+    for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
+      TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
+      for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
+       AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
+       TLorentzVector p12  = *ph1  + *ph2;
+       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
+       
+       snprintf(key,55,"hMiMassPtAll_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt()) ;
+       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+         snprintf(key,55,"hMiMassPtCPV_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt()) ;
+       }
+       if(ph1->IsDispOK() && ph2->IsDispOK()){
+         snprintf(key,55,"hMiMassPtDisp_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt()) ;
+         
+         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+           snprintf(key,55,"hMiMassPtBoth_cen%d",fCenBin) ;
+           FillHistogram(key,p12.M() ,p12.Pt()) ;
+         }
+       }
+      } // end of loop i2
+    }
+  } // end of loop i1
+  
+  
+  //Now we either add current events to stack or remove
+  //If no photons in current event - no need to add it to mixed
+  if(fPHOSEvent->GetEntriesFast()>0){
+    prevPHOS->AddFirst(fPHOSEvent) ;
+    fPHOSEvent=0;
+    if(prevPHOS->GetSize()>100){//Remove redundant events
+      TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
+      prevPHOS->RemoveLast() ;
+      delete tmp ;
+    }
+  }
+  // Post output data.
+  PostData(1, fOutputContainer);
+  fEventCounter++;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPi0Efficiency::Terminate(Option_t *)
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+  
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskPi0Efficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
+{
+  //Check if this channel belogs to the good ones
+
+  if(strcmp(det,"PHOS")==0){
+    if(mod>5 || mod<1){
+      AliError(Form("No bad map for PHOS module %d ",mod)) ;
+      return kTRUE ;
+    }
+    if(!fPHOSBadMap[mod]){
+      AliError(Form("No Bad map for PHOS module %d",mod)) ;
+      return kTRUE ;
+    }
+    if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
+      return kFALSE ;
+    else
+      return kTRUE ;
+  }
+  else{
+    AliError(Form("Can not find bad channels for detector %s ",det)) ;
+  }
+  return kTRUE ;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x)const{
+  //FillHistogram
+  TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
+  if(tmpI){
+    tmpI->Fill(x) ;
+    return ;
+  }
+  TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
+  if(tmpF){
+    tmpF->Fill(x) ;
+    return ;
+  }
+  TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
+  if(tmpD){
+    tmpD->Fill(x) ;
+    return ;
+  }
+  AliInfo(Form("can not find histogram <%s> ",key)) ;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
+  //FillHistogram
+  TObject * tmp = fOutputContainer->FindObject(key) ;
+  if(!tmp){
+    AliInfo(Form("can not find histogram <%s> ",key)) ;
+    return ;
+  }
+  if(tmp->IsA() == TClass::GetClass("TH1F")){
+    ((TH1F*)tmp)->Fill(x,y) ;
+    return ;
+  }
+  if(tmp->IsA() == TClass::GetClass("TH2F")){
+    ((TH2F*)tmp)->Fill(x,y) ;
+    return ;
+  }
+  AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
+  //Fills 1D histograms with key
+  TObject * tmp = fOutputContainer->FindObject(key) ;
+  if(!tmp){
+    AliInfo(Form("can not find histogram <%s> ",key)) ;
+    return ;
+  }
+  if(tmp->IsA() == TClass::GetClass("TH2F")){
+    ((TH2F*)tmp)->Fill(x,y,z) ;
+    return ;
+  }
+  if(tmp->IsA() == TClass::GetClass("TH3F")){
+    ((TH3F*)tmp)->Fill(x,y,z) ;
+    return ;
+  }
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskPi0Efficiency::TestLambda(Double_t l1,Double_t l2){
+  Double_t l1Mean=1.22 ;
+  Double_t l2Mean=2.0 ;
+  Double_t l1Sigma=0.42 ;
+  Double_t l2Sigma=0.71 ;
+  Double_t c=-0.59 ;
+  Double_t R2=(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma+(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma-c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
+  return (R2<9.) ;
+  
+}
+//___________________________________________________________________________
+void AliAnalysisTaskPi0Efficiency::ProcessMC(){
+  //fill histograms for efficiensy etc. calculation
+  const Double_t rcut = 1. ; //cut for primary particles
+  //---------First pi0/eta-----------------------------
+  char partName[10] ;
+  char hkey[55] ;
+
+  AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
+  TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
+  for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
+     AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(i);
+    if(particle->GetPdgCode() == 111)
+      sprintf(partName,"pi0") ;
+    else
+      if(particle->GetPdgCode() == 221)
+        sprintf(partName,"eta") ;
+      else
+        if(particle->GetPdgCode() == 22)
+           sprintf(partName,"gamma") ;
+       else
+           continue ;
+
+    //Primary particle
+    Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
+    if(r >rcut)
+      continue ;
+
+    Double_t pt = particle->Pt() ;
+    //Total number of pi0 with creation radius <1 cm
+    sprintf(hkey,"hMC_all_%s_cen%d",partName,fCenBin) ;
+    FillHistogram(hkey,pt) ;
+    if(TMath::Abs(particle->Y())<0.12){
+      sprintf(hkey,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
+      FillHistogram(hkey,pt) ;
+    }
+
+    sprintf(hkey,"hMC_rap_%s_cen%d",partName,fCenBin) ;
+    FillHistogram(hkey,particle->Y()) ;
+    
+    Double_t phi=particle->Phi() ;
+    while(phi<0.)phi+=TMath::TwoPi() ;
+    while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
+    sprintf(hkey,"hMC_phi_%s_cen%d",partName,fCenBin) ;
+    FillHistogram(hkey,phi) ;
+   
+
+    //Check if one of photons converted
+    if(particle->GetNDaughters()!=2)
+     continue ; //Do not account Dalitz decays
+
+/*
+    TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
+    TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
+    //Number of pi0s decayed into acceptance
+    Int_t mod1,mod2 ;
+    Double_t x=0.,z=0. ;
+    Bool_t hitPHOS1 = fPHOSGeo->ImpactOnEmc(gamma1, mod1, z,x) ;
+    Bool_t hitPHOS2 = fPHOSGeo->ImpactOnEmc(gamma2, mod2, z,x) ;
+
+    Bool_t goodPair=kFALSE ;
+    if( hitPHOS1 && hitPHOS2){
+      sprintf(hkey,"hMC_PHOSacc_%s",partName) ;
+      FillHistogram(hkey,pt) ;
+      goodPair=kTRUE ;
+    }
+
+*/
+  }
+  //Now calculate "Real" distribution of clusters with primary
+  TClonesArray cluPrim("AliCaloPhoton",200) ; //clusters with primary
+  Int_t multClust = event->GetNumberOfCaloClusters();
+  Int_t inPHOS=0 ;
+  Double_t vtx0[3] = {0,0,0}; 
+  for (Int_t i=0; i<multClust; i++) {
+    AliVCluster *clu = event->GetCaloCluster(i);
+    if ( !clu->IsPHOS() || clu->E()<0.3) continue;
+    if(clu->GetLabel()<0) continue ;
+
+    Float_t  position[3];
+    clu->GetPosition(position);
+    TVector3 global(position) ;
+    Int_t relId[4] ;
+    fPHOSGeo->GlobalPos2RelId(global,relId) ;
+    Int_t mod  = relId[0] ;
+    Int_t cellX = relId[2];
+    Int_t cellZ = relId[3] ;
+    if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) 
+      continue ;
+    if(clu->GetNCells()<3)
+      continue ;
+
+    TLorentzVector pv1 ;
+    clu->GetMomentum(pv1 ,vtx0);
+    
+    if(inPHOS>=cluPrim.GetSize()){
+      cluPrim.Expand(inPHOS+50) ;
+    }
+    AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
+    //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
+    ph->SetModule(mod) ;
+    ph->SetMomV2(&pv1) ;
+    ph->SetNCells(clu->GetNCells());
+    ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
+    ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
+
+    inPHOS++ ;
+
+  }
+  
+  //Single photon
+  char key[55] ;
+  for (Int_t i1=0; i1<inPHOS; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
+    snprintf(key,55,"hMCPhotAll_cen%d",fCenBin) ;
+    FillHistogram(key,ph1->Pt()) ;
+    if(ph1->IsCPVOK() ){
+      snprintf(key,55,"hMCPhotCPV_cen%d",fCenBin) ;
+      FillHistogram(key,ph1->Pt()) ;
+      
+    }
+    if(ph1->IsDispOK()){
+      snprintf(key,55,"hMCPhotDisp_cen%d",fCenBin) ;
+      FillHistogram(key,ph1->Pt()) ;
+      if(ph1->IsCPVOK()){
+       snprintf(key,55,"hMCPhotBoth_cen%d",fCenBin) ;
+       FillHistogram(key,ph1->Pt()) ;
+      }
+    } // end of loop i2
+  } // end of loop i1 
+
+  // Fill Real disribution
+  for (Int_t i1=0; i1<inPHOS-1; i1++) {
+    AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
+    for (Int_t i2=i1+1; i2<inPHOS; i2++) {
+      AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
+      TLorentzVector p12  = *ph1  + *ph2;
+      TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());      
+  
+      snprintf(key,55,"hMCMassPtAll_cen%d",fCenBin) ;
+      FillHistogram(key,p12.M() ,p12.Pt()) ;
+      if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+       snprintf(key,55,"hMCMassPtCPV_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt()) ;
+      }
+      if(ph1->IsDispOK() && ph2->IsDispOK()){
+       snprintf(key,55,"hMCMassPtDisp_cen%d",fCenBin) ;
+       FillHistogram(key,p12.M() ,p12.Pt()) ;
+
+       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
+         snprintf(key,55,"hMCMassPtBoth_cen%d",fCenBin) ;
+         FillHistogram(key,p12.M() ,p12.Pt()) ;
+        }
+      }
+    } // end of loop i2
+  } // end of loop i1
+
+  
+}
+
+
diff --git a/PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0Efficiency.h b/PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0Efficiency.h
new file mode 100644 (file)
index 0000000..2203dd0
--- /dev/null
@@ -0,0 +1,76 @@
+#ifndef ALIANALYSISTASKPI0EFFICIENCY_H\r
+#define ALIANALYSISTASKPI0EFFICIENCY_H\r
+\r
+// example of an analysis task creating a p_t spectrum\r
+// Authors: Panos Cristakoglou, Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing\r
+\r
+class TObjArray;\r
+class TH1F;\r
+class TH2I;\r
+class TH2F;\r
+class TH3F;\r
+class TF1 ;\r
+class AliStack ;\r
+class AliESDtrackCuts;\r
+class AliPHOSGeometry;\r
+class AliAODEvent ;\r
+class AliPHOSCalibData;\r
+class AliAODTrack ;\r
+\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class AliAnalysisTaskPi0Efficiency : public AliAnalysisTaskSE {\r
+public:\r
+  AliAnalysisTaskPi0Efficiency(const char *name = "AliAnalysisTaskPi0Efficiency");\r
+  virtual ~AliAnalysisTaskPi0Efficiency() {}\r
+  \r
+  virtual void   UserCreateOutputObjects();\r
+  virtual void   UserExec(Option_t *option);\r
+  virtual void   Terminate(Option_t *);\r
+  void SetPHOSBadMap(Int_t mod,TH2I * h)\r
+  {\r
+    if(fPHOSBadMap[mod]) delete fPHOSBadMap[mod] ;\r
+    fPHOSBadMap[mod]=new TH2I(*h) ;\r
+    printf("Set %s \n",fPHOSBadMap[mod]->GetName());\r
+  }\r
+  \r
+private:\r
+  AliAnalysisTaskPi0Efficiency(const AliAnalysisTaskPi0Efficiency&); // not implemented\r
+  AliAnalysisTaskPi0Efficiency& operator=(const AliAnalysisTaskPi0Efficiency&); // not implemented\r
+  Bool_t IsGoodChannel(const char * det, Int_t mod,Int_t ix, Int_t iz); //Use addisional bad map for PHOS\r
+  void FillHistogram(const char * key,Double_t x) const ; //Fill 1D histogram witn name key\r
+  void FillHistogram(const char * key,Double_t x, Double_t y) const ; //Fill 2D histogram witn name key\r
+  void FillHistogram(const char * key,Double_t x, Double_t y, Double_t z) const ; //Fill 3D histogram witn name key\r
+  Bool_t TestLambda(Double_t l1,Double_t l2) ;  //Evaluate Dispersion cuts for photons\r
+  void ProcessMC() ;\r
\r
+private:\r
+  AliStack * fStack ;\r
+  TList * fOutputContainer;        //final histogram container\r
+  TList * fPHOSEvents[1][10][11] ; //Containers for events with PHOS photons\r
+  TClonesArray * fPHOSEvent ;      //PHOS photons in current event\r
+  AliPHOSCalibData *fPHOSCalibData; // PHOS calibration object\r
+  TF1 *fNonLinCorr;          // Non-linearity correction\r
\r
+  //Reaction plain for v2\r
+  Float_t fRPfull ; //!Reaction plain calculated with full TPC \r
+  Float_t fRPA ;    //!Reaction plain calculated with A-side TPC: eta>0.15 \r
+  Float_t fRPC ;    //!Reaction plain calculated with C-side TPC: eta<-0.15\r
+  Float_t fRPFar ;  //!Reaction plain calculated with TPC: eta>0.6 \r
+  Float_t fRPAFar ; //!Reaction plain calculated with A-side TPC: eta>0.6 \r
+  Float_t fRPCFar ; //!Reaction plain calculated with C-side TPC: eta<-0.6\r
+\r
+  Float_t fCentrality ; //!Centrality of the currecnt event\r
+\r
+  Int_t fCenBin ;       //! Current centrality bin\r
+\r
+  TH2I *fPHOSBadMap[6] ;    //Container for PHOS bad channels map\r
+\r
+  AliPHOSGeometry  *fPHOSGeo;  //! PHOS geometry\r
+  Int_t fEventCounter;         // number of analyzed events\r
+\r
+  ClassDef(AliAnalysisTaskPi0Efficiency, 1); // PHOS analysis task\r
+};\r
+\r
+#endif\r
diff --git a/PWG4/UserTasks/PHOS_embedding/AliPHOSDigitDecalibrate.cxx b/PWG4/UserTasks/PHOS_embedding/AliPHOSDigitDecalibrate.cxx
new file mode 100644 (file)
index 0000000..6db616c
--- /dev/null
@@ -0,0 +1,141 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  PHOS tender, recalibrate PHOS clusters                                   //
+//  and do track matching                                                    //
+//  Author : Dmitri Peressounko (RRC KI)                                     //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+
+#include <AliLog.h>
+#include <AliESDEvent.h>
+#include <AliAnalysisManager.h>
+#include <AliTender.h>
+#include <TH2.h>
+#include <TROOT.h>
+
+#include "AliPHOSDigitDecalibrate.h"
+#include "AliPHOSGeometry.h"
+
+AliPHOSDigitDecalibrate::AliPHOSDigitDecalibrate() :
+  AliTenderSupply()
+  ,fPHOSGeo(0x0)
+  ,fPHOSCalibData(0x0)
+{
+       //
+       // default ctor
+       //
+}
+
+//_____________________________________________________
+AliPHOSDigitDecalibrate::AliPHOSDigitDecalibrate(const char *name, const AliTender *tender) :
+  AliTenderSupply(name,tender)
+  ,fPHOSGeo(0x0)
+  ,fPHOSCalibData(0x0)
+{
+       //
+       // named ctor
+       //
+}
+
+//_____________________________________________________
+AliPHOSDigitDecalibrate::~AliPHOSDigitDecalibrate()
+{
+  //Destructor
+  for(int m=0; m<5; m++){
+    if(hDec[m]){
+      delete hDec[m];
+      hDec[m]=0;
+    }
+  }
+}
+
+//_____________________________________________________
+void AliPHOSDigitDecalibrate::Init()
+{
+  //
+  // Initialise PHOS tender
+  //
+    
+
+  
+}
+
+//_____________________________________________________
+void AliPHOSDigitDecalibrate::ProcessEvent()
+{
+  //Choose PHOS clusters and recalibrate them
+  //that it recalculate energy, position and distance 
+  //to closest track extrapolation     
+
+  AliESDEvent *event=fTender->GetEvent();
+  if (!event) return;
+
+  // Init goemetry
+  if(!fPHOSGeo){
+    fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP") ;
+  }
+
+
+  AliESDCaloCells * cells = event->GetPHOSCells() ;
+
+  for (Short_t icell = 0; icell < cells->GetNumberOfCells(); icell++) {
+    Short_t id=0;
+    Double_t time=0., amp=0. ;
+    if (cells->GetCell(icell, id, amp, time) != kTRUE)
+      break;
+
+    Int_t relId[4] ;
+    fPHOSGeo->AbsToRelNumbering(id,relId);
+    Int_t   module = relId[0];
+    Int_t   column = relId[3];
+    Int_t   row    = relId[2];
+   
+//    amp=amp*fPHOSCalibData->GetADCchannelEmc(module,column,row);
+    amp=amp*hDec[module-1]->GetBinContent(row,column);
+
+    cells->SetCell(icell, id, amp, time);     
+
+  }
+
+
+
+
+}
+//_____________________________________________________
+void  AliPHOSDigitDecalibrate::SetDecalibration(Int_t mod, TH2F * dec){
+
+  if(mod<0 || mod>4){
+    printf("Don't know module %d \n",mod) ;
+    return ;
+  }
+  if(dec==0)
+    return ;
+  if(hDec[mod])
+    delete hDec[mod] ;
+  gROOT->cd() ;
+  hDec[mod] = new TH2F(*dec) ;
+  char key[55] ;
+  sprintf(key,"DecalibrationModule%d",mod) ;
+  hDec[mod]->SetName(key) ;
+
+
+
+}
+
diff --git a/PWG4/UserTasks/PHOS_embedding/AliPHOSDigitDecalibrate.h b/PWG4/UserTasks/PHOS_embedding/AliPHOSDigitDecalibrate.h
new file mode 100644 (file)
index 0000000..006aaca
--- /dev/null
@@ -0,0 +1,47 @@
+#ifndef ALIPHOSDIGITDECALIBRATE_H
+#define ALIPHOSDIGITDECALIBRATE_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+////////////////////////////////////////////////////////////////////////
+//                                                                    //
+//  PHOS tender, apply corrections to PHOS clusters                   //
+//  and do track matching                                             //
+//  Author : Dmitri Peressounko (RRC KI)                              //
+//                                                                    //
+////////////////////////////////////////////////////////////////////////
+
+#include <AliTenderSupply.h>
+
+class TVector3;
+class AliPHOSGeometry;
+class AliPHOSCalibData ; 
+class AliPHOSDigitDecalibrate: public AliTenderSupply {
+  
+public:
+  AliPHOSDigitDecalibrate();
+  AliPHOSDigitDecalibrate(const char *name, const AliTender *tender=NULL);
+  virtual ~AliPHOSDigitDecalibrate();
+
+  virtual void   Init();
+  virtual void   ProcessEvent();
+  
+  void  SetDecalibration(Int_t mod, TH2F * dec) ;
+
+protected:
+  AliPHOSDigitDecalibrate(const AliPHOSDigitDecalibrate&c);
+  AliPHOSDigitDecalibrate& operator= (const AliPHOSDigitDecalibrate&c);
+private:
+
+  TH2F * hDec[5] ;                           //! Decalibration coeff.
+  AliPHOSGeometry   *fPHOSGeo;               //! PHOS geometry
+  AliPHOSCalibData *fPHOSCalibData;          //! PHOS calibration object
+
+  ClassDef(AliPHOSDigitDecalibrate, 1); // PHOS tender task
+};
+
+
+#endif
+
diff --git a/PWG4/UserTasks/PHOS_embedding/AliPHOSEmbedding.cxx b/PWG4/UserTasks/PHOS_embedding/AliPHOSEmbedding.cxx
new file mode 100644 (file)
index 0000000..443fddd
--- /dev/null
@@ -0,0 +1,1052 @@
+// Analysis task to perform embedding of the simulated AOD
+// into real data (ESD)
+// Output conntainer contain AOD with
+//   0. Data before embedding
+//   1. Data+with embedded signal
+//   2. MC information from signal
+// Data for embedding are obtained using standard Manager
+// Signal is read from TChain with aodTree
+//
+// Authors: Dmitri Peressounko
+// Date   : 28.05.2011
+
+#include "TChain.h"
+#include "TFile.h"
+#include "TROOT.h"
+#include "TClonesArray.h"
+#include "TH2.h"
+#include "TGeoManager.h"
+#include "TGeoGlobalMagField.h"
+
+#include "AliGeomManager.h"
+#include "AliGRPObject.h"
+#include "AliCDBPath.h"
+#include "AliCDBEntry.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliVEventHandler.h"
+#include "AliInputEventHandler.h"
+#include "AliMultiplicity.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
+#include "AliCDBManager.h"
+#include "AliAODHandler.h"
+#include "AliMagF.h"
+
+#include "AliPHOSReconstructor.h"
+#include "AliPHOSClusterizerv1.h"
+#include "AliPHOSDigit.h"
+#include "AliPHOSGeometry.h"
+
+#include "AliPHOSEmbedding.h"
+
+ClassImp(AliPHOSEmbedding)
+
+
+//________________________________________________________________________
+AliPHOSEmbedding::AliPHOSEmbedding(const char *name) 
+  : AliAnalysisTaskSE(name),
+    fAODChain(0x0),
+    fDigitsTree(0x0),
+    fClustersTree(0x0),
+    fTreeOut(0x0),
+    fDigitsArr(0x0),
+    fEmbeddedClusters(0x0),
+    fEmbeddedCells(0x0),
+    fCellsPHOS(0x0),
+    fClusterizer(0x0),
+    fPHOSReconstructor(0x0),
+    fNSignal(0),
+    fNCaloClustersOld(0),
+    fInitialized(0)   
+{
+  // Constructor
+  for(int i=0;i<5;i++)fOldPHOSCalibration[i]=0x0 ;
+
+}
+
+//________________________________________________________________________
+void AliPHOSEmbedding::UserCreateOutputObjects()
+{
+ //prepare output containers
+  
+  //Create branch for output MC particles  
+  TClonesArray *mcparticles = new TClonesArray("AliAODMCParticle",0);
+  mcparticles->SetName(AliAODMCParticle::StdBranchName());
+  AddAODBranch("TClonesArray", &mcparticles);
+
+  //Branch with clusters after embedding
+  fEmbeddedClusters = new TClonesArray("AliAODCaloCluster",0);
+  fEmbeddedClusters->SetName("EmbeddedCaloClusters");
+  AddAODBranch("TClonesArray", &fEmbeddedClusters);
+  
+
+  //Branch with cells after embedding
+  fEmbeddedCells = new AliAODCaloCells();
+  fEmbeddedCells->SetName("EmbeddedPHOScells");
+  AddAODBranch("AliAODCaloCells", &fEmbeddedCells);
+
+  
+  fDigitsArr = new TClonesArray("AliPHOSDigit",3*56*64) ;
+
+  AliAODHandler* handler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+  if (handler) {
+      fTreeOut = handler->GetTree();
+  }
+  else {
+       AliWarning("No output AOD Event Handler connected.") ;
+  }
+
+}
+
+//________________________________________________________________________
+void AliPHOSEmbedding::UserExec(Option_t *) {
+  // Main loop, called for each event
+  // Perform embedding
+  
+  AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
+  if (!event) {
+    Printf("ERROR: Could not retrieve event");
+    PostData(0, fTreeOut);
+    return;
+  }
+
+  AliCentrality *centrality = event->GetCentrality(); 
+  Float_t acentrality=centrality->GetCentralityPercentileUnchecked("V0M");
+  if( acentrality <= 0. || acentrality>80.){
+    return;
+  }
+
+  
+  //Read next AOD event
+  //If necesary method checks if AOD event is good to embed
+  //e.g. there are PHOS clusters etc.
+  AliAODEvent * signal = GetNextSignalEvent() ;
+  if (!signal) {
+    Printf("ERROR: Could not retrieve signal");
+    PostData(0, fTreeOut);
+    return;
+  }
+  
+  Init() ;
+
+  //Remember PHOS digits before any modification
+  if(fCellsPHOS)
+    delete fCellsPHOS ;
+  fCellsPHOS = new AliESDCaloCells(*(event->GetPHOSCells())) ;
+
+  //First re-reconstruct existing digits to ensure identical reconstruction before and 
+  //after embeding
+  AliAODEvent * nullSignal=0x0 ;
+  MakeEmbedding(event,nullSignal) ;
+
+  //Convert ESD + MC information from signal to AOD
+  ConvertESDtoAOD(event);
+  
+  //Now perform real embedding  
+  //Embed signal clusters into event
+  MakeEmbedding(event,signal) ;
+
+  //Convert ESD to AOD
+  ConvertEmbeddedClusters(event) ;
+  //Fill MC information
+  ConvertMCParticles(signal) ;
+
+  delete signal ;
+
+  PostData(0, fTreeOut);
+}
+//______________________________________________________________________________
+void AliPHOSEmbedding::Init(){
+  if(fInitialized)
+    return ;
+  
+  AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
+  Int_t runNum = event->GetRunNumber();
+  AliCDBManager::Instance()->SetRun(runNum);
+
+  fPHOSReconstructor = new AliPHOSReconstructor() ;
+  AliCDBPath path("PHOS","Calib","RecoParam");
+  AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
+  TObjArray* recoParamArray = (TObjArray*)entry->GetObject();
+  AliPHOSRecoParam* recoParam = (AliPHOSRecoParam*)recoParamArray->At(2); 
+  fPHOSReconstructor->SetRecoParam(recoParam) ;
+  
+  InitMF() ;  
+  InitGeometry() ;
+
+  fInitialized=kTRUE ;
+}
+
+//________________________________________________________________________
+AliAODEvent* AliPHOSEmbedding::GetNextSignalEvent(){
+  //Read AOD event from the chain
+  
+  if(fAODChain==0){
+    AliError("No chain to read signal events") ;
+    return 0 ;
+  }
+  
+  AliAODEvent* aod = new AliAODEvent;
+  aod->ReadFromTree(fAODChain);
+  
+  if(fNSignal>=fAODChain->GetEntries()){
+    delete aod ;
+    return 0 ;
+  }
+  fAODChain->GetEvent(fNSignal) ;
+  fNSignal++ ;
+
+  //check if event is read and there is at least 2 PHOS clusters
+  //Otherwise fill histogram and read next event
+  //TODO
+
+  return aod ;
+  
+}
+
+//________________________________________________________________________
+void AliPHOSEmbedding::MakeEmbedding(AliESDEvent *event,AliAODEvent * signal){
+  //Perform embedding of the signal to the event
+
+  gROOT->cd(); //make sure that the digits and RecPoints Trees are memory resident
+  if(fDigitsTree){
+    delete fDigitsTree ;
+  }
+  fDigitsTree = new TTree("digitstree","digitstree");
+  fDigitsTree->Branch("PHOS","TClonesArray", &fDigitsArr, 32000);
+  
+  if(fClustersTree){
+    delete fClustersTree ;
+  }
+  fClustersTree = new TTree("clustertree","clustertree");
+
+  //Remember number of Clusters before we added new ones...
+  fNCaloClustersOld=event->GetNumberOfCaloClusters() ;
+
+
+/*
+  printf("---------------Before---------\n") ;
+  Int_t n=event->GetNumberOfCaloClusters() ;
+  for (Int_t iClust=0; iClust<event->GetNumberOfCaloClusters(); ++iClust) {
+    AliESDCaloCluster * cluster = event->GetCaloCluster(iClust);
+     if(!cluster->IsPHOS())
+      continue;
+   Float_t pos[3] = { 0.};
+    cluster->GetPosition(pos);
+    printf("Cluster(%d): E=%f, x=%f, z=%f, CPV=%f \n",iClust,cluster->E(),pos[0],pos[2],cluster->GetEmcCpvDistance()) ;
+    UShort_t * index    = cluster->GetCellsAbsId() ;
+    for(Int_t ic=0; ic < cluster->GetNCells(); ic++ )
+    printf("Dig(%d)=%d, ",ic,index[ic]) ;
+  printf("\n") ;
+       
+       
+  }
+  printf("---------------END Before---------\n") ;
+*/
+
+  
+  //create digits
+  MakeDigits(signal);    
+  
+  //clusterize and make tracking
+  fPHOSReconstructor->Reconstruct(fDigitsTree,fClustersTree) ;
+
+  //Note that the cirrent ESDEvent will be modified!
+  fPHOSReconstructor->FillESD(fDigitsTree, fClustersTree, event) ;
+
+
+/*
+   printf("---------------After---------\n") ;
+  for (Int_t iClust=n; iClust<event->GetNumberOfCaloClusters(); ++iClust) {
+    AliESDCaloCluster * cluster = event->GetCaloCluster(iClust);
+    if(!cluster->IsPHOS())
+      continue;
+    Float_t pos[3] = { 0.};
+    cluster->GetPosition(pos);
+    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()) ;
+    UShort_t * index    = cluster->GetCellsAbsId() ;
+    for(Int_t ic=0; ic < cluster->GetNCells(); ic++ )
+      printf("Dig(%d)=%d, ",ic,index[ic]) ;
+  printf("\n") ;
+  }
+  printf("---------------END After---------\n") ;
+*/
+
+
+}
+//______________________________________________________________________________
+void AliPHOSEmbedding::ConvertESDtoAOD(AliESDEvent* esd) 
+{
+  // ESD Filter analysis task executed for each event
+  
+  if(!esd)return;
+    
+  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
+  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
+      
+  ConvertHeader(*esd);
+  
+  // Fetch Stack for debuggging if available 
+  
+  // Update the header
+  Int_t nTracks    = esd->GetNumberOfTracks();
+  Int_t nV0s      = esd->GetNumberOfV0s();
+  Int_t nCascades = esd->GetNumberOfCascades();
+  Int_t nKinks    = esd->GetNumberOfKinks();
+  Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;
+  Int_t nPileSPDVertices=1+esd->GetNumberOfPileupVerticesSPD(); // also SPD main vertex
+  Int_t nPileTrkVertices=esd->GetNumberOfPileupVerticesTracks();
+  nVertices+=nPileSPDVertices;
+  nVertices+=nPileTrkVertices;
+  Int_t nJets     = 0;
+  Int_t nCaloClus = esd->GetNumberOfCaloClusters();
+  Int_t nFmdClus  = 0;
+  Int_t nPmdClus  = esd->GetNumberOfPmdTracks();
+    
+       
+  AODEvent()->ResetStd(nTracks, nVertices, nV0s, nCascades, nJets, nCaloClus, nFmdClus, nPmdClus);
+      
+  ConvertPrimaryVertices(*esd);
+  
+  ConvertCaloClusters(*esd);
+  
+  ConvertEMCALCells(*esd);
+  
+  ConvertPHOSCells(*esd);
+  
+}
+//______________________________________________________________________________
+void AliPHOSEmbedding::ConvertPHOSCells(const AliESDEvent& esd)
+{
+  // fill PHOS cell info
+  if (esd.GetPHOSCells()) { // protection against missing ESD information
+    AliESDCaloCells &esdPHcells = *(esd.GetPHOSCells());
+    Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
+    
+    AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
+    aodPHcells.CreateContainer(nPHcell);
+    aodPHcells.SetType(AliAODCaloCells::kPHOSCell);
+    for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      
+      aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell));
+    }
+    aodPHcells.Sort();
+  }
+}
+
+//______________________________________________________________________________
+void AliPHOSEmbedding::ConvertEMCALCells(const AliESDEvent& esd)
+{
+  // fill EMCAL cell info
+  if (esd.GetEMCALCells()) { // protection against missing ESD information
+    AliESDCaloCells &esdEMcells = *(esd.GetEMCALCells());
+    Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
+    
+    AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
+    aodEMcells.CreateContainer(nEMcell);
+    aodEMcells.SetType(AliAODCaloCells::kEMCALCell);
+    for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
+      aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell));
+    }
+    aodEMcells.Sort();
+  }
+}
+//______________________________________________________________________________
+void AliPHOSEmbedding::ConvertCaloClusters(const AliESDEvent& esd)
+{
+  
+  // Access to the AOD container of clusters
+  TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
+  Int_t jClusters(0);
+  
+  for (Int_t iClust=fNCaloClustersOld; iClust<esd.GetNumberOfCaloClusters(); ++iClust) {
+    
+    AliESDCaloCluster * cluster = esd.GetCaloCluster(iClust);
+    
+    Int_t  id        = cluster->GetID();
+    Int_t  nLabel    = cluster->GetNLabels();
+    Int_t *labels    = cluster->GetLabels();
+    
+    Float_t energy = cluster->E();
+    Float_t posF[3] = { 0.};
+    cluster->GetPosition(posF);
+    
+    AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,
+                                                                                      nLabel,
+                                                                                      labels,
+                                                                                      energy,
+                                                                                      posF,
+                                                                                      NULL,
+                                                                                      cluster->GetType(),0);
+    
+    Double_t dx=cluster->GetTrackDx() ;
+    Double_t dz=cluster->GetTrackDz() ;
+    Float_t cpv=99. ; //No track matched by default
+    TArrayI * itracks = cluster->GetTracksMatched() ;  
+    if(itracks->GetSize()>0){
+      Int_t iTr = itracks->At(0);
+      if(iTr>=0 && iTr<esd.GetNumberOfTracks()){
+        AliESDtrack *track = esd.GetTrack(iTr) ;
+        Double_t pt = track->Pt() ;
+        Short_t charge = track->Charge() ;
+        cpv=TestCPV(dx, dz, pt,charge) ;
+      }
+    }
+    
+    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
+                                cluster->GetDispersion(),
+                                cluster->GetM20(), cluster->GetM02(),
+                                cpv,  
+                                cluster->GetNExMax(),cluster->GetTOF()) ;
+    
+    caloCluster->SetPIDFromESD(cluster->GetPID());
+    caloCluster->SetNCells(cluster->GetNCells());
+    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
+    caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
+
+    /*    
+    TArrayI* matchedT =        cluster->GetTracksMatched();
+    if (fNumberOfTracks>0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {       
+      for (Int_t im = 0; im < matchedT->GetSize(); im++) {
+        Int_t iESDtrack = matchedT->At(im);;
+        if (fAODTrackRefs->At(iESDtrack) != 0) {
+          caloCluster->AddTrackMatched((AliAODTrack*)fAODTrackRefs->At(iESDtrack));
+        }
+      }
+    }
+    */
+    
+  } 
+  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters     
+}
+//______________________________________________________________________________
+void AliPHOSEmbedding::ConvertEmbeddedClusters(const AliESDEvent* esd)
+{
+   //Copy PHOS clusters and cells after embedding
+  
+  // Access to the AOD container of clusters
+  Int_t jClusters(0);
+  fEmbeddedClusters->Clear();
+  fEmbeddedClusters->Expand(esd->GetNumberOfCaloClusters()-fNCaloClustersOld) ;
+  for (Int_t iClust=fNCaloClustersOld; iClust<esd->GetNumberOfCaloClusters(); ++iClust) {
+    
+    AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
+    
+    Int_t  id        = cluster->GetID();
+    Int_t  nLabel    = cluster->GetNLabels();
+    Int_t *labels    = cluster->GetLabels();
+    
+    Float_t energy = cluster->E();
+    Float_t posF[3] = { 0.};
+    cluster->GetPosition(posF);
+    
+    AliAODCaloCluster *caloCluster = new((*fEmbeddedClusters)[jClusters++]) AliAODCaloCluster(id,
+                                                                                      nLabel,
+                                                                                      labels,
+                                                                                      energy,
+                                                                                      posF,
+                                                                                      NULL,
+                                                                                      cluster->GetType(),0);
+
+    
+    Double_t dx=cluster->GetTrackDx() ;
+    Double_t dz=cluster->GetTrackDz() ;
+    Float_t cpv=99. ; //No track matched by default
+    TArrayI * itracks = cluster->GetTracksMatched() ;  
+    if(itracks->GetSize()>0){
+      Int_t iTr = itracks->At(0);
+      if(iTr>=0 && iTr<esd->GetNumberOfTracks()){
+        AliESDtrack *track = esd->GetTrack(iTr) ;
+        Double_t pt = track->Pt() ;
+        Short_t charge = track->Charge() ;
+        cpv=TestCPV(dx, dz, pt,charge) ;
+      }
+    }
+
+    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
+                                cluster->GetDispersion(),
+                                cluster->GetM20(), cluster->GetM02(),
+                                cpv,  
+                                cluster->GetNExMax(),cluster->GetTOF()) ;
+    
+    caloCluster->SetPIDFromESD(cluster->GetPID());
+    caloCluster->SetNCells(cluster->GetNCells());    
+    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
+    caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());            
+  } 
+  
+  for(Int_t i=0;i<fEmbeddedClusters->GetEntriesFast();i++){
+    AliAODCaloCluster *caloCluster =(AliAODCaloCluster *)fEmbeddedClusters->At(i) ;
+    caloCluster->GetID() ;
+    
+  }
+  
+  //Now cells
+  if (esd->GetPHOSCells()) { // protection against missing ESD information
+    AliESDCaloCells &esdPHcells = *(esd->GetPHOSCells());
+    Int_t nPHcell = esdPHcells.GetNumberOfCells() ;    
+    fEmbeddedCells->CreateContainer(nPHcell);
+    fEmbeddedCells->SetType(AliAODCaloCells::kPHOSCell);
+    for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      
+      fEmbeddedCells->SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell));
+    }
+    fEmbeddedCells->Sort();
+  }
+
+  
+}
+//______________________________________________________________________________
+void AliPHOSEmbedding::ConvertMCParticles(const AliAODEvent* aod)
+{
+  //Copy MC branches to new AOD
+  
+  TClonesArray *mcArray = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
+  TClonesArray * aodMcParticles = (TClonesArray*)AODEvent()->FindListObject(AliAODMCParticle::StdBranchName());
+  for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
+    
+    AliAODMCParticle* aodpart =  (AliAODMCParticle*) mcArray->At(i);
+//    AliAODMCParticle * mcp =new ((*aodMcParticles)[i]) AliAODMCParticle(*aodpart);
+    new ((*aodMcParticles)[i]) AliAODMCParticle(*aodpart);
+    
+  }
+
+  /*  
+  // MC header
+  AliAODMCHeader * mcHeader = (AliAODMCHeader *)aod.FindListObject(AliAODMCHeader::StdBranchName());
+  AliAODMCHeader * aodMcHeader = new AliAODMCHeader(*mcHeader);
+  aodMcHeader->SetName(AliAODMCHeader::StdBranchName());
+  AddAODBranch("AliAODMCHeader",&aodMcHeader);    
+  */
+
+}
+//______________________________________________________________________________
+void AliPHOSEmbedding::ConvertHeader(AliESDEvent & esd){
+  
+  AliAODHeader* header = AODEvent()->GetHeader();
+  
+  header->SetRunNumber(esd.GetRunNumber());
+  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
+
+  TTree* tree = fInputHandler->GetTree();
+  if (tree) {
+    TFile* file = tree->GetCurrentFile();
+    if (file) header->SetESDFileName(file->GetName());
+  }
+  
+  header->SetBunchCrossNumber(esd.GetBunchCrossNumber());
+  header->SetOrbitNumber(esd.GetOrbitNumber());
+  header->SetPeriodNumber(esd.GetPeriodNumber());
+  header->SetEventType(esd.GetEventType());
+  
+  header->SetEventNumberESDFile(esd.GetHeader()->GetEventNumberInFile());
+  if(const_cast<AliESDEvent&>(esd).GetCentrality()){
+    header->SetCentrality(const_cast<AliESDEvent&>(esd).GetCentrality());
+  }
+  else{
+    header->SetCentrality(0);
+  }
+  if(const_cast<AliESDEvent&>(esd).GetEventplane()){
+    header->SetEventplane(const_cast<AliESDEvent&>(esd).GetEventplane());
+  }
+  else{
+    header->SetEventplane(0);
+  }
+
+  // Trigger
+  header->SetFiredTriggerClasses(esd.GetFiredTriggerClasses());
+  header->SetTriggerMask(esd.GetTriggerMask()); 
+  header->SetTriggerCluster(esd.GetTriggerCluster());
+  header->SetL0TriggerInputs(esd.GetHeader()->GetL0TriggerInputs());    
+  header->SetL1TriggerInputs(esd.GetHeader()->GetL1TriggerInputs());    
+  header->SetL2TriggerInputs(esd.GetHeader()->GetL2TriggerInputs());    
+  
+  header->SetMagneticField(esd.GetMagneticField());
+  header->SetMuonMagFieldScale(esd.GetCurrentDip()/6000.);
+  header->SetZDCN1Energy(esd.GetZDCN1Energy());
+  header->SetZDCP1Energy(esd.GetZDCP1Energy());
+  header->SetZDCN2Energy(esd.GetZDCN2Energy());
+  header->SetZDCP2Energy(esd.GetZDCP2Energy());
+  header->SetZDCEMEnergy(esd.GetZDCEMEnergy(0),esd.GetZDCEMEnergy(1));
+  
+  // ITS Cluster Multiplicty
+  const AliMultiplicity *mult = esd.GetMultiplicity();
+  for (Int_t ilay = 0; ilay < 6; ilay++) header->SetITSClusters(ilay, mult->GetNumberOfITSClusters(ilay));
+  
+  // TPC only Reference Multiplicty
+  //  Int_t refMult  = fTPCaloneTrackCuts ? (Short_t)fTPCaloneTrackCuts->GetReferenceMultiplicity(&esd, kTRUE) : -1;
+  header->SetTPConlyRefMultiplicity(-1);
+  
+  //
+  Float_t diamxy[2]={esd.GetDiamondX(),esd.GetDiamondY()};
+  Float_t diamcov[3]; 
+  esd.GetDiamondCovXY(diamcov);
+  header->SetDiamond(diamxy,diamcov);
+  header->SetDiamondZ(esd.GetDiamondZ(),esd.GetSigma2DiamondZ());
+
+  
+  //PHOS matrixes
+  char path[255] ;
+  TGeoHMatrix * m ;
+  for(Int_t mod=0; mod<5; mod++){
+    snprintf(path,255,"/ALIC_1/PHOS_%d",mod+1) ; //In Geometry modules numbered 1,2,.,5
+    if (gGeoManager->cd(path)){
+      m = gGeoManager->GetCurrentMatrix() ;
+      header->SetPHOSMatrix(new TGeoHMatrix(*m),mod) ;
+    }
+    else{
+      header->SetPHOSMatrix(NULL,mod) ;
+    }
+  } 
+  
+  //EventPlane
+  Double_t epQ=TPCrp(&esd);
+  //Standard Eventplane setter is too complicated...
+  //use durty hack
+  header->SetZDCN1Energy(epQ) ;
+  
+  AliCentrality *centrality = esd.GetCentrality(); 
+  Double_t acentrality=centrality->GetCentralityPercentileUnchecked("V0M");
+  header->SetZDCN2Energy(acentrality) ;
+
+  
+}//______________________________________________________________________________
+void AliPHOSEmbedding::ConvertPrimaryVertices(AliESDEvent const&esd){
+  // Access to the AOD container of vertices
+
+  // Access to the AOD container of vertices
+  Int_t fNumberOfVertices = 0;
+    
+  Double_t pos[3] = { 0. };
+  Double_t covVtx[6] = { 0. };
+
+  // Add primary vertex. The primary tracks will be defined
+  // after the loops on the composite objects (V0, cascades, kinks)
+  const AliESDVertex *vtx = esd.GetPrimaryVertex();
+  
+  vtx->GetXYZ(pos); // position
+  vtx->GetCovMatrix(covVtx); //covariance matrix
+  
+  TClonesArray * vertexes=AODEvent()->GetVertices();
+
+  AliAODVertex* primaryVertex = new((*vertexes)[fNumberOfVertices++])
+  AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
+  primaryVertex->SetName(vtx->GetName());
+  primaryVertex->SetTitle(vtx->GetTitle());
+  
+  TString vtitle = vtx->GetTitle();
+  if (!vtitle.Contains("VertexerTracks")) 
+    primaryVertex->SetNContributors(vtx->GetNContributors());
+  
+  if (fDebug > 0) primaryVertex->Print();  
+  
+  // Add SPD "main" vertex 
+  const AliESDVertex *vtxS = esd.GetPrimaryVertexSPD();
+  vtxS->GetXYZ(pos); // position
+  vtxS->GetCovMatrix(covVtx); //covariance matrix
+  AliAODVertex * mVSPD = new((*vertexes)[fNumberOfVertices++])
+  AliAODVertex(pos, covVtx, vtxS->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainSPD);
+  mVSPD->SetName(vtxS->GetName());
+  mVSPD->SetTitle(vtxS->GetTitle());
+  mVSPD->SetNContributors(vtxS->GetNContributors()); 
+  
+  // Add SPD pileup vertices
+  for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesSPD(); ++iV)
+  {
+    const AliESDVertex *vtxP = esd.GetPileupVertexSPD(iV);
+    vtxP->GetXYZ(pos); // position
+    vtxP->GetCovMatrix(covVtx); //covariance matrix
+    AliAODVertex * pVSPD =  new((*vertexes)[fNumberOfVertices++])
+    AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupSPD);
+    pVSPD->SetName(vtxP->GetName());
+    pVSPD->SetTitle(vtxP->GetTitle());
+    pVSPD->SetNContributors(vtxP->GetNContributors()); 
+  }
+  
+  // Add TRK pileup vertices
+  for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesTracks(); ++iV)
+  {
+    const AliESDVertex *vtxP = esd.GetPileupVertexTracks(iV);
+    vtxP->GetXYZ(pos); // position
+    vtxP->GetCovMatrix(covVtx); //covariance matrix
+    AliAODVertex * pVTRK = new((*vertexes)[fNumberOfVertices++])
+    AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupTracks);
+    pVTRK->SetName(vtxP->GetName());
+    pVTRK->SetTitle(vtxP->GetTitle());
+    pVTRK->SetNContributors(vtxP->GetNContributors());
+  }
+
+  
+}
+//__________________________________________________________________________________
+void AliPHOSEmbedding::MakeDigits(AliAODEvent * signal){
+  
+    //-------------------------------------------------------------------------------------
+  //Transform CaloCells into Digits
+  //-------------------------------------------------------------------------------------
+  
+  fDigitsArr->Clear() ;
+  fDigitsTree->Branch("PHOS","TClonesArray", &fDigitsArr, 32000);
+
+  //First copy data digits
+  Int_t ndigit=0 ;
+  for (Short_t icell = 0; icell < fCellsPHOS->GetNumberOfCells(); icell++) {
+    Short_t id=0;
+    Double_t time=0., amp=0. ;
+    if (fCellsPHOS->GetCell(icell, id, amp, time) != kTRUE)
+      break;
+        
+    new((*fDigitsArr)[ndigit]) AliPHOSDigit(-1,id,float(amp),float(time),ndigit);
+    ndigit++;
+  }
+
+
+  //Add Digits from Signal
+  TClonesArray sdigits("AliPHOSDigit",0) ;
+  Int_t isdigit=0 ;
+  if(signal){
+    AliAODCaloCells* cellsS = signal->GetPHOSCells();
+    Int_t cellLabels[cellsS->GetNumberOfCells()] ;
+    Int_t cellSecondLabels[cellsS->GetNumberOfCells()] ;
+    for(Int_t i=0;i<cellsS->GetNumberOfCells();i++){
+      cellLabels[i]=-1 ;
+      cellSecondLabels[i]=-1;
+   }
+//  printf("===========Signal clusters==============\n") ;
+    //------------------------------------------------------------------------------------
+    //Ancestry information
+    //------------------------------------------------------------------------------------
+    sdigits.Expand(cellsS->GetNumberOfCells());
+    for(Int_t i=0; i<signal->GetNumberOfCaloClusters(); i++) {    
+      //cluster from embedded signal
+      AliVCluster *clus = signal->GetCaloCluster(i);  
+    
+      if(!clus->IsPHOS())
+        continue;
+/*    
+    printf("Signal clu(%d): E=%f \n",i,clus->E()) ;
+    UShort_t * ind    = clus->GetCellsAbsId() ;
+    for(Int_t ic=0; ic < clus->GetNCells(); ic++ )
+      printf("Dig(%d)=%d, ",ic,ind[ic]) ;
+  printf("\n") ;
+*/   
+    
+      Int_t label = clus->GetLabel();
+      Int_t label2 = -1 ;
+      if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
+    
+      UShort_t * index    = clus->GetCellsAbsId() ;
+    
+      for(Int_t ic=0; ic < clus->GetNCells(); ic++ ){
+        for (Int_t icell = 0; icell < cellsS->GetNumberOfCells(); icell++){
+          Short_t cellNumber;
+          Double_t cellAmplitude=0., cellTime=0. ;
+          cellsS->GetCell(icell, cellNumber, cellAmplitude, cellTime) ;
+          if(cellNumber==index[ic]){
+             cellLabels[icell]=label;
+              cellSecondLabels[icell]=label2;
+             break ;
+         }
+        }
+      }
+    }
+// printf("================End Signal==================\n") ;
+
+    for (Int_t icell = 0; icell < cellsS->GetNumberOfCells(); icell++) {
+      Short_t cellNumber;
+      Double_t cellAmplitude=0., cellTime=0. ;
+      if (cellsS->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
+        break;
+      //Add only digits related to the cluster, no noisy digits...
+      if(cellLabels[icell]==-1)
+        continue ;
+    
+      new(sdigits[isdigit]) AliPHOSDigit(cellLabels[icell],cellNumber,float(cellAmplitude),float(cellTime),isdigit);    
+      isdigit++;
+    }
+  }
+  
+  //If necessary, take into account difference between calibration used to produce ESD and
+  //final calibration
+  if(fOldPHOSCalibration[0]){ //there is a difference
+    for(Int_t i=0; i<isdigit;i++){
+      AliPHOSDigit * sdigit=static_cast<AliPHOSDigit*>(sdigits.At(i)) ;
+      Int_t relId[4] ;
+      AliPHOSGeometry::GetInstance()->AbsToRelNumbering(sdigit->GetId(),relId);
+      Float_t calibESD=fOldPHOSCalibration[relId[0]-1]->GetBinContent(relId[2],relId[3]) ;
+      Float_t calibNew=fPHOSReconstructor->Calibrate(1.,sdigit->GetId()) ;
+      if(calibNew>0.)
+        sdigit->SetEnergy(sdigit->GetEnergy()*calibESD/calibNew) ;
+    }
+  }
+  
+  //Merge digits  
+  Int_t icurrent = 0 ; //index of the last used digit in underlying event
+  fDigitsArr->Expand(ndigit+isdigit) ;
+  for(Int_t i=0; i<isdigit;i++){
+    AliPHOSDigit * sdigit=static_cast<AliPHOSDigit*>(sdigits.At(i)) ;
+    Bool_t added=kFALSE ;
+    for(Int_t id=icurrent;id<ndigit;id++){
+      AliPHOSDigit * digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(id)) ;
+      if(sdigit->GetId() ==  digit->GetId() ){
+        *digit += *sdigit ;  //add energies
+       icurrent=id+1 ;
+        added=kTRUE ;
+       break ; //no more digits with same ID in the list
+      }
+      if(sdigit->GetId() < digit->GetId() ){
+        icurrent=id ;
+        break ; //no more digits with same ID in the list
+      }
+    }
+    if(!added){
+      new((*fDigitsArr)[ndigit]) AliPHOSDigit(*sdigit) ;
+      ndigit++ ;
+    }
+  }
+  //Change Amp back from Energy to ADC counts
+  for(Int_t i=0; i<ndigit;i++){
+     AliPHOSDigit * digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
+     Float_t calib = 0 ; 
+     if(fOldPHOSCalibration[0]){ //Use old calibration
+        Int_t relId[4] ;
+       AliPHOSGeometry::GetInstance()->AbsToRelNumbering(digit->GetId(),relId);
+        calib=fOldPHOSCalibration[relId[0]-1]->GetBinContent(relId[2],relId[3]) ;
+     }
+     else{
+       calib=fPHOSReconstructor->Calibrate(1.,digit->GetId()) ;        
+     }
+     if(calib>0.)
+       digit->SetEnergy(digit->GetEnergy()/calib) ;
+  }  
+  
+   
+  fDigitsArr->Sort() ;
+  for (Int_t i = 0 ; i < ndigit ; i++) { 
+    AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(i) ) ; 
+    digit->SetIndexInList(i) ;     
+  }
+  fDigitsTree->Fill() ;
+
+}
+//__________________________________________________________________________________
+void AliPHOSEmbedding::SetOldCalibration(TH2F **calibH){ 
+  //Copy histograms with calibration coeff
+  gROOT->cd() ;
+  char name[55] ;
+  for(Int_t i=0;i<5;i++){
+    fOldPHOSCalibration[i] = new TH2F(*((TH2F*)calibH[i])) ;
+    sprintf(name,"OldPHOSCalibrationMod%d",i) ;
+    fOldPHOSCalibration[i]->SetName(name) ;
+  }
+}  
+//___________________________________________________________________________
+Double_t AliPHOSEmbedding::TPCrp(const AliESDEvent * event){
+  //calculate reaction plain in TPC
+  Double_t rp=999 ; //!Reaction plain calculated with full TPC 
+  
+  Double_t cosF=0.,sinF=0.,nF=0.;
+  for (Int_t i=0;i<event->GetNumberOfTracks();i++) {
+    AliESDtrack *track = event->GetTrack(i) ;
+    if(!SelectTrack(track))
+      continue ;
+    Double_t eta= track->Eta() ;
+    if(TMath::Abs(eta)> 0.9)
+      continue ;
+    Double_t phi=track->Phi();
+    Double_t cos2phi=TMath::Cos(2.*phi) ;
+    Double_t sin2phi=TMath::Sin(2.*phi) ;
+    cosF+=cos2phi ;
+    sinF+=sin2phi ;
+    nF+=1. ;
+  }
+  if(nF>0){
+    rp=0.5*TMath::ATan2(sinF,cosF) ;
+  }
+  return rp ;
+  
+}
+//___________________________________________________________________________
+Bool_t AliPHOSEmbedding::SelectTrack(AliESDtrack * t){
+  //estimate if this track can be used for the RP calculation
+  Float_t pt=t->Pt();
+  if(pt<0.15 || pt>5.)
+    return kFALSE ;
+//  if(!fESDtrackCuts->AcceptTrack(t))
+//    return kFALSE ;
+  if(t->GetTPCNcls()<70) 
+    return kFALSE;
+/*
+  Float_t dcaXY=t->DCA();
+  Float_t dcaZ=t->ZAtDCA() ;
+  t->GetImpactParametersTPC(&dcaXY,&dcaZ);
+  if (TMath::Abs(dcaZ)>3)
+    return kFALSE;
+  if (TMath::Abs(dcaXY)>3.) 
+    return kFALSE;
+*/  
+  return kTRUE ;
+
+
+}
+//____________________________________________________________________________
+void AliPHOSEmbedding::InitMF(){
+  
+  printf("............Init MF \n") ;
+  //------------------------------------
+  // Initialization of the Mag.Fiels from GRP entry
+  // Copied from AliReconstruction::InitGRP()
+  //------------------------------------
+  AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
+  AliGRPObject * aGRPData=0 ;
+  if (entry) {
+
+    TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
+
+    if (m) {
+       AliInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
+       m->Print();
+       aGRPData = new AliGRPObject();
+       aGRPData->ReadValuesFromMap(m);
+    }
+
+    else {
+       AliInfo("Found an AliGRPObject in GRP/GRP/Data, reading it");
+       aGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
+       entry->SetOwner(0);
+    }
+  }
+
+  if (!aGRPData) {
+     AliError("No GRP entry found in OCDB!");
+     return ;
+  }
+  //*** Dealing with the magnetic field map
+    
+  TString lhcState = aGRPData->GetLHCState();
+  if (lhcState==AliGRPObject::GetInvalidString()) {
+    AliError("GRP/GRP/Data entry:  missing value for the LHC state ! Using UNKNOWN");
+    lhcState = "UNKNOWN";
+  }
+
+  TString beamType = aGRPData->GetBeamType();
+  if (beamType==AliGRPObject::GetInvalidString()) {
+    AliError("GRP/GRP/Data entry:  missing value for the beam type ! Using UNKNOWN");
+    beamType = "UNKNOWN";
+  }
+
+  Float_t beamEnergy = aGRPData->GetBeamEnergy();
+  if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
+    AliError("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
+    beamEnergy = 0;
+  }
+
+  TString runType = aGRPData->GetRunType();
+  if (runType==AliGRPObject::GetInvalidString()) {
+    AliError("GRP/GRP/Data entry:  missing value for the run type ! Using UNKNOWN");
+    runType = "UNKNOWN";
+  }
+  
+  if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
+    if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
+      AliInfo("PHOSEmbedding: MF is locked - GRP information will be ignored !");
+      AliInfo("Running with the externally locked B field !");
+    }
+    else {
+      AliInfo("Destroying existing B field instance!");
+      delete TGeoGlobalMagField::Instance();
+    }   
+  }
+  if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
+    // Construct the field map out of the information retrieved from GRP.
+    Bool_t ok = kTRUE;
+    // L3
+    Float_t l3Current = aGRPData->GetL3Current((AliGRPObject::Stats)0);
+    if (l3Current == AliGRPObject::GetInvalidFloat()) {
+      AliError("GRP/GRP/Data entry:  missing value for the L3 current !");
+      ok = kFALSE;
+    }
+
+    Char_t l3Polarity = aGRPData->GetL3Polarity();
+    if (l3Polarity == AliGRPObject::GetInvalidChar()) {
+      AliError("GRP/GRP/Data entry:  missing value for the L3 polarity !");
+      ok = kFALSE;
+    }
+
+    // Dipole
+    Float_t diCurrent = aGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
+    if (diCurrent == AliGRPObject::GetInvalidFloat()) {
+      AliError("GRP/GRP/Data entry:  missing value for the dipole current !");
+      ok = kFALSE;
+    }
+
+    Char_t diPolarity = aGRPData->GetDipolePolarity();
+    if (diPolarity == AliGRPObject::GetInvalidChar()) {
+      AliError("GRP/GRP/Data entry:  missing value for the dipole polarity !");
+      ok = kFALSE;
+    }
+
+    // read special bits for the polarity convention and map type
+    Int_t  polConvention = aGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
+    Bool_t uniformB = aGRPData->IsUniformBMap();
+
+    if (ok) {
+      AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1),
+                                             TMath::Abs(diCurrent) * (diPolarity ? -1:1),
+                                             polConvention,uniformB,beamEnergy, beamType.Data());
+      if (fld) {
+        TGeoGlobalMagField::Instance()->SetField( fld );
+        TGeoGlobalMagField::Instance()->Lock();
+        AliInfo("Running with the B field constructed out of GRP !");
+      }
+      else AliFatal("Failed to create a B field map !");
+    }
+    else AliFatal("B field is neither set nor constructed from GRP ! Exitig...");
+  }
+  
+}
+//____________________________________________________________________________
+void AliPHOSEmbedding::InitGeometry(){
+
+  // Import ideal TGeo geometry and apply misalignment
+  if (!gGeoManager) {
+    AliGeomManager::LoadGeometry("geometry.root");
+    if (!gGeoManager) {
+      AliFatal("Can not load geometry");
+    }
+    if(!AliGeomManager::CheckSymNamesLUT("PHOS")) {
+      AliFatal("CheckSymNamesLUT");
+    }
+  }
+    
+
+  TString detStr = "PHOS";
+  TString loadAlObjsListOfDets = "PHOS";
+        
+  if(AliGeomManager::GetNalignable("GRP") != 0)
+      loadAlObjsListOfDets.Prepend("GRP "); //add alignment objects for non-sensitive modules
+  AliGeomManager::ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
+
+  AliCDBManager::Instance()->UnloadFromCache("*/Align/*");
+  AliCDBManager::Instance()->UnloadFromCache("GRP/Geometry/Data");
+}
+//____________________________________________________________________________
+Float_t AliPHOSEmbedding::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
+  //Parameterization of LHC10h period
+  //_true if neutral_
+  
+  Double_t z0=0.84 ; //Both charges
+  Double_t sz=TMath::Min(2.75,2.*3.*3./(pt*pt+3.*3.)+1.);
+  Double_t x0=TMath::Min(7.6, 27.4628*TMath::Exp(-(pt+13.3572)*(pt+13.3572)/2./6.29214/6.29214)+
+                                    243.020*0.0526626*0.0526626/(pt*pt+0.0526626*0.0526626)) ;
+  if(charge>0)
+    x0=-x0 ;
+  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);
+  Double_t rz=(dz-z0)/sz ;
+  Double_t rx=(dx-x0)/sx ;
+  return TMath::Sqrt(rx*rx+rz*rz)/2. ;
+}
diff --git a/PWG4/UserTasks/PHOS_embedding/AliPHOSEmbedding.h b/PWG4/UserTasks/PHOS_embedding/AliPHOSEmbedding.h
new file mode 100644 (file)
index 0000000..b67bb8f
--- /dev/null
@@ -0,0 +1,87 @@
+#ifndef AliPHOSEmbedding_h\r
+#define AliPHOSEmbedding_h\r
+\r
+// Class to perform embedding on the AOD level\r
+// Author: D.Peressounko\r
+\r
+class TChain ;\r
+class TClonesArray ;\r
+class TH2F ;\r
+\r
+class AliPHOSClusterizerv1 ;\r
+class AliPHOSReconstructor ;\r
+class AliAODEvent ;\r
+class AliESDEvent ;\r
+class AliESDtrack ;\r
+class AliESDCaloCells ;\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class AliPHOSEmbedding : public AliAnalysisTaskSE {\r
+public:\r
+  AliPHOSEmbedding(const char *name = "AliPHOSEmbedding");\r
+  virtual ~AliPHOSEmbedding() {}\r
+  \r
+  virtual void   UserCreateOutputObjects();\r
+  virtual void   UserExec(Option_t *option);\r
+  virtual void   Terminate(Option_t *){}\r
+  \r
+\r
+  void SetSignalChain(TChain * signal){fAODChain =signal;}\r
+\r
+  void SetOldCalibration(TH2F **calib) ; \r
+  //Calibration used in reconstruction of real data (ESDs)\r
+  //If not set, assume same calibration as set by default\r
+\r
+private:\r
+  AliPHOSEmbedding(const AliPHOSEmbedding&); // not implemented\r
+  AliPHOSEmbedding& operator=(const AliPHOSEmbedding&); // not implemented\r
+\r
+  void Init() ;\r
+  void InitMF() ; //Mag.Field initialization for track matching\r
+  void InitGeometry() ;\r
+  \r
+  AliAODEvent * GetNextSignalEvent(void) ;\r
+\r
+  void MakeEmbedding(AliESDEvent * data, AliAODEvent * signal) ;\r
+  void MakeDigits(AliAODEvent* signal) ;  \r
+\r
+  void ConvertESDtoAOD(AliESDEvent *esd) ;\r
+  void ConvertHeader(AliESDEvent &esd) ;\r
+  void ConvertPrimaryVertices(const AliESDEvent &esd) ;\r
+  void ConvertCaloClusters(const AliESDEvent &esd) ;\r
+  void ConvertEMCALCells(const AliESDEvent &esd) ;\r
+  void ConvertPHOSCells(const AliESDEvent &esd) ;\r
+  \r
+  void ConvertEmbeddedClusters(const AliESDEvent *esd) ;\r
+  void ConvertEmbeddedCells(const AliESDEvent *esd) ;\r
+  void ConvertMCParticles(const AliAODEvent *aod) ;\r
+\r
+  Double_t TPCrp(const AliESDEvent * event) ;\r
+  Bool_t SelectTrack(AliESDtrack * t) ;  \r
+  Float_t TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge) ;\r
+  \r
+\r
+  TChain * fAODChain ; //Signal\r
+\r
+  TTree * fDigitsTree ;  //! Digits\r
+  TTree * fClustersTree; //! Clusters\r
+  TTree * fTreeOut; //Output AOD\r
+  TClonesArray * fDigitsArr ; //!\r
+\r
+  TClonesArray * fEmbeddedClusters ; //!\r
+  AliAODCaloCells * fEmbeddedCells ; //!\r
+  AliESDCaloCells * fCellsPHOS ; //! Old PHOS cells\r
+\r
+  AliPHOSClusterizerv1 * fClusterizer ; //!\r
+  AliPHOSReconstructor * fPHOSReconstructor ; //!\r
+  \r
+  TH2F * fOldPHOSCalibration[5] ; //! Calibration coeff. used in ESD production\r
+\r
+  Int_t fNSignal ; // Number of signal evetns processed  \r
+  Int_t fNCaloClustersOld ; //Number of CaloClusters already in ESD\r
+  Bool_t fInitialized ; //!\r
+  ClassDef(AliPHOSEmbedding, 1); // PHOS analysis task\r
+};\r
+\r
+#endif\r
diff --git a/PWG4/UserTasks/PHOS_embedding/Analyze.C b/PWG4/UserTasks/PHOS_embedding/Analyze.C
new file mode 100644 (file)
index 0000000..b681a92
--- /dev/null
@@ -0,0 +1,68 @@
+void Analyze()
+{
+    
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libPhysics.so");
+  
+  //load analysis framework
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice"); //AliAnalysisTaskSE
+  
+  gSystem->AddIncludePath("-I$ALICE_ROOT/include -I$ALICE_ROOT/PHOS");
+
+  // A task can be compiled dynamically with AClic
+  gROOT->LoadMacro("AliCaloPhoton.cxx+g");
+  gROOT->LoadMacro("AliAnalysisTaskPi0Efficiency.cxx+g");
+  
+  // Create the chain
+  TChain* chain = new TChain("aodTree");
+  chain->AddFile("AliAODout.root") ;
+
+  // Make the analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("Pi0EmbeddingManager");
+  
+  // ESD input handler
+  AliAODInputHandler* aodH = new AliAODInputHandler();
+  mgr->SetInputEventHandler(aodH);
+
+  // Output
+  
+  // Debug level
+  mgr->SetDebugLevel(0);
+
+  AliAnalysisTaskPi0Efficiency * task = new AliAnalysisTaskPi0Efficiency("Pi0Efficiency") ;
+  
+  TFile *fBadMap = TFile::Open("BadMap_LHC10h.root");
+  if(fBadMap->IsOpen()){
+    printf("\n\n...Adding PHOS bad channel map \n") ;
+    gROOT->cd();
+    char key[55] ;
+    for(Int_t mod=1;mod<4; mod++){
+      sprintf(key,"PHOS_BadMap_mod%d",mod) ;
+      TH2I * h = (TH2I*)fBadMap->Get(key) ;
+      if(h)
+       task->SetPHOSBadMap(mod,h) ;
+    }
+    fBadMap->Close() ;
+  }
+
+
+  mgr->AddTask(task);
+
+
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput   = mgr->GetCommonInputContainer(); 
+  AliAnalysisDataContainer *coutput = mgr->CreateContainer("histESD",TList::Class(),AliAnalysisManager::kOutputContainer,"histos.root");
+  
+  // Connect input/output
+  mgr->ConnectInput(task , 0, cinput);
+  mgr->ConnectOutput(task, 1,coutput);
+
+  if (mgr->InitAnalysis()) {
+    mgr->PrintStatus();
+    mgr->StartAnalysis("local", chain);
+  }
+  
+}
diff --git a/PWG4/UserTasks/PHOS_embedding/AnalyzeDiff.C b/PWG4/UserTasks/PHOS_embedding/AnalyzeDiff.C
new file mode 100644 (file)
index 0000000..6d7acef
--- /dev/null
@@ -0,0 +1,68 @@
+void AnalyzeDiff()
+{
+    
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libPhysics.so");
+  
+  //load analysis framework
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice"); //AliAnalysisTaskSE
+  
+  gSystem->AddIncludePath("-I$ALICE_ROOT/include -I$ALICE_ROOT/PHOS");
+
+  // A task can be compiled dynamically with AClic
+  gROOT->LoadMacro("AliCaloPhoton.cxx+g");
+  gROOT->LoadMacro("AliAnalysisTaskPi0DiffEfficiency.cxx+g");
+  
+  // Create the chain
+  TChain* chain = new TChain("aodTree");
+  chain->AddFile("AliAODout.root") ;
+
+  // Make the analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("Pi0EmbeddingManager");
+  
+  // ESD input handler
+  AliAODInputHandler* aodH = new AliAODInputHandler();
+  mgr->SetInputEventHandler(aodH);
+
+  // Output
+  
+  // Debug level
+  mgr->SetDebugLevel(0);
+
+  AliAnalysisTaskPi0DiffEfficiency * task = new AliAnalysisTaskPi0DiffEfficiency("Pi0Efficiency") ;
+  
+  TFile *fBadMap = TFile::Open("BadMap_LHC10h.root");
+  if(fBadMap->IsOpen()){
+    printf("\n\n...Adding PHOS bad channel map \n") ;
+    gROOT->cd();
+    char key[55] ;
+    for(Int_t mod=1;mod<4; mod++){
+      sprintf(key,"PHOS_BadMap_mod%d",mod) ;
+      TH2I * h = (TH2I*)fBadMap->Get(key) ;
+      if(h)
+       task->SetPHOSBadMap(mod,h) ;
+    }
+    fBadMap->Close() ;
+  }
+
+
+  mgr->AddTask(task);
+
+
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput   = mgr->GetCommonInputContainer(); 
+  AliAnalysisDataContainer *coutput = mgr->CreateContainer("histESD",TList::Class(),AliAnalysisManager::kOutputContainer,"histosDiff.root");
+  
+  // Connect input/output
+  mgr->ConnectInput(task , 0, cinput);
+  mgr->ConnectOutput(task, 1,coutput);
+
+  if (mgr->InitAnalysis()) {
+    mgr->PrintStatus();
+    mgr->StartAnalysis("local", chain);
+  }
+  
+}
diff --git a/PWG4/UserTasks/PHOS_embedding/Config.C b/PWG4/UserTasks/PHOS_embedding/Config.C
new file mode 100644 (file)
index 0000000..f0eb2f5
--- /dev/null
@@ -0,0 +1,697 @@
+//
+// Configuration for the first physics production 2008
+//
+
+// One can use the configuration macro in compiled mode by
+// root [0] gSystem->Load("libgeant321");
+// root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\
+//                   -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
+// root [0] .x grun.C(1,"Config.C++")
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include <TRandom.h>
+#include <TDatime.h>
+#include <TSystem.h>
+#include <TVirtualMC.h>
+#include <TGeant3TGeo.h>
+#include "STEER/AliRunLoader.h"
+#include "STEER/AliRun.h"
+#include "STEER/AliConfig.h"
+#include "PYTHIA6/AliDecayerPythia.h"
+#include "PYTHIA6/AliGenPythia.h"
+#include "TDPMjet/AliGenDPMjet.h"
+#include "STEER/AliMagFCheb.h"
+#include "STRUCT/AliBODY.h"
+#include "STRUCT/AliMAG.h"
+#include "STRUCT/AliABSOv3.h"
+#include "STRUCT/AliDIPOv3.h"
+#include "STRUCT/AliHALLv3.h"
+#include "STRUCT/AliFRAMEv2.h"
+#include "STRUCT/AliSHILv3.h"
+#include "STRUCT/AliPIPEv3.h"
+#include "ITS/AliITSv11Hybrid.h"
+#include "TPC/AliTPCv2.h"
+#include "TOF/AliTOFv6T0.h"
+#include "HMPID/AliHMPIDv3.h"
+#include "ZDC/AliZDCv3.h"
+#include "TRD/AliTRDv1.h"
+#include "TRD/AliTRDgeometry.h"
+#include "FMD/AliFMDv1.h"
+#include "MUON/AliMUONv1.h"
+#include "PHOS/AliPHOSv1.h"
+#include "PHOS/AliPHOSSimParam.h"
+#include "PMD/AliPMDv1.h"
+#include "T0/AliT0v1.h"
+#include "EMCAL/AliEMCALv2.h"
+#include "ACORDE/AliACORDEv1.h"
+#include "VZERO/AliVZEROv7.h"
+#endif
+
+
+enum PDC06Proc_t 
+{
+  kPythia6, kPythia6D6T, kPythia6ATLAS, kPythia6ATLAS_Flat, kPythiaPerugia0, kPhojet, kRunMax
+};
+
+const char * pprRunName[] = {
+  "kPythia6", "kPythia6D6T", "kPythia6ATLAS", "kPythia6ATLAS_Flat", "kPythiaPerugia0", "kPhojet" 
+};
+
+enum Mag_t
+{
+  kNoField, k5kG, kFieldMax
+};
+
+const char * pprField[] = {
+  "kNoField", "k5kG"
+};
+
+//--- Functions ---
+class AliGenPythia;
+AliGenerator *MbPythia();
+AliGenerator *MbPythiaTuneD6T();
+AliGenerator *MbPhojet();
+void ProcessEnvironmentVars();
+
+// Geterator, field, beam energy
+static PDC06Proc_t   proc     = kPhojet;
+static Mag_t         mag      = k5kG;
+static Float_t       energy   = 10000; // energy in CMS
+static Int_t         runNumber = 0;
+
+//========================//
+// Set Random Number seed //
+//========================//
+TDatime dt;
+static UInt_t seed    = dt.Get();
+
+// Comment line
+static TString comment;
+
+void Config()
+{
+    
+
+  // Get settings from environment variables
+  ProcessEnvironmentVars();
+
+  gRandom->SetSeed(seed);
+  cerr<<"Seed for random number generation= "<<seed<<endl; 
+
+  // Libraries required by geant321
+#if defined(__CINT__)
+  gSystem->Load("liblhapdf");      // Parton density functions
+  gSystem->Load("libEGPythia6");   // TGenerator interface
+  if (proc == kPythia6 || proc == kPhojet) {
+    gSystem->Load("libpythia6");        // Pythia 6.2
+  } else {
+    gSystem->Load("libpythia6.4.21");   // Pythia 6.4
+  }
+  gSystem->Load("libAliPythia6");  // ALICE specific implementations
+  gSystem->Load("libgeant321");
+#endif
+
+  new TGeant3TGeo("C++ Interface to Geant3");
+
+  //=======================================================================
+  //  Create the output file
+
+   
+  AliRunLoader* rl=0x0;
+
+  cout<<"Config.C: Creating Run Loader ..."<<endl;
+  rl = AliRunLoader::Open("galice.root",
+                         AliConfig::GetDefaultEventFolderName(),
+                         "recreate");
+  if (rl == 0x0)
+    {
+      gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
+      return;
+    }
+  rl->SetCompressionLevel(2);
+  rl->SetNumberOfEventsPerFile(3000);
+  gAlice->SetRunLoader(rl);
+  // gAlice->SetGeometryFromFile("geometry.root");
+  // gAlice->SetGeometryFromCDB();
+  
+  // Set the trigger configuration: proton-proton
+ // gAlice->SetTriggerDescriptor("p-p");
+
+  //
+  //=======================================================================
+  // ************* STEERING parameters FOR ALICE SIMULATION **************
+  // --- Specify event type to be tracked through the ALICE setup
+  // --- All positions are in cm, angles in degrees, and P and E in GeV
+
+
+    gMC->SetProcess("DCAY",1);
+    gMC->SetProcess("PAIR",1);
+    gMC->SetProcess("COMP",1);
+    gMC->SetProcess("PHOT",1);
+    gMC->SetProcess("PFIS",0);
+    gMC->SetProcess("DRAY",0);
+    gMC->SetProcess("ANNI",1);
+    gMC->SetProcess("BREM",1);
+    gMC->SetProcess("MUNU",1);
+    gMC->SetProcess("CKOV",1);
+    gMC->SetProcess("HADR",1);
+    gMC->SetProcess("LOSS",2);
+    gMC->SetProcess("MULS",1);
+    gMC->SetProcess("RAYL",1);
+
+    Float_t cut = 1.e-3;        // 1MeV cut by default
+    Float_t tofmax = 1.e10;
+
+    gMC->SetCut("CUTGAM", cut);
+    gMC->SetCut("CUTELE", cut);
+    gMC->SetCut("CUTNEU", cut);
+    gMC->SetCut("CUTHAD", cut);
+    gMC->SetCut("CUTMUO", cut);
+    gMC->SetCut("BCUTE",  cut); 
+    gMC->SetCut("BCUTM",  cut); 
+    gMC->SetCut("DCUTE",  cut); 
+    gMC->SetCut("DCUTM",  cut); 
+    gMC->SetCut("PPCUTM", cut);
+    gMC->SetCut("TOFMAX", tofmax); 
+
+
+
+
+  //======================//
+  // Set External decayer //
+  //======================//
+  TVirtualMCDecayer* decayer = new AliDecayerPythia();
+  decayer->SetForceDecay(kAll);
+  decayer->Init();
+  gMC->SetExternalDecayer(decayer);
+
+  //=========================//
+  // Generator Configuration //
+  //=========================//
+/*
+  AliGenerator* gener = 0x0;
+  
+  if (proc == kPythia6) {
+      gener = MbPythia();
+  } else if (proc == kPythia6D6T) {
+      gener = MbPythiaTuneD6T();
+  } else if (proc == kPythia6ATLAS) {
+      gener = MbPythiaTuneATLAS();
+  } else if (proc == kPythiaPerugia0) {
+      gener = MbPythiaTunePerugia0();
+  } else if (proc == kPythia6ATLAS_Flat) {
+      gener = MbPythiaTuneATLAS_Flat();
+  } else if (proc == kPhojet) {
+      gener = MbPhojet();
+  }
+*/  
+/*
+   AliGenCocktail *gener  = new AliGenCocktail();
+
+   AliGenPHOSlib * lib = new AliGenPHOSlib() ;
+   //4pi0
+   AliGenParam *genPHOS = new AliGenParam(4,AliGenPHOSlib::kPion,lib->GetPt(AliGenPHOSlib::kPi0),lib->GetY(AliGenPHOSlib::kPi0Flat),lib->GetIp(AliGenPHOSlib::kPi0)) ;
+   genPHOS->SetPhiRange(250.,330.) ;
+   genPHOS->SetYRange(-0.15.,0.15) ;
+   genPHOS->SetPtRange(0.5,30.) ;
+   gener->AddGenerator(genPHOS,"PHOS",1.) ;
+*/
+
+printf("Creating FLAT generator \n") ;
+
+   gener = new AliGenBox(2) ;
+   gener->SetPhiRange(250.,330.) ;
+   gener->SetYRange(-0.15.,0.15) ;
+   gener->SetEtaRange(-0.15,0.15) ;
+   gener->SetPtRange(0.5,30.) ;
+   gener->SetPart(111) ;
+
+  //
+  //
+  // Size of the interaction diamond
+  // Longitudinal
+  Float_t sigmaz  = 5.4 / TMath::Sqrt(2.); // [cm]
+  if (energy == 900)
+    //sigmaz  = 10.5 / TMath::Sqrt(2.); // [cm]
+    //sigmaz = 3.7;
+  if (energy == 7000)
+    sigmaz  = 6.3 / TMath::Sqrt(2.); // [cm]
+  
+  //
+  // Transverse
+  // beta*
+  Float_t betast  = 10.;                 // beta* [m]
+  if (runNumber >= 117048) betast = 2.;
+  printf("beta* for run# %8d is %13.3f", runNumber, betast);
+  //
+  Float_t eps     = 3.75e-6;            // emittance [m]
+  Float_t gamma   = energy / 2.0 / 0.938272;  // relativistic gamma [1]
+  Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.;  // [cm]
+  printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
+    
+  gener->SetSigma(sigmaxy, sigmaxy, sigmaz);      // Sigma in (X,Y,Z) (cm) on IP position
+  gener->SetVertexSmear(kPerEvent);
+  gener->Init();
+
+  printf("\n \n Comment: %s \n \n", comment.Data());
+
+  rl->CdGAFile();
+  
+  Int_t iABSO  = 1;
+  Int_t iACORDE= 0;
+  Int_t iDIPO  = 1;
+  Int_t iEMCAL = 1;
+  Int_t iFMD   = 1;
+  Int_t iFRAME = 1;
+  Int_t iHALL  = 1;
+  Int_t iITS   = 1;
+  Int_t iMAG   = 1;
+  Int_t iMUON  = 1;
+  Int_t iPHOS  = 1;
+  Int_t iPIPE  = 1;
+  Int_t iPMD   = 1;
+  Int_t iHMPID = 1;
+  Int_t iSHIL  = 1;
+  Int_t iT0    = 1;
+  Int_t iTOF   = 1;
+  Int_t iTPC   = 1;
+  Int_t iTRD   = 1;
+  Int_t iVZERO = 1;
+  Int_t iZDC   = 1;
+  
+
+    //=================== Alice BODY parameters =============================
+    AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
+
+
+    if (iMAG)
+    {
+        //=================== MAG parameters ============================
+        // --- Start with Magnet since detector layouts may be depending ---
+        // --- on the selected Magnet dimensions ---
+        AliMAG *MAG = new AliMAG("MAG", "Magnet");
+    }
+
+
+    if (iABSO)
+    {
+        //=================== ABSO parameters ============================
+        AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");
+    }
+
+    if (iDIPO)
+    {
+        //=================== DIPO parameters ============================
+
+        AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");
+    }
+
+    if (iHALL)
+    {
+        //=================== HALL parameters ============================
+
+        AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");
+    }
+
+
+    if (iFRAME)
+    {
+        //=================== FRAME parameters ============================
+
+        AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
+       FRAME->SetHoles(1);
+    }
+
+    if (iSHIL)
+    {
+        //=================== SHIL parameters ============================
+
+        AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");
+    }
+
+
+    if (iPIPE)
+    {
+        //=================== PIPE parameters ============================
+
+        AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
+    }
+    if (iITS)
+    {
+        //=================== ITS parameters ============================
+
+       AliITS *ITS  = new AliITSv11Hybrid("ITS","ITS v11Hybrid");
+    }
+
+    if (iTPC)
+    {
+      //============================ TPC parameters =====================
+
+        AliTPC *TPC = new AliTPCv2("TPC", "Default");
+    }
+
+
+    if (iTOF) {
+        //=================== TOF parameters ============================
+
+       AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
+    }
+
+
+    if (iHMPID)
+    {
+        //=================== HMPID parameters ===========================
+
+        AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");
+
+    }
+
+
+    if (iZDC)
+    {
+        //=================== ZDC parameters ============================
+
+        AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
+    }
+
+    if (iTRD)
+    {
+        //=================== TRD parameters ============================
+
+        AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
+        AliTRDgeometry *geoTRD = TRD->GetGeometry();
+       // Partial geometry: modules at 0,1,7,8,9,16,17
+       // starting at 3h in positive direction
+       geoTRD->SetSMstatus(2,0);
+       geoTRD->SetSMstatus(3,0);
+       geoTRD->SetSMstatus(4,0);
+        geoTRD->SetSMstatus(5,0);
+       geoTRD->SetSMstatus(6,0);
+        geoTRD->SetSMstatus(11,0);
+        geoTRD->SetSMstatus(12,0);
+        geoTRD->SetSMstatus(13,0);
+        geoTRD->SetSMstatus(14,0);
+        geoTRD->SetSMstatus(15,0);
+        geoTRD->SetSMstatus(16,0);
+    }
+
+    if (iFMD)
+    {
+        //=================== FMD parameters ============================
+
+       AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
+   }
+
+    if (iMUON)
+    {
+        //=================== MUON parameters ===========================
+        // New MUONv1 version (geometry defined via builders)
+
+        AliMUON *MUON = new AliMUONv1("MUON", "default");
+        
+        // activate trigger efficiency by cells
+
+        MUON->SetTriggerEffCells(1);     
+
+     }
+
+    if (iPHOS)
+    {
+        //=================== PHOS parameters ===========================
+
+     AliPHOS *PHOS = new AliPHOSv1("PHOS", "noCPV_Modules123");
+
+    }
+
+
+    if (iPMD)
+    {
+        //=================== PMD parameters ============================
+
+        AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
+    }
+
+    if (iT0)
+    {
+        //=================== T0 parameters ============================
+        AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
+    }
+
+    if (iEMCAL)
+    {
+        //=================== EMCAL parameters ============================
+
+        AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEAR");
+    }
+
+     if (iACORDE)
+    {
+        //=================== ACORDE parameters ============================
+
+        AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
+    }
+
+     if (iVZERO)
+    {
+        //=================== ACORDE parameters ============================
+
+        AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
+    }
+}
+//
+//           PYTHIA
+//
+
+AliGenerator* MbPythia()
+{
+      comment = comment.Append(" pp: Pythia low-pt");
+//
+//    Pythia
+      AliGenPythia* pythia = new AliGenPythia(-1); 
+      pythia->SetMomentumRange(0, 999999.);
+      pythia->SetThetaRange(0., 180.);
+      pythia->SetYRange(-12.,12.);
+      pythia->SetPtRange(0,1000.);
+      pythia->SetProcess(kPyMb);
+      pythia->SetEnergyCMS(energy);
+      
+      return pythia;
+}
+
+AliGenerator* MbPythiaTuneD6T()
+{
+      comment = comment.Append(" pp: Pythia low-pt");
+//
+//    Pythia
+      AliGenPythia* pythia = new AliGenPythia(-1); 
+      pythia->SetMomentumRange(0, 999999.);
+      pythia->SetThetaRange(0., 180.);
+      pythia->SetYRange(-12.,12.);
+      pythia->SetPtRange(0,1000.);
+      pythia->SetProcess(kPyMb);
+      pythia->SetEnergyCMS(energy);
+//    Tune
+//    109     D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally)
+      pythia->SetTune(109); // F I X 
+      pythia->SetStrucFunc(kCTEQ6l);
+//
+      return pythia;
+}
+
+AliGenerator* MbPythiaTunePerugia0()
+{
+      comment = comment.Append(" pp: Pythia low-pt (Perugia0)");
+//
+//    Pythia
+      AliGenPythia* pythia = new AliGenPythia(-1); 
+      pythia->SetMomentumRange(0, 999999.);
+      pythia->SetThetaRange(0., 180.);
+      pythia->SetYRange(-12.,12.);
+      pythia->SetPtRange(0,1000.);
+      pythia->SetProcess(kPyMb);
+      pythia->SetEnergyCMS(energy);
+//    Tune
+//    320     Perugia 0
+      pythia->SetTune(320); 
+      pythia->UseNewMultipleInteractionsScenario();
+//
+      return pythia;
+}
+
+
+AliGenerator* MbPythiaTuneATLAS()
+{
+      comment = comment.Append(" pp: Pythia low-pt");
+//
+//    Pythia
+      AliGenPythia* pythia = new AliGenPythia(-1); 
+      pythia->SetMomentumRange(0, 999999.);
+      pythia->SetThetaRange(0., 180.);
+      pythia->SetYRange(-12.,12.);
+      pythia->SetPtRange(0,1000.);
+      pythia->SetProcess(kPyMb);
+      pythia->SetEnergyCMS(energy);
+//    Tune
+//    C   306 ATLAS-CSC: Arthur Moraes' (new) ATLAS tune (needs CTEQ6L externally)
+      pythia->SetTune(306);
+      pythia->SetStrucFunc(kCTEQ6l);
+//
+      return pythia;
+}
+
+AliGenerator* MbPythiaTuneATLAS_Flat()
+{
+      AliGenPythia* pythia = MbPythiaTuneATLAS();
+      
+      comment = comment.Append("; flat multiplicity distribution");
+      
+      // set high multiplicity trigger
+      // this weight achieves a flat multiplicity distribution
+      TH1 *weight = new TH1D("weight","weight",201,-0.5,200.5);
+      weight->SetBinContent(1,5.49443);
+      weight->SetBinContent(2,8.770816);
+      weight->SetBinContent(6,0.4568624);
+      weight->SetBinContent(7,0.2919915);
+      weight->SetBinContent(8,0.6674189);
+      weight->SetBinContent(9,0.364737);
+      weight->SetBinContent(10,0.8818444);
+      weight->SetBinContent(11,0.531885);
+      weight->SetBinContent(12,1.035197);
+      weight->SetBinContent(13,0.9394057);
+      weight->SetBinContent(14,0.9643193);
+      weight->SetBinContent(15,0.94543);
+      weight->SetBinContent(16,0.9426507);
+      weight->SetBinContent(17,0.9423649);
+      weight->SetBinContent(18,0.789456);
+      weight->SetBinContent(19,1.149026);
+      weight->SetBinContent(20,1.100491);
+      weight->SetBinContent(21,0.6350525);
+      weight->SetBinContent(22,1.351941);
+      weight->SetBinContent(23,0.03233504);
+      weight->SetBinContent(24,0.9574557);
+      weight->SetBinContent(25,0.868133);
+      weight->SetBinContent(26,1.030998);
+      weight->SetBinContent(27,1.08897);
+      weight->SetBinContent(28,1.251382);
+      weight->SetBinContent(29,0.1391099);
+      weight->SetBinContent(30,1.192876);
+      weight->SetBinContent(31,0.448944);
+      weight->SetBinContent(32,1);
+      weight->SetBinContent(33,1);
+      weight->SetBinContent(34,1);
+      weight->SetBinContent(35,1);
+      weight->SetBinContent(36,0.9999997);
+      weight->SetBinContent(37,0.9999997);
+      weight->SetBinContent(38,0.9999996);
+      weight->SetBinContent(39,0.9999996);
+      weight->SetBinContent(40,0.9999995);
+      weight->SetBinContent(41,0.9999993);
+      weight->SetBinContent(42,1);
+      weight->SetBinContent(43,1);
+      weight->SetBinContent(44,1);
+      weight->SetBinContent(45,1);
+      weight->SetBinContent(46,1);
+      weight->SetBinContent(47,0.9999999);
+      weight->SetBinContent(48,0.9999998);
+      weight->SetBinContent(49,0.9999998);
+      weight->SetBinContent(50,0.9999999);
+      weight->SetBinContent(51,0.9999999);
+      weight->SetBinContent(52,0.9999999);
+      weight->SetBinContent(53,0.9999999);
+      weight->SetBinContent(54,0.9999998);
+      weight->SetBinContent(55,0.9999998);
+      weight->SetBinContent(56,0.9999998);
+      weight->SetBinContent(57,0.9999997);
+      weight->SetBinContent(58,0.9999996);
+      weight->SetBinContent(59,0.9999995);
+      weight->SetBinContent(60,1);
+      weight->SetBinContent(61,1);
+      weight->SetBinContent(62,1);
+      weight->SetBinContent(63,1);
+      weight->SetBinContent(64,1);
+      weight->SetBinContent(65,0.9999999);
+      weight->SetBinContent(66,0.9999998);
+      weight->SetBinContent(67,0.9999998);
+      weight->SetBinContent(68,0.9999999);
+      weight->SetBinContent(69,1);
+      weight->SetBinContent(70,1);
+      weight->SetBinContent(71,0.9999997);
+      weight->SetBinContent(72,0.9999995);
+      weight->SetBinContent(73,0.9999994);
+      weight->SetBinContent(74,1);
+      weight->SetBinContent(75,1);
+      weight->SetBinContent(76,1);
+      weight->SetBinContent(77,1);
+      weight->SetBinContent(78,0.9999999);
+      weight->SetBinContent(79,1);
+      weight->SetBinContent(80,1);
+      weight->SetEntries(526);
+        
+      Int_t limit = weight->GetRandom();
+      pythia->SetTriggerChargedMultiplicity(limit, 1.4);
+      
+      comment = comment.Append(Form("; multiplicity threshold set to %d in |eta| < 1.4", limit));
+
+      return pythia;
+}
+
+AliGenerator* MbPhojet()
+{
+      comment = comment.Append(" pp: Pythia low-pt");
+//
+//    DPMJET
+#if defined(__CINT__)
+  gSystem->Load("libdpmjet");      // Parton density functions
+  gSystem->Load("libTDPMjet");      // Parton density functions
+#endif
+      AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); 
+      dpmjet->SetMomentumRange(0, 999999.);
+      dpmjet->SetThetaRange(0., 180.);
+      dpmjet->SetYRange(-12.,12.);
+      dpmjet->SetPtRange(0,1000.);
+      dpmjet->SetProcess(kDpmMb);
+      dpmjet->SetEnergyCMS(energy);
+
+      return dpmjet;
+}
+
+void ProcessEnvironmentVars()
+{
+    // Run type
+    if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
+      for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
+       if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) {
+         proc = (PDC06Proc_t)iRun;
+         cout<<"Run type set to "<<pprRunName[iRun]<<endl;
+       }
+      }
+    }
+
+    // Field
+    if (gSystem->Getenv("CONFIG_FIELD")) {
+      for (Int_t iField = 0; iField < kFieldMax; iField++) {
+       if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) {
+         mag = (Mag_t)iField;
+         cout<<"Field set to "<<pprField[iField]<<endl;
+       }
+      }
+    }
+
+    // Energy
+    if (gSystem->Getenv("CONFIG_ENERGY")) {
+      energy = atoi(gSystem->Getenv("CONFIG_ENERGY"));
+      cout<<"Energy set to "<<energy<<" GeV"<<endl;
+    }
+
+    // Random Number seed
+    if (gSystem->Getenv("CONFIG_SEED")) {
+      seed = atoi(gSystem->Getenv("CONFIG_SEED"));
+    }
+
+    // Run number
+    if (gSystem->Getenv("DC_RUN")) {
+      runNumber = atoi(gSystem->Getenv("DC_RUN"));
+    }
+}
diff --git a/PWG4/UserTasks/PHOS_embedding/CreateAOD.C b/PWG4/UserTasks/PHOS_embedding/CreateAOD.C
new file mode 100644 (file)
index 0000000..e292050
--- /dev/null
@@ -0,0 +1,115 @@
+void CreateAOD()
+{
+  // Main
+  
+  // LoadLibraries
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libXMLIO.so");
+  gSystem->Load("libMatrix.so");
+  gSystem->Load("libPhysics.so");
+  gSystem->Load("libSTEERBase.so");
+  gSystem->Load("libESD.so");
+  gSystem->Load("libAOD.so");
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libCORRFW.so");
+  gSystem->Load("libPWG3base.so");
+  gSystem->Load("libPWG3muon.so");
+  gSystem->Load("libTENDER.so");
+//gSystem->Load("libTENDERSupplies.so"); 
+
+  gSystem->AddIncludePath("-I$ALICE_ROOT/include -I$ALICE_ROOT/PHOS");
+  gROOT->LoadMacro("AliPHOSDigitDecalibrate.cxx++") ;
+
+  //Data
+  TChain *chain       = new TChain("esdTree") ;
+  chain->Add("AliESDs.root");
+  
+  if(chain){
+    AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
+    
+    //--------------------------------------
+    // Make the analysis manager
+    //-------------------------------------
+    AliAnalysisManager *mgr  = new AliAnalysisManager("Manager", "Manager");
+    
+    AliMCEventHandler* mcHandler = new AliMCEventHandler();
+    mcHandler->SetReadTR(kFALSE);//Do not search TrackRef file
+    mgr->SetMCtruthEventHandler(mcHandler);
+    
+    AliAODHandler* aodoutHandler   = new AliAODHandler();
+    aodoutHandler->SetOutputFileName("AliAOD.root");
+    ////aodoutHandler->SetCreateNonStandardAOD();
+    mgr->SetOutputEventHandler(aodoutHandler);
+    
+    AliESDInputHandler *esdHandler = new AliESDInputHandler();
+    esdHandler->SetReadFriends(kFALSE);
+    mgr->SetInputEventHandler(esdHandler);
+    
+    //mgr->SetDebugLevel(-1); // For debugging, do not uncomment if you want no messages.
+    
+    //-------------------------------------------------------------------------
+    //Define task, put here any other task that you want to use.
+    //-------------------------------------------------------------------------
+    AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
+    AliAnalysisDataContainer *coutput1;
+    
+    coutput1 = mgr->GetCommonOutputContainer();
+
+
+
+    //========= Add tender to the ANALYSIS manager and set default storage =====
+
+    AliTender *tender=new AliTender("AnalysisTender");
+    tender->SetDefaultCDBStorage("local://$ALICE_ROOT/OCDB"); //comment this you acess grid
+    mgr->AddTask(tender);
+
+    //========= Attach PHOS supply ======
+    AliPHOSDigitDecalibrate *PHOSSupply=new AliPHOSDigitDecalibrate("PHOStender");
+    tender->AddSupply(PHOSSupply);
+    TFile * fDecalib = TFile::Open("Decalibraiton.root") ;
+    char key[55] ;
+    for(Int_t mod=0; mod<5; mod++){
+      sprintf(key,"DecalibrMod%d",mod) ;
+      TH2F * h = (TH2F*)fDecalib->Get(key) ;
+      PHOSSupply->SetDecalibration(mod,h) ;
+    }
+    fDecalib->Close() ;
+
+    //================================================
+    //              data containers
+    //================================================
+
+    //            define output containers, please use 'username'_'somename'
+    AliAnalysisDataContainer *coutput1 =
+      mgr->CreateContainer("tender_event", AliESDEvent::Class(),
+                           AliAnalysisManager::kExchangeContainer,"default_tender");
+    //           connect containers
+    mgr->ConnectInput  (tender,  0, mgr->GetCommonInputContainer() );
+    mgr->ConnectOutput (tender,  1, coutput1);
+
+
+
+
+
+
+    
+    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C");
+    AliAnalysisTaskESDfilter *taskesdfilter = AddTaskESDFilter(kTRUE);
+    
+    //-----------------------
+    // Run the analysis
+    //-----------------------    
+    mgr->InitAnalysis();
+    mgr->PrintStatus();
+    mgr->StartAnalysis("local",chain);
+    
+    cout <<" Analysis ended sucessfully "<< endl ;
+    
+  }
+  else cout << "Chain was not produced ! "<<endl;
+  
+}
diff --git a/PWG4/UserTasks/PHOS_embedding/Embedding.C b/PWG4/UserTasks/PHOS_embedding/Embedding.C
new file mode 100644 (file)
index 0000000..5a5a786
--- /dev/null
@@ -0,0 +1,180 @@
+void Embedding(const char* dataset="collection.xml")
+{
+    
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libPhysics.so");
+  
+  //load analysis framework
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice"); //AliAnalysisTaskSE
+  
+  gSystem->AddIncludePath("-I$ALICE_ROOT/include -I$ALICE_ROOT/PHOS");
+
+  // A task can be compiled dynamically with AClic
+  gROOT->LoadMacro("AliPHOSEmbedding.cxx+g");
+  
+  // Connect to alien
+  TString token = gSystem->Getenv("GRID_TOKEN") ;
+  if (1) // token == "OK" ) 
+    TGrid::Connect("alien://");
+  else 
+    AliInfo("You are not connected to the GRID") ; 
+  cout << "Pi0Analysis: processing collection " << dataset << endl;
+  
+  // Create the chain
+  TChain* chain = new TChain("esdTree");
+
+  chain->Add("/home/prsnko/PbPb/Embedding/Data/AliESDs.root") ;
+/*
+
+  TGridCollection * collection = dynamic_cast<TGridCollection*>(TAlienCollection::Open(dataset));
+  
+  TAlienResult* result = collection->GetGridResult("",0 ,0);
+  TList* rawFileList = result->GetFileInfoList();
+  for (Int_t counter=0 ; counter < rawFileList->GetEntries() ; counter++) {
+    TFileInfo * fi =  static_cast<TFileInfo*>(rawFileList->At(counter)) ; 
+    const char * rawFile = fi->GetCurrentUrl()->GetUrl() ;  
+    printf("Processing %s\n", rawFile) ;
+    chain->Add(rawFile);
+    printf("Chain: %d entries.\n",chain->GetEntries()); 
+  }
+  TFileInfo * fi =  static_cast<TFileInfo*>(rawFileList->At(0));
+  const char * fn = fi->GetCurrentUrl()->GetUrl() ;
+*/
+  char runNum[7]={"139438"}; 
+//  for(Int_t i=0;i<6;i++)runNum[i]=fn[35+i] ;
+  runNum[6]=0 ;
+  Int_t iRunNum=atoi(runNum) ;
+  printf("Run number=%d \n",iRunNum) ;
+
+  //Run AOD simulation
+  int nrun = atoi(runNum);
+  int nevent = 0;
+  int seed = 0;
+
+  char sseed[1024];
+  char sevent[1024];
+  char sprocess[1024];
+  char sfield[1024];
+  char senergy[1024];
+
+  sprintf(sevent,"");
+  sprintf(sprocess,"");
+  sprintf(sfield,"");
+  sprintf(senergy,"");
+
+  seed = 0;
+  sprintf(sseed,"%d",seed);
+
+  if (seed==0) {
+    fprintf(stderr,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+    fprintf(stderr,"!!!!  WARNING! Seeding variable for MC is 0          !!!!\n");
+    fprintf(stderr,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+  } else {
+    fprintf(stdout,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+    fprintf(stdout,"!!!  MC Seed is %d \n",seed);
+    fprintf(stdout,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+  }
+  
+// set the seed environment variable
+  gSystem->Setenv("CONFIG_SEED",sseed);
+  gSystem->Setenv("CONFIG_RUN_TYPE","kPythia6"); // kPythia6 or kPhojet^M
+  gSystem->Setenv("CONFIG_FIELD","k5kG");      // kNoField or k5kG^M
+  gSystem->Setenv("CONFIG_ENERGY","2760");    // 900 or 10000 (GeV)
+  gSystem->Setenv("DC_RUN",runNum);    //run number 
+  
+
+  char nSimEvents[55] ;
+  sprintf(nSimEvents,"%d",chain->GetEntries());
+  gSystem->Setenv("SIM_EVENTS",nSimEvents); 
+  gSystem->Exec("mv geometry.root geometry_PHOS.root") ;
+//  gSystem->Exec("aliroot -b -q simrun.C > simrun.log 2>&1");
+  gSystem->Exec("mv geometry_PHOS.root geometry.root") ;
+
+  // Make the analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("Pi0EmbeddingManager");
+  
+  // ESD input handler
+  AliESDInputHandler* esdH = new AliESDInputHandler();
+  esdH->SetReadFriends(kFALSE);
+  mgr->SetInputEventHandler(esdH);
+
+  // Output
+  AliAODHandler* aodHandler   = new AliAODHandler();
+  aodHandler->SetOutputFileName("AliAODout.root");
+  mgr->SetOutputEventHandler(aodHandler);
+
+  
+  // Debug level
+  mgr->SetDebugLevel(0);
+
+  // Add physics selection
+//  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+//  AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(kFALSE,kFALSE);
+
+  //Add centrality
+  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
+  AliCentralitySelectionTask *taskCentrality = AddTaskCentrality() ;
+  taskCentrality->SetPass(2); // remember to set the pass you are processing!!! 
+
+  // Add my task
+  AliPHOSEmbedding *task1 = new AliPHOSEmbedding("Embedding");
+
+  TChain* chainAOD = new TChain("aodTree");
+  chainAOD->AddFile("AliAOD.root") ;
+  task1->SetSignalChain(chainAOD) ;
+ // task1->SelectCollisionCandidates();
+
+
+  TFile *fOldCalib = TFile::Open("OldCalibration.root");
+  if(fOldCalib->IsOpen()){
+    printf("\n\n...Adding PHOS calibration used in ESD production \n") ;
+    char key[55] ;
+    TH2F * hCalib[5] ;
+    for(Int_t mod=0;mod<5; mod++){
+      sprintf(key,"calibrationMod%d",mod) ;
+      hCalib[mod] = (TH2F*)fOldCalib->Get(key) ;
+    }
+    task1->SetOldCalibration(hCalib) ;
+    fOldCalib->Close() ;
+  }
+
+  mgr->AddTask(task1);
+  
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput   = mgr->GetCommonInputContainer(); 
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("output0",
+                      TTree::Class(), AliAnalysisManager::kOutputContainer);
+
+
+  // Connect input/output
+  mgr->ConnectInput(task1 , 0, cinput);
+  mgr->ConnectOutput(task1, 0,coutput1);
+
+  AliLog::SetClassDebugLevel("AliGeomManager", 10) ;
+
+  AliCDBManager::Instance()->SetRun(iRunNum) ;
+//  AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB") ;
+//  AliCDBManager::Instance()->SetDefaultStorage("raw://") ;
+  AliCDBManager::Instance()->SetDefaultStorage("alien://folder=/alice/data/2010/OCDB") ;
+  AliCDBManager::Instance()->SetSpecificStorage("PHOS/*/*",
+                                                "local://OCDB");
+printf("RunNunm===%d \n",iRunNum) ;
+  
+
+  if (mgr->InitAnalysis()) {
+    mgr->PrintStatus();
+    mgr->StartAnalysis("local", chain);
+  }
+  
+  if(iRunNum<=137848){ //period 1
+    gSystem->Exec("mv BadMap_LHC10h_period1.root BadMap_LHC10h.root") ;
+  }
+  else{
+    gSystem->Exec("mv BadMap_LHC10h_period234.root BadMap_LHC10h.root") ;
+  }
+//  gSystem->Exec("aliroot -b -q Analyze.C> analyze.log 2>&1");
+
+}
diff --git a/PWG4/UserTasks/PHOS_embedding/Readme b/PWG4/UserTasks/PHOS_embedding/Readme
new file mode 100644 (file)
index 0000000..4145cd5
--- /dev/null
@@ -0,0 +1,24 @@
+These are tasks to perform PHOS embedding. Steps are following:
+1. Steering macro Embeding.C reads list of data ESDs to embed and calculates number of events
+2. Performs simulation of the same number of events, using enclosed Config.C 
+   (single pi0 per event with flat or Hagedorn spectrum)
+   then simulated signals are reconstructed clusters and converted to AOD
+3. Then class AliPHOSEmbedding is called
+   - read data ESDs, extract PHOS calocells, convert them to PHOSDigits
+   - run default reconstruction to avoid possible problems with different 
+     OCDBs used now and in official reconstruction. Apply final calibration if necessary
+   - merge digits with (re-calibrated if necessary) signal digits
+   - reconstruct digits with embedding signal (with final calibration if necessary)
+4. As a result of all these procedures AOD is produced with
+   - MC information about embedded pi0s
+   - standard CaloClosters, CaloCells with final calibration
+   - Embedded CaloClusters, embedded CaloCells with final calibration
+5. Run analysis with
+   - AliAnalysisTaskPi0Efficiency: Signal (nominator in efficiency calculation) is any cluster with 
+     contribution from simulated pi0s. Only for cross-check, this appoach is not quite correct at high accupancy.
+   - AliAnalysisTaskPi0DiffEfficiency: correct approach, Signal (nominator in efficiency calculation) is
+     difference between Real inv. mass distribution after and before embedding.
+
+Note that to perform MC simulation one needs either access to official OCDB (Why do I need TPC drift velocity to simulate PHOS????)
+or local copy of OCDB. Also histograms with new/old calibrations should be provided if necessary.
+
diff --git a/PWG4/UserTasks/PHOS_embedding/rec.C b/PWG4/UserTasks/PHOS_embedding/rec.C
new file mode 100644 (file)
index 0000000..86bef89
--- /dev/null
@@ -0,0 +1,46 @@
+void rec() {
+
+  AliReconstruction reco;
+
+//
+// switch off cleanESD, write ESDfriends and Alignment data
+  
+//  reco.SetCleanESD(kFALSE);
+//  reco.SetWriteESDfriend();
+//  reco.SetWriteAlignmentData();
+    reco.SetRunVertexFinder(kFALSE) ;
+    reco.SetRunV0Finder(kFALSE) ;
+    reco.SetRunMultFinder(kFALSE) ;
+
+   AliQAManager::QAManager(AliQAv1::kSIMMODE) ;
+
+//
+// ITS Efficiency and tracking errors
+
+//
+// Residual OCDB
+
+  reco.SetDefaultStorage("alien://folder=/alice/data/2010/OCDB");
+//  reco.SetDefaultStorage("local://./OCDB");
+  reco.SetSpecificStorage("PHOS/*/*","local://./OCDB");
+
+
+//  reco.SetRunReconstruction("PHOS EMCAL") ;
+  reco.SetRunReconstruction("PHOS") ;
+  reco.SetRunQA(":") ;
+  reco.SetRunGlobalQA(kFALSE) ;
+//
+// GRP from local OCDB
+
+ reco.SetSpecificStorage("GRP/GRP/Data",
+                          Form("local://%s",gSystem->pwd()));
+//  reco.SetSpecificStorage("GRP/GRP/Data", "alien://Folder=/alice/data/2010/OCDB");
+
+
+//
+  TStopwatch timer;
+  timer.Start();
+  reco.Run();
+  timer.Stop();
+  timer.Print();
+}
diff --git a/PWG4/UserTasks/PHOS_embedding/sim.C b/PWG4/UserTasks/PHOS_embedding/sim.C
new file mode 100644 (file)
index 0000000..9782229
--- /dev/null
@@ -0,0 +1,47 @@
+void sim(Int_t nev=20) {
+
+  if (gSystem->Getenv("SIM_EVENTS"))
+    nev = atoi(gSystem->Getenv("SIM_EVENTS"));
+
+  printf("GENERATE << %d >> events \n",nev);
+
+
+  AliSimulation simulator;
+  simulator.SetMakeSDigits("PHOS");
+  simulator.SetMakeDigits("PHOS");
+//
+// Ideal OCDB
+//  simulator.SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+  simulator.SetDefaultStorage("local://./OCDB");
+//  simulator.SetSpecificStorage("GRP/GRP/Data",
+//                               Form("local://%s",gSystem->pwd()));
+
+//  simulator.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/");
+
+ //simulator.SetSpecificStorage("GRP/Calib/MeanVertexSPD", "alien://folder=/alice/data/2010/OCDB");
+
+  //PHOS bad map from RAW OCDB
+  simulator.SetSpecificStorage("PHOS/*/*/","local://./OCDB");
+//  simulator.SetSpecificStorage("PHOS/Calib/EmcBadChannels/","local://./OCDB");
+//  simulator.SetSpecificStorage("PHOS/Calib/EmcGainPedestals/","local://./OCDB");
+
+  simulator.SetRunHLT("");
+//
+
+  simulator.SetSpecificStorage("GRP/GRP/Data", "alien://Folder=/alice/data/2010/OCDB");
+
+// Vertex and Mag.field from OCDB
+
+//  simulator.UseVertexFromCDB();
+  simulator.UseMagFieldFromGRP();
+  simulator.SetRunQA(":") ;
+
+//
+// The rest
+
+  TStopwatch timer;
+  timer.Start();
+  simulator.Run(nev);
+  timer.Stop();
+  timer.Print();
+}
diff --git a/PWG4/UserTasks/PHOS_embedding/simrun.C b/PWG4/UserTasks/PHOS_embedding/simrun.C
new file mode 100644 (file)
index 0000000..eccadac
--- /dev/null
@@ -0,0 +1,16 @@
+// simrun.C\r
+{\r
+// set job and simulation variables as :\r
+// root.exe -b -q simrun.C  --run <x> --event <y> --process <kPythia6/kPhojet/kPythia6ATLAS_Flat/kPythia6D6T> --field <kNoField/k5kG> --energy <900/2360/10000>\r
+\r
+  cout<<">>>>> SIMULATION <<<<<"<<endl;\r
+  gSystem->Exec("aliroot -b -q sim.C > sim.log 2>&1");\r
+  gSystem->Exec("mv syswatch.log simwatch.log");\r
+  gSystem->Exec("rm *ITS*.root  *TPC*.root *TOF*.root *TRD*.root *PMD*.root ") ;\r
+  cout<<">>>>> RECONSTRUCTION <<<<<"<<endl;\r
+  gSystem->Exec("aliroot -b -q rec.C > rec.log 2>&1");\r
+  gSystem->Exec("mv syswatch.log recwatch.log");\r
+  cout<<">>>>> AOD <<<<<"<<endl;\r
+  gSystem->Exec("aliroot -b -q CreateAOD.C > aod.log 2>&1");\r
+\r
+}\r