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
#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+;
--- /dev/null
+#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 ;
+
+}
+
+
+
--- /dev/null
+#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
--- /dev/null
+#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
+
+
+}
+
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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) ;
+
+
+
+}
+
--- /dev/null
+#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
+
--- /dev/null
+// 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. ;
+}
--- /dev/null
+#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
--- /dev/null
+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);
+ }
+
+}
--- /dev/null
+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);
+ }
+
+}
--- /dev/null
+//
+// 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"));
+ }
+}
--- /dev/null
+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;
+
+}
--- /dev/null
+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");
+
+}
--- /dev/null
+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.
+
--- /dev/null
+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();
+}
--- /dev/null
+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();
+}
--- /dev/null
+// 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
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff