new functionality and new class added
[u/mrichter/AliRoot.git] / PWG0 / dNdPt / AlidNdPtAnalysis.cxx
index 4297b2ecfad2b6f5e3f58121e8db61a3f085268a..08a1c15c20a0a2291ff47cd4effd5f65bc834a86 100644 (file)
@@ -39,6 +39,8 @@
 #include "AliMCEvent.h"  \r
 #include "AliESDtrackCuts.h"  \r
 #include "AliLog.h" \r
+#include "AliMultiplicity.h"\r
+#include "AliTracker.h"\r
 \r
 #include "AlidNdPtEventCuts.h"\r
 #include "AlidNdPtAcceptanceCuts.h"\r
@@ -136,6 +138,7 @@ ClassImp(AlidNdPtAnalysis)
   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) { \r
     fMCTrackHist1[i]=0;     \r
     fMCPrimTrackHist1[i]=0;     \r
+    fMCPrimTrackHist2[i]=0;     \r
     fMCSecTrackHist1[i]=0;     \r
     fRecTrackHist1[i]=0;     \r
     fRecTrackMultHist1[i]=0;     \r
@@ -229,6 +232,7 @@ AlidNdPtAnalysis::AlidNdPtAnalysis(Char_t* name, Char_t* title): AlidNdPt(name,t
   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) { \r
     fMCTrackHist1[i]=0;     \r
     fMCPrimTrackHist1[i]=0;     \r
+    fMCPrimTrackHist2[i]=0;     \r
     fMCSecTrackHist1[i]=0;     \r
     fRecTrackHist1[i]=0;     \r
     fRecTrackMultHist1[i]=0; \r
@@ -305,6 +309,7 @@ AlidNdPtAnalysis::~AlidNdPtAnalysis() {
   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) { \r
     if(fMCTrackHist1[i]) delete fMCTrackHist1[i]; fMCTrackHist1[i]=0;\r
     if(fMCPrimTrackHist1[i]) delete fMCPrimTrackHist1[i]; fMCPrimTrackHist1[i]=0;\r
+    if(fMCPrimTrackHist2[i]) delete fMCPrimTrackHist2[i]; fMCPrimTrackHist2[i]=0;\r
     if(fMCSecTrackHist1[i]) delete fMCSecTrackHist1[i]; fMCSecTrackHist1[i]=0;\r
     if(fRecTrackHist1[i]) delete fRecTrackHist1[i]; fRecTrackHist1[i]=0;\r
     if(fRecTrackMultHist1[i]) delete fRecTrackMultHist1[i]; fRecTrackMultHist1[i]=0;\r
@@ -344,7 +349,7 @@ void AlidNdPtAnalysis::Init(){
   Double_t maxMultTrueEventMatrix[2]={149.5,149.5}; \r
   fEventMultCorrelationMatrix = new THnSparseF("fEventMultCorrelationMatrix","mult:true_mult",2,binsMultTrueEventMatrix,minMultTrueEventMatrix,maxMultTrueEventMatrix);\r
   fEventMultCorrelationMatrix->GetAxis(0)->SetTitle("track multiplicity");\r
-  fEventMultCorrelationMatrix->GetAxis(1)->SetTitle("multiplicity");\r
+  fEventMultCorrelationMatrix->GetAxis(1)->SetTitle("true multiplicity");\r
   fEventMultCorrelationMatrix->Sumw2();\r
   \r
   Int_t binsTrackPtCorrelationMatrix[3]={ptNbins,ptNbins,etaNbins};\r
@@ -745,13 +750,13 @@ void AlidNdPtAnalysis::Init(){
   fMCTrackHist1[i]->GetAxis(2)->SetTitle("mcPhi (rad)");\r
   fMCTrackHist1[i]->Sumw2();\r
 \r
-  Int_t binsMCPrimTrackHist2[5]=  {ptNbins,etaNbins,6,20,4000};\r
-  Double_t minMCPrimTrackHist2[5]={0.,-1.,0.,0.,0.}; \r
-  Double_t maxMCPrimTrackHist2[5]={10.,1.,6.,20.,4000.}; \r
+  Int_t binsMCPrimTrackHist1[5]=  {ptNbins,etaNbins,6,20,4000};\r
+  Double_t minMCPrimTrackHist1[5]={0.,-1.,0.,0.,0.}; \r
+  Double_t maxMCPrimTrackHist1[5]={10.,1.,6.,20.,4000.}; \r
   sprintf(name,"fMCPrimTrackHist1_%d",i);\r
   sprintf(title,"mcPt:mcEta:pid:mech:mother");\r
   \r
-  fMCPrimTrackHist1[i] = new THnSparseF(name,title,5,binsMCPrimTrackHist2,minMCPrimTrackHist2,maxMCPrimTrackHist2);\r
+  fMCPrimTrackHist1[i] = new THnSparseF(name,title,5,binsMCPrimTrackHist1,minMCPrimTrackHist1,maxMCPrimTrackHist1);\r
   fMCPrimTrackHist1[i]->SetBinEdges(0,binsPt);\r
   fMCPrimTrackHist1[i]->SetBinEdges(1,binsEta);\r
   fMCPrimTrackHist1[i]->GetAxis(0)->SetTitle("mcPt (GeV/c)");\r
@@ -761,6 +766,18 @@ void AlidNdPtAnalysis::Init(){
   fMCPrimTrackHist1[i]->GetAxis(4)->SetTitle("mother");\r
   fMCPrimTrackHist1[i]->Sumw2();\r
 \r
+  Int_t binsMCPrimTrackHist2[5]=  {4000,20,4000};\r
+  Double_t minMCPrimTrackHist2[5]={0.,0.,0.}; \r
+  Double_t maxMCPrimTrackHist2[5]={4000.,20.,4000.}; \r
+  sprintf(name,"fMCPrimTrackHist2_%d",i);\r
+  sprintf(title,"pdg:mech:mother");\r
+  \r
+  fMCPrimTrackHist2[i] = new THnSparseF(name,title,5,binsMCPrimTrackHist2,minMCPrimTrackHist2,maxMCPrimTrackHist2);\r
+  fMCPrimTrackHist2[i]->GetAxis(0)->SetTitle("pdg");\r
+  fMCPrimTrackHist2[i]->GetAxis(1)->SetTitle("mech");\r
+  fMCPrimTrackHist2[i]->GetAxis(2)->SetTitle("mother");\r
+  fMCPrimTrackHist2[i]->Sumw2();\r
+\r
   Int_t binsMCSecTrackHist1[5]=  {ptNbins,etaNbins,6,20,4000};\r
   Double_t minMCSecTrackHist1[5]={0.,-1.,0.,0.,0.}; \r
   Double_t maxMCSecTrackHist1[5]={10.,1.,6.,20.,4000.}; \r
@@ -872,16 +889,16 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
     return;\r
   }\r
 \r
-  // get physics trigger selection \r
-  AliPhysicsSelection *trigSel = GetPhysicsTriggerSelection();\r
-  if(!trigSel) {\r
-    AliDebug(AliLog::kError, "cannot get trigSel");\r
-    return;\r
-  }\r
-\r
   // trigger selection\r
   Bool_t isEventTriggered = kTRUE;\r
-  if(evtCuts->IsTriggerRequired())  {\r
+  if(evtCuts->IsTriggerRequired())  \r
+  {\r
+    AliPhysicsSelection *trigSel = GetPhysicsTriggerSelection();\r
+    if(!trigSel) {\r
+      AliDebug(AliLog::kError, "cannot get trigSel");\r
+      return;\r
+    }\r
+\r
     if(IsUseMCInfo()) { \r
       //static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;\r
       //isEventTriggered = triggerAnalysis->IsTriggerFired(esdEvent, GetTrigger());\r
@@ -889,7 +906,6 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
       isEventTriggered = trigSel->IsCollisionCandidate(esdEvent);\r
     }\r
     else {\r
-      //isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());\r
       isEventTriggered = trigSel->IsCollisionCandidate(esdEvent);\r
     }\r
   }\r
@@ -968,11 +984,37 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
   {  \r
      multMBTracks = AlidNdPtHelper::GetTPCMBTrackMult(esdEvent,evtCuts,accCuts,esdTrackCuts);\r
   } \r
-  else if(GetAnalysisMode() == AlidNdPtHelper::kTPCSPDvtx || GetAnalysisMode()==AlidNdPtHelper::kTPCSPDvtxUpdate) \r
+  else if( GetAnalysisMode() == AlidNdPtHelper::kTPCITS ||  GetAnalysisMode() == AlidNdPtHelper::kTPCSPDvtx || \r
+           GetAnalysisMode()==AlidNdPtHelper::kTPCSPDvtxUpdate || GetAnalysisMode()==AlidNdPtHelper::kTPCITSHybrid ) \r
   {\r
-     //multMBTracks = AlidNdPtHelper::GetSPDMBTrackMult(esdEvent,0.0);\r
-     if(vtxESD->GetStatus())\r
-       multMBTracks = vtxESD->GetNContributors();\r
+     //if(vtxESD->GetStatus())\r
+     //  multMBTracks = vtxESD->GetNContributors();\r
+\r
+     // origin Jan Fiete GO\r
+     const AliMultiplicity* mult = esdEvent->GetMultiplicity();\r
+     if (mult) {\r
+       Int_t trackletCount = 0;\r
+       for(Int_t i=0; i<mult->GetNumberOfTracklets(); ++i) {\r
+         Float_t deltaPhi = mult->GetDeltaPhi(i);\r
+         // prevent values to be shifted by 2 Pi()\r
+         if (deltaPhi < -TMath::Pi())\r
+           deltaPhi += TMath::Pi() * 2;\r
+         if (deltaPhi > TMath::Pi())\r
+          deltaPhi -= TMath::Pi() * 2;\r
+\r
+         //if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)\r
+         // continue;\r
+\r
+       trackletCount++;\r
+       }\r
+       //multMBTracks = mult->GetNumberOfTracklets();\r
+       multMBTracks = trackletCount;\r
+       //printf("trackletCount %d \n", trackletCount);\r
+     }\r
+     else {\r
+       AliDebug(AliLog::kError, Form("Multiplicty %p", mult));\r
+       return; \r
+     }\r
   } \r
   else {\r
     AliDebug(AliLog::kError, Form("Found analysis type %s", GetAnalysisMode()));\r
@@ -983,17 +1025,20 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
   Int_t multAll=0, multAcc=0, multRec=0;\r
   Int_t *labelsAll=0, *labelsAcc=0, *labelsRec=0;\r
 \r
+  // cosmics analysis\r
+  Int_t cosmicCount = 0;\r
+  // high-pt tracks\r
+  Int_t highPtCount = 0;\r
+\r
   // check event cuts\r
   if(isEventOK && isEventTriggered)\r
   {\r
     // get all charged tracks\r
-    //allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,vtxESD,GetAnalysisMode());\r
     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());\r
     if(!allChargedTracks) return;\r
 \r
     Int_t entries = allChargedTracks->GetEntries();\r
     //printf("entries %d \n",entries);\r
-     Bool_t isCosmic = kFALSE;\r
 \r
     labelsAll = new Int_t[entries];\r
     labelsAcc = new Int_t[entries];\r
@@ -1016,6 +1061,29 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
       if(!esdTrackCuts->AcceptTrack(track))\r
        continue;\r
 \r
+      //\r
+      Bool_t isOK = kFALSE;\r
+      Double_t x[3]; track->GetXYZ(x);\r
+      Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
+\r
+      //\r
+      // if TPC-ITS hybrid tracking (kTPCITSHybrid)\r
+      // replace track parameters with TPC-ony track parameters\r
+      //\r
+      if( GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybrid ) \r
+      {\r
+        // Relate TPC-only tracks to SPD vertex\r
+        isOK = track->RelateToVertexTPCBxByBz(esdEvent->GetPrimaryVertexSPD(), b, kVeryBig);\r
+        if(!isOK) continue;\r
+\r
+       // replace esd track parameters with TPCinner\r
+        AliExternalTrackParam  *tpcTrack  = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
+       if (!tpcTrack) return;\r
+        track->Set(tpcTrack->GetX(),tpcTrack->GetAlpha(),tpcTrack->GetParameter(),tpcTrack->GetCovariance());\r
+\r
+        if(tpcTrack) delete tpcTrack; \r
+      } \r
+\r
       FillHistograms(track,stack,AlidNdPtHelper::kAllTracks); \r
       labelsAll[multAll] = TMath::Abs(track->GetLabel());\r
       multAll++;\r
@@ -1023,13 +1091,18 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
          //FillHistograms(track,stack,AlidNdPtHelper::kAccTracks); \r
         //labelsAcc[multAcc] = TMath::Abs(track->GetLabel());\r
         //multAcc++;\r
+      \r
+         // check high-pt tracks\r
+         if(accCuts->AcceptTrack(track) && track->Pt() > 6.)\r
+        {\r
+            //printf(" high pt: pt %f \n",track->Pt());\r
+            highPtCount++;\r
+        }\r
 \r
-         // cosmics analysis\r
-         isCosmic = kFALSE;\r
-         if( GetParticleMode()==AlidNdPtHelper::kCosmics )\r
+         // check cosmics tracks\r
+         if( GetParticleMode()==AlidNdPtHelper::kCosmic )\r
          {\r
-           Int_t mult = 0;\r
-           if(accCuts->AcceptTrack(track) ) \r
+           if(accCuts->AcceptTrack(track)) \r
           { \r
              for(Int_t j=0; j<entries;++j) \r
              {\r
@@ -1039,21 +1112,22 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
 \r
                if( esdTrackCuts->AcceptTrack(track1) && accCuts->AcceptTrack(track1) ) \r
                { \r
-                 mult++;\r
-                 isCosmic = AlidNdPtHelper::IsCosmicTrack(track,track1);\r
+                if ( AlidNdPtHelper::IsCosmicTrack(track,track1) ) { \r
+                  cosmicCount++;\r
+                  break;\r
+                }\r
               }\r
             }\r
           }\r
-          if(isCosmic) \r
-             printf("evt. number %d , mult %d \n", esdEvent->GetEventNumberInFile(), mult);\r
-\r
-           if(!isCosmic) continue;\r
+          // if(!isCosmic) continue;\r
          }\r
 \r
+         // update track parameters using vertex point \r
         if(GetAnalysisMode() == AlidNdPtHelper::kTPCSPDvtxUpdate) {\r
           // update track parameters\r
              AliExternalTrackParam cParam;\r
-            track->RelateToVertexTPC(esdEvent->GetPrimaryVertexSPD(),esdEvent->GetMagneticField(),kVeryBig,&cParam);\r
+             isOK = track->RelateToVertexTPCBxByBz(esdEvent->GetPrimaryVertexSPD(), b, kVeryBig, &cParam);\r
+             if(!isOK) continue;\r
             track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());\r
 \r
              if(accCuts->AcceptTrack(track)) {\r
@@ -1068,8 +1142,16 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
               labelsRec[multRec] = TMath::Abs(track->GetLabel());\r
               multRec++;\r
             }\r
-          }\r
-         }\r
+        }\r
+     }\r
+     if(cosmicCount) \r
+       printf("COSMIC EVENT: number %d , mult %d \n", esdEvent->GetEventNumberInFile(), multRec);\r
+\r
+     if(highPtCount) \r
+       printf("HIGH PT EVENT: number %d , mult %d \n", esdEvent->GetEventNumberInFile(), multRec);\r
+\r
+     if(multRec > 30) \r
+       printf("HIGH MULT EVENT: number %d , mult %d \n", esdEvent->GetEventNumberInFile(), multRec);\r
 \r
      // fill track multiplicity histograms\r
      FillHistograms(allChargedTracks,labelsAll,multAll,labelsAcc,multAcc,labelsRec,multRec);\r
@@ -1090,7 +1172,7 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
      // multiplicity correlation matrix\r
      //\r
      Double_t vMultTrueEventMatrix[2] = {multRec,multMCTrueTracks};\r
-     fEventMultCorrelationMatrix->Fill(vMultTrueEventMatrix);\r
+     if(isEventOK && isEventTriggered) fEventMultCorrelationMatrix->Fill(vMultTrueEventMatrix);\r
 \r
      // \r
      // event level corrections (zv,N_MB)\r
@@ -1247,7 +1329,7 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
 \r
          // only charged particles\r
          Double_t charge = particle->GetPDG()->Charge()/3.;\r
-         if (charge == 0.0)\r
+         if (TMath::Abs(charge) < 0.001)\r
          continue;\r
 \r
          // only postive charged \r
@@ -1265,7 +1347,8 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
          if(accCuts->AcceptTrack(particle)) \r
         {\r
 \r
-           if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) ) fGenPrimTrackMatrix->Fill(vTrackMatrix);\r
+           if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) ) \r
+            fGenPrimTrackMatrix->Fill(vTrackMatrix);\r
 \r
           // fill control histograms\r
            if(fHistogramsOn) \r
@@ -1298,7 +1381,10 @@ void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mc
              if(iMc == labelsRec[iRec]) \r
             {\r
                fRecTrackMatrix->Fill(vTrackMatrix);\r
-               if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) ) fRecPrimTrackMatrix->Fill(vTrackMatrix);\r
+\r
+               if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) ) \r
+                fRecPrimTrackMatrix->Fill(vTrackMatrix);\r
+\r
                if(!prim) fRecSecTrackMatrix->Fill(vTrackMatrix);\r
 \r
               // fill control histograms\r
@@ -1388,12 +1474,14 @@ void AlidNdPtAnalysis::FillHistograms(AliESDtrack *const esdTrack, AliStack *con
   Double_t values[3] = {pt,eta,phi};     \r
   fRecTrackHist1[trackObj]->Fill(values);\r
 \r
+  /*\r
   Double_t values1[5] = {nClust,chi2PerCluster,pt,eta,phi};      \r
   if(trackObj == AlidNdPtHelper::kRecTracks)  \r
   {\r
     if(fHistogramsOn)\r
       fRecTrackHist2->Fill(values1);\r
   }\r
+  */\r
  \r
   //\r
   // Fill rec vs MC information\r
@@ -1461,6 +1549,7 @@ void AlidNdPtAnalysis::FillHistograms(AliStack *const stack, Int_t label, AlidNd
 \r
   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3 \r
   if(TMath::Abs(gq) < 0.001) return;\r
+\r
   Float_t gpt = particle->Pt();\r
   //Float_t qgpt = particle->Pt() * gq;\r
   Float_t geta = particle->Eta();\r
@@ -1482,10 +1571,16 @@ void AlidNdPtAnalysis::FillHistograms(AliStack *const stack, Int_t label, AlidNd
   fMCTrackHist1[trackObj]->Fill(vMCTrackHist1);\r
 \r
   Double_t vMCPrimTrackHist1[5] = {gpt,geta,pid,mech,motherPdg};\r
-  if(prim) fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);\r
-  else     { \r
-         fMCSecTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);\r
+  Double_t vMCPrimTrackHist2[5] = {TMath::Abs(particle->GetPdgCode()),mech,motherPdg};\r
+  //if(prim && AliPWG0Helper::IsPrimaryCharged(particle, label)) fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);\r
+  if(prim) { \r
+    fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);\r
+    if(pid == 5) fMCPrimTrackHist2[trackObj]->Fill(vMCPrimTrackHist2);\r
   }\r
+  else { \r
+    fMCSecTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);\r
+  }\r
+\r
 }\r
 \r
 //_____________________________________________________________________________\r
@@ -1575,6 +1670,7 @@ Long64_t AlidNdPtAnalysis::Merge(TCollection* const list)
       fMCTrackHist1[i]->Add(entry->fMCTrackHist1[i]);\r
 \r
       fMCPrimTrackHist1[i]->Add(entry->fMCPrimTrackHist1[i]);\r
+      fMCPrimTrackHist2[i]->Add(entry->fMCPrimTrackHist2[i]);\r
       fMCSecTrackHist1[i]->Add(entry->fMCSecTrackHist1[i]);\r
 \r
       fRecTrackHist1[i]->Add(entry->fRecTrackHist1[i]);\r
@@ -2105,6 +2201,15 @@ void AlidNdPtAnalysis::Analyse()
   h2c->SetName("eff_pt_protons");\r
   aFolderObj->Add(h2c);\r
 \r
+  fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,5); \r
+  fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,5); \r
+  h1 = fMCPrimTrackHist1[1]->Projection(0);\r
+  h2 = fMCPrimTrackHist1[2]->Projection(0);\r
+  h2c = (TH1D *)h2->Clone();\r
+  h2c->Divide(h1);\r
+  h2c->SetName("eff_pt_selected");\r
+  aFolderObj->Add(h2c);\r
+\r
   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,6); \r
   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,6); \r
   h1 = fMCPrimTrackHist1[1]->Projection(0);\r