//______________________________________________________________________________________________
AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
-AliAnalysisTaskSE(name),fEMCALGeo(0x0),
+AliAnalysisTaskSE(name),
+fEMCALGeo(0x0), fLoadMatrices(0),
+fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
+fTriggerName("EMC"),
+fRecoUtils(new AliEMCALRecoUtils),
+fOADBFilePath(""), fCorrectClusters(kFALSE),
+fCaloClustersArr(0x0), fEMCALCells(0x0),
+fCuts(0x0), fOutputContainer(0x0),
+fVertex(), fFilteredInput(kFALSE),
fEmin(0.5), fEmax(15.),
fL0min(0.01), fL0max(0.5),
fDTimeCut(100.), fTimeMax(1000000), fTimeMin(-1000000),
fAsyCut(1.), fMinNCells(2), fGroupNCells(0),
-fLogWeight(4.5), fSameSM(kFALSE), fFilteredInput(kFALSE),
-fCorrectClusters(kFALSE), fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
-fTriggerName("EMC"), fOADBFilePath(""),
-fRecoUtils(new AliEMCALRecoUtils),
-fCuts(0x0), fLoadMatrices(0),
+fLogWeight(4.5), fSameSM(kFALSE),
fNMaskCellColumns(11), fMaskCellColumns(0x0),
fInvMassCutMin(110.), fInvMassCutMax(160.),
//Histograms
-fOutputContainer(0x0), fNbins(300),
+fNbins(300),
fMinBin(0.), fMaxBin(300.),
fNTimeBins(1000), fMinTimeBin(0.), fMaxTimeBin(1000.),
fHmgg(0x0), fHmggDifferentSM(0x0),
}
}
+ fVertex[0]=fVertex[1]=fVertex[2]=-1000;
+
fHTpi0[0]= 0 ;
fHTpi0[1]= 0 ;
fHTpi0[2]= 0 ;
}
-
-//________________________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices()
+//____________________________________________________________
+void AliAnalysisTaskEMCALPi0CalibSelection::CorrectClusters()
{
- // Init geometry and set the geometry matrix, for the first event, skip the rest
- // Also set once the run dependent calibrations
+ // loop over EMCAL clusters
+ //----------------------------------------------------------
+ // First recalibrate and recalculate energy and position
- if(fhNEvents->GetEntries()!=1) return;
-
- Int_t runnumber = InputEvent()->GetRunNumber() ;
- if(fLoadMatrices)
+ if(fCorrectClusters)
{
- printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - Load user defined EMCAL geometry matrices\n");
- // OADB if available
- AliOADBContainer emcGeoMat("AliEMCALgeo");
-
- if(fOADBFilePath=="") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;
+ if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton)
+ {
+ AliFatal(Form("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType()));
+ }
- emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
+ if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
- TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
+ Float_t pos[]={0,0,0};
- for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
+ for(Int_t iClu=0; iClu < fCaloClustersArr->GetEntriesFast(); iClu++)
{
+ AliVCluster *c1 = (AliVCluster *) fCaloClustersArr->At(iClu);
- if (!fMatrix[mod]) // Get it from OADB
- {
- if(fDebug > 1 )
- printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - EMCAL matrices SM %d, %p\n",
- mod,((TGeoHMatrix*) matEMCAL->At(mod)));
- //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
-
- fMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
- }
+ Float_t e1i = c1->E(); // cluster energy before correction
+ if (e1i < fEmin) continue;
+ else if (e1i > fEmax) continue;
- if(fMatrix[mod])
- {
- if(DebugLevel() > 1)
- fMatrix[mod]->Print();
-
- fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
+ else if (c1->GetNCells() < fMinNCells) continue;
+
+ else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
+
+ if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
+
+ if(DebugLevel() > 2)
+ {
+ printf("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
+ c1->GetPosition(pos);
+ printf("Std : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
}
-
- }//SM loop
- }//Load matrices
- else if(!gGeoManager)
- {
- printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - Get geo matrices from data");
- //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
- if(!strcmp(InputEvent()->GetName(),"AliAODEvent"))
- {
- if(DebugLevel() > 1)
- Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
- }//AOD
- else
- {
- if(DebugLevel() > 1)
- printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
- for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
+ //Correct cluster energy and position if requested, and not corrected previously, by default Off
+ if(fRecoUtils->IsRecalibrationOn())
{
- if(DebugLevel() > 1)
- InputEvent()->GetEMCALMatrix(mod)->Print();
-
- if(InputEvent()->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(InputEvent()->GetEMCALMatrix(mod),mod) ;
-
- }
-
- }//ESD
- }//Load matrices from Data
-
+ fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, fEMCALCells);
+ fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, fEMCALCells,c1);
+ fRecoUtils->RecalculateClusterPID(c1);
+ }
+
+ if(DebugLevel() > 2)
+ printf("Energy: after recalibration %f; \n",c1->E());
+
+ // Recalculate cluster position
+ fRecoUtils->RecalculateClusterPosition(fEMCALGeo, fEMCALCells,c1);
+
+ // Correct Non-Linearity
+ c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
+
+ if(DebugLevel() > 2)
+ printf("\t after linearity correction %f\n",c1->E());
+
+ //In case of MC analysis, to match resolution/calibration in real data
+ c1->SetE(fRecoUtils->SmearClusterEnergy(c1));
+
+ if(DebugLevel() > 2)
+ printf("\t after smearing %f\n",c1->E());
+
+ if(DebugLevel() > 2)
+ {
+ printf("Cor : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
+ c1->GetPosition(pos);
+ printf("Cor : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
+ }
+ }
+ }
}
-//___________________________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
+//__________________________________________________________
+void AliAnalysisTaskEMCALPi0CalibSelection::FillHistograms()
{
- //Create output container, init geometry
-
- fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
- Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
-
- fOutputContainer = new TList();
- const Int_t buffersize = 255;
- char hname[buffersize], htitl[buffersize];
-
- fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
- fOutputContainer->Add(fhNEvents);
-
- fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
- fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
- fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHmgg);
-
- fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
- fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
- fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHmggDifferentSM);
-
- fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
- fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
- fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHOpeningAngle);
-
- fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
- fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
- fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHOpeningAngleDifferentSM);
-
- fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
- fHAsymmetry->SetXTitle("a");
- fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHAsymmetry);
-
- fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
- fHAsymmetryDifferentSM->SetXTitle("a");
- fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHAsymmetryDifferentSM);
-
-
- //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
-
- fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
- fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
- fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHmggMaskFrame);
+ // Now fill the invariant mass analysis with the corrected clusters, and other general histograms
+
+ Int_t absId1 = -1;
+ Int_t iSupMod1 = -1;
+ Int_t iphi1 = -1;
+ Int_t ieta1 = -1;
+ Int_t absId2 = -1;
+ Int_t iSupMod2 = -1;
+ Int_t iphi2 = -1;
+ Int_t ieta2 = -1;
+ Bool_t shared = kFALSE;
- fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
- fNbins,fMinBin,fMaxBin,100,0,10);
- fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
- fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHmggDifferentSMMaskFrame);
+ TLorentzVector p1;
+ TLorentzVector p2;
+ TLorentzVector p12;
+ Float_t pos[]={0,0,0};
- for(Int_t iSM = 0; iSM < nSM; iSM++)
+ Int_t bc = InputEvent()->GetBunchCrossNumber();
+ Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
+
+ for(Int_t iClu=0; iClu<fCaloClustersArr->GetEntriesFast()-1; iClu++)
{
- snprintf(hname, buffersize, "hmgg_SM%d",iSM);
- snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
- fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
- fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
- fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHmggSM[iSM]);
+ AliVCluster *c1 = (AliVCluster *) fCaloClustersArr->At(iClu);
- snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
- snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
- fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
- fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
- fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHmggSMMaskFrame[iSM]);
+ if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
+ Float_t e1i = c1->E(); // cluster energy before correction
- if(iSM < nSM/2)
+ if (e1i < fEmin) continue;
+ else if (e1i > fEmax) continue;
+
+ else if (!fRecoUtils->IsGoodCluster(c1,fEMCALGeo,fEMCALCells,bc)) continue;
+
+ else if (c1->GetNCells() < fMinNCells) continue;
+
+ else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
+
+ if(DebugLevel() > 2)
+ {
+ printf("IMA : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
+ c1->GetPosition(pos);
+ printf("IMA : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
+ }
+
+ fRecoUtils->GetMaxEnergyCell(fEMCALGeo, fEMCALCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
+ c1->GetMomentum(p1,fVertex);
+
+ //Check if cluster is in fidutial region, not too close to borders
+ Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, fEMCALCells);
+
+ // Clusters not facing frame structures
+ Bool_t mask1 = MaskFrameCluster(iSupMod1, ieta1);
+ //if(mask1) printf("Reject eta %d SM %d\n",ieta1, iSupMod1);
+
+ Double_t time1 = c1->GetTOF()*1.e9;
+
+ if(time1 > fTimeMax || time1 < fTimeMin) continue;
+
+ fhClusterTime ->Fill(c1->E(),time1);
+ fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
+
+ // Combine cluster with other clusters and get the invariant mass
+ for (Int_t jClu=iClu+1; jClu < fCaloClustersArr->GetEntriesFast(); jClu++)
{
- snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
- snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
- fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
- fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
- fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHmggPairSameSectorSM[iSM]);
+ AliAODCaloCluster *c2 = (AliAODCaloCluster *) fCaloClustersArr->At(jClu);
- snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
- snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
- fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
- fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
- fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
- fOutputContainer->Add(fHmggPairSameSectorSMMaskFrame[iSM]);
+ Float_t e2i = c2->E();
+ if (e2i < fEmin) continue;
+ else if (e2i > fEmax) continue;
- fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
- Form("cluster pair time difference vs E, Sector %d",iSM),
- 100,0,10, 200,-100,100);
- fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
- fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
- fOutputContainer->Add(fhClusterPairDiffTimeSameSector[iSM]);
+ else if (!fRecoUtils->IsGoodCluster(c2,fEMCALGeo,fEMCALCells,bc))continue;
+ else if (c2->GetNCells() < fMinNCells) continue;
- }
-
- if(iSM < nSM-2)
- {
- snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
+ else if (c2->GetM02() < fL0min || c2->GetM02() > fL0max) continue;
+
+
+ fRecoUtils->GetMaxEnergyCell(fEMCALGeo, fEMCALCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
+ c2->GetMomentum(p2,fVertex);
+
+ p12 = p1+p2;
+ Float_t invmass = p12.M()*1000;
+
+ //Asimetry cut
+ Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
+
+ if(asym > fAsyCut) continue;
+
+ //Time cut
+ Double_t time2 = c2->GetTOF()*1.e9;
+
+ if(time2 > fTimeMax || time2 < fTimeMin) continue;
+
+ fhClusterPairDiffTime->Fill(p12.E(),time1-time2);
+ if(TMath::Abs(time1-time2) > fDTimeCut) continue;
+
+ if(invmass < fMaxBin && invmass > fMinBin )
+ {
+ //Check if cluster is in fidutial region, not too close to borders
+ Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, fEMCALCells);
+
+ // Clusters not facing frame structures
+ Bool_t mask2 = MaskFrameCluster(iSupMod2, ieta2);
+ //if(mask2) printf("Reject eta %d SM %d\n",ieta2, iSupMod2);
+
+ if(in1 && in2)
+ {
+ fHmgg->Fill(invmass,p12.Pt());
+
+ if(iSupMod1==iSupMod2)
+ {
+ fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
+ fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(p12.E(),time1-time2);
+ }
+ else
+ fHmggDifferentSM ->Fill(invmass,p12.Pt());
+
+ // Same sector
+ Int_t j=0;
+ for(Int_t i = 0; i < nSM/2; i++)
+ {
+ j=2*i;
+ if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j))
+ {
+ fHmggPairSameSectorSM[i]->Fill(invmass,p12.Pt());
+ fhClusterPairDiffTimeSameSector[i]->Fill(p12.E(),time1-time2);
+ }
+ }
+
+ // Same side
+ for(Int_t i = 0; i < nSM-2; i++)
+ {
+ if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i))
+ {
+ fHmggPairSameSideSM[i]->Fill(invmass,p12.Pt());
+ fhClusterPairDiffTimeSameSide[i]->Fill(p12.E(),time1-time2);
+ }
+ }
+
+
+ if(!mask1 && !mask2)
+ {
+ fHmggMaskFrame->Fill(invmass,p12.Pt());
+
+ if(iSupMod1==iSupMod2) fHmggSMMaskFrame[iSupMod1]->Fill(invmass,p12.Pt());
+ else fHmggDifferentSMMaskFrame ->Fill(invmass,p12.Pt());
+
+ // Same sector
+ j=0;
+ for(Int_t i = 0; i < nSM/2; i++)
+ {
+ j=2*i;
+ if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) fHmggPairSameSectorSMMaskFrame[i]->Fill(invmass,p12.Pt());
+ }
+
+ // Same side
+ for(Int_t i = 0; i < nSM-2; i++)
+ {
+ if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) fHmggPairSameSideSMMaskFrame[i]->Fill(invmass,p12.Pt());
+ }
+
+ }// Pair not facing frame
+
+
+ if(invmass > fInvMassCutMin && invmass < fInvMassCutMax) //restrict to clusters really close to pi0 peak
+ {
+
+ // Check time of cells in both clusters, and fill time histogram
+ for(Int_t icell = 0; icell < c1->GetNCells(); icell++)
+ {
+ Int_t absID = c1->GetCellAbsId(icell);
+ fHTpi0[bc%4]->Fill(absID, fEMCALCells->GetCellTime(absID)*1.e9);
+ }
+
+ for(Int_t icell = 0; icell < c2->GetNCells(); icell++)
+ {
+ Int_t absID = c2->GetCellAbsId(icell);
+ fHTpi0[bc%4]->Fill(absID, fEMCALCells->GetCellTime(absID)*1.e9);
+ }
+
+ //Opening angle of 2 photons
+ Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
+ //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
+
+
+ fHOpeningAngle ->Fill(opangle,p12.Pt());
+ fHAsymmetry ->Fill(asym,p12.Pt());
+
+ if(iSupMod1==iSupMod2)
+ {
+ fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
+ fHAsymmetrySM[iSupMod1] ->Fill(asym,p12.Pt());
+ }
+ else
+ {
+ fHOpeningAngleDifferentSM ->Fill(opangle,p12.Pt());
+ fHAsymmetryDifferentSM ->Fill(asym,p12.Pt());
+ }
+
+ if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0))
+ {
+ fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt());
+ fHAsymmetryPairSM[0] ->Fill(asym,p12.Pt());
+
+ }
+ if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1))
+ {
+ fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt());
+ fHAsymmetryPairSM[1] ->Fill(asym,p12.Pt());
+ }
+
+ if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0))
+ {
+ fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt());
+ fHAsymmetryPairSM[2] ->Fill(asym,p12.Pt());
+ }
+ if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2))
+ {
+ fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt());
+ fHAsymmetryPairSM[3] ->Fill(asym,p12.Pt());
+ }
+
+ }// pair in 100 < mass < 160
+
+ }//in acceptance cuts
+
+ //In case of filling only channels with second cluster in same SM
+ if(fSameSM && iSupMod1!=iSupMod2) continue;
+
+ if (fGroupNCells == 0)
+ {
+ fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
+ fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
+
+ if(invmass > fInvMassCutMin && invmass < fInvMassCutMax)//restrict to clusters really close to pi0 peak
+ {
+ fhTowerDecayPhotonHit [iSupMod1]->Fill(ieta1,iphi1);
+ fhTowerDecayPhotonEnergy [iSupMod1]->Fill(ieta1,iphi1,p1.E());
+ fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
+
+ fhTowerDecayPhotonHit [iSupMod2]->Fill(ieta2,iphi2);
+ fhTowerDecayPhotonEnergy [iSupMod2]->Fill(ieta2,iphi2,p2.E());
+ fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
+
+ if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
+ if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
+
+ }// pair in mass of pi0
+ }
+ else {
+ //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
+ for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++)
+ {
+ for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++)
+ {
+ Int_t absId11 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod1, iphi1+j, ieta1+i);
+ Int_t absId22 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod2, iphi2+j, ieta2+i);
+ Bool_t ok1 = kFALSE;
+ Bool_t ok2 = kFALSE;
+ for(Int_t icell = 0; icell < c1->GetNCells(); icell++){
+ if(c1->GetCellsAbsId()[icell] == absId11) ok1=kTRUE;
+ }
+ for(Int_t icell = 0; icell < c2->GetNCells(); icell++){
+ if(c2->GetCellsAbsId()[icell] == absId22) ok2=kTRUE;
+ }
+
+ if(ok1 && (ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24))
+ {
+ fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
+ }
+ if(ok2 && (ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24))
+ {
+ fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
+ }
+ }// j loop
+ }//i loop
+ }//group cells
+
+ if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Mass in (SM%d,%d,%d) and (SM%d,%d,%d): %.3f GeV E1_i=%f E1_ii=%f E2_i=%f E2_ii=%f\n",
+ iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
+ }
+
+ }
+
+ } // end of loop over EMCAL clusters
+}
+
+//________________________________________________________________
+void AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices()
+{
+ // Init geometry and set the geometry matrix, for the first event, skip the rest
+ // Also set once the run dependent calibrations
+
+
+ Int_t runnumber = InputEvent()->GetRunNumber() ;
+
+ if(fLoadMatrices)
+ {
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - Load user defined EMCAL geometry matrices\n");
+
+ // OADB if available
+ AliOADBContainer emcGeoMat("AliEMCALgeo");
+
+ if(fOADBFilePath=="") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;
+
+ emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
+
+ TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
+
+ for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
+ {
+
+ if (!fMatrix[mod]) // Get it from OADB
+ {
+ if(fDebug > 1 )
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - EMCAL matrices SM %d, %p\n",
+ mod,((TGeoHMatrix*) matEMCAL->At(mod)));
+ //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
+
+ fMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
+ }
+
+ if(fMatrix[mod])
+ {
+ if(DebugLevel() > 1)
+ fMatrix[mod]->Print();
+
+ fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
+ }
+
+ }//SM loop
+ }//Load matrices
+ else if(!gGeoManager)
+ {
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - Get geo matrices from data");
+ //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
+ if(!strcmp(InputEvent()->GetName(),"AliAODEvent"))
+ {
+ if(DebugLevel() > 1)
+ Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
+ }//AOD
+ else
+ {
+ if(DebugLevel() > 1)
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
+
+ for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
+ {
+ if(DebugLevel() > 1)
+ InputEvent()->GetEMCALMatrix(mod)->Print();
+
+ if(InputEvent()->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(InputEvent()->GetEMCALMatrix(mod),mod) ;
+
+ }
+
+ }//ESD
+ }//Load matrices from Data
+
+}
+
+//______________________________________________________________________
+void AliAnalysisTaskEMCALPi0CalibSelection::InitTemperatureCorrections()
+{
+ // Apply run dependent calibration correction
+
+ if(!fRecoUtils->IsRunDepRecalibrationOn()) return;
+
+ AliOADBContainer *contRFTD=new AliOADBContainer("");
+
+ contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
+
+ Int_t runnumber = InputEvent()->GetRunNumber() ;
+
+ TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
+
+ //If it did not exist for this run, get closes one
+ if (!htd)
+ {
+ AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
+ // let's get the closest runnumber instead then..
+ Int_t lower = 0;
+ Int_t ic = 0;
+ Int_t maxEntry = contRFTD->GetNumberOfEntries();
+
+ while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
+ lower = ic;
+ ic++;
+ }
+
+ Int_t closest = lower;
+ if ( (ic<maxEntry) &&
+ (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
+ closest = ic;
+ }
+
+ AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
+ htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
+ }
+
+ // Fill parameters
+ if(htd)
+ {
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
+
+ Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
+
+ for (Int_t ism = 0; ism < nSM; ++ism)
+ {
+ for (Int_t icol = 0; icol < 48; ++icol)
+ {
+ for (Int_t irow = 0; irow < 24; ++irow)
+ {
+ Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
+
+ Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
+
+ if(DebugLevel() > 3)
+ printf(" ism %d, icol %d, irow %d,absID %d - Calib factor %1.5f - ",ism, icol, irow, absID, factor);
+
+ factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
+
+ fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
+
+ if(DebugLevel() > 3)
+ printf(" T factor %1.5f - final factor %1.5f \n",htd->GetBinContent(absID) / 10000.,
+ fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow));
+
+ } // columns
+ } // rows
+ } // SM loop
+ }else printf("AliAnalysisTaskEMCALPi0CalibSelection::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n");
+
+}
+
+//___________________________________________________________________
+void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
+{
+ //Create output container, init geometry
+
+ fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
+ Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
+
+ fOutputContainer = new TList();
+ const Int_t buffersize = 255;
+ char hname[buffersize], htitl[buffersize];
+
+ fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
+ fOutputContainer->Add(fhNEvents);
+
+ fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmgg);
+
+ fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmggDifferentSM);
+
+ fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
+ fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
+ fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHOpeningAngle);
+
+ fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
+ fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
+ fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHOpeningAngleDifferentSM);
+
+ fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
+ fHAsymmetry->SetXTitle("a");
+ fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHAsymmetry);
+
+ fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
+ fHAsymmetryDifferentSM->SetXTitle("a");
+ fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHAsymmetryDifferentSM);
+
+
+ //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
+
+ fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmggMaskFrame);
+
+ fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
+ fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmggDifferentSMMaskFrame);
+
+
+ for(Int_t iSM = 0; iSM < nSM; iSM++)
+ {
+ snprintf(hname, buffersize, "hmgg_SM%d",iSM);
+ snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
+ fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmggSM[iSM]);
+
+ snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
+ snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
+ fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmggSMMaskFrame[iSM]);
+
+
+ if(iSM < nSM/2)
+ {
+ snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
+ snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
+ fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmggPairSameSectorSM[iSM]);
+
+ snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
+ snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
+ fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmggPairSameSectorSMMaskFrame[iSM]);
+
+ fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
+ Form("cluster pair time difference vs E, Sector %d",iSM),
+ 100,0,10, 200,-100,100);
+ fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
+ fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
+ fOutputContainer->Add(fhClusterPairDiffTimeSameSector[iSM]);
+
+
+ }
+
+ if(iSM < nSM-2)
+ {
+ snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
fHmggPairSameSideSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
fHmggPairSameSideSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
}
//______________________________________________________________________________________________________
-Bool_t AliAnalysisTaskEMCALPi0CalibSelection::MaskFrameCluster(const Int_t iSM, const Int_t ieta) const
+Bool_t AliAnalysisTaskEMCALPi0CalibSelection::MaskFrameCluster(Int_t iSM, Int_t ieta) const
{
//Check if cell is in one of the regions where we have significant amount of material in front of EMCAL
//__________________________________________________________________________
void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
{
- //Analysis per event.
-
- if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton)
- {
- printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
- abort();
- }
+ // Do analysis, first select the events, then correct the clusters if needed
+ // and finally fill the histograms per channel after recalibration
+ //Event selection
if(fTriggerName!="")
{
AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (InputEvent());
AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (InputEvent());
- TString triggerClass = "";
- if (esdevent) triggerClass = esdevent->GetFiredTriggerClasses();
- else if(aodevent) triggerClass = aodevent->GetFiredTriggerClasses();
-
- if(triggerClass.Contains(fTriggerName))
- {
- //printf("Reject Event %d, FiredClass %s\n",(Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data());
- return;
- }
- }
-
- fhNEvents->Fill(0); //Event analyzed
-
- //Get the input event
- AliVEvent* event = 0;
- if(fFilteredInput) event = AODEvent();
- else event = InputEvent();
-
- if(!event)
- {
- printf("Input event not available!\n");
- return;
- }
-
- if(DebugLevel() > 1)
- printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
-
- //Get the primary vertex
- Double_t v[3];
- event->GetPrimaryVertex()->GetXYZ(v) ;
-
- if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
-
- //Int_t runNum = aod->GetRunNumber();
- //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
-
- InitGeometryMatrices();
-
- Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
-
- if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
- Int_t absId1 = -1;
- Int_t iSupMod1 = -1;
- Int_t iphi1 = -1;
- Int_t ieta1 = -1;
- Int_t absId2 = -1;
- Int_t iSupMod2 = -1;
- Int_t iphi2 = -1;
- Int_t ieta2 = -1;
- Bool_t shared = kFALSE;
-
- TLorentzVector p1;
- TLorentzVector p2;
- TLorentzVector p12;
-
- //Get the list of clusters
- TRefArray * caloClustersArr = new TRefArray();
-
- event->GetEMCALClusters(caloClustersArr);
-
- const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
-
- if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
-
- // Get EMCAL cells
- AliVCaloCells *emCells = event->GetEMCALCells();
-
- // loop over EMCAL clusters
- //----------------------------------------------------------
- // First recalibrate and recalculate energy and position
- Float_t pos[]={0,0,0};
-
- if(fCorrectClusters)
- {
- for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++)
- {
- AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
-
- Float_t e1i = c1->E(); // cluster energy before correction
- if (e1i < fEmin) continue;
- else if (e1i > fEmax) continue;
-
- else if (c1->GetNCells() < fMinNCells) continue;
-
- else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
-
- if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
-
- if(DebugLevel() > 2)
- {
- printf("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
- c1->GetPosition(pos);
- printf("Std : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
- }
-
- //Correct cluster energy and position if requested, and not corrected previously, by default Off
- if(fRecoUtils->IsRecalibrationOn())
- {
- fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
- fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, emCells,c1);
- fRecoUtils->RecalculateClusterPID(c1);
- }
-
- if(DebugLevel() > 2)
- printf("Energy: after recalibration %f; \n",c1->E());
-
- // Recalculate cluster position
- fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
-
- // Correct Non-Linearity
- c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
-
- if(DebugLevel() > 2)
- printf("\t after linearity correction %f\n",c1->E());
-
- //In case of MC analysis, to match resolution/calibration in real data
- c1->SetE(fRecoUtils->SmearClusterEnergy(c1));
-
- if(DebugLevel() > 2)
- printf("\t after smearing %f\n",c1->E());
-
- if(DebugLevel() > 2)
- {
- printf("Cor : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
- c1->GetPosition(pos);
- printf("Cor : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
- }
- }
- }
-
- //----------------------------------------------------------
- //Now the invariant mass analysis with the corrected clusters
- Int_t bc = event->GetBunchCrossNumber();
-
- for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++)
- {
- AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
-
- if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
-
- Float_t e1i = c1->E(); // cluster energy before correction
-
- if (e1i < fEmin) continue;
- else if (e1i > fEmax) continue;
-
- else if (!fRecoUtils->IsGoodCluster(c1,fEMCALGeo,emCells,bc)) continue;
-
- else if (c1->GetNCells() < fMinNCells) continue;
-
- else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
-
- if(DebugLevel() > 2)
- {
- printf("IMA : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
- c1->GetPosition(pos);
- printf("IMA : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
- }
-
- fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
- c1->GetMomentum(p1,v);
-
- //Check if cluster is in fidutial region, not too close to borders
- Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, emCells);
-
- // Clusters not facing frame structures
- Bool_t mask1 = MaskFrameCluster(iSupMod1, ieta1);
- //if(mask1) printf("Reject eta %d SM %d\n",ieta1, iSupMod1);
-
- Double_t time1 = c1->GetTOF()*1.e9;
+ TString triggerClass = "";
+ if (esdevent) triggerClass = esdevent->GetFiredTriggerClasses();
+ else if(aodevent) triggerClass = aodevent->GetFiredTriggerClasses();
- if(time1 > fTimeMax || time1 < fTimeMin) continue;
-
- fhClusterTime ->Fill(c1->E(),time1);
- fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
+ if(DebugLevel() > 1)
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Event %d, FiredClass %s",
+ (Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data());
- // Combine cluster with other clusters and get the invariant mass
- for (Int_t jClu=iClu+1; jClu<kNumberOfEMCALClusters; jClu++)
+ if(!triggerClass.Contains(fTriggerName))
{
- AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
-
- Float_t e2i = c2->E();
- if (e2i < fEmin) continue;
- else if (e2i > fEmax) continue;
-
- else if (!fRecoUtils->IsGoodCluster(c2,fEMCALGeo,emCells,bc))continue;
-
- else if (c2->GetNCells() < fMinNCells) continue;
-
- else if (c2->GetM02() < fL0min || c2->GetM02() > fL0max) continue;
-
-
- fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
- c2->GetMomentum(p2,v);
-
- p12 = p1+p2;
- Float_t invmass = p12.M()*1000;
-
- //Asimetry cut
- Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
-
- if(asym > fAsyCut) continue;
-
- //Time cut
- Double_t time2 = c2->GetTOF()*1.e9;
-
- if(time2 > fTimeMax || time2 < fTimeMin) continue;
-
- fhClusterPairDiffTime->Fill(p12.E(),time1-time2);
- if(TMath::Abs(time1-time2) > fDTimeCut) continue;
-
- if(invmass < fMaxBin && invmass > fMinBin )
- {
- //Check if cluster is in fidutial region, not too close to borders
- Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, emCells);
-
- // Clusters not facing frame structures
- Bool_t mask2 = MaskFrameCluster(iSupMod2, ieta2);
- //if(mask2) printf("Reject eta %d SM %d\n",ieta2, iSupMod2);
-
- if(in1 && in2)
- {
- fHmgg->Fill(invmass,p12.Pt());
-
- if(iSupMod1==iSupMod2)
- {
- fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
- fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(p12.E(),time1-time2);
- }
- else
- fHmggDifferentSM ->Fill(invmass,p12.Pt());
-
- // Same sector
- Int_t j=0;
- for(Int_t i = 0; i < nSM/2; i++)
- {
- j=2*i;
- if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j))
- {
- fHmggPairSameSectorSM[i]->Fill(invmass,p12.Pt());
- fhClusterPairDiffTimeSameSector[i]->Fill(p12.E(),time1-time2);
- }
- }
-
- // Same side
- for(Int_t i = 0; i < nSM-2; i++)
- {
- if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i))
- {
- fHmggPairSameSideSM[i]->Fill(invmass,p12.Pt());
- fhClusterPairDiffTimeSameSide[i]->Fill(p12.E(),time1-time2);
- }
- }
-
-
- if(!mask1 && !mask2)
- {
- fHmggMaskFrame->Fill(invmass,p12.Pt());
-
- if(iSupMod1==iSupMod2) fHmggSMMaskFrame[iSupMod1]->Fill(invmass,p12.Pt());
- else fHmggDifferentSMMaskFrame ->Fill(invmass,p12.Pt());
-
- // Same sector
- j=0;
- for(Int_t i = 0; i < nSM/2; i++)
- {
- j=2*i;
- if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) fHmggPairSameSectorSMMaskFrame[i]->Fill(invmass,p12.Pt());
- }
-
- // Same side
- for(Int_t i = 0; i < nSM-2; i++)
- {
- if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) fHmggPairSameSideSMMaskFrame[i]->Fill(invmass,p12.Pt());
- }
-
- }// Pair not facing frame
-
-
- if(invmass > fInvMassCutMin && invmass < fInvMassCutMax) //restrict to clusters really close to pi0 peak
- {
-
- // Check time of cells in both clusters, and fill time histogram
- for(Int_t icell = 0; icell < c1->GetNCells(); icell++)
- {
- Int_t absID = c1->GetCellAbsId(icell);
- fHTpi0[bc%4]->Fill(absID, emCells->GetCellTime(absID)*1.e9);
- }
-
- for(Int_t icell = 0; icell < c2->GetNCells(); icell++)
- {
- Int_t absID = c2->GetCellAbsId(icell);
- fHTpi0[bc%4]->Fill(absID, emCells->GetCellTime(absID)*1.e9);
- }
-
- //Opening angle of 2 photons
- Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
- //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
-
-
- fHOpeningAngle ->Fill(opangle,p12.Pt());
- fHAsymmetry ->Fill(asym,p12.Pt());
-
- if(iSupMod1==iSupMod2)
- {
- fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
- fHAsymmetrySM[iSupMod1] ->Fill(asym,p12.Pt());
- }
- else
- {
- fHOpeningAngleDifferentSM ->Fill(opangle,p12.Pt());
- fHAsymmetryDifferentSM ->Fill(asym,p12.Pt());
- }
-
- if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0))
- {
- fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt());
- fHAsymmetryPairSM[0] ->Fill(asym,p12.Pt());
-
- }
- if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1))
- {
- fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt());
- fHAsymmetryPairSM[1] ->Fill(asym,p12.Pt());
- }
-
- if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0))
- {
- fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt());
- fHAsymmetryPairSM[2] ->Fill(asym,p12.Pt());
- }
- if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2))
- {
- fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt());
- fHAsymmetryPairSM[3] ->Fill(asym,p12.Pt());
- }
-
- }// pair in 100 < mass < 160
-
- }//in acceptance cuts
-
- //In case of filling only channels with second cluster in same SM
- if(fSameSM && iSupMod1!=iSupMod2) continue;
-
- if (fGroupNCells == 0)
- {
- fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
- fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
-
- if(invmass > fInvMassCutMin && invmass < fInvMassCutMax)//restrict to clusters really close to pi0 peak
- {
- fhTowerDecayPhotonHit [iSupMod1]->Fill(ieta1,iphi1);
- fhTowerDecayPhotonEnergy [iSupMod1]->Fill(ieta1,iphi1,p1.E());
- fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
-
- fhTowerDecayPhotonHit [iSupMod2]->Fill(ieta2,iphi2);
- fhTowerDecayPhotonEnergy [iSupMod2]->Fill(ieta2,iphi2,p2.E());
- fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
-
- if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
- if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
-
- }// pair in mass of pi0
- }
- else {
- //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
- for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++)
- {
- for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++)
- {
- Int_t absId11 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod1, iphi1+j, ieta1+i);
- Int_t absId22 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod2, iphi2+j, ieta2+i);
- Bool_t ok1 = kFALSE;
- Bool_t ok2 = kFALSE;
- for(Int_t icell = 0; icell < c1->GetNCells(); icell++){
- if(c1->GetCellsAbsId()[icell] == absId11) ok1=kTRUE;
- }
- for(Int_t icell = 0; icell < c2->GetNCells(); icell++){
- if(c2->GetCellsAbsId()[icell] == absId22) ok2=kTRUE;
- }
-
- if(ok1 && (ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24))
- {
- fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
- }
- if(ok2 && (ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24))
- {
- fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
- }
- }// j loop
- }//i loop
- }//group cells
-
- if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Mass in (SM%d,%d,%d) and (SM%d,%d,%d): %.3f GeV E1_i=%f E1_ii=%f E2_i=%f E2_ii=%f\n",
- iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
- }
-
- }
-
- } // end of loop over EMCAL clusters
+ if(DebugLevel() > 1) printf("Reject event! \n");
+ return;
+ }
+ else
+ if(DebugLevel() > 1) printf("Accept Event! \n");
+ }
+
+ //Get the input event
+ AliVEvent* event = 0;
+ if(fFilteredInput) event = AODEvent();
+ else event = InputEvent();
+
+ if(!event)
+ {
+ printf("Input event not available!\n");
+ return;
+ }
+
+ if(DebugLevel() > 1)
+ printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
+
+ //Get the primary vertex
+ event->GetPrimaryVertex()->GetXYZ(fVertex) ;
+
+ if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",fVertex[0],fVertex[1],fVertex[2]);
+
+ //Int_t runNum = aod->GetRunNumber();
+ //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
+
+ fhNEvents->Fill(0); //Count the events to be analyzed
+
+ // Acccess once the geometry matrix and temperature corrections
+ if(fhNEvents->GetEntries()==1)
+ {
+ InitGeometryMatrices();
+
+ InitTemperatureCorrections();
+ }
+
+ //Get the list of clusters and cells
+ fEMCALCells = event->GetEMCALCells();
+
+ fCaloClustersArr = new TRefArray();
+ event->GetEMCALClusters(fCaloClustersArr);
- delete caloClustersArr;
+ if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d - N CaloCells %d \n",
+ fCaloClustersArr->GetEntriesFast(), fEMCALCells->GetNumberOfCells());
+
+ CorrectClusters(); // Non linearity, new calibration, T calibration
+
+ FillHistograms();
+
+ delete fCaloClustersArr;
PostData(1,fOutputContainer);
printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
if(fLoadMatrices) {for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) if(fMatrix[ism]) fMatrix[ism]->Print() ; }
-
}
//____________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff