10 #include "TParticle.h"
15 #include "AliAODMCParticle.h"
16 #include "AliAnalysisManager.h"
17 #include "AliMCEventHandler.h"
18 #include "AliMCEvent.h"
20 #include "AliAnalysisTaskSE.h"
21 #include "AliAnalysisTaskPi0Efficiency.h"
22 #include "AliCaloPhoton.h"
23 #include "AliPHOSGeometry.h"
24 #include "AliPHOSEsdCluster.h"
25 #include "AliPHOSCalibData.h"
26 #include "AliAODEvent.h"
27 #include "AliAODCaloCluster.h"
28 #include "AliAODVertex.h"
29 #include "AliESDtrackCuts.h"
32 #include "AliCDBManager.h"
33 #include "AliCentrality.h"
35 // Analysis task to fill histograms with PHOS ESD clusters and cells
36 // Authors: Yuri Kharlov
39 ClassImp(AliAnalysisTaskPi0Efficiency)
41 //________________________________________________________________________
42 AliAnalysisTaskPi0Efficiency::AliAnalysisTaskPi0Efficiency(const char *name)
43 : AliAnalysisTaskSE(name),
61 for(Int_t i=0;i<1;i++){
62 for(Int_t j=0;j<10;j++)
63 for(Int_t k=0;k<11;k++)
64 fPHOSEvents[i][j][k]=0 ;
67 // Output slots #0 write into a TH1 container
68 DefineOutput(1,TList::Class());
70 // Set bad channel map
72 for(Int_t i=0; i<6; i++){
73 snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
74 fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
76 // Initialize the PHOS geometry
77 fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
82 //________________________________________________________________________
83 void AliAnalysisTaskPi0Efficiency::UserCreateOutputObjects()
89 if(fOutputContainer != NULL){
90 delete fOutputContainer;
92 fOutputContainer = new TList();
93 fOutputContainer->SetOwner(kTRUE);
96 fOutputContainer->Add(new TH1F("hSelEvents","Event celection", 10,0.,10.)) ;
99 fOutputContainer->Add(new TH1F("hZvertex","Z vertex position", 50,-25.,25.)) ;
102 fOutputContainer->Add(new TH1F("hCentrality","Event centrality", 100,0.,100.)) ;
105 fOutputContainer->Add(new TH2F("hCluM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
106 fOutputContainer->Add(new TH2F("hCluM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
107 fOutputContainer->Add(new TH2F("hCluM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
117 for(Int_t cent=0; cent<6; cent++){
119 snprintf(key,55,"hPhotAll_cen%d",cent) ;
120 fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
121 snprintf(key,55,"hPhotCPV_cen%d",cent) ;
122 fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
123 snprintf(key,55,"hPhotDisp_cen%d",cent) ;
124 fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
125 snprintf(key,55,"hPhotBoth_cen%d",cent) ;
126 fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
128 snprintf(key,55,"hMassPtAll_cen%d",cent) ;
129 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
130 snprintf(key,55,"hMassPtCPV_cen%d",cent) ;
131 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
132 snprintf(key,55,"hMassPtDisp_cen%d",cent) ;
133 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
134 snprintf(key,55,"hMassPtBoth_cen%d",cent) ;
135 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
138 snprintf(key,55,"hMiMassPtAll_cen%d",cent) ;
139 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
140 snprintf(key,55,"hMiMassPtCPV_cen%d",cent) ;
141 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
142 snprintf(key,55,"hMiMassPtDisp_cen%d",cent) ;
143 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
144 snprintf(key,55,"hMiMassPtBoth_cen%d",cent) ;
145 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
147 snprintf(key,55,"hMCMassPtAll_cen%d",cent) ;
148 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
149 snprintf(key,55,"hMCMassPtCPV_cen%d",cent) ;
150 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
151 snprintf(key,55,"hMCMassPtDisp_cen%d",cent) ;
152 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
153 snprintf(key,55,"hMCMassPtBoth_cen%d",cent) ;
154 fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
157 snprintf(key,55,"hMCPhotAll_cen%d",cent) ;
158 fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
159 snprintf(key,55,"hMCPhotCPV_cen%d",cent) ;
160 fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
161 snprintf(key,55,"hMCPhotDisp_cen%d",cent) ;
162 fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
163 snprintf(key,55,"hMCPhotBoth_cen%d",cent) ;
164 fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
172 for(Int_t cent=0; cent<6; cent++){
173 snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
174 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
175 snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
176 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
177 snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
178 fOutputContainer->Add(new TH1F("hMC_rap_eta","Rapidity eta",200,-1.,1.)) ;
179 snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
180 fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
181 snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
182 fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
183 snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
184 fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
185 snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
186 fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
187 snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
188 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
189 snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
190 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
191 snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
192 fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
193 snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
194 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
195 snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
196 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
199 PostData(1, fOutputContainer);
203 //________________________________________________________________________
204 void AliAnalysisTaskPi0Efficiency::UserExec(Option_t *)
206 // Main loop, called for each event
209 FillHistogram("hSelEvents",0.5) ;
211 AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
213 Printf("ERROR: Could not retrieve event");
214 PostData(1, fOutputContainer);
218 FillHistogram("hSelEvents",1.5) ;
219 AliAODHeader *header = event->GetHeader() ;
221 // Checks if we have a primary vertex
222 // Get primary vertices form ESD
223 const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
225 // don't rely on ESD vertex, assume (0,0,0)
226 Double_t vtx0[3] ={0.,0.,0.};
227 Double_t vtx5[3] ={0.,0.,0.};
229 vtx5[0] = esdVertex5->GetX();
230 vtx5[1] = esdVertex5->GetY();
231 vtx5[2] = esdVertex5->GetZ();
234 FillHistogram("hZvertex",esdVertex5->GetZ());
235 if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
236 PostData(1, fOutputContainer);
239 FillHistogram("hSelEvents",2.5) ;
242 // Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
243 // if(zvtx<0)zvtx=0 ;
244 // if(zvtx>9)zvtx=9 ;
247 // fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile,
248 // //a float from 0 to 100 (or to the trigger efficiency)
249 fCentrality=header->GetZDCN2Energy() ;
251 if( fCentrality < 0. || fCentrality>80.){
252 PostData(1, fOutputContainer);
255 FillHistogram("hSelEvents",3.5) ;
256 Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
258 while(fCenBin<6 && fCentrality > bins[fCenBin+1])
263 fRPfull= header->GetZDCN1Energy() ;
264 if(fRPfull==999){ //reaction plain was not defined
265 PostData(1, fOutputContainer);
269 FillHistogram("hSelEvents",4.5) ;
270 //All event selections done
271 FillHistogram("hCentrality",fCentrality) ;
272 //Reaction plain is defined in the range (-pi/2;pi/2)
274 Int_t irp=Int_t(10.*(fRPfull+TMath::PiOver2())/TMath::Pi());
277 if(!fPHOSEvents[zvtx][fCenBin][irp])
278 fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
279 TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
281 // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
282 if(fEventCounter == 0) {
283 for(Int_t mod=0; mod<5; mod++) {
284 const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
285 fPHOSGeo->SetMisalMatrix(m,mod) ;
286 Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
294 fPHOSEvent->Clear() ;
296 fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
301 TVector3 vertex(vtx0);
302 TClonesArray * clusters = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
303 Int_t multClust = clusters->GetEntriesFast();
304 for (Int_t i=0; i<multClust; i++) {
305 AliVCluster *clu = (AliVCluster*) clusters->At(i);
306 if ( !clu->IsPHOS() || clu->E()<0.3) continue;
309 clu->GetPosition(position);
310 TVector3 global(position) ;
312 fPHOSGeo->GlobalPos2RelId(global,relId) ;
313 Int_t mod = relId[0] ;
314 Int_t cellX = relId[2];
315 Int_t cellZ = relId[3] ;
316 if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
318 if(clu->GetNCells()<3)
321 snprintf(key,55,"hCluM%d",mod) ;
322 FillHistogram(key,cellX,cellZ,1.);
325 clu->GetMomentum(pv1 ,vtx0);
327 if(inPHOS>=fPHOSEvent->GetSize()){
328 fPHOSEvent->Expand(inPHOS+50) ;
330 new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
331 AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
334 ph->SetNCells(clu->GetNCells());
335 ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
336 ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
341 for (Int_t i1=0; i1<inPHOS; i1++) {
342 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
343 snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
344 FillHistogram(key,ph1->Pt()) ;
346 snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
347 FillHistogram(key,ph1->Pt()) ;
350 snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
351 FillHistogram(key,ph1->Pt()) ;
353 snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
354 FillHistogram(key,ph1->Pt()) ;
359 // Fill Real disribution
360 for (Int_t i1=0; i1<inPHOS-1; i1++) {
361 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
362 for (Int_t i2=i1+1; i2<inPHOS; i2++) {
363 AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
364 TLorentzVector p12 = *ph1 + *ph2;
365 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
367 snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
368 FillHistogram(key,p12.M() ,p12.Pt()) ;
369 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
370 snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
371 FillHistogram(key,p12.M() ,p12.Pt()) ;
373 if(ph1->IsDispOK() && ph2->IsDispOK()){
374 snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
375 FillHistogram(key,p12.M() ,p12.Pt()) ;
377 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
378 snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
379 FillHistogram(key,p12.M() ,p12.Pt()) ;
386 for (Int_t i1=0; i1<inPHOS; i1++) {
387 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
388 for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
389 TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
390 for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
391 AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
392 TLorentzVector p12 = *ph1 + *ph2;
393 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
395 snprintf(key,55,"hMiMassPtAll_cen%d",fCenBin) ;
396 FillHistogram(key,p12.M() ,p12.Pt()) ;
397 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
398 snprintf(key,55,"hMiMassPtCPV_cen%d",fCenBin) ;
399 FillHistogram(key,p12.M() ,p12.Pt()) ;
401 if(ph1->IsDispOK() && ph2->IsDispOK()){
402 snprintf(key,55,"hMiMassPtDisp_cen%d",fCenBin) ;
403 FillHistogram(key,p12.M() ,p12.Pt()) ;
405 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
406 snprintf(key,55,"hMiMassPtBoth_cen%d",fCenBin) ;
407 FillHistogram(key,p12.M() ,p12.Pt()) ;
415 //Now we either add current events to stack or remove
416 //If no photons in current event - no need to add it to mixed
417 if(fPHOSEvent->GetEntriesFast()>0){
418 prevPHOS->AddFirst(fPHOSEvent) ;
420 if(prevPHOS->GetSize()>100){//Remove redundant events
421 TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
422 prevPHOS->RemoveLast() ;
427 PostData(1, fOutputContainer);
431 //________________________________________________________________________
432 void AliAnalysisTaskPi0Efficiency::Terminate(Option_t *)
434 // Draw result to the screen
435 // Called once at the end of the query
439 //________________________________________________________________________
440 Bool_t AliAnalysisTaskPi0Efficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
442 //Check if this channel belogs to the good ones
444 if(strcmp(det,"PHOS")==0){
446 AliError(Form("No bad map for PHOS module %d ",mod)) ;
449 if(!fPHOSBadMap[mod]){
450 AliError(Form("No Bad map for PHOS module %d",mod)) ;
453 if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
459 AliError(Form("Can not find bad channels for detector %s ",det)) ;
463 //_____________________________________________________________________________
464 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x)const{
466 TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
471 TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
476 TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
481 AliInfo(Form("can not find histogram <%s> ",key)) ;
483 //_____________________________________________________________________________
484 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
486 TObject * tmp = fOutputContainer->FindObject(key) ;
488 AliInfo(Form("can not find histogram <%s> ",key)) ;
491 if(tmp->IsA() == TClass::GetClass("TH1F")){
492 ((TH1F*)tmp)->Fill(x,y) ;
495 if(tmp->IsA() == TClass::GetClass("TH2F")){
496 ((TH2F*)tmp)->Fill(x,y) ;
499 AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
502 //_____________________________________________________________________________
503 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
504 //Fills 1D histograms with key
505 TObject * tmp = fOutputContainer->FindObject(key) ;
507 AliInfo(Form("can not find histogram <%s> ",key)) ;
510 if(tmp->IsA() == TClass::GetClass("TH2F")){
511 ((TH2F*)tmp)->Fill(x,y,z) ;
514 if(tmp->IsA() == TClass::GetClass("TH3F")){
515 ((TH3F*)tmp)->Fill(x,y,z) ;
519 //_____________________________________________________________________________
520 Bool_t AliAnalysisTaskPi0Efficiency::TestLambda(Double_t l1,Double_t l2){
521 Double_t l1Mean=1.22 ;
522 Double_t l2Mean=2.0 ;
523 Double_t l1Sigma=0.42 ;
524 Double_t l2Sigma=0.71 ;
526 Double_t R2=(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma+(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma-c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
530 //___________________________________________________________________________
531 void AliAnalysisTaskPi0Efficiency::ProcessMC(){
532 //fill histograms for efficiensy etc. calculation
533 const Double_t rcut = 1. ; //cut for primary particles
534 //---------First pi0/eta-----------------------------
538 AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
540 TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
541 for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
542 AliAODMCParticle* particle = (AliAODMCParticle*) mcArray->At(i);
543 if(particle->GetPdgCode() == 111)
544 snprintf(partName,10,"pi0") ;
546 if(particle->GetPdgCode() == 221)
547 snprintf(partName,10,"eta") ;
549 if(particle->GetPdgCode() == 22)
550 snprintf(partName,10,"gamma") ;
555 Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
559 Double_t pt = particle->Pt() ;
560 //Total number of pi0 with creation radius <1 cm
561 snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCenBin) ;
562 FillHistogram(hkey,pt) ;
563 if(TMath::Abs(particle->Y())<0.12){
564 snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
565 FillHistogram(hkey,pt) ;
568 snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCenBin) ;
569 FillHistogram(hkey,particle->Y()) ;
571 Double_t phi=particle->Phi() ;
572 while(phi<0.)phi+=TMath::TwoPi() ;
573 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
574 snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCenBin) ;
575 FillHistogram(hkey,phi) ;
578 //Check if one of photons converted
579 if(particle->GetNDaughters()!=2)
580 continue ; //Do not account Dalitz decays
583 TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
584 TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
585 //Number of pi0s decayed into acceptance
588 Bool_t hitPHOS1 = fPHOSGeo->ImpactOnEmc(gamma1, mod1, z,x) ;
589 Bool_t hitPHOS2 = fPHOSGeo->ImpactOnEmc(gamma2, mod2, z,x) ;
591 Bool_t goodPair=kFALSE ;
592 if( hitPHOS1 && hitPHOS2){
593 sprintf(hkey,"hMC_PHOSacc_%s",partName) ;
594 FillHistogram(hkey,pt) ;
601 //Now calculate "Real" distribution of clusters with primary
602 TClonesArray cluPrim("AliCaloPhoton",200) ; //clusters with primary
603 Int_t multClust = event->GetNumberOfCaloClusters();
605 Double_t vtx0[3] = {0,0,0};
606 for (Int_t i=0; i<multClust; i++) {
607 AliVCluster *clu = event->GetCaloCluster(i);
608 if ( !clu->IsPHOS() || clu->E()<0.3) continue;
609 if(clu->GetLabel()<0) continue ;
612 clu->GetPosition(position);
613 TVector3 global(position) ;
615 fPHOSGeo->GlobalPos2RelId(global,relId) ;
616 Int_t mod = relId[0] ;
617 Int_t cellX = relId[2];
618 Int_t cellZ = relId[3] ;
619 if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
621 if(clu->GetNCells()<3)
625 clu->GetMomentum(pv1 ,vtx0);
627 if(inPHOS>=cluPrim.GetSize()){
628 cluPrim.Expand(inPHOS+50) ;
630 AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
631 //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
634 ph->SetNCells(clu->GetNCells());
635 ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
636 ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
644 for (Int_t i1=0; i1<inPHOS; i1++) {
645 AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
646 snprintf(key,55,"hMCPhotAll_cen%d",fCenBin) ;
647 FillHistogram(key,ph1->Pt()) ;
649 snprintf(key,55,"hMCPhotCPV_cen%d",fCenBin) ;
650 FillHistogram(key,ph1->Pt()) ;
654 snprintf(key,55,"hMCPhotDisp_cen%d",fCenBin) ;
655 FillHistogram(key,ph1->Pt()) ;
657 snprintf(key,55,"hMCPhotBoth_cen%d",fCenBin) ;
658 FillHistogram(key,ph1->Pt()) ;
663 // Fill Real disribution
664 for (Int_t i1=0; i1<inPHOS-1; i1++) {
665 AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
666 for (Int_t i2=i1+1; i2<inPHOS; i2++) {
667 AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
668 TLorentzVector p12 = *ph1 + *ph2;
669 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
671 snprintf(key,55,"hMCMassPtAll_cen%d",fCenBin) ;
672 FillHistogram(key,p12.M() ,p12.Pt()) ;
673 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
674 snprintf(key,55,"hMCMassPtCPV_cen%d",fCenBin) ;
675 FillHistogram(key,p12.M() ,p12.Pt()) ;
677 if(ph1->IsDispOK() && ph2->IsDispOK()){
678 snprintf(key,55,"hMCMassPtDisp_cen%d",fCenBin) ;
679 FillHistogram(key,p12.M() ,p12.Pt()) ;
681 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
682 snprintf(key,55,"hMCMassPtBoth_cen%d",fCenBin) ;
683 FillHistogram(key,p12.M() ,p12.Pt()) ;