]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/Correlations/DPhi/AliDptDptInMC.cxx
Update from Prabhat - Dpt Dpt corr
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / AliDptDptInMC.cxx
index 8df44d39928389f005f897b6fc5a1b48c85dfa0c..a7d419a18d02940fdfbc0369c80504c8295270b4 100644 (file)
@@ -52,7 +52,7 @@
 #include "AliDptDptInMC.h"
 
 #include "AliPID.h"
-//#include "AliPIDResponse.h"
+#include "AliPIDResponse.h"
 
 #include "AliESDVertex.h"
 #include "AliESDEvent.h"
@@ -75,7 +75,7 @@
 #include "AliGenEventHeader.h"
 #include "AliPID.h"
 #include "AliAODPid.h"
-//#include "AliPIDResponse.h"
+#include "AliPIDResponse.h"
 #include "AliAODpidUtil.h"
 #include "AliPIDCombined.h"
 
@@ -89,7 +89,7 @@ AliDptDptInMC::AliDptDptInMC()
 fAODEvent(0), 
 fESDEvent(0),             //! ESD Event
 fInputHandler(0),
-//fPIDResponse(0x0),
+fPIDResponse(0x0),
 _outputHistoList(0),
 fArrayMC (0),
 _twoPi         ( 2.0 * 3.1415927),
@@ -384,7 +384,7 @@ AliDptDptInMC::AliDptDptInMC(const TString & name)
 fAODEvent(0), 
 fESDEvent(0),  
 fInputHandler(0),
-  //fPIDResponse(0),
+fPIDResponse(0),
 _outputHistoList(0),
 fArrayMC (0),
 _twoPi         ( 2.0 * 3.1415927),
@@ -781,6 +781,7 @@ void AliDptDptInMC::UserCreateOutputObjects()
   //_dedx_1     = new float[arraySize];
   
   __n1_1_vsPt              = getDoubleArray(_nBins_pt_1,        0.);
+
   __n1_1_vsEtaPhi          = getDoubleArray(_nBins_etaPhi_1,    0.);
   __s1pt_1_vsEtaPhi        = getDoubleArray(_nBins_etaPhi_1,    0.);
   __n1_1_vsZEtaPhiPt       = getFloatArray(_nBins_zEtaPhiPt_1,  0.);
@@ -1001,10 +1002,14 @@ void  AliDptDptInMC::createHistograms()
   name = "zV"; _vertexZ = createHisto1D(name,name,_nBins_vertexZ, _min_vertexZ, _max_vertexZ, "z-Vertex (cm)", _title_counts);
   
 
-  name = "Eta";     _etadis   = createHisto1F(name,name, 200, -1.0, 1.0, "#eta","counts");
-  name = "Phi";     _phidis   = createHisto1F(name,name, 360, 0.0, 6.4, "#phi","counts");
-  name = "DCAz";    _dcaz     = createHisto1F(name,name, 340, -3.3, 3.3, "dcaZ","counts");
-  name = "DCAxy";   _dcaxy    = createHisto1F(name,name, 100, -0.1, 2.5, "dcaXY","counts");
+  //name = "Eta";     _etadis   = createHisto1F(name,name, 200, -1.0, 1.0, "#eta","counts");
+  //name = "Phi";     _phidis   = createHisto1F(name,name, 360, 0.0, 6.4, "#phi","counts");
+  //name = "DCAz";    _dcaz     = createHisto1F(name,name, 340, -3.3, 3.3, "dcaZ","counts");
+  //name = "DCAxy";   _dcaxy    = createHisto1F(name,name, 100, -0.1, 2.5, "dcaXY","counts");
+
+  name = "Eta";     _etadis   = createHisto1F(name,name, 250, 0.0, 2.5, "#eta","counts"); //temporaryly
+  name = "Phi";     _phidis   = createHisto1F(name,name, 250, 0.0, 2.5, "#phi","counts");
+  name = "DCAz";    _dcaz     = createHisto1F(name,name, 250, 0.0, 2.5, "dcaZ","counts");
 
   //name = "Nclus1";   _Ncluster1    = createHisto1F(name,name, 200, 0, 200, "Ncluster1","counts");
   //name = "Nclus2";   _Ncluster2    = createHisto1F(name,name, 200, 0, 200, "Ncluster2","counts");
@@ -1012,6 +1017,7 @@ void  AliDptDptInMC::createHistograms()
   if (_singlesOnly)
     {
     name = n1Name+part_1_Name+vsPt;              _n1_1_vsPt              = createHisto1F(name,name, _nBins_pt_1,  _min_pt_1,  _max_pt_1,   _title_pt_1,  _title_AvgN_1);
+
     name = n1Name+part_1_Name+vsZ+vsEtaPhi+vsPt; _n1_1_vsZVsEtaVsPhiVsPt = createHisto3F(name,name, _nBins_vertexZ,_min_vertexZ,_max_vertexZ, _nBins_etaPhi_1, 0., double(_nBins_etaPhi_1), _nBins_pt_1, _min_pt_1, _max_pt_1, "zVertex", _title_etaPhi_1,  _title_pt_1);
     //name = "dedxVsP_1";                          _dedxVsP_1              = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
     //name = "corrDedxVsP_1";                      _corrDedxVsP_1          = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
@@ -1150,11 +1156,11 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
     {
       return;
     }
-  //  fPIDResponse =inputHandler->GetPIDResponse();
-  //if (!fPIDResponse){
-  //AliFatal("This Task needs the PID response attached to the inputHandler");
-  //return;
-  //}
+  fPIDResponse =inputHandler->GetPIDResponse();
+  if (!fPIDResponse){
+    AliFatal("This Task needs the PID response attached to the inputHandler");
+    return;
+  }
   // count all events looked at here                                                                                                        
   _eventCount++;
 
@@ -1259,7 +1265,7 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
                   continue;
                 }
              
-             //if(!aodTrack->IsPhysicalPrimary()) continue;
+             if(!aodTrack->IsPhysicalPrimary()) continue;
 
              q      = aodTrack->Charge();
               charge = int(q);
@@ -1272,23 +1278,47 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
              // Kinematics cuts from ESD track cuts                                                           
              if( pt < 0.2 || pt > 2.0)      continue;
              if( eta < -0.8 || eta > 0.8)  continue;
-             
-             // Remove neutral tracks                                                                                         
-              if(q == 0) continue;  
+
+             if(q == 0) continue;  
+             _dcaz->Fill(pt); //AllCh 
+
+             if(TMath::Abs(aodTrack->GetPdgCode()) == 211)
+               {
+                 _etadis->Fill(pt); //pion
+               } 
+             if(TMath::Abs(aodTrack->GetPdgCode()) == 321)
+               {
+                 _phidis->Fill(pt); //kaon
+               } 
+
+             //_etadis->Fill(eta);
+              //_phidis->Fill(phi); 
 
              if(fExcludeResonancesInMC)
                 {
-                 if (aodTrack->IsSecondaryFromWeakDecay()) continue;
-               }
+                  //cout<<"***************Prabhat on Weak Decay Particles ************"<<endl;                               
+                 Int_t gMotherIndex = aodTrack->GetMother();
+                  if(gMotherIndex != -1) {
+                    AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
+                    if(motherTrack) {
+                      Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
+                     
+                      if(pdgCodeOfMother == 311  ||
+                         pdgCodeOfMother == -311 ||
+                         pdgCodeOfMother == 310  ||
+                         pdgCodeOfMother == 3122 ||
+                         pdgCodeOfMother == -3122 ||
+                        pdgCodeOfMother == 111 ||
+                        pdgCodeOfMother == 22 ) continue;
+                    }
+                  }
+                }
              //Exclude electrons with PDG                                                                                   
                                                                                                                                
               if(fExcludeElectronsInMC) {
                 if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
               }
 
-             _etadis->Fill(eta);
-              _phidis->Fill(phi);
              //Particle loop 1st particle
               if (q > 0)
                 {
@@ -1331,6 +1361,7 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
                     {
 
                       __n1_1_vsPt[iPt]               += corr;          //cout << "step 15" << endl;                                            
+                      
                       __n1_1_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;                                               
 
                     }
@@ -1451,24 +1482,22 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
            AliError("ERROR: Could not retrieve MC event");
          }
 
-         TExMap *trackMap = new TExMap();//Mapping matrix----
-          //1st loop track for Global tracks                                                              
-                                                                                                           
-          for(Int_t i = 0; i < fAODEvent->GetNumberOfTracks(); i++)
-            {
-              AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(i));
-              if(!aodTrack) {
-                AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
-                continue;
-              }
-              Int_t gID = aodTrack->GetID();
-              if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks                       
-                                                                                                           
-            }
-
-          AliAODTrack* newAodTrack;
-
-         for (int iTrack=0; iTrack < fAODEvent->GetNumberOfTracks(); iTrack++)
+         TExMap *trackMap = new TExMap();//Mapping matrix---- 
+         //1st loop track for Global tracks                                                                                                    
+         for(Int_t i = 0; i < fAODEvent->GetNumberOfTracks(); i++)
+           {
+             AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(i));
+             if(!aodTrack) {
+               AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
+               continue;
+             }
+             Int_t gID = aodTrack->GetID();
+             if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks                                                             
+           }
+         
+         AliAODTrack* newAodTrack;
+         
+         for (int iTrack=0; iTrack< fAODEvent->GetNumberOfTracks(); iTrack++)
            {
              
              AliAODTrack *t = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(iTrack));
@@ -1480,6 +1509,8 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
              
              bitOK  = t->TestFilterBit(_trackFilterBit);
              if (!bitOK) continue;
+             Int_t gID = t->GetID();
+             newAodTrack = gID >= 0 ?t : fAODEvent->GetTrack(trackMap->GetValue(-1-gID));
              
              q      = t->Charge();
              charge = int(q);
@@ -1492,43 +1523,92 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
              //Float_t dcaXY = t->DCA();  
              //Float_t dcaZ  = t->ZAtDCA(); 
              
-             Int_t gID = t->GetID();
-              newAodTrack = gID >= 0 ?t : fAODEvent->GetTrack(trackMap->GetValue(-1-gID));
-
+             // get the electron nsigma                                                                                                
+             //Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kElectron));
+             Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kPion));
+              Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kKaon));
+              //Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kProton));
+             
+             //Make the decision based on the n-sigma of electrons exclusively 
+             /* if(nSigma < fNSigmaCut
+                 && nSigmaPions   > fNSigmaCut
+                 && nSigmaKaons   > fNSigmaCut
+                 && nSigmaProtons > fNSigmaCut ) continue;
+             */
+
+             if(q == 0) continue;            
              // Kinematics cuts                                                                                 
              if( pt < 0.2 || pt > 2.0)      continue;
              if( eta < _min_eta_1 || eta > _max_eta_1)  continue;
              
-             //W/Wo Secondaries
-             //if (!AODmcTrack->IsPhysicalPrimary()) continue;
+             _dcaz->Fill(pt);
+          
+             if (nSigmaPions < fNSigmaCut)
+               {
+                 _etadis->Fill(pt);
+               }
+             if (nSigmaKaons < fNSigmaCut)
+               {
+                 _phidis->Fill(pt);
+               }
 
-             
+             /*
              Double_t pos[3];
-              newAodTrack->GetXYZ(pos);
-
-              Double_t DCAX = pos[0] - vertexX;
-              Double_t DCAY = pos[1] - vertexY;
-              Double_t DCAZ = pos[2] - vertexZ;
+             newAodTrack->GetXYZ(pos);
+             
+             Double_t DCAX = pos[0] - vertexX;
+             Double_t DCAY = pos[1] - vertexY;
+             Double_t DCAZ = pos[2] - vertexZ;
+             
+             Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
+             
+             if (DCAZ     <  _dcaZMin ||
+                 DCAZ     >  _dcaZMax ||
+                 DCAXY    >  _dcaXYMax ) continue;
+             */
 
-              Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
 
-             Int_t label = TMath::Abs(t->GetLabel());
-             AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
-             if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;       
+      
+             //==== QA ===========================                                          
+             //_dcaz->Fill(DCAZ);                                                           
+             //_dcaxy->Fill(DCAXY);                                                         
+             //_etadis->Fill(eta);                                                          
+             //_phidis->Fill(phi); 
+             //===================================   
              
+             //W/Wo Secondaries
+             //if (!AODmcTrack->IsPhysicalPrimary()) continue;
+             
+             //cout<<"***************Prabhat on Weak Decay Particles ************"<<endl;
              if(fExcludeResonancesInMC)
                {
-                 if (AODmcTrack->IsSecondaryFromWeakDecay()) continue;
-                 //if (AODmcTrack->IsSecondaryFromMaterial()) continue;
+                 Int_t label = TMath::Abs(t->GetLabel());
+                 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
+                 
+                 Int_t gMotherIndex = AODmcTrack->GetMother();
+                 if(gMotherIndex != -1) {
+                   AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
+                   if(motherTrack) {
+                     Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
+                     
+                     if(pdgCodeOfMother == 311  ||
+                        pdgCodeOfMother == -311 ||
+                        pdgCodeOfMother == 310  ||
+                        pdgCodeOfMother == 3122 ||
+                        pdgCodeOfMother == -3122 ||
+                        pdgCodeOfMother == 111 ||
+                        pdgCodeOfMother == 22 ) continue;
+                   }
+                 }
                }
-                     
-             //==== QA ===========================                                          
-             _dcaz->Fill(DCAZ);                                                           
-             _dcaxy->Fill(DCAXY);
-             _etadis->Fill(eta);                                                          
-             _phidis->Fill(phi); 
-             //===================================   
-
+             
+             Int_t label = TMath::Abs(t->GetLabel());
+             AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
+             if (AODmcTrack)
+               {
+                 if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
+               }
+             
              //Particle 1                                                                                                                  
              if (t->Charge() > 0)
                {
@@ -1571,6 +1651,7 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
                    {
 
                      __n1_1_vsPt[iPt]               += corr;          //cout << "step 15" << endl;                                           
+
                      __n1_1_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;                                              
 
                    }
@@ -1684,13 +1765,13 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
     } //AOD events             
 
 
-  _m0->Fill(_mult0);
-  _m1->Fill(_mult1);
-  _m2->Fill(_mult2);
-  _m3->Fill(_mult3);
-  _m4->Fill(_mult4);
-  _m5->Fill(_mult5);
-  _m6->Fill(_mult6);
+  //  _m0->Fill(_mult0);
+  //_m1->Fill(_mult1);
+  //_m2->Fill(_mult2);
+  //_m3->Fill(_mult3);
+  //_m4->Fill(_mult4);
+  //_m5->Fill(_mult5);
+  //_m6->Fill(_mult6);
   _vertexZ->Fill(vertexZ);
 
   if (_singlesOnly)