from Carlos Perez:

adds more MC QA and PID for Lambda

diff --git a/PWG/FLOW/Tasks/AliAnalysisTaskFlowStrange.cxx b/PWG/FLOW/Tasks/AliAnalysisTaskFlowStrange.cxx
index 8230dcf0a4953ba5954e64a38a98811c9a62c8a5..3817c6d5ab99d8fc9fa9e32b11ce89c838c2ddaf 100644 (file)
--- a/PWG/FLOW/Tasks/AliAnalysisTaskFlowStrange.cxx
+++ b/PWG/FLOW/Tasks/AliAnalysisTaskFlowStrange.cxx
@@ -457,10 +457,21 @@ void AliAnalysisTaskFlowStrange::AddQACandidates() {
   if(fReadMC) {\r
     tList=new TList(); tList->SetName("RecMth"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);\r
     tH1D = new TH1D("MCOrigin", "MCOrigin;Rad2",1000,0,100); tList->Add(tH1D);\r
+    tH2D = new TH2D("PTRes", "PTRes;MCPt;|DAT-MC|/MC",100,0,20,50,0,1); tList->Add(tH2D);\r
+    tH2D = new TH2D("ETARes","ETARes;MCETA;|DAT-MC|/MC",16,0,0.8,50,0,1); tList->Add(tH2D);\r
+    tH2D = new TH2D("RXYRes","RXYRes;MCRXY;|DAT-MC|/MC",100,0,20,50,0,1); tList->Add(tH2D);\r
+    tH2D = new TH2D("DLERes","DLERes;MCDLE;|DAT-MC|/MC",100,0,20,50,0,1); tList->Add(tH2D);\r
+    tH2D = new TH2D("RAPRes","RAPRes;MCRAP;|DAT-MC|/MC",10,0,0.5,50,0,1); tList->Add(tH2D);\r
+    tH2D = new TH2D("APARes","APARes;MCAPA;|DAT-MC|/MC",24,-1.2,1.2,50,0,1); tList->Add(tH2D);\r
+    tH2D = new TH2D("APQRes","APQRes;MCAPQ;|DAT-MC|/MC",25,0,0.25,50,0,1); tList->Add(tH2D);\r
+\r
     tList=new TList(); tList->SetName("TrkMth"); tList->SetOwner(); AddTracksSpy(tList); fList->Add(tList);\r
+    tList=new TList(); tList->SetName("MthFDW"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);\r
+    tH1D = new TH1D("MCOrigin", "MCOrigin;Rad2",1000,0,100); tList->Add(tH1D);\r
   }\r
   //stack\r
   if(fReadMC) {\r
+    tList=new TList(); tList->SetName("GenTru"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);\r
     tList=new TList(); tList->SetName("MCTK0sGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);\r
     tList=new TList(); tList->SetName("MCTLdaGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);\r
     tList=new TList(); tList->SetName("MCTPhiGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);\r
@@ -852,12 +863,57 @@ void AliAnalysisTaskFlowStrange::ReadFromESD(AliESDEvent *tESD) {
 //=======================================================================\r
 void AliAnalysisTaskFlowStrange::ReadStack(TClonesArray* mcArray) {\r
   if(!mcArray) return;\r
-  AliAODMCParticle *myMCTrack;\r
+  AliAODMCParticle *myMCTrack, *iMCDau, *jMCDau;\r
   for(int i=0; i!=mcArray->GetEntriesFast(); ++i) {\r
     myMCTrack = dynamic_cast<AliAODMCParticle*>(mcArray->At( i ));\r
     if(!myMCTrack) continue;\r
+    int tPDG=310;\r
+    if(fSpecie>0) tPDG = 3122;\r
+    if( TMath::Abs(myMCTrack->PdgCode())==tPDG )\r
+      if( myMCTrack->GetNDaughters() == 2 ) {\r
+        Int_t iDau = myMCTrack->GetDaughter(0);\r
+        Int_t jDau = myMCTrack->GetDaughter(1);\r
+        AliAODMCParticle *posDau=NULL;\r
+        AliAODMCParticle *negDau=NULL;\r
+        if(iDau>0&&jDau>0) {\r
+          iMCDau = dynamic_cast<AliAODMCParticle*>(mcArray->At( iDau ));\r
+          jMCDau = dynamic_cast<AliAODMCParticle*>(mcArray->At( jDau ));\r
+          if(iMCDau) {\r
+            if(iMCDau->Charge()>0) posDau=iMCDau;\r
+            else negDau=iMCDau;\r
+          }\r
+          if(jMCDau) {\r
+            if(jMCDau->Charge()>0) posDau=jMCDau;\r
+            else negDau=jMCDau;\r
+          }\r
+        } //got two daughters\r
+        if(posDau&&negDau) {\r
+          Double_t dx = myMCTrack->Xv() - posDau->Xv();\r
+          Double_t dy = myMCTrack->Yv() - posDau->Yv();\r
+          Double_t dz = myMCTrack->Zv() - posDau->Zv();\r
+          fDecayRadXY = TMath::Sqrt( dx*dx + dy*dy );\r
+          TVector3 momPos(posDau->Px(),posDau->Py(),posDau->Pz());\r
+          TVector3 momNeg(negDau->Px(),negDau->Py(),negDau->Pz());\r
+          TVector3 momTot(myMCTrack->Px(),myMCTrack->Py(),myMCTrack->Pz());\r
+          Double_t qlpos = momPos.Dot(momTot)/momTot.Mag();\r
+          Double_t qlneg = momNeg.Dot(momTot)/momTot.Mag();\r
+          fDecayQt = momPos.Perp(momTot);\r
+          fDecayAlpha = 1.-2./(1.+qlpos/qlneg);\r
+          fDecayMass = myMCTrack->GetCalcMass();\r
+          Double_t energy = myMCTrack->E();\r
+          Double_t gamma = energy/fDecayMass;\r
+          fDecayDecayLength = TMath::Sqrt(dx*dx+dy*dy+dz*dz)/gamma;\r
+          fDecayPt = myMCTrack->Pt();\r
+          fDecayPhi = myMCTrack->Phi();\r
+          fDecayEta = myMCTrack->Eta();\r
+          fDecayRapidity = myMCTrack->Y();\r
+          fDecayDCAdaughters = 0;\r
+          fDecayCosinePointingAngleXY = 1;\r
+          fDecayProductIPXY = -1;\r
+          if(AcceptCandidate()) FillCandidateSpy("GenTru");\r
+        }\r
+      } // k0/lda with two daughters\r
     //==== BEGIN TRACK CUTS\r
-    if(myMCTrack->Pt()<0.2)   continue;\r
     if(myMCTrack->Eta()<-0.8) continue;\r
     if(myMCTrack->Eta()>+0.8) continue;\r
     //==== END TRACK CUTS\r
@@ -865,22 +921,23 @@ void AliAnalysisTaskFlowStrange::ReadStack(TClonesArray* mcArray) {
     case (310): //k0s\r
       FillMCParticleSpy( "MCTK0s", myMCTrack );\r
       if( myMCTrack->IsPrimary() )\r
-       FillMCParticleSpy( "MCTK0sGenAcc", myMCTrack );\r
+        FillMCParticleSpy( "MCTK0sGenAcc", myMCTrack );\r
       break;\r
     case (3122): //lda\r
       FillMCParticleSpy( "MCTLda", myMCTrack );\r
       if( myMCTrack->IsPrimary() )\r
-       FillMCParticleSpy( "MCTLdaGenAcc", myMCTrack );\r
+        FillMCParticleSpy( "MCTLdaGenAcc", myMCTrack );\r
       break;\r
     case (333): //phi\r
       if( myMCTrack->IsPrimary() )\r
-       FillMCParticleSpy( "MCTPhiGenAcc", myMCTrack );\r
+        FillMCParticleSpy( "MCTPhiGenAcc", myMCTrack );\r
       break;\r
     case (3312): //xi\r
       if( myMCTrack->IsPrimary() )\r
-       FillMCParticleSpy( "MCTXiGenAcc", myMCTrack );\r
+        FillMCParticleSpy( "MCTXiGenAcc", myMCTrack );\r
       break;\r
     }\r
+    \r
   }\r
 }\r
 //=======================================================================\r
@@ -1042,14 +1099,15 @@ void AliAnalysisTaskFlowStrange::ReadFromAODv0(AliAODEvent *tAOD) {
     ((TH2D*)((TList*)fList->FindObject("RecAll"))->FindObject("D0PD0N"))->Fill( dD0P, dD0N );\r
     ((TH2D*)((TList*)fList->FindObject("RecAll"))->FindObject("XPOSXNEG"))->Fill( xa, xb );\r
     if(!AcceptCandidate()) continue;\r
-    // PID for lambda::proton\r
-    //if( fSpecie>0 ) {\r
-    //  if(fDecayAlpha>0) {\r
-    //    if( !PassesPIDCuts(&ieT) ) continue; //positive track\r
-    //  } else {\r
-    //    if( !PassesPIDCuts(&jeT) ) continue; //negative track\r
-    //  }\r
-    //}\r
+    //PID for lambda::proton\r
+    if( fSpecie>0 )\r
+      if(fDecayPt<1.2) {\r
+        if(fDecayAlpha>0) {\r
+          if( !PassesPIDCuts(&ieT) ) continue; //positive track\r
+        } else {\r
+          if( !PassesPIDCuts(&jeT) ) continue; //negative track\r
+        }\r
+      }\r
     if(fQAlevel>1) {\r
       if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("RecSelOP");\r
       else FillCandidateSpy("RecSelIP");\r
@@ -1073,29 +1131,80 @@ void AliAnalysisTaskFlowStrange::ReadFromAODv0(AliAODEvent *tAOD) {
     //===== BEGIN OF MCMATCH                                                                                                                         \r
     if(mcArray) {\r
       bool matched = false;\r
+      bool feeddown = false;\r
       Int_t labelpos = iT->GetLabel();\r
       Int_t labelneg = jT->GetLabel();\r
       Double_t rOri=-1;\r
+      Double_t mcPt=-1, mcEta=-1, mcRap=-1, mcRxy=-1, mcDle=-1, mcApa=-100, mcApq=-1;\r
+      Double_t resPt=-1, resEta=-1, resRap=-1, resRxy=-1, resDle=-1, resApa=-1, resApq=-1;\r
       if( labelpos>0 && labelneg>0 ) {\r
         AliAODMCParticle *mcpos = (AliAODMCParticle*) mcArray->At( labelpos );\r
         AliAODMCParticle *mcneg = (AliAODMCParticle*) mcArray->At( labelneg );\r
         Int_t pdgRecPos = mcpos->GetPdgCode();\r
         Int_t pdgRecNeg = mcneg->GetPdgCode();\r
-        if( pdgRecPos==211&&pdgRecNeg==-211 ) if(mcpos->GetMother()>0) {\r
-          if( mcpos->GetMother()==mcneg->GetMother() ) {\r
-            AliAODMCParticle *mcmot = (AliAODMCParticle*) mcArray->At( mcpos->GetMother() );\r
-            rOri = TMath::Sqrt( mcmot->Xv()*mcmot->Xv() + mcmot->Yv()*mcmot->Yv() );\r
-            if( TMath::Abs(mcmot->GetPdgCode())==310) {\r
-              if(mcmot->GetNDaughters()==2) matched=true;\r
-            }\r
+        int pospdg=211, negpdg=211;\r
+        int mompdg=310, fdwpdg=333;\r
+        if(fSpecie>0) {\r
+          mompdg=3122;\r
+          fdwpdg=3312;\r
+          if(fDecayAlpha) {\r
+            pospdg=2212; negpdg=211;\r
+          } else {\r
+            negpdg=2212; pospdg=211;\r
           }\r
         }\r
-      }\r
+        if( TMath::Abs(pdgRecPos)==pospdg&&TMath::Abs(pdgRecNeg)==negpdg )\r
+          if(mcpos->GetMother()>0)\r
+            if( mcpos->GetMother()==mcneg->GetMother() ) {\r
+              AliAODMCParticle *mcmot = (AliAODMCParticle*) mcArray->At( mcpos->GetMother() );\r
+              rOri = TMath::Sqrt( mcmot->Xv()*mcmot->Xv() + mcmot->Yv()*mcmot->Yv() );\r
+              mcPt = mcmot->Pt();\r
+              mcEta = TMath::Abs( mcmot->Eta() );\r
+              mcRap = TMath::Abs( mcmot->Y() );\r
+              if(!TMath::AreEqualAbs(mcPt,0,1e-6))  resPt = TMath::Abs(fDecayPt - mcPt) / mcPt;\r
+              if(!TMath::AreEqualAbs(mcEta,0,1e-6)) resEta = TMath::Abs(fDecayEta - mcEta) / mcEta;\r
+              if(!TMath::AreEqualAbs(mcRap,0,1e-6)) resRap = TMath::Abs(fDecayRapidity - mcRap) / mcRap;\r
+              if( TMath::Abs(mcmot->GetPdgCode())==mompdg) {\r
+                if(mcmot->GetNDaughters()==2) {\r
+                  matched=true;\r
+                  Double_t dx = mcmot->Xv() - mcpos->Xv();\r
+                  Double_t dy = mcmot->Yv() - mcpos->Yv();\r
+                  Double_t dz = mcmot->Zv() - mcpos->Zv();\r
+                  Double_t mcGamma = mcmot->E() / mcmot->GetCalcMass();\r
+                  mcRxy = TMath::Sqrt( dx*dx + dy*dy );\r
+                  mcDle = TMath::Sqrt(dx*dx+dy*dy+dz*dz)/mcGamma;\r
+                  if(!TMath::AreEqualAbs(mcRxy,0,1e-6)) resRxy = TMath::Abs(fDecayRadXY - mcRxy) / mcRxy;\r
+                  if(!TMath::AreEqualAbs(mcDle,0,1e-6)) resDle = TMath::Abs(fDecayDecayLength - mcDle) / mcDle;\r
+                  TVector3 momPos(mcpos->Px(),mcpos->Py(),mcpos->Pz());\r
+                  TVector3 momNeg(mcneg->Px(),mcneg->Py(),mcneg->Pz());\r
+                  TVector3 momTot(mcmot->Px(),mcmot->Py(),mcmot->Pz());\r
+                  Double_t qlpos = momPos.Dot(momTot)/momTot.Mag();\r
+                  Double_t qlneg = momNeg.Dot(momTot)/momTot.Mag();\r
+                  mcApq = momPos.Perp(momTot);\r
+                  mcApa = 1.-2./(1.+qlpos/qlneg);\r
+                  if(!TMath::AreEqualAbs(mcApq,0,1e-6)) resApq = TMath::Abs(fDecayQt - mcApq) / mcApq;\r
+                  if(!TMath::AreEqualAbs(mcApa,0,1e-6)) resApa = TMath::Abs(fDecayAlpha - mcApa) / mcApa;\r
+                }\r
+                if(mcmot->GetMother()>0) {\r
+                  AliAODMCParticle *mcfdw = (AliAODMCParticle*) mcArray->At( mcmot->GetMother() );\r
+                  if( TMath::Abs(mcfdw->GetPdgCode())==fdwpdg)\r
+                    feeddown=true;\r
+                } // k0/lda have mother\r
+              } // mother matches k0/lda\r
+            } // both have same mother\r
+      } // match to MCparticles\r
       if(matched) {\r
         FillCandidateSpy("RecMth");\r
         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("D0PD0N"))->Fill( dD0P,dD0N );\r
         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("XPOSXNEG"))->Fill( xa, xb );\r
         ((TH1D*)((TList*)fList->FindObject("RecMth"))->FindObject("MCOrigin"))->Fill( rOri );\r
+        ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("PTRes"))->Fill( mcPt, resPt );\r
+        ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("ETARes"))->Fill( mcEta, resEta );\r
+        ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("RAPRes"))->Fill( mcRap, resRap );\r
+        ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("RXYRes"))->Fill( mcRxy, resRxy );\r
+        ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("DLERes"))->Fill( mcDle, resDle );\r
+        ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("APARes"))->Fill( mcApa, resApa );\r
+        ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("APQRes"))->Fill( mcApq, resApq );\r
         LoadTrack(&ieT,iT->Chi2perNDF());\r
         ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);\r
         fDaughterImpactParameterXY = ip[0];\r
@@ -1107,6 +1216,11 @@ void AliAnalysisTaskFlowStrange::ReadFromAODv0(AliAODEvent *tAOD) {
         fDaughterImpactParameterZ = ip[1];\r
         FillTrackSpy("TrkMth");\r
       }\r
+      if(feeddown) {\r
+        FillCandidateSpy("MthFDW");\r
+        ((TH2D*)((TList*)fList->FindObject("MthFDW"))->FindObject("D0PD0N"))->Fill( dD0P,dD0N );\r
+        ((TH1D*)((TList*)fList->FindObject("MthFDW"))->FindObject("MCOrigin"))->Fill( rOri );\r
+      }\r
     }\r
     //===== END OF MCMATCH\r
   }\r