Int_t label=label1[0];
if(label < 0) return -1;
- while(label > -1 && counter1 < 99){
+ while(label > -1 && counter1 < 99)
+ {
counter1++;
AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
- if(mom){
+ if(mom)
+ {
label = mom->GetMother() ;
label1[counter1]=label;
}
//printf("\t counter %d, label %d\n", counter1,label);
}
+
//printf("Org label2=%d,\n",label2[0]);
label=label2[0];
if(label < 0) return -1;
- while(label > -1 && counter2 < 99){
+
+ while(label > -1 && counter2 < 99)
+ {
counter2++;
AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
- if(mom){
+ if(mom)
+ {
label = mom->GetMother() ;
label2[counter2]=label;
}
//printf("\t counter %d, label %d\n", counter2,label);
}
}//AOD MC
- else { //Kine stack from ESDs
+ else
+ { //Kine stack from ESDs
AliStack * stack = reader->GetStack();
+
Int_t label=label1[0];
- while(label > -1 && counter1 < 99){
+ while(label > -1 && counter1 < 99)
+ {
counter1++;
TParticle * mom = stack->Particle(label);
- if(mom){
+ if(mom)
+ {
label = mom->GetFirstMother() ;
label1[counter1]=label;
}
//printf("\t counter %d, label %d\n", counter1,label);
}
+
//printf("Org label2=%d,\n",label2[0]);
+
label=label2[0];
- while(label > -1 && counter2 < 99){
+ while(label > -1 && counter2 < 99)
+ {
counter2++;
TParticle * mom = stack->Particle(label);
- if(mom){
+ if(mom)
+ {
label = mom->GetFirstMother() ;
label2[counter2]=label;
}
}// Kine stack from ESDs
}//First labels not the same
- if((counter1==99 || counter2==99) && fDebug >=0) printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2);
+ if((counter1==99 || counter2==99) && fDebug >=0)
+ printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2);
//printf("CheckAncestor:\n");
+
Int_t commonparents = 0;
Int_t ancLabel = -1;
//printf("counters %d %d \n",counter1, counter2);
- for (Int_t c1 = 0; c1 < counter1; c1++) {
- for (Int_t c2 = 0; c2 < counter2; c2++) {
- if(label1[c1]==label2[c2] && label1[c1]>-1) {
+ for (Int_t c1 = 0; c1 < counter1; c1++)
+ {
+ for (Int_t c2 = 0; c2 < counter2; c2++)
+ {
+ if(label1[c1]==label2[c2] && label1[c1]>-1)
+ {
ancLabel = label1[c1];
commonparents++;
- if(reader->ReadAODMCParticles()){
+
+ if(reader->ReadAODMCParticles())
+ {
AliAODMCParticle * mom = (AliAODMCParticle *) reader->GetAODMCParticles()->At(label1[c1]);
- if (mom) {
+
+ if (mom)
+ {
ancPDG = mom->GetPdgCode();
ancStatus = mom->GetStatus();
momentum.SetPxPyPzE(mom->Px(),mom->Py(),mom->Pz(),mom->E());
prodVertex.SetXYZ(mom->Xv(),mom->Yv(),mom->Zv());
}
}
- else {
+ else
+ {
TParticle * mom = (reader->GetStack())->Particle(label1[c1]);
- if (mom) {
+
+ if (mom)
+ {
ancPDG = mom->GetPdgCode();
ancStatus = mom->GetStatusCode();
mom->Momentum(momentum);
prodVertex.SetXYZ(mom->Vx(),mom->Vy(),mom->Vz());
}
}
+
//First ancestor found, end the loops
counter1=0;
counter2=0;
}//second cluster loop
}//first cluster loop
+ if(ancLabel < 0)
+ {
+ ancPDG = -10000;
+ ancStatus = -10000;
+ momentum.SetXYZT(0,0,0,0);
+ prodVertex.SetXYZ(-10,-10,-10);
+ }
+
return ancLabel;
}
//about its heritage, but one can also use it directly with stack
//particles not connected to reconstructed entities
- if(!stack) {
+ if(!stack)
+ {
if (fDebug >=0)
printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
return -1;
Int_t tag = 0;
Int_t label=labels[0];//Most significant particle contributing to the cluster
- if(label >= 0 && label < stack->GetNtrack()){
+ if(label >= 0 && label < stack->GetNtrack())
+ {
//MC particle of interest is the "mom" of the entity
TParticle * mom = stack->Particle(label);
Int_t iMom = label;
Int_t mPdg = TMath::Abs(mPdgSign);
Int_t mStatus = mom->GetStatusCode() ;
Int_t iParent = mom->GetFirstMother() ;
+
if(fDebug > 0 && label < 8 && fMCGenerator!="") printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
//GrandParent of the entity
TParticle * parent = NULL;
Int_t pPdg = -1;
Int_t pStatus =-1;
- if(iParent >= 0){
+ if(iParent >= 0)
+ {
parent = stack->Particle(iParent);
- if(parent){
+
+ if(parent)
+ {
pPdg = TMath::Abs(parent->GetPdgCode());
pStatus = parent->GetStatusCode();
}
}
else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
- if(fDebug > 2 ) {
+ if(fDebug > 2 )
+ {
printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n");
printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
}
//Check if "mother" of entity is converted, if not, get the first non converted mother
- if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0){
+ if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0)
+ {
SetTagBit(tag,kMCConversion);
+
//Check if the mother is photon or electron with status not stable
- while ((pPdg == 22 || pPdg == 11) && mStatus != 1) {
+ while ((pPdg == 22 || pPdg == 11) && mStatus != 1)
+ {
//Mother
iMom = mom->GetFirstMother();
mom = stack->Particle(iMom);
mPdg = TMath::Abs(mPdgSign);
mStatus = mom->GetStatusCode() ;
iParent = mom->GetFirstMother() ;
+
if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
//GrandParent
- if(iParent >= 0){
+ if(iParent >= 0)
+ {
parent = stack->Particle(iParent);
- if(parent){
+ if(parent)
+ {
pPdg = TMath::Abs(parent->GetPdgCode());
pStatus = parent->GetStatusCode();
}
pStatus = 0;
break;
}
- }//while
- if(fDebug > 2 ) {
+ }//while
+
+ if(fDebug > 2 )
+ {
printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n");
printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
}
}//mother and parent are electron or photon and have status 0
- else if((mPdg == 22 || mPdg == 11) && mStatus == 0){
+ else if((mPdg == 22 || mPdg == 11) && mStatus == 0)
+ {
//Still a conversion but only one electron/photon generated. Just from hadrons but not decays.
if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
- pPdg == 2212 || pPdg == 130 || pPdg == 13 ) {
+ pPdg == 2212 || pPdg == 130 || pPdg == 13 )
+ {
SetTagBit(tag,kMCConversion);
iMom = mom->GetFirstMother();
mom = stack->Particle(iMom);
mPdgSign = mom->GetPdgCode();
mPdg = TMath::Abs(mPdgSign);
- if(fDebug > 2 ) {
+ if(fDebug > 2 )
+ {
printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n");
printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
}
else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
//check for pi0 and eta (shouldn't happen unless their decays were turned off)
- else if(mPdg == 111) {
+ else if(mPdg == 111)
+ {
SetTagBit(tag,kMCPi0Decay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
}
- else if(mPdg == 221) {
+ else if(mPdg == 221)
+ {
SetTagBit(tag,kMCEtaDecay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster
}
//Photons
- else if(mPdg == 22){
+ else if(mPdg == 22)
+ {
SetTagBit(tag,kMCPhoton);
- if(mStatus == 1){ //undecayed particle
- if(fMCGenerator == "PYTHIA"){
- if(iParent < 8 && iParent > 5) {//outgoing partons
+ if(mStatus == 1)
+ { //undecayed particle
+ if(fMCGenerator == "PYTHIA")
+ {
+ if(iParent < 8 && iParent > 5)
+ {//outgoing partons
if(pPdg == 22) SetTagBit(tag,kMCPrompt);
- else SetTagBit(tag,kMCFragmentation);
+ else SetTagBit(tag,kMCFragmentation);
}//Outgoing partons
- else if(iParent <= 5) {
+ else if(iParent <= 5)
+ {
SetTagBit(tag, kMCISR); //Initial state radiation
}
- else if(pStatus == 11){//Decay
- if(pPdg == 111) {
+ else if(pStatus == 11)
+ {//Decay
+ if(pPdg == 111)
+ {
SetTagBit(tag,kMCPi0Decay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon, parent pi0 with status 11 \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
}
- else if (pPdg == 221) {
+ else if (pPdg == 221)
+ {
SetTagBit(tag, kMCEtaDecay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon, parent pi0 with status 11 \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
}
else SetTagBit(tag,kMCOtherDecay);
}//Decay
- else {
+ else
+ {
if(fDebug > 1 && parent) printf("AliMCAnalysisUtils::CheckOrigingInStack() - what is it in PYTHIA? Wrong generator setting? Mother mPdg %d, status %d \n Parent iParent %d, pPdg %d %s, status %d\n",
mPdg, mStatus,iParent, pPdg, parent->GetName(),pStatus);
- if(pPdg == 111) {
+
+ if(pPdg == 111)
+ {
SetTagBit(tag,kMCPi0Decay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon, parent pi0 with status 11 \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
}
- else if (pPdg == 221) {
+ else if (pPdg == 221)
+ {
SetTagBit(tag, kMCEtaDecay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon, parent pi0 with status 11 \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
}
else SetTagBit(tag,kMCOtherDecay);
}
}//PYTHIA
- else if(fMCGenerator == "HERWIG"){
- if(pStatus < 197){//Not decay
- while(1){
- if(parent){
+ else if(fMCGenerator == "HERWIG")
+ {
+ if(pStatus < 197)
+ {//Not decay
+ while(1)
+ {
+ if(parent)
+ {
if(parent->GetFirstMother()<=5) break;
iParent = parent->GetFirstMother();
parent=stack->Particle(iParent);
} else break;
}//Look for the parton
- if(iParent < 8 && iParent > 5) {
+ if(iParent < 8 && iParent > 5)
+ {
if(pPdg == 22) SetTagBit(tag,kMCPrompt);
- else SetTagBit(tag,kMCFragmentation);
+ else SetTagBit(tag,kMCFragmentation);
}
else SetTagBit(tag,kMCISR);//Initial state radiation
}//Not decay
else{//Decay
- if(pPdg == 111) {
+ if(pPdg == 111)
+ {
SetTagBit(tag,kMCPi0Decay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG pi0 decay photon \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
}
- else if (pPdg == 221) {
+ else if (pPdg == 221)
+ {
SetTagBit(tag,kMCEtaDecay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG eta decay photon \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
}
else SetTagBit(tag,kMCOtherDecay);
}//Status 1 : created by event generator
- else if(mStatus == 0){ // geant
- if(pPdg == 111) {
+ else if(mStatus == 0)
+ { // geant
+ if(pPdg == 111)
+ {
SetTagBit(tag,kMCPi0Decay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC pi0 decay photon \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
}
- else if (pPdg == 221) {
+ else if (pPdg == 221)
+ {
SetTagBit(tag,kMCEtaDecay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC eta decay photon \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
}
else SetTagBit(tag,kMCOtherDecay);
//Electron check. Where did that electron come from?
else if(mPdg == 11){ //electron
- if(pPdg == 11 && parent){
+ if(pPdg == 11 && parent)
+ {
Int_t iGrandma = parent->GetFirstMother();
- if(iGrandma >= 0) {
+ if(iGrandma >= 0)
+ {
TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
Int_t gPdg = TMath::Abs(gma->GetPdgCode());
- if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
+ if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
}
}
- SetTagBit(tag,kMCElectron);
+
+ SetTagBit(tag,kMCElectron);
+
if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
+
if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
- if(parent){
+ if(parent)
+ {
Int_t iGrandma = parent->GetFirstMother();
- if(iGrandma >= 0) {
+ if(iGrandma >= 0)
+ {
TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
Int_t gPdg = TMath::Abs(gma->GetPdgCode());
if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
else SetTagBit(tag,kMCEFromC); //c-->e
- } else SetTagBit(tag,kMCEFromC); //c-->e
+ }
+ else SetTagBit(tag,kMCEFromC); //c-->e
}//parent
- } else {
+ }
+ else
+ {
//if it is not from any of the above, where is it from?
if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
+
else SetTagBit(tag,kMCOtherDecay);
+
if(fDebug > 0 && parent) printf("AliMCAnalysisUtils::CheckOriginInStack() - Status %d Electron from other origin: %s (pPdg = %d) %s (mpdg = %d)\n",mStatus,parent->GetName(),pPdg,mom->GetName(),mPdg);
}
}//electron check
//Cluster was made by something else
- else {
- if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
+ else
+ {
+ if(fDebug > 0)
+ printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",
+ mom->GetName(),mPdg,pPdg);
+
SetTagBit(tag,kMCUnknown);
}
}//Good label value
- else{// Bad label
-
+ else
+ {// Bad label
if(label < 0 && (fDebug >= 0))
printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***: label %d \n", label);
- if(label >= stack->GetNtrack() && (fDebug >= 0))
+
+ if(label >= stack->GetNtrack() && (fDebug >= 0))
printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
+
SetTagBit(tag,kMCUnknown);
}//Bad label
{
// Play with the MCParticles in AOD if available. Tag particles depending on their origin.
// Do same things as in CheckOriginInStack but different input.
- if(!mcparticles) {
+ if(!mcparticles)
+ {
if(fDebug >= 0)
printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
return -1;
Int_t label=labels[0];//Most significant particle contributing to the cluster
Int_t nprimaries = mcparticles->GetEntriesFast();
- if(label >= 0 && label < nprimaries){
+ if(label >= 0 && label < nprimaries)
+ {
//Mother
AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
Int_t iMom = label;
Int_t mPdgSign = mom->GetPdgCode();
Int_t mPdg = TMath::Abs(mPdgSign);
Int_t iParent = mom->GetMother() ;
+
if(fDebug > 0 && label < 8 && fMCGenerator!="") printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
//GrandParent
AliAODMCParticle * parent = NULL ;
Int_t pPdg = -1;
- if(iParent >= 0){
+ if(iParent >= 0)
+ {
parent = (AliAODMCParticle *) mcparticles->At(iParent);
pPdg = TMath::Abs(parent->GetPdgCode());
}
else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
- if(fDebug > 2 ) {
+ if(fDebug > 2 )
+ {
printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n");
- printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
+
+ printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
+ iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
+
if(parent)
- printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
+ printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
+ iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
}
//Check if mother is converted, if not, get the first non converted mother
- if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()){
+ if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary())
+ {
SetTagBit(tag,kMCConversion);
+
//Check if the mother is photon or electron with status not stable
- while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) {
+ while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary())
+ {
//Mother
iMom = mom->GetMother();
mom = (AliAODMCParticle *) mcparticles->At(iMom);
if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
//GrandParent
- if(iParent >= 0 && parent){
+ if(iParent >= 0 && parent)
+ {
parent = (AliAODMCParticle *) mcparticles->At(iParent);
pPdg = TMath::Abs(parent->GetPdgCode());
}
}//while
- if(fDebug > 2 ) {
+ if(fDebug > 2 )
+ {
printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n");
- printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
+
+ printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
+ iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
+
if(parent)
- printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
+ printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
+ iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
}
}//mother and parent are electron or photon and have status 0 and parent is photon or electron
- else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary()){
+ else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary())
+ {
//Still a conversion but only one electron/photon generated. Just from hadrons
if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
- pPdg == 2212 || pPdg == 130 || pPdg == 13 ) {
+ pPdg == 2212 || pPdg == 130 || pPdg == 13 )
+ {
SetTagBit(tag,kMCConversion);
iMom = mom->GetMother();
mom = (AliAODMCParticle *) mcparticles->At(iMom);
mPdgSign = mom->GetPdgCode();
mPdg = TMath::Abs(mPdgSign);
- if(fDebug > 2 ) {
+ if(fDebug > 2 )
+ {
printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n");
printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
}
else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
//check for pi0 and eta (shouldn't happen unless their decays were turned off)
- else if(mPdg == 111) {
+ else if(mPdg == 111)
+ {
SetTagBit(tag,kMCPi0Decay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
}
- else if(mPdg == 221) {
+ else if(mPdg == 221)
+ {
SetTagBit(tag,kMCEtaDecay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
}
//Photons
- else if(mPdg == 22){
+ else if(mPdg == 22)
+ {
SetTagBit(tag,kMCPhoton);
if(mom->IsPhysicalPrimary() && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG")) //undecayed particle
{
- if(iParent < 8 && iParent > 5 ) {//outgoing partons
+ if(iParent < 8 && iParent > 5 )
+ {//outgoing partons
if(pPdg == 22) SetTagBit(tag,kMCPrompt);
else SetTagBit(tag,kMCFragmentation);
}//Outgoing partons
- else if(iParent <= 5 && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG")) {
+ else if(iParent <= 5 && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG"))
+ {
SetTagBit(tag, kMCISR); //Initial state radiation
}
- else if(parent && parent->IsPrimary() && !parent->IsPhysicalPrimary()){//Decay
- if(pPdg == 111){
+ else if(parent && parent->IsPrimary() && !parent->IsPhysicalPrimary())
+ {//Decay
+ if(pPdg == 111)
+ {
SetTagBit(tag,kMCPi0Decay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
}
- else if (pPdg == 221) {
+ else if (pPdg == 221)
+ {
SetTagBit(tag, kMCEtaDecay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
}
else SetTagBit(tag,kMCOtherDecay);
}//Decay
- else {
- if(parent)printf("AliMCAnalysisUtils::CheckOriginInAOD() - what is it? Mother mPdg %d, is primary? %d, is physical %d \n Parent iParent %d, pPdg %d, is primary? %d, is physical? %d\n",
- mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(),iParent, pPdg,parent->IsPrimary(), parent->IsPhysicalPrimary());
+ else
+ {
+ if(parent)
+ printf("AliMCAnalysisUtils::CheckOriginInAOD() - what is it? Mother mPdg %d, is primary? %d, is physical %d \n Parent iParent %d, pPdg %d, is primary? %d, is physical? %d\n",
+ mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(),iParent, pPdg,parent->IsPrimary(), parent->IsPhysicalPrimary());
+
SetTagBit(tag,kMCOtherDecay);//Check
}
}//Physical primary
- else if(!mom->IsPrimary()){ //Decays
- if(pPdg == 111){
+ else if(!mom->IsPrimary())
+ { //Decays
+ if(pPdg == 111)
+ {
SetTagBit(tag,kMCPi0Decay);
- if(fDebug > 2 )
- printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC pi0 decay photon \n");
- CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
+
+ if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC pi0 decay photon \n");
+
+ CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
}
- else if (pPdg == 221) {
+ else if (pPdg == 221)
+ {
SetTagBit(tag,kMCEtaDecay);
+
if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC eta decay photon \n");
+
CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
}
else SetTagBit(tag,kMCOtherDecay);
}//not primary : geant generated, decays
- else {
+ else
+ {
//printf("UNKNOWN 1, mom pdg %d, primary %d, physical primary %d; parent %d, pdg %d, primary %d, physical primary %d \n",
//mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(), iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
SetTagBit(tag,kMCUnknown);
}//Mother Photon
//Electron check. Where did that electron come from?
- else if(mPdg == 11){ //electron
- if(pPdg == 11 && parent){
+ else if(mPdg == 11)
+ { //electron
+ if(pPdg == 11 && parent)
+ {
Int_t iGrandma = parent->GetMother();
- if(iGrandma >= 0) {
+ if(iGrandma >= 0)
+ {
AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma);
Int_t gPdg = TMath::Abs(gma->GetPdgCode());
- if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
+ if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
}
}
- SetTagBit(tag,kMCElectron);
- if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
+
+ SetTagBit(tag,kMCElectron);
+
+ if (fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
- else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //c-hadron decay check
- if(parent){
+ else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000))
+ { //c-hadron decay check
+ if(parent)
+ {
Int_t iGrandma = parent->GetMother();
- if(iGrandma >= 0) {
+ if(iGrandma >= 0)
+ {
AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother
Int_t gPdg = TMath::Abs(gma->GetPdgCode());
if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
else SetTagBit(tag,kMCEFromC); //c-hadron decay
- } else SetTagBit(tag,kMCEFromC); //c-hadron decay
+ }
+ else SetTagBit(tag,kMCEFromC); //c-hadron decay
}//parent
- } else { //prompt or other decay
+ } else
+ { //prompt or other decay
TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
- if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",foo->GetName(), pPdg,foo1->GetName(),mPdg);
+
+ if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",
+ foo->GetName(), pPdg,foo1->GetName(),mPdg);
+
if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
- else SetTagBit(tag,kMCOtherDecay);
+ else SetTagBit(tag,kMCOtherDecay);
}
}//electron check
//cluster was made by something else
- else {
- if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",mPdg,pPdg);
+ else
+ {
+ if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",
+ mPdg,pPdg);
SetTagBit(tag,kMCUnknown);
}
}//Good label value
- else{//Bad label
-
+ else
+ {//Bad label
if(label < 0 && (fDebug >= 0) )
printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***: label %d \n", label);
- if(label >= mcparticles->GetEntriesFast() && (fDebug >= 0) )
+
+ if(label >= mcparticles->GetEntriesFast() && (fDebug >= 0) )
printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***: label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
+
SetTagBit(tag,kMCUnknown);
}//Bad label
//Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack
if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) {
- if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
+ if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
labels[0],stack->GetNtrack(), nlabels);
return;
}
TParticle * meson = stack->Particle(mesonIndex);
Int_t mesonPdg = meson->GetPdgCode();
- if(mesonPdg!=111 && mesonPdg!=221){
- printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Wrong pi0/eta PDG : %d \n",mesonPdg);
+ if(mesonPdg!=111 && mesonPdg!=221)
+ {
+ printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
return;
}
- if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s, label %d\n",meson->GetName(), mesonIndex);
+ if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %s, label %d\n",meson->GetName(), mesonIndex);
//Check if meson decayed into 2 daughters or if both were kept.
- if(meson->GetNDaughters() != 2){
+ if(meson->GetNDaughters() != 2)
+ {
if(fDebug > 2)
- printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
+ printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
return;
}
TParticle *photon1 = stack->Particle(iPhoton1);
//Check if both daughters are photons
- if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22){
+ if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22)
+ {
if(fDebug > 2)
- printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
+ printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
return;
}
if(fDebug > 2)
- printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
+ printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
//Check if both photons contribute to the cluster
Bool_t okPhoton0 = kFALSE;
Bool_t okPhoton1 = kFALSE;
- if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Labels loop:\n");
+ if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Labels loop:\n");
- for(Int_t i = 0; i < nlabels; i++){
+ Bool_t conversion = kFALSE;
+
+ for(Int_t i = 0; i < nlabels; i++)
+ {
if(fDebug > 3) printf("\t at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1);
//If we already found both, break the loop
if(okPhoton0 && okPhoton1) break;
- Int_t index = labels[i];
- if (iPhoton0 == index) {
+ Int_t index = labels[i];
+ if (iPhoton0 == index)
+ {
okPhoton0 = kTRUE;
continue;
}
- else if (iPhoton1 == index) {
+ else if (iPhoton1 == index)
+ {
okPhoton1 = kTRUE;
continue;
}
//Trace back the mother in case it was a conversion
- if(index >= stack->GetNtrack()){
- printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) Particle index %d larger than size of list %d !!\n",index,stack->GetNtrack());
+ if(index >= stack->GetNtrack())
+ {
+ printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) Particle index %d larger than size of list %d !!\n",
+ index,stack->GetNtrack());
continue;
}
TParticle * daught = stack->Particle(index);
Int_t tmpindex = daught->GetFirstMother();
if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex);
- while(tmpindex>=0){
+ while(tmpindex>=0)
+ {
//MC particle of interest is the mother
if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex);
daught = stack->Particle(tmpindex);
- if (iPhoton0 == tmpindex) {
- okPhoton0 = kTRUE;
+ if (iPhoton0 == tmpindex)
+ {
+ conversion = kTRUE;
+ okPhoton0 = kTRUE;
break;
}
- else if (iPhoton1 == tmpindex) {
- okPhoton1 = kTRUE;
+ else if (iPhoton1 == tmpindex)
+ {
+ conversion = kTRUE;
+ okPhoton1 = kTRUE;
break;
}
+
tmpindex = daught->GetFirstMother();
+
}//While to check if pi0/eta daughter was one of these contributors to the cluster
if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0)
- printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Something happens, first label should be from a photon decay!\n");
+ printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Something happens, first label should be from a photon decay!\n");
}//loop on list of labels
//If both photons contribute tag as the corresponding meson.
- if(okPhoton0 && okPhoton1){
+ if(okPhoton0 && okPhoton1)
+ {
if(fDebug > 2)
- printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s OVERLAPPED DECAY \n", meson->GetName());
+ printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %s OVERLAPPED DECAY \n", meson->GetName());
+
+ if(!CheckTagBit(tag,kMCConversion) && conversion) SetTagBit(tag,kMCConversion) ;
if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
- else SetTagBit(tag,kMCEta);
+ else SetTagBit(tag,kMCEta);
}
}
if(meson->GetNDaughters() != 2)
{
if(fDebug > 2)
- printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
+ printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
return;
}
if(fDebug > 3)
printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n");
+ Bool_t conversion = kFALSE;
+
for(Int_t i = 0; i < nlabels; i++)
{
if(fDebug > 3)
if(okPhoton0 && okPhoton1) break;
Int_t index = labels[i];
- if (iPhoton0 == index) {
+ if (iPhoton0 == index)
+ {
okPhoton0 = kTRUE;
continue;
}
- else if (iPhoton1 == index) {
+ else if (iPhoton1 == index)
+ {
okPhoton1 = kTRUE;
continue;
}
//Trace back the mother in case it was a conversion
- if(index >= mcparticles->GetEntriesFast()){
+ if(index >= mcparticles->GetEntriesFast())
+ {
printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Particle index %d larger than size of list %d !!\n",index,mcparticles->GetEntriesFast());
continue;
}
printf("\t parent index %d\n",tmpindex);
daught = (AliAODMCParticle*) mcparticles->At(tmpindex);
//printf("tmpindex %d\n",tmpindex);
- if (iPhoton0 == tmpindex) {
- okPhoton0 = kTRUE;
+ if (iPhoton0 == tmpindex)
+ {
+ conversion = kTRUE;
+ okPhoton0 = kTRUE;
break;
}
- else if (iPhoton1 == tmpindex) {
- okPhoton1 = kTRUE;
+ else if (iPhoton1 == tmpindex)
+ {
+ conversion = kTRUE;
+ okPhoton1 = kTRUE;
break;
- }
+ }
+
tmpindex = daught->GetMother();
+
}//While to check if pi0/eta daughter was one of these contributors to the cluster
if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 )
}//loop on list of labels
//If both photons contribute tag as the corresponding meson.
- if(okPhoton0 && okPhoton1){
+ if(okPhoton0 && okPhoton1)
+ {
if(fDebug > 2)
- printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
+ printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",
+ (TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
+
+ if(!CheckTagBit(tag,kMCConversion) && conversion)
+ {
+ if(fDebug > 2)
+ printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Second decay photon produced a conversion \n");
+ SetTagBit(tag,kMCConversion) ;
+ }
if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
else SetTagBit(tag,kMCEta);
}
//_____________________________________________________________________________________
-TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int_t pdg, const AliCaloTrackReader* reader, Bool_t & ok)
+TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int_t pdg, const AliCaloTrackReader* reader,
+ Bool_t & ok, Int_t & momlabel)
{
//Return the kinematics of the particle that generated the signal
ok = kFALSE;
return grandmom;
}
+
if(label >= 0 && label < reader->GetStack()->GetNtrack())
{
TParticle * momP = reader->GetStack()->Particle(label);
grandmomPDG = grandmomP->GetPdgCode();
if(grandmomPDG==pdg)
{
+ momlabel = grandmomLabel;
grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
break;
}
if(grandmomPDG==pdg)
{
//printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
-
+ momlabel = grandmomLabel;
grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
break;
}
fhEtaEtaPhiNLocMax1(0), fhEtaEtaPhiNLocMax2(0), fhEtaEtaPhiNLocMaxN(0),
fhPi0EPairDiffTimeNLM1(0), fhPi0EPairDiffTimeNLM2(0), fhPi0EPairDiffTimeNLMN(0),
fhEtaEPairDiffTimeNLM1(0), fhEtaEPairDiffTimeNLM2(0), fhEtaEPairDiffTimeNLMN(0),
- fhMCPi0HighNLMPair(0), fhMCPi0LowNLMPair(0),
- fhMCPi0AnyNLMPair(0), fhMCPi0NoneNLMPair(0),
- fhMCPi0HighNLMPairNoMCMatch(0), fhMCPi0LowNLMPairNoMCMatch(0),
- fhMCPi0AnyNLMPairNoMCMatch(0), fhMCPi0NoneNLMPairNoMCMatch(0),
- fhMCPi0DecayPhotonHitHighLM(0), fhMCPi0DecayPhotonAdjHighLM(0),
- fhMCPi0DecayPhotonHitOtherLM(0), fhMCPi0DecayPhotonAdjOtherLM(0),
- fhMCPi0DecayPhotonHitNoLM(0),
- fhMCEOverlapType(0), fhMCEOverlapTypeMatch(0)
+ fhMCPi0HighNLMPair(0), fhMCPi0LowNLMPair(0),
+ fhMCPi0AnyNLMPair(0), fhMCPi0NoneNLMPair(0),
+ fhMCPi0HighNLMPairNoMCMatch(0), fhMCPi0LowNLMPairNoMCMatch(0),
+ fhMCPi0AnyNLMPairNoMCMatch(0), fhMCPi0NoneNLMPairNoMCMatch(0),
+ fhMCPi0HighNLMPairOverlap(0), fhMCPi0LowNLMPairOverlap(0),
+ fhMCPi0AnyNLMPairOverlap(0), fhMCPi0NoneNLMPairOverlap(0),
+ fhMCPi0HighNLMPairNoMCMatchOverlap(0), fhMCPi0LowNLMPairNoMCMatchOverlap(0),
+ fhMCPi0AnyNLMPairNoMCMatchOverlap(0), fhMCPi0NoneNLMPairNoMCMatchOverlap(0),
+ fhMCPi0DecayPhotonHitHighLM(0), fhMCPi0DecayPhotonAdjHighLM(0),
+ fhMCPi0DecayPhotonHitOtherLM(0), fhMCPi0DecayPhotonAdjOtherLM(0),
+ fhMCPi0DecayPhotonAdjacent(0), fhMCPi0DecayPhotonHitNoLM(0),
+ fhMCPi0DecayPhotonHitHighLMOverlap(0), fhMCPi0DecayPhotonAdjHighLMOverlap(0),
+ fhMCPi0DecayPhotonHitOtherLMOverlap(0), fhMCPi0DecayPhotonAdjOtherLMOverlap(0),
+ fhMCPi0DecayPhotonAdjacentOverlap(0), fhMCPi0DecayPhotonHitNoLMOverlap(0),
+ fhMCEOverlapType(0), fhMCEOverlapTypeMatch(0)
{
//default ctor
}
//_______________________________________________________________________________________________________
-void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* cluster, const Int_t mcindex)
+void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* cluster,
+ const Int_t mcindex, const Int_t noverlaps)
{
// Check origin NLM tower of the cluster, when MC gives merged pi0
if(!IsDataMC()) return;
if(mcindex != kmcPi0 && mcindex != kmcPi0Conv) return;
-
+
const UInt_t nc = cluster->GetNCells();
Int_t list[nc];
Float_t elist[nc];
Int_t nMax = GetCaloUtils()->GetNumberOfLocalMaxima(cluster, GetEMCALCells(),list, elist);
-
-// printf("AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin() - Cluster E %2.2f; NLM = %d, cluster MC labels:\n",cluster->E(),nMax);
-//
+
+// if(mcindex==kmcPi0) printf("** Normal Pi0 **\n");
+// if(mcindex==kmcPi0Conv) printf("** Converted Pi0 **\n");
+// printf("** N max %d - Overlaps = %d **\n",nMax, noverlaps);
+//
+// // Study the mothers of cluster
+// printf("Cluster MC labels %d \n", cluster->GetNLabels());
// for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
// {
-// Bool_t ok =kFALSE,gok = kFALSE;
-// Int_t pdg = -22222, status = -1;
-// Int_t gpdg = -22222, gstatus = -1;
-// Int_t ggpdg = -22222, ggstatus = -1;
-// Int_t gLabel = -1, ggLabel = -1;
-//
-// Int_t label = cluster->GetLabels()[ilab];
-// TLorentzVector primary =GetMCAnalysisUtils()->GetMother (label,GetReader(), pdg, status, ok);
-// TLorentzVector gprimary =GetMCAnalysisUtils()->GetGrandMother(label,GetReader(), gpdg, gstatus,gok, gLabel,ggLabel);
-// TLorentzVector ggprimary =GetMCAnalysisUtils()->GetMother(ggLabel ,GetReader(),ggpdg,ggstatus,gok);
-// printf("\t %d; mother: Label %d; PDG %d; E %2.2f - grand mother label %d; PDG %d; E %2.2f- great grand mother label %d; PDG %d; E %2.2f\n",
-// ilab,label,pdg,primary.E(), gLabel,gpdg,gprimary.E(), ggLabel,ggpdg,ggprimary.E());
-//
+// Int_t mclabel = cluster->GetLabels()[ilab];
+//
+// Bool_t mOK = 0;
+// Int_t mpdg = -999999;
+// Int_t mstatus = -1;
+// Int_t grandLabel = -1;
+// TLorentzVector mother = GetMCAnalysisUtils()->GetMother(mclabel,GetReader(),mpdg,mstatus,mOK,grandLabel);
+//
+// printf("******** mother %d : Label %d, pdg %d; status %d, E %2.2f, Eta %2.2f, Phi %2.2f, ok %d, mother label %d\n",
+// ilab, mclabel, mpdg, mstatus,mother.E(), mother.Eta(),mother.Phi()*TMath::RadToDeg(),mOK,grandLabel);
+//
+// if( ( mpdg == 22 || TMath::Abs(mpdg)==11 ) && grandLabel >=0 )
+// {
+// while( ( mpdg == 22 || TMath::Abs(mpdg)==11 ) && grandLabel >=0 )
+// {
+// Int_t newLabel = -1;
+// TLorentzVector grandmother = GetMCAnalysisUtils()->GetMother(grandLabel,GetReader(),mpdg,mstatus,mOK,newLabel);
+// printf("\t grandmother %d : Label %d, pdg %d; status %d, E %2.2f, Eta %2.2f, Phi %2.2f, ok %d, mother label %d\n",
+// ilab, grandLabel, mpdg, mstatus,grandmother.E(), grandmother.Eta(), grandmother.Phi()*TMath::RadToDeg(),mOK,newLabel);
+// grandLabel = newLabel;
+//
+// }
+// }
// }
//
-// printf("AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin() - Cluster Cells MC labels:\n");
-//
-// for (UInt_t icell = 0; icell < nc; icell++ )
+// printf("Cells in cluster %d\n",cluster->GetNCells() );
+// for(Int_t icell = 0; icell < cluster->GetNCells(); icell++)
// {
-// Bool_t ok =kFALSE,gok = kFALSE;
-// Int_t pdg = -22222, status = -1;
-// Int_t gpdg = -22222, gstatus = -1;
-// Int_t ggpdg = -22222, ggstatus = -1;
-// Int_t gLabel = -1, ggLabel = -1;
-//
-// Int_t absId = cluster->GetCellAbsId(icell);
-//
-// printf("cell abs Id %d, amplitude %f\n",absId,GetEMCALCells()->GetCellAmplitude(absId));
+// Int_t absIdCell = cluster->GetCellAbsId(icell);
+// Int_t mcLabel = GetEMCALCells()->GetCellMCLabel(absIdCell);
+// GetReader()->RemapMCLabelForAODs(mcLabel);
//
-// Int_t label = GetEMCALCells()->GetCellMCLabel(absId);
-// TLorentzVector primary =GetMCAnalysisUtils()->GetMother (label,GetReader(), pdg, status, ok);
-// TLorentzVector gprimary =GetMCAnalysisUtils()->GetGrandMother(label,GetReader(), gpdg, gstatus,gok, gLabel,ggLabel);
-// TLorentzVector ggprimary =GetMCAnalysisUtils()->GetMother(ggLabel ,GetReader(),ggpdg,ggstatus,gok);
-// printf(" %d; mother: Label %d; PDG %d; E %2.2f - grand mother label %d; PDG %d; E %2.2f- great grand mother label %d; PDG %d; E %2.2f\n",
-// icell,label,pdg,primary.E(), gLabel,gpdg,gprimary.E(), ggLabel,ggpdg,ggprimary.E());
-
+// printf(" \t cell i %d, abs %d, amp %2.3f, mclabel %d\n",icell,absIdCell,GetEMCALCells()->GetCellAmplitude(absIdCell),mcLabel);
// }
- //printf("en %2.3f, nc %d, nMax %d \n",cluster->E(),nc,nMax);
-
//If only one maxima, consider all the towers in the cluster
if(nMax==1)
{
- for (UInt_t icell = 0; icell < nc; icell++ )
- {
- list [icell] = cluster->GetCellAbsId(icell);
- elist[icell] = GetEMCALCells()->GetCellAmplitude(list[icell]);
- }
+ for (UInt_t icell = 0; icell < nc; icell++ )
+ {
+ list [icell] = cluster->GetCellAbsId(icell);
+ elist[icell] = GetEMCALCells()->GetCellAmplitude(list[icell]);
+ }
}
Int_t nmaxima = nMax;
if(i==imax) continue;
//printf("j %d: AbsId %d; E %2.3f\n",i,list[i],elist[i]);
-
+
if(elist[i] > emax2)
{
//printf("Highest : %d and %d\n",imax,imax2);
+ //---------------------------------------------------------
//---------------------------------------------------------
// Compare ancestors of all local maxima at cell MC level
//---------------------------------------------------------
+ //---------------------------------------------------------
// Check that the highest mc label and the max cluster label are the same
Int_t mcLabelMax = -1 ;
Int_t mcLabelclusterMax = cluster->GetLabels()[0];
Bool_t matchHighLMAndHighMC = kFALSE;
+ //printf("MC label: LM1 %d, LM2 %d, cluster %d\n",mcLabelMax,mcLabelMax2,mcLabelclusterMax);
+
if(mcLabelclusterMax == mcLabelMax && mcLabelclusterMax >= 0)
{
matchHighLMAndHighMC = kTRUE;
- //printf("*** MATCH cluster and LM maximum ***\n");
+ //printf("\t *** MATCH cluster and LM maximum ***\n");
}
else
{
- //printf("*** NO MATCH cluster and LM maximum, check second ***\n");
+ //printf("\t *** NO MATCH cluster and LM maximum, check second ***\n");
if(mcLabelclusterMax == mcLabelMax2 && mcLabelclusterMax >= 0)
{
- //printf("\t *** MATCH cluster and 2nd LM maximum ***\n");
+ //printf("\t \t *** MATCH cluster and 2nd LM maximum ***\n");
matchHighLMAndHighMC = kTRUE;
}
else
{
- //printf("\t *** NO MATCH***\n");
+ //printf("\t \t *** NO MATCH***\n");
matchHighLMAndHighMC = kFALSE;
}
}
// printf("Max index %d; mother: Label %d; PDG %d; E %2.2f - grand mother label %d; PDG %d; E %2.2f- great grand mother label %d; PDG %d; E %2.2f\n",
// i,mcLabel1,pdg,primary.E(), gLabel,gpdg,gprimary.E(), ggLabel,ggpdg,ggprimary.E());
// }
-
-
for(Int_t i = 0; i < nmaxima-1; i++)
{
Bool_t ok =kFALSE;
Int_t pdg = -22222, status = -1;
TLorentzVector primary =GetMCAnalysisUtils()->GetMother(ancLabel,GetReader(), pdg, status, ok);
-
//printf("\t i %d label %d - j %d label %d; ancestor label %d, PDG %d-%d; E %2.2f; high %d, any %d \n",i,mcLabel1,j,mcLabel2, ancLabel, ancPDG,pdg, primary.E(), high, low);
}
//printf("nMax %d; Match MC? %d; high %d; low %d\n",nMax,matchHighLMAndHighMC,high,low);
- if(matchHighLMAndHighMC)
+ if(!noverlaps)
{
- if (high && !low) fhMCPi0HighNLMPair->Fill(en,nMax);
- else if(low && !high) fhMCPi0LowNLMPair ->Fill(en,nMax);
- else if(low && high) fhMCPi0AnyNLMPair ->Fill(en,nMax);
- else fhMCPi0NoneNLMPair->Fill(en,nMax);
+ if(matchHighLMAndHighMC)
+ {
+ if (high && !low) fhMCPi0HighNLMPair->Fill(en,nMax);
+ else if(low && !high) fhMCPi0LowNLMPair ->Fill(en,nMax);
+ else if(low && high) fhMCPi0AnyNLMPair ->Fill(en,nMax);
+ else fhMCPi0NoneNLMPair->Fill(en,nMax);
+ }
+ else
+ {
+ if (high && !low) fhMCPi0HighNLMPairNoMCMatch->Fill(en,nMax);
+ else if(low && !high) fhMCPi0LowNLMPairNoMCMatch ->Fill(en,nMax);
+ else if(low && high) fhMCPi0AnyNLMPairNoMCMatch ->Fill(en,nMax);
+ else fhMCPi0NoneNLMPairNoMCMatch->Fill(en,nMax);
+ }
}
else
{
- if (high && !low) fhMCPi0HighNLMPairNoMCMatch->Fill(en,nMax);
- else if(low && !high) fhMCPi0LowNLMPairNoMCMatch ->Fill(en,nMax);
- else if(low && high) fhMCPi0AnyNLMPairNoMCMatch ->Fill(en,nMax);
- else fhMCPi0NoneNLMPairNoMCMatch->Fill(en,nMax);
+ if(matchHighLMAndHighMC)
+ {
+ if (high && !low) fhMCPi0HighNLMPairOverlap->Fill(en,nMax);
+ else if(low && !high) fhMCPi0LowNLMPairOverlap->Fill(en,nMax);
+ else if(low && high) fhMCPi0AnyNLMPairOverlap->Fill(en,nMax);
+ else fhMCPi0NoneNLMPairOverlap->Fill(en,nMax);
+ }
+ else
+ {
+ if (high && !low) fhMCPi0HighNLMPairNoMCMatchOverlap->Fill(en,nMax);
+ else if(low && !high) fhMCPi0LowNLMPairNoMCMatchOverlap->Fill(en,nMax);
+ else if(low && high) fhMCPi0AnyNLMPairNoMCMatchOverlap->Fill(en,nMax);
+ else fhMCPi0NoneNLMPairNoMCMatchOverlap->Fill(en,nMax);
+ }
}
-
+ //----------------------------------------------------------------------
//----------------------------------------------------------------------
// Compare MC decay photon projection to cell location and Local Maxima
//----------------------------------------------------------------------
+ //----------------------------------------------------------------------
// Get the mother pi0
Float_t phi0 = photon0Kine.Phi();
Float_t phi1 = photon1Kine.Phi();
+// printf("MC pi0 label %d E %2.2f, eta %2.2f, phi %2.2f: \n \t photon0 label %d E %2.2f, eta %2.2f, phi %2.2f \n \t photon1 label %d E %2.2f eta %2.2f, phi %2.2f\n",
+// label , pi0Kine.E() , pi0Kine.Eta(),pi0Kine.Phi()*TMath::RadToDeg(),
+// label0, photon0Kine.E(), eta0, phi0*TMath::RadToDeg(),
+// label1, photon1Kine.E(), eta1, phi1*TMath::RadToDeg());
+//
+// TLorentzVector momclus;
+// cluster->GetMomentum(momclus,GetVertex(0));
+// printf("Cluster E %2.2F eta %2.2f, phi %f\n",momclus.E(),momclus.Eta(),momclus.Phi()*TMath::RadToDeg());
+
if(phi0 < 0 ) phi0+=TMath::TwoPi();
if(phi1 < 0 ) phi1+=TMath::TwoPi();
return;
}
+ //-----------------------------------------------
// Check that the 2 photons hit the Local Maxima
+ //-----------------------------------------------
- //printf("Photons AbsId (%d,%d); Local Maxima AbsId(%d,%d)\n",absId0,absId1,list[imax],list[imax2]);
+ // printf("Photons AbsId (%d,%d); Local Maxima AbsId(%d,%d)\n",absId0,absId1,list[imax],list[imax2]);
//printf("Photon1 (eta,phi)=(%f,%f); Photon2 (eta,phi)=(%f,%f);\n",eta0,phi0*TMath::RadToDeg(),eta1,phi1*TMath::RadToDeg());
- Bool_t matchMCHitLM = kFALSE;
- if(imax >= 0 && imax2 >=0 )
+ Bool_t match0 = kFALSE;
+ Bool_t match1 = kFALSE;
+ Int_t imatch0 = -1;
+ Int_t imatch1 = -1;
+ if(imax >= 0 && imax2 >=0 && absId0 > 0 && absId1 > 0 )
{
- if (absId0 == list[imax] && absId1 == list[imax2]) matchMCHitLM = kTRUE;
- else if(absId1 == list[imax] && absId0 == list[imax2]) matchMCHitLM = kTRUE;
+ if (absId0 == list[imax] ) { match0 = kTRUE ; imatch0 = imax ; }
+ else if(absId0 == list[imax2]) { match0 = kTRUE ; imatch0 = imax2 ; }
+
+ if (absId1 == list[imax] ) { match1 = kTRUE ; imatch1 = imax ; }
+ else if(absId1 == list[imax2]) { match1 = kTRUE ; imatch1 = imax2 ; }
}
- //Check the adjacent cells
- Bool_t adjacent = kFALSE;
- Int_t ieta0=-1; Int_t iphi0 = 0; Int_t rcu0 = 0;
- GetModuleNumberCellIndexes(absId0,fCalorimeter, ieta0, iphi0, rcu0);
- Int_t ieta1=-1; Int_t iphi1 = 0; Int_t rcu1 = 0;
- GetModuleNumberCellIndexes(absId1,fCalorimeter, ieta1, iphi1, rcu1);
+ //printf("primary imatch0 %d, imatch1 %d\n",imatch0,imatch1);
+
+ // If one or the 2 not matched, check with the other MC labels
+ // only in case there was a conversion
+
+ Int_t absId0second = -1;
+ Int_t absId1second = -1;
+ Int_t secLabel0 = -1;
+ Int_t secLabel1 = -1;
+ Int_t mcLabel0 = -1;
+ Int_t mcLabel1 = -1;
+ Bool_t secOK = 0;
+ Int_t secpdg = -999999;
+ Int_t secstatus = -1;
+ Int_t secgrandLabel = -1;
+
+ if(match0) { secLabel0 = label0 ; mcLabel0 = label0 ; }
+ if(match1) { secLabel1 = label1 ; mcLabel1 = label1 ; }
- if(!matchMCHitLM)
+ if((!match0 || !match1) && mcindex == kmcPi0Conv)
{
- Int_t ietam0=-1; Int_t iphim0 = 0; Int_t rcum0 = 0 ;
- if(imax >= 0) GetModuleNumberCellIndexes(list[imax] ,fCalorimeter, ietam0, iphim0, rcum0);
- Int_t ietam1=-1; Int_t iphim1 = 0; Int_t rcum1 = 0 ;
- if(imax2 >= 0) GetModuleNumberCellIndexes(list[imax2],fCalorimeter, ietam1, iphim1, rcum1);
+ for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
+ {
+ Int_t mclabel = cluster->GetLabels()[ilab];
+
+ //printf("Check label %d - %d\n",ilab,mclabel);
+
+ if(mclabel == label0 || mclabel == label1)
+ {
+ //printf("continue: secLabel %d, label0 %d, label1 %d\n",mclabel,label0,label1);
+ if(mclabel == label0 && secLabel0 < 0) { secLabel0 = label0 ; mcLabel0 = label0 ; }
+ if(mclabel == label1 && secLabel1 < 0) { secLabel1 = label1 ; mcLabel1 = label1 ; }
+ continue ;
+ }
+
+ //printf("Before while: secLabel0 %d, secLabel1 %d\n",secLabel0,secLabel1);
+
+ // match mc label and parent photon
+ Int_t tmplabel = mclabel;
+ while((secLabel0 < 0 || secLabel1 < 0) && tmplabel > 0 )
+ {
+ TLorentzVector mother = GetMCAnalysisUtils()->GetMother(tmplabel,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
+
+ //printf("\t \t while secLabel %d, mom %d, granmom %d\n",mclabel,tmplabel,secgrandLabel);
+
+ if((secgrandLabel == label0) || (secgrandLabel == label1 ))
+ {
+ //printf("mcMatch! grand label %d, secLabel %d\n",secgrandLabel, mclabel);
+ if(!match0 && mcLabel1 != secgrandLabel) { secLabel0 = mclabel; mcLabel0 = secgrandLabel; }
+ if(!match1 && mcLabel0 != secgrandLabel) { secLabel1 = mclabel; mcLabel1 = secgrandLabel; }
+ }
+
+ //printf("\t GrandMother %d, secLabel0 %d, secLabel1 %d \n",secgrandLabel, secLabel0,secLabel1);
- Bool_t adjacent0m0 = kFALSE;
- Bool_t adjacent1m0 = kFALSE;
- Bool_t adjacent0m1 = kFALSE;
- Bool_t adjacent1m1 = kFALSE;
+ tmplabel = secgrandLabel;
+ }
+ }
- if( iphim0>=0 && ietam0 >= 0 )
+ // Get the position of the found secondaries mother
+ if(!match0 && secLabel0 > 0)
{
- if(TMath::Abs(ieta0-ietam0) == 1 && TMath::Abs(iphi0-iphim0) == 0 ) adjacent0m0 = kTRUE;
- if(TMath::Abs(ieta0-ietam0) == 0 && TMath::Abs(iphi0-iphim0) == 1 ) adjacent0m0 = kTRUE;
- if(TMath::Abs(ieta0-ietam0) == 1 && TMath::Abs(iphi0-iphim0) == 1 ) adjacent0m0 = kTRUE;
-
- if(TMath::Abs(ieta1-ietam0) == 1 && TMath::Abs(iphi1-iphim0) == 0 ) adjacent1m0 = kTRUE;
- if(TMath::Abs(ieta1-ietam0) == 0 && TMath::Abs(iphi1-iphim0) == 1 ) adjacent1m0 = kTRUE;
- if(TMath::Abs(ieta1-ietam0) == 1 && TMath::Abs(iphi1-iphim0) == 1 ) adjacent1m0 = kTRUE;
+ TLorentzVector mother = GetMCAnalysisUtils()->GetMother(secLabel0,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
+
+ Float_t eta = mother.Eta();
+ Float_t phi = mother.Phi();
+ if(phi < 0 ) phi+=TMath::TwoPi();
+ GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta, phi, absId0second);
+
+ //printf("Secondary MC0 label %d, absId %d E %2.2F eta %2.2f, phi %f\n", secLabel0,absId0second, mother.E(),mother.Eta(),mother.Phi()*TMath::RadToDeg());
+
+ if(absId0second == list[imax] ) { match0 = kTRUE ; imatch0 = imax ; }
+ if(absId0second == list[imax2]) { match0 = kTRUE ; imatch0 = imax2 ; }
}
-
- if( iphim1>=0 && ietam1 >= 0 )
+
+ if(!match1 && secLabel1 > 0)
{
- if(TMath::Abs(ieta0-ietam1) == 1 && TMath::Abs(iphi0-iphim1) == 0 ) adjacent0m1 = kTRUE;
- if(TMath::Abs(ieta0-ietam1) == 0 && TMath::Abs(iphi0-iphim1) == 1 ) adjacent0m1 = kTRUE;
- if(TMath::Abs(ieta0-ietam1) == 1 && TMath::Abs(iphi0-iphim1) == 1 ) adjacent0m1 = kTRUE;
-
- if(TMath::Abs(ieta1-ietam1) == 1 && TMath::Abs(iphi1-iphim1) == 0 ) adjacent1m1 = kTRUE;
- if(TMath::Abs(ieta1-ietam1) == 0 && TMath::Abs(iphi1-iphim1) == 1 ) adjacent1m1 = kTRUE;
- if(TMath::Abs(ieta1-ietam1) == 1 && TMath::Abs(iphi1-iphim1) == 1 ) adjacent1m1 = kTRUE;
+ TLorentzVector mother = GetMCAnalysisUtils()->GetMother(secLabel1,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
+
+ Float_t eta = mother.Eta();
+ Float_t phi = mother.Phi();
+ if(phi < 0 ) phi+=TMath::TwoPi();
+ GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta, phi, absId1second);
+
+ //printf("Secondary MC1 label %d absId %d E %2.2F eta %2.2f, phi %f\n",secLabel1, absId1second, mother.E(),mother.Eta(),mother.Phi()*TMath::RadToDeg());
+
+ if(absId1second == list[imax] ) { match1 = kTRUE ; imatch1 = imax ; }
+ if(absId1second == list[imax2]) { match1 = kTRUE ; imatch1 = imax2 ; }
}
+
+ //printf("secondary label mc0 %d, mc1 %d, imatch0 %d, imatch1 %d\n",secLabel0,secLabel1,imatch0,imatch1);
+
+ }
+
+ //printf("imatch0 %d, imatch1 %d\n",imatch0,imatch1);
+ if( match0 && match1 )
+ {
+ //printf("a) Both Photons hit local maxima \n");
+
+ if(!noverlaps)fhMCPi0DecayPhotonHitHighLM ->Fill(en,nMax);
+ else fhMCPi0DecayPhotonHitHighLMOverlap->Fill(en,nMax);
+
+ return ;
+ }
+
+ //printf("Any match? photon0 %d, photon1 %d\n",match0,match1);
+ //if(!match0 && !match1) printf("WARNING, LM not matched to any photon decay!\n");
+
+ //---------------------------------------------
+ // Check the adjacent cells to the local maxima
+ //---------------------------------------------
+
+// Int_t ieta0=-1; Int_t iphi0 = 0; Int_t rcu0 = 0;
+// Int_t sm0 = GetModuleNumberCellIndexes(absId0,fCalorimeter, ieta0, iphi0, rcu0);
+// Int_t ieta1=-1; Int_t iphi1 = 0; Int_t rcu1 = 0;
+// Int_t sm1 = GetModuleNumberCellIndexes(absId1,fCalorimeter, ieta1, iphi1, rcu1);
+//
+// printf("Photon1 (id,sm,eta,phi)=(%d,%d,%d,%d), Photon2 (id,sm,eta,phi)=(%d,%d,%d,%d)\n",
+// absId0,sm0,ieta0,iphi0,absId1,sm1,ieta1,iphi1);
+//
+// Int_t ietam0=-1; Int_t iphim0 = 0; Int_t rcum0 = 0; Int_t smm0 = -1 ;
+// if(imax >= 0) smm0 = GetModuleNumberCellIndexes(list[imax] ,fCalorimeter, ietam0, iphim0, rcum0);
+// Int_t ietam1=-1; Int_t iphim1 = 0; Int_t rcum1 = 0; Int_t smm1 = -1 ;
+// if(imax2 >= 0) smm1 = GetModuleNumberCellIndexes(list[imax2],fCalorimeter, ietam1, iphim1, rcum1);
+//
+// printf("Max (id, sm,eta,phi)=(%d,%d,%d,%d), Max2 (id, sm,eta,phi)=(%d,%d,%d,%d), imatch0 %d, imatch1 %d\n",
+// list[imax],smm0,ietam0,iphim0,list[imax2],smm1,ietam1,iphim1,imatch0,imatch1);
+
+ if(!match0)
+ {
+ if(imatch1!=imax && GetCaloUtils()->AreNeighbours(fCalorimeter,absId0,list[imax])) match0 = kTRUE;
+ //printf("imax - match0? (%d-%d)=%d, (%d-%d)=%d\n",ieta0,ietam0,ieta0-ietam0, iphi0,iphim0,iphi0-iphim0);
+ if(imatch1!=imax2 && GetCaloUtils()->AreNeighbours(fCalorimeter,absId0,list[imax2]) ) match0 = kTRUE;
+ //printf("imax2 - match0? (%d-%d)=%d, (%d-%d)=%d\n",ieta0,ietam1,ieta0-ietam1, iphi0,iphim1,iphi0-iphim1);
+ }
+
+ if(!match1)
+ {
+ if(imatch0!=imax && GetCaloUtils()->AreNeighbours(fCalorimeter,absId1,list[imax]) ) match1 = kTRUE;
+ //printf("imax - match1? (%d-%d)=%d, (%d-%d)=%d\n",ieta1,ietam0,ieta1-ietam0, iphi1,iphim0,iphi1-iphim0);
+
+ if(imatch0!=imax2 && GetCaloUtils()->AreNeighbours(fCalorimeter,absId1,list[imax2])) match1 = kTRUE;
+ //printf("imax2 - match1? (%d-%d)=%d, (%d-%d)=%d\n",ieta1,ietam1,ieta1-ietam1, iphi1,iphim1,iphi1-iphim1);
+ }
- if(adjacent1m1 || adjacent1m0 || adjacent0m1 || adjacent0m0 ) adjacent = kTRUE;
+ //printf("Local Maxima: adjacent0 %d,adjacent1 %d \n",match0,match1);
+
+ if(match0 && match1)
+ {
+ //printf("b) Both Photons hit local maxima or cell adjacent or 2 cells adjacent \n");
+ if(!noverlaps) fhMCPi0DecayPhotonAdjHighLM ->Fill(en,nMax);
+ else fhMCPi0DecayPhotonAdjHighLMOverlap->Fill(en,nMax);
+
+ return;
}
+
+
+ // Decay photon cells are adjacent?
+
+ if( (match0 || match1) && GetCaloUtils()->AreNeighbours(fCalorimeter,absId0,absId1) )
+ {
+ //printf("c) Both Photons hit a local maxima and in adjacent cells \n");
+ if(!noverlaps) fhMCPi0DecayPhotonAdjacent ->Fill(en,nMax);
+ else fhMCPi0DecayPhotonAdjacentOverlap ->Fill(en,nMax);
+
+ return;
+ }
+
+ //--------------------
+ // Other Local maxima
+ //--------------------
Bool_t matchMCHitOtherLM = kFALSE;
- if(!adjacent && !matchMCHitLM)
+ if(!match1)
+ {
+ for(Int_t i = 0; i < nmaxima; i++)
+ {
+ if(imax!=i && imax2!=i && absId1 == list[i]) { match1 = kTRUE; matchMCHitOtherLM = kTRUE; }
+ }
+ }
+
+ if(!match0)
{
- for(Int_t i = 0; i < nmaxima-1; i++)
+ for(Int_t i = 0; i < nmaxima; i++)
{
- Int_t mcLabel1 = GetEMCALCells()->GetCellMCLabel(list[i]);
- GetReader()->RemapMCLabelForAODs(mcLabel1);
+ if(imax!=i && imax2!=i && absId0 == list[i]) { match0 = kTRUE; matchMCHitOtherLM = kTRUE; }
+ }
+ }
+
+ if(matchMCHitOtherLM)
+ {
+ //printf("d) One Photon hits a local maxima, the other another not high \n");
+
+ if(!noverlaps) fhMCPi0DecayPhotonHitOtherLM ->Fill(en,nMax);
+ else fhMCPi0DecayPhotonHitOtherLMOverlap->Fill(en,nMax);
+
+ return ;
+ }
+
+ // Adjacent to other maxima
+
+ Bool_t adjacentOther1 = kFALSE;
+ if(match0)
+ {
+ for(Int_t i = 0; i < nmaxima; i++)
+ {
+ Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
+ GetModuleNumberCellIndexes(list[i] ,fCalorimeter, ieta, iphi, rcu);
- for(Int_t j = i+1; j < nmaxima; j++)
- {
- if (absId0==list[i] && absId1 == list[j]) matchMCHitOtherLM = kTRUE;
- else if(absId1==list[i] && absId0 == list[j]) matchMCHitOtherLM = kTRUE;
- }
+ //printf(" Other Max (eta,phi)=(%d,%d)\n",ieta,iphi);
+
+ if(GetCaloUtils()->AreNeighbours(fCalorimeter,absId1,list[i]) ) adjacentOther1 = kTRUE;
+
+ //printf("Other Maxima: adjacentOther1 %d\n",adjacentOther1);
}
}
- Bool_t adjacentOther = kFALSE;
- if(!adjacent && ! matchMCHitLM && !matchMCHitOtherLM)
+ Bool_t adjacentOther0 = kFALSE;
+ if(match1)
{
- for(Int_t i = 0; i < nmaxima-1; i++)
+ for(Int_t i = 0; i < nmaxima; i++)
{
- Int_t mcLabel1 = GetEMCALCells()->GetCellMCLabel(list[i]);
- GetReader()->RemapMCLabelForAODs(mcLabel1);
+ Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
+ GetModuleNumberCellIndexes(list[i] ,fCalorimeter, ieta, iphi, rcu);
- for(Int_t j = i+1; j < nmaxima; j++)
- {
- Int_t ietam0=-1; Int_t iphim0 = 0; Int_t rcum0 = 0;
- GetModuleNumberCellIndexes(list[i] ,fCalorimeter, ietam0, iphim0, rcum0);
- Int_t ietam1=-1; Int_t iphim1 = 0; Int_t rcum1 = 0;
- GetModuleNumberCellIndexes(list[j],fCalorimeter, ietam1, iphim1, rcum1);
-
- Bool_t adjacent0m0 = kFALSE;
- Bool_t adjacent1m0 = kFALSE;
- Bool_t adjacent0m1 = kFALSE;
- Bool_t adjacent1m1 = kFALSE;
-
- if(TMath::Abs(ieta0-ietam0) == 1 && TMath::Abs(iphi0-iphim0) == 0 ) adjacent0m0 = kTRUE;
- if(TMath::Abs(ieta0-ietam0) == 0 && TMath::Abs(iphi0-iphim0) == 1 ) adjacent0m0 = kTRUE;
- if(TMath::Abs(ieta0-ietam0) == 1 && TMath::Abs(iphi0-iphim0) == 1 ) adjacent0m0 = kTRUE;
-
- if(TMath::Abs(ieta1-ietam0) == 1 && TMath::Abs(iphi1-iphim0) == 0 ) adjacent1m0 = kTRUE;
- if(TMath::Abs(ieta1-ietam0) == 0 && TMath::Abs(iphi1-iphim0) == 1 ) adjacent1m0 = kTRUE;
- if(TMath::Abs(ieta1-ietam0) == 1 && TMath::Abs(iphi1-iphim0) == 1 ) adjacent1m0 = kTRUE;
-
- if(TMath::Abs(ieta0-ietam1) == 1 && TMath::Abs(iphi0-iphim1) == 0 ) adjacent0m1 = kTRUE;
- if(TMath::Abs(ieta0-ietam1) == 0 && TMath::Abs(iphi0-iphim1) == 1 ) adjacent0m1 = kTRUE;
- if(TMath::Abs(ieta0-ietam1) == 1 && TMath::Abs(iphi0-iphim1) == 1 ) adjacent0m1 = kTRUE;
-
- if(TMath::Abs(ieta1-ietam1) == 1 && TMath::Abs(iphi1-iphim1) == 0 ) adjacent1m1 = kTRUE;
- if(TMath::Abs(ieta1-ietam1) == 0 && TMath::Abs(iphi1-iphim1) == 1 ) adjacent1m1 = kTRUE;
- if(TMath::Abs(ieta1-ietam1) == 1 && TMath::Abs(iphi1-iphim1) == 1 ) adjacent1m1 = kTRUE;
-
- if(adjacent1m1 || adjacent1m0 || adjacent0m1 || adjacent0m0 ) adjacentOther = kTRUE;
-
- }
+ //printf(" Other Max (eta,phi)=(%d,%d)\n",ieta,iphi);
+
+ if(GetCaloUtils()->AreNeighbours(fCalorimeter,absId0,list[i]) ) adjacentOther0 = kTRUE;
+
+ //printf("Other Maxima: adjacentOther0 %d\n",adjacentOther0);
}
}
- //printf("nMax %d; Match HitHigh? %d; Hit Adj High? %d; Hit Other %d, Hit Other Ajd %d\n",nMax,matchMCHitLM,adjacent,matchMCHitOtherLM,adjacentOther);
-
- if (matchMCHitLM) fhMCPi0DecayPhotonHitHighLM ->Fill(en,nMax);
- else if(adjacent) fhMCPi0DecayPhotonAdjHighLM ->Fill(en,nMax);
- else if(matchMCHitOtherLM) fhMCPi0DecayPhotonHitOtherLM->Fill(en,nMax);
- else if(adjacentOther) fhMCPi0DecayPhotonAdjOtherLM->Fill(en,nMax);
- else fhMCPi0DecayPhotonHitNoLM ->Fill(en,nMax);
+ if((match0 && adjacentOther1) || (match1 && adjacentOther0))
+ {
+ //printf("e) One Photon hits a local maxima, the other another not high, adjacent \n");
+
+ if(!noverlaps) fhMCPi0DecayPhotonAdjOtherLM ->Fill(en,nMax);
+ else fhMCPi0DecayPhotonAdjOtherLMOverlap->Fill(en,nMax);
+
+ return;
+ }
+
+ //printf("f) No hit found \n");
+ if(!noverlaps) fhMCPi0DecayPhotonHitNoLM ->Fill(en,nMax);
+ else fhMCPi0DecayPhotonHitNoLMOverlap->Fill(en,nMax);
}
outputContainer->Add(fhMCPi0NoneNLMPairNoMCMatch) ;
- fhMCPi0DecayPhotonHitHighLM = new TH2F("hMCPi0DecayPhotonHitHighLM ","NLM vs E for merged pi0 cluster, decay photon hit High Local Maxima",
+
+
+
+ fhMCPi0HighNLMPairOverlap = new TH2F("hMCPi0HighNLMPairOverlap","NLM vs E for merged pi0 cluster, high energy NLM pair are decays",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0HighNLMPairOverlap ->SetYTitle("N maxima");
+ fhMCPi0HighNLMPairOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0HighNLMPairOverlap) ;
+
+ fhMCPi0LowNLMPairOverlap = new TH2F("hMCPi0LowNLMPairOverlap","NLM vs E for merged pi0 cluster, lower energy NLM pair are decays",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0LowNLMPairOverlap ->SetYTitle("N maxima");
+ fhMCPi0LowNLMPairOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0LowNLMPairOverlap) ;
+
+ fhMCPi0AnyNLMPairOverlap = new TH2F("hMCPi0AnyNLMPairOverlap","NLM vs E for merged pi0 cluster, both high and other energy NLM pair are decays",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0AnyNLMPairOverlap ->SetYTitle("N maxima");
+ fhMCPi0AnyNLMPairOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0AnyNLMPairOverlap) ;
+
+ fhMCPi0NoneNLMPairOverlap = new TH2F("hMCPi0NoneNLMPairOverlap","NLM vs E for merged pi0 cluster, no NLM pair are decays",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0NoneNLMPairOverlap ->SetYTitle("N maxima");
+ fhMCPi0NoneNLMPairOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0NoneNLMPairOverlap) ;
+
+ fhMCPi0HighNLMPairNoMCMatchOverlap = new TH2F("hMCPi0HighNLMPairNoMCMatchOverlap","NLM vs E for merged pi0 cluster, high energy NLM pair are decays",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0HighNLMPairNoMCMatchOverlap ->SetYTitle("N maxima");
+ fhMCPi0HighNLMPairNoMCMatchOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0HighNLMPairNoMCMatchOverlap) ;
+
+ fhMCPi0LowNLMPairNoMCMatchOverlap = new TH2F("hMCPi0LowNLMPairNoMCMatchOverlap","NLM vs E for merged pi0 cluster, lower energy NLM pair are decays",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0LowNLMPairNoMCMatchOverlap ->SetYTitle("N maxima");
+ fhMCPi0LowNLMPairNoMCMatchOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0LowNLMPairNoMCMatchOverlap) ;
+
+ fhMCPi0AnyNLMPairNoMCMatchOverlap = new TH2F("hMCPi0AnyNLMPairNoMCMatchOverlap","NLM vs E for merged pi0 cluster, both high and other energy NLM pair are decays",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0AnyNLMPairNoMCMatchOverlap ->SetYTitle("N maxima");
+ fhMCPi0AnyNLMPairNoMCMatchOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0AnyNLMPairNoMCMatchOverlap) ;
+
+ fhMCPi0NoneNLMPairNoMCMatchOverlap = new TH2F("hMCPi0NoneNLMPairNoMCMatchOverlap","NLM vs E for merged pi0 cluster, no NLM pair are decays",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0NoneNLMPairNoMCMatchOverlap ->SetYTitle("N maxima");
+ fhMCPi0NoneNLMPairNoMCMatchOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0NoneNLMPairNoMCMatchOverlap) ;
+
+ fhMCPi0DecayPhotonHitHighLM = new TH2F("hMCPi0DecayPhotonHitHighLM","NLM vs E for merged pi0 cluster, decay photon hit High Local Maxima",
nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
fhMCPi0DecayPhotonHitHighLM ->SetYTitle("N maxima");
fhMCPi0DecayPhotonHitHighLM ->SetXTitle("E (GeV)");
outputContainer->Add(fhMCPi0DecayPhotonHitHighLM ) ;
- fhMCPi0DecayPhotonAdjHighLM = new TH2F("hMCPi0DecayPhotonAdjHighLM ","NLM vs E for merged pi0 cluster, decay photon hit cells adjacent to High Local Maxima",
+ fhMCPi0DecayPhotonAdjHighLM = new TH2F("hMCPi0DecayPhotonAdjHighLM","NLM vs E for merged pi0 cluster, decay photon hit cells adjacent to High Local Maxima",
nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
fhMCPi0DecayPhotonAdjHighLM ->SetYTitle("N maxima");
fhMCPi0DecayPhotonAdjHighLM ->SetXTitle("E (GeV)");
outputContainer->Add(fhMCPi0DecayPhotonAdjHighLM ) ;
- fhMCPi0DecayPhotonHitOtherLM = new TH2F("hMCPi0DecayPhotonHitOtherLM ","NLM vs E for merged pi0 cluster, decay photon hit Other Local Maxima",
+ fhMCPi0DecayPhotonHitOtherLM = new TH2F("hMCPi0DecayPhotonHitOtherLM","NLM vs E for merged pi0 cluster, decay photon hit Other Local Maxima",
nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
fhMCPi0DecayPhotonHitOtherLM ->SetYTitle("N maxima");
fhMCPi0DecayPhotonHitOtherLM ->SetXTitle("E (GeV)");
outputContainer->Add(fhMCPi0DecayPhotonHitOtherLM ) ;
- fhMCPi0DecayPhotonAdjOtherLM = new TH2F("hMCPi0DecayPhotonAdjOtherLM ","NLM vs E for merged pi0 cluster, decay photon hit cells adjacent to Other Local Maxima",
+ fhMCPi0DecayPhotonAdjOtherLM = new TH2F("hMCPi0DecayPhotonAdjOtherLM","NLM vs E for merged pi0 cluster, decay photon hit cells adjacent to Other Local Maxima",
nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
fhMCPi0DecayPhotonAdjOtherLM ->SetYTitle("N maxima");
fhMCPi0DecayPhotonAdjOtherLM ->SetXTitle("E (GeV)");
outputContainer->Add(fhMCPi0DecayPhotonAdjOtherLM ) ;
- fhMCPi0DecayPhotonHitNoLM = new TH2F("hMCPi0DecayPhotonHitNoLM ","NLM vs E for merged pi0 cluster, decay photon do not hit Local Maxima",
- nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0DecayPhotonAdjacent = new TH2F("hMCPi0DecayPhotonAdjacent","NLM vs E for merged pi0 cluster, decay photon hit adjacent cells",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0DecayPhotonAdjacent ->SetYTitle("N maxima");
+ fhMCPi0DecayPhotonAdjacent ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0DecayPhotonAdjacent ) ;
+
+ fhMCPi0DecayPhotonHitNoLM = new TH2F("hMCPi0DecayPhotonHitNoLM","NLM vs E for merged pi0 cluster, decay photon do not hit Local Maxima",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
fhMCPi0DecayPhotonHitNoLM ->SetYTitle("N maxima");
fhMCPi0DecayPhotonHitNoLM ->SetXTitle("E (GeV)");
outputContainer->Add(fhMCPi0DecayPhotonHitNoLM ) ;
+
+ fhMCPi0DecayPhotonHitHighLMOverlap = new TH2F("hMCPi0DecayPhotonHitHighLMOverlap","NLM vs E for merged pi0 cluster, decay photon hit High Local Maxima, there was an overlap",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0DecayPhotonHitHighLMOverlap ->SetYTitle("N maxima");
+ fhMCPi0DecayPhotonHitHighLMOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0DecayPhotonHitHighLMOverlap ) ;
+
+ fhMCPi0DecayPhotonAdjHighLMOverlap = new TH2F("hMCPi0DecayPhotonAdjHighLMOverlap","NLM vs E for merged pi0 cluster, decay photon hit cells adjacent to High Local Maxima, there was an overlap",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0DecayPhotonAdjHighLMOverlap ->SetYTitle("N maxima");
+ fhMCPi0DecayPhotonAdjHighLMOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0DecayPhotonAdjHighLMOverlap ) ;
+
+ fhMCPi0DecayPhotonHitOtherLMOverlap = new TH2F("hMCPi0DecayPhotonHitOtherLMOverlap","NLM vs E for merged pi0 cluster, decay photon hit Other Local Maxima, there was an overlap",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0DecayPhotonHitOtherLMOverlap ->SetYTitle("N maxima");
+ fhMCPi0DecayPhotonHitOtherLMOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0DecayPhotonHitOtherLMOverlap ) ;
+
+ fhMCPi0DecayPhotonAdjOtherLMOverlap = new TH2F("hMCPi0DecayPhotonAdjOtherLMOverlap","NLM vs E for merged pi0 cluster, decay photon hit cells adjacent to Other Local Maxima, there was an overlap",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0DecayPhotonAdjOtherLMOverlap ->SetYTitle("N maxima");
+ fhMCPi0DecayPhotonAdjOtherLMOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0DecayPhotonAdjOtherLMOverlap ) ;
+
+ fhMCPi0DecayPhotonAdjacentOverlap = new TH2F("hMCPi0DecayPhotonAdjacentOverlap","NLM vs E for merged pi0 cluster, decay photon hit adjacent cells, there was an overlap",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0DecayPhotonAdjacentOverlap ->SetYTitle("N maxima");
+ fhMCPi0DecayPhotonAdjacentOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0DecayPhotonAdjacentOverlap ) ;
+
+ fhMCPi0DecayPhotonHitNoLMOverlap = new TH2F("hMCPi0DecayPhotonHitNoLMOverlap","NLM vs E for merged pi0 cluster, decay photon do not hit Local Maxima, there was an overlap",
+ nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
+ fhMCPi0DecayPhotonHitNoLMOverlap ->SetYTitle("N maxima");
+ fhMCPi0DecayPhotonHitNoLMOverlap ->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPi0DecayPhotonHitNoLMOverlap ) ;
+
+
fhMCEOverlapType = new TH2F("hMCEOverlapType","Kind of overlap particle, neutral clusters",
nptbins,ptmin,ptmax,5,0,5);
//fhMCEOverlapType ->SetYTitle("Overlap Type");
fhMCEOverlapType->GetYaxis()->SetBinLabel(3 ,"hadron^{#pm}");
fhMCEOverlapType->GetYaxis()->SetBinLabel(4 ,"hadron^{0}");
fhMCEOverlapType->GetYaxis()->SetBinLabel(5 ,"??");
- fhMCEOverlapType ->SetXTitle("Cluster E (GeV)");
+ fhMCEOverlapType->SetXTitle("Cluster E (GeV)");
outputContainer->Add(fhMCEOverlapType) ;
fhMCEOverlapTypeMatch = new TH2F("hMCEOverlapTypeMatched","Kind of overlap particle, charged clusters",
TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
eprim = primary.E();
+ Int_t mesonLabel = -1;
+
if(mcindex == kmcPi0 || mcindex == kmcEta || mcindex == kmcPi0Conv)
{
if(mcindex == kmcPi0 || mcindex == kmcPi0Conv)
{
asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,111,GetReader(),ok));
- TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok);
+ TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
if(grandmom.E() > 0 && ok) eprim = grandmom.E();
}
else
{
asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,221,GetReader(),ok));
- TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok);
+ TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
if(grandmom.E() > 0 && ok) eprim = grandmom.E();
}
}
TLorentzVector momentum; TVector3 prodVertex;
Int_t ancLabel = 0;
noverlaps = 0;
+
for (UInt_t ilab = 1; ilab < cluster->GetNLabels(); ilab++ )
{
ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab],
GetReader(),ancPDG,ancStatus,momentum,prodVertex);
- if(ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG!=111 && ancPDG!=221)
+ //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
+ // ilab,cluster->GetLabels()[0],cluster->GetLabels()[ilab],ancLabel,ancPDG, mcindex);
+
+ Bool_t overlap = kFALSE;
+
+ //if(mcindex==kmcPi0 || mcindex==kmcPi0Conv || mcindex == kmcEta) printf("\t Meson MC : Label %d\n",mesonLabel);
+
+ if ( ancLabel < 0 )
{
- noverlaps++;
-
- // What is the origin of the overlap?
- Bool_t mOK = 0, gOK = 0;
- Int_t mpdg = -999999, gpdg = -1;
- Int_t mstatus = -1, gstatus = -1;
- Int_t gLabel = -1, ggLabel = -1;
- TLorentzVector mother = GetMCAnalysisUtils()->GetMother (cluster->GetLabels()[ilab],GetReader(),mpdg,mstatus,mOK);
- TLorentzVector grandmother = GetMCAnalysisUtils()->GetGrandMother(cluster->GetLabels()[ilab],GetReader(),gpdg,gstatus,gOK, gLabel,ggLabel);
-
- //printf("Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
-
- if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
- ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
- gLabel >=0 )
- {
- Int_t label = gLabel;
- while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
- {
- mpdg=gpdg;
- grandmother = GetMCAnalysisUtils()->GetGrandMother(label,GetReader(),gpdg,gstatus,ok, gLabel,ggLabel);
- label=gLabel;
- }
- }
-
- //printf("; Final PDG %d\n",mpdg);
- Float_t histobin = -1;
- if (mpdg==22) histobin = 0.5;
- else if(TMath::Abs(mpdg)==11) histobin = 1.5;
- else if(mpdg==-999999) histobin = 4.5;
- else {
- Double_t charge = TDatabasePDG::Instance()->GetParticle(mpdg)->Charge();
- if(TMath::Abs(charge) > 0 ) histobin = 2.5;
- else histobin = 3.5;
- //printf("charge %f\n",charge);
- }
-
- //printf("pdg = %d, histobin %2.1f\n",mpdg,histobin);
- if(histobin > 0)
+ overlap = kTRUE;
+ //printf("\t \t \t No Label = %d\n",ancLabel);
+ }
+ else if( ( ancPDG==111 || ancPDG==221 ) && (mcindex == kmcPi0 || mcindex == kmcPi0Conv || mcindex == kmcEta) && mesonLabel != ancLabel)
+ {
+ //printf("\t \t meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
+ overlap = kTRUE;
+ }
+ else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
+ {
+ //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
+ overlap = kTRUE ;
+ }
+
+ if( !overlap ) continue ;
+
+ // We have at least one overlap
+
+ //printf("Overlap!!!!!!!!!!!!!!\n");
+
+ noverlaps++;
+
+ // What is the origin of the overlap?
+ Bool_t mOK = 0, gOK = 0;
+ Int_t mpdg = -999999, gpdg = -1;
+ Int_t mstatus = -1, gstatus = -1;
+ Int_t gLabel = -1, ggLabel = -1;
+ TLorentzVector mother = GetMCAnalysisUtils()->GetMother (cluster->GetLabels()[ilab],GetReader(),mpdg,mstatus,mOK);
+ TLorentzVector grandmother = GetMCAnalysisUtils()->GetGrandMother(cluster->GetLabels()[ilab],GetReader(),gpdg,gstatus,gOK, gLabel,ggLabel);
+
+ //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
+
+ if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
+ ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
+ gLabel >=0 )
+ {
+ Int_t label = gLabel;
+ while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
{
- if(matched)fhMCEOverlapType ->Fill(cluster->E(),histobin);
- else fhMCEOverlapTypeMatch->Fill(cluster->E(),histobin);
+ mpdg=gpdg;
+ grandmother = GetMCAnalysisUtils()->GetGrandMother(label,GetReader(),gpdg,gstatus,ok, gLabel,ggLabel);
+ label=gLabel;
}
}
+
+ //printf("; Final PDG %d\n",mpdg);
+
+ Float_t histobin = -1;
+ if (mpdg==22) histobin = 0.5;
+ else if(TMath::Abs(mpdg)==11) histobin = 1.5;
+ else if(mpdg==-999999) histobin = 4.5;
+ else
+ {
+ Double_t charge = TDatabasePDG::Instance()->GetParticle(mpdg)->Charge();
+ if(TMath::Abs(charge) > 0 ) histobin = 2.5;
+ else histobin = 3.5;
+ //printf("charge %f\n",charge);
+ }
+
+ //printf("\t pdg = %d, histobin %2.1f\n",mpdg,histobin);
+ if(histobin > 0)
+ {
+ if(matched)fhMCEOverlapType ->Fill(cluster->E(),histobin);
+ else fhMCEOverlapTypeMatch->Fill(cluster->E(),histobin);
+ }
}
+
}
//___________________________________________
//
// For cluster with MC pi0 and more than 1 maxima
- CheckLocalMaximaMCOrigin(cluster, mcindex);
+ CheckLocalMaximaMCOrigin(cluster, mcindex,noverlaps);
//----------------
// Fill histograms
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff