]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.cxx
bugfix in event class filling of AliTHn for Psi BF
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskEventMixingBF.cxx
index ecfd846fe696e16ffb88f4ce8c8f16786a7646f6..fa8d441d34d126295b2f59e57177a332d36b088e 100755 (executable)
@@ -387,12 +387,16 @@ void AliAnalysisTaskEventMixingBF::Terminate(Option_t *) {
 \r
 void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)\r
 {\r
+  // Main loop for event mixing\r
 \r
   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 \r
   AliMixInputEventHandler *mixIEH = SetupEventsForMixing();\r
 \r
   Float_t fCentrality           = 0.;\r
+\r
+  // for HBT like cuts need magnetic field sign\r
+  Float_t bSign = 0; // only used in AOD so far\r
   \r
   // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
   vector<Double_t> *chargeVector[9];          // original charge\r
@@ -400,13 +404,13 @@ void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
     chargeVector[i]        = new vector<Double_t>;\r
   }\r
   \r
-  Double_t v_charge;\r
-  Double_t v_y;\r
-  Double_t v_eta;\r
-  Double_t v_phi;\r
-  Double_t v_p[3];\r
-  Double_t v_pt;\r
-  Double_t v_E;\r
+  Double_t vCharge;\r
+  Double_t vY;\r
+  Double_t vEta;\r
+  Double_t vPhi;\r
+  Double_t vP[3];\r
+  Double_t vPt;\r
+  Double_t vE;\r
 \r
   Int_t iMainTrackUsed = -1;\r
 \r
@@ -426,77 +430,81 @@ void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
        Printf("ERROR: aodEventMix not available");\r
        return;\r
       }\r
+\r
+     // for HBT like cuts need magnetic field sign\r
+     bSign = (aodEventMain->GetMagneticField() > 0) ? 1 : -1;\r
       \r
-     //AliAODHeader *aodHeaderMain = aodEventMain->GetHeader();\r
-     //AliAODHeader *aodHeaderMix  = aodEventMix->GetHeader();    \r
-  \r
+     AliAODHeader *aodHeaderMain = aodEventMain->GetHeader();  \r
 \r
       // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
       fHistEventStats->Fill(1); //all events\r
 \r
-      // Bool_t isSelectedMain = kTRUE;\r
-      // Bool_t isSelectedMix = kTRUE;\r
-\r
-      // if(fUseOfflineTrigger)\r
-      //       isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
-      // isSelectedMix = ((AliInputEventHandler*)((AliMultiInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetFirstMultiInputHandler())->IsEventSelected();\r
+      // this is not needed (checked in mixing handler!)\r
+      Bool_t isSelectedMain = kTRUE;\r
+      Bool_t isSelectedMix = kTRUE;\r
+      \r
+      if(fUseOfflineTrigger){\r
+               isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
+       isSelectedMix = ((AliInputEventHandler*)((AliMultiInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetFirstMultiInputHandler())->IsEventSelected();\r
+      }\r
       \r
-      // if(isSelectedMain && isSelectedMix) {\r
-      //       fHistEventStats->Fill(2); //triggered events\r
+      if(isSelectedMain && isSelectedMix) {\r
+               fHistEventStats->Fill(2); //triggered events\r
        \r
-      //       //Centrality stuff (centrality in AOD header)\r
-      //       if(fUseCentrality) {\r
-      //         fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
+       //Centrality stuff (centrality in AOD header)\r
+       if(fUseCentrality) {\r
+         fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
          \r
-      //         // QA for centrality estimators\r
-      //         fHistCentStats->Fill(0.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0M"));\r
-      //         fHistCentStats->Fill(1.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("FMD"));\r
-      //         fHistCentStats->Fill(2.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TRK"));\r
-      //         fHistCentStats->Fill(3.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKL"));\r
-      //         fHistCentStats->Fill(4.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL0"));\r
-      //         fHistCentStats->Fill(5.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL1"));\r
-      //         fHistCentStats->Fill(6.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
-      //         fHistCentStats->Fill(7.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
-      //         fHistCentStats->Fill(8.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
+         // QA for centrality estimators\r
+         fHistCentStats->Fill(0.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0M"));\r
+         fHistCentStats->Fill(1.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("FMD"));\r
+         fHistCentStats->Fill(2.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TRK"));\r
+         fHistCentStats->Fill(3.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKL"));\r
+         fHistCentStats->Fill(4.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL0"));\r
+         fHistCentStats->Fill(5.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL1"));\r
+         fHistCentStats->Fill(6.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
+         fHistCentStats->Fill(7.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
+         fHistCentStats->Fill(8.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
          \r
-      //         // take only events inside centrality class\r
-      //         if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
-      //           return;\r
+         // take only events inside centrality class\r
+         if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
+           return;\r
          \r
-      //         // centrality QA (V0M)\r
-      //         fHistV0M->Fill(aodEventMain->GetVZEROData()->GetMTotV0A(), aodEventMain->GetVZEROData()->GetMTotV0C());\r
+         // centrality QA (V0M)\r
+         fHistV0M->Fill(aodEventMain->GetVZEROData()->GetMTotV0A(), aodEventMain->GetVZEROData()->GetMTotV0C());\r
          \r
-      //         // centrality QA (reference tracks)\r
-      //         fHistRefTracks->Fill(0.,aodHeaderMain->GetRefMultiplicity());\r
-      //         fHistRefTracks->Fill(1.,aodHeaderMain->GetRefMultiplicityPos());\r
-      //         fHistRefTracks->Fill(2.,aodHeaderMain->GetRefMultiplicityNeg());\r
-      //         fHistRefTracks->Fill(3.,aodHeaderMain->GetTPConlyRefMultiplicity());\r
-      //         fHistRefTracks->Fill(4.,aodHeaderMain->GetNumberOfITSClusters(0));\r
-      //         fHistRefTracks->Fill(5.,aodHeaderMain->GetNumberOfITSClusters(1));\r
-      //         fHistRefTracks->Fill(6.,aodHeaderMain->GetNumberOfITSClusters(2));\r
-      //         fHistRefTracks->Fill(7.,aodHeaderMain->GetNumberOfITSClusters(3));\r
-      //         fHistRefTracks->Fill(8.,aodHeaderMain->GetNumberOfITSClusters(4));\r
-      //       }\r
+         // centrality QA (reference tracks)\r
+         fHistRefTracks->Fill(0.,aodHeaderMain->GetRefMultiplicity());\r
+         fHistRefTracks->Fill(1.,aodHeaderMain->GetRefMultiplicityPos());\r
+         fHistRefTracks->Fill(2.,aodHeaderMain->GetRefMultiplicityNeg());\r
+         fHistRefTracks->Fill(3.,aodHeaderMain->GetTPConlyRefMultiplicity());\r
+         fHistRefTracks->Fill(4.,aodHeaderMain->GetNumberOfITSClusters(0));\r
+         fHistRefTracks->Fill(5.,aodHeaderMain->GetNumberOfITSClusters(1));\r
+         fHistRefTracks->Fill(6.,aodHeaderMain->GetNumberOfITSClusters(2));\r
+         fHistRefTracks->Fill(7.,aodHeaderMain->GetNumberOfITSClusters(3));\r
+         fHistRefTracks->Fill(8.,aodHeaderMain->GetNumberOfITSClusters(4));\r
+       }\r
        \r
-      //       const AliAODVertex *vertexMain = aodEventMain->GetPrimaryVertex();\r
-      //       const AliAODVertex *vertexMix  = aodEventMix->GetPrimaryVertex();\r
-      //       \r
-      //       if(vertexMain && vertexMix) {\r
-      //         Double32_t fCovMain[6];\r
-      //         Double32_t fCovMix[6];\r
-      //         vertexMain->GetCovarianceMatrix(fCovMain);\r
-      //         vertexMix->GetCovarianceMatrix(fCovMix);\r
+       // // this is crashing (Bug in ROOT (to be solved)) but not needed (checked in mixing handler)\r
+       // const AliAODVertex *vertexMain = aodEventMain->GetPrimaryVertex();\r
+       // const AliAODVertex *vertexMix  = aodEventMix->GetPrimaryVertex();\r
+       \r
+       // if(vertexMain && vertexMix) {\r
+       //    Double32_t fCovMain[6];\r
+       //    Double32_t fCovMix[6];\r
+       //    vertexMain->GetCovarianceMatrix(fCovMain);\r
+       //    vertexMix->GetCovarianceMatrix(fCovMix);\r
          \r
-      //         if(vertexMain->GetNContributors() > 0 && vertexMix->GetNContributors() > 0) {\r
-      //           if(fCovMain[5] != 0 && fCovMix[5] != 0) {\r
-      //             fHistEventStats->Fill(3); //events with a proper vertex\r
-      //             if(TMath::Abs(vertexMain->GetX()) < fVxMax && TMath::Abs(vertexMix->GetX()) < fVxMax ) {\r
-      //               if(TMath::Abs(vertexMain->GetY()) < fVyMax && TMath::Abs(vertexMix->GetY()) < fVyMax) {\r
-      //                 if(TMath::Abs(vertexMain->GetZ()) < fVzMax && TMath::Abs(vertexMix->GetZ()) < fVzMax) {\r
-      //                   fHistEventStats->Fill(4); //analyzed events\r
-      //                   fHistVx->Fill(vertexMain->GetX());\r
-      //                   fHistVy->Fill(vertexMain->GetY());\r
-      //                   fHistVz->Fill(vertexMain->GetZ());\r
+       //    if(vertexMain->GetNContributors() > 0 && vertexMix->GetNContributors() > 0) {\r
+       //     if(fCovMain[5] != 0 && fCovMix[5] != 0) {\r
+       //       fHistEventStats->Fill(3); //events with a proper vertex\r
+       //       if(TMath::Abs(vertexMain->GetX()) < fVxMax && TMath::Abs(vertexMix->GetX()) < fVxMax ) {\r
+       //      if(TMath::Abs(vertexMain->GetY()) < fVyMax && TMath::Abs(vertexMix->GetY()) < fVyMax) {\r
+       //        if(TMath::Abs(vertexMain->GetZ()) < fVzMax && TMath::Abs(vertexMix->GetZ()) < fVzMax) {\r
+       //          fHistEventStats->Fill(4); //analyzed events\r
+       //          fHistVx->Fill(vertexMain->GetX());\r
+       //          fHistVy->Fill(vertexMain->GetY());\r
+       //          fHistVz->Fill(vertexMain->GetZ());\r
 \r
                    // Loop over tracks in main event\r
                    for (Int_t iTracksMain = 0; iTracksMain < aodEventMain->GetNumberOfTracks(); iTracksMain++) {\r
@@ -513,25 +521,25 @@ void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
                      fHistTrackStats->Fill(aodTrackMain->GetFilterMap());\r
                      if(!aodTrackMain->TestFilterBit(nAODtrackCutBit)) continue;\r
                      \r
-                     v_charge = aodTrackMain->Charge();\r
-                     v_y      = aodTrackMain->Y();\r
-                     v_eta    = aodTrackMain->Eta();\r
-                     v_phi    = aodTrackMain->Phi() * TMath::RadToDeg();\r
-                     v_E      = aodTrackMain->E();\r
-                     v_pt     = aodTrackMain->Pt();\r
-                     aodTrackMain->PxPyPz(v_p);\r
+                     vCharge = aodTrackMain->Charge();\r
+                     vY      = aodTrackMain->Y();\r
+                     vEta    = aodTrackMain->Eta();\r
+                     vPhi    = aodTrackMain->Phi() * TMath::RadToDeg();\r
+                     vE      = aodTrackMain->E();\r
+                     vPt     = aodTrackMain->Pt();\r
+                     aodTrackMain->PxPyPz(vP);\r
                      \r
-                     Float_t DCAxyMain = aodTrackMain->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
-                     Float_t DCAzMain  = aodTrackMain->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
+                     Float_t dcaXYMain = aodTrackMain->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
+                     Float_t dcaZMain  = aodTrackMain->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
                      \r
                      \r
                      // Kinematics cuts from ESD track cuts\r
-                     if( v_pt < fPtMin || v_pt > fPtMax)      continue;\r
-                     if( v_eta < fEtaMin || v_eta > fEtaMax)  continue;\r
+                     if( vPt < fPtMin || vPt > fPtMax)      continue;\r
+                     if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
                      \r
                      // Extra DCA cuts (for systematic studies [!= -1])\r
                      if( fDCAxyCut != -1 && fDCAzCut != -1){\r
-                       if(TMath::Sqrt((DCAxyMain*DCAxyMain)/(fDCAxyCut*fDCAxyCut)+(DCAzMain*DCAzMain)/(fDCAzCut*fDCAzCut)) > 1 ){\r
+                       if(TMath::Sqrt((dcaXYMain*dcaXYMain)/(fDCAxyCut*fDCAxyCut)+(dcaZMain*dcaZMain)/(fDCAzCut*fDCAzCut)) > 1 ){\r
                          continue;  // 2D cut\r
                        }\r
                      }\r
@@ -546,22 +554,22 @@ void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
                      \r
                      // fill QA histograms\r
                      fHistClus->Fill(aodTrackMain->GetITSNcls(),aodTrackMain->GetTPCNcls());\r
-                     fHistDCA->Fill(DCAzMain,DCAxyMain);\r
+                     fHistDCA->Fill(dcaZMain,dcaXYMain);\r
                      fHistChi2->Fill(aodTrackMain->Chi2perNDF());\r
-                     fHistPt->Fill(v_pt);\r
-                     fHistEta->Fill(v_eta);\r
-                     fHistPhi->Fill(v_phi);\r
+                     fHistPt->Fill(vPt);\r
+                     fHistEta->Fill(vEta);\r
+                     fHistPhi->Fill(vPhi);\r
                      \r
                      // fill charge vector\r
-                     chargeVector[0]->push_back(v_charge);\r
-                     chargeVector[1]->push_back(v_y);\r
-                     chargeVector[2]->push_back(v_eta);\r
-                     chargeVector[3]->push_back(v_phi);\r
-                     chargeVector[4]->push_back(v_p[0]);\r
-                     chargeVector[5]->push_back(v_p[1]);\r
-                     chargeVector[6]->push_back(v_p[2]);\r
-                     chargeVector[7]->push_back(v_pt);\r
-                     chargeVector[8]->push_back(v_E);\r
+                     chargeVector[0]->push_back(vCharge);\r
+                     chargeVector[1]->push_back(vY);\r
+                     chargeVector[2]->push_back(vEta);\r
+                     chargeVector[3]->push_back(vPhi);\r
+                     chargeVector[4]->push_back(vP[0]);\r
+                     chargeVector[5]->push_back(vP[1]);\r
+                     chargeVector[6]->push_back(vP[2]);\r
+                     chargeVector[7]->push_back(vPt);\r
+                     chargeVector[8]->push_back(vE);\r
 \r
                      // -------------------------------------------------------------               \r
                      // for each track in main event loop over all tracks in mix event\r
@@ -580,25 +588,25 @@ void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
                        fHistTrackStats->Fill(aodTrackMix->GetFilterMap());\r
                        if(!aodTrackMix->TestFilterBit(nAODtrackCutBit)) continue;\r
                        \r
-                       v_charge = aodTrackMix->Charge();\r
-                       v_y      = aodTrackMix->Y();\r
-                       v_eta    = aodTrackMix->Eta();\r
-                       v_phi    = aodTrackMix->Phi() * TMath::RadToDeg();\r
-                       v_E      = aodTrackMix->E();\r
-                       v_pt     = aodTrackMix->Pt();\r
-                       aodTrackMix->PxPyPz(v_p);\r
+                       vCharge = aodTrackMix->Charge();\r
+                       vY      = aodTrackMix->Y();\r
+                       vEta    = aodTrackMix->Eta();\r
+                       vPhi    = aodTrackMix->Phi() * TMath::RadToDeg();\r
+                       vE      = aodTrackMix->E();\r
+                       vPt     = aodTrackMix->Pt();\r
+                       aodTrackMix->PxPyPz(vP);\r
                      \r
-                       Float_t DCAxyMix = aodTrackMix->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
-                       Float_t DCAzMix  = aodTrackMix->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
+                       Float_t dcaXYMix = aodTrackMix->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
+                       Float_t dcaZMix  = aodTrackMix->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
                        \r
                        \r
                        // Kinematics cuts from ESD track cuts\r
-                       if( v_pt < fPtMin || v_pt > fPtMax)      continue;\r
-                       if( v_eta < fEtaMin || v_eta > fEtaMax)  continue;\r
+                       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
+                       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
                        \r
                        // Extra DCA cuts (for systematic studies [!= -1])\r
                        if( fDCAxyCut != -1 && fDCAxyCut != -1){\r
-                         if(TMath::Sqrt((DCAxyMix*DCAxyMix)/(fDCAxyCut*fDCAxyCut)+(DCAzMix*DCAzMix)/(fDCAzCut*fDCAzCut)) > 1 ){\r
+                         if(TMath::Sqrt((dcaXYMix*dcaXYMix)/(fDCAxyCut*fDCAxyCut)+(dcaZMix*dcaZMix)/(fDCAzCut*fDCAzCut)) > 1 ){\r
                            continue;  // 2D cut\r
                          }\r
                        }\r
@@ -613,22 +621,22 @@ void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
                        \r
                        // fill QA histograms\r
                        fHistClus->Fill(aodTrackMix->GetITSNcls(),aodTrackMix->GetTPCNcls());\r
-                       fHistDCA->Fill(DCAzMix,DCAxyMix);\r
+                       fHistDCA->Fill(dcaZMix,dcaXYMix);\r
                        fHistChi2->Fill(aodTrackMix->Chi2perNDF());\r
-                       fHistPt->Fill(v_pt);\r
-                       fHistEta->Fill(v_eta);\r
-                       fHistPhi->Fill(v_phi);\r
+                       fHistPt->Fill(vPt);\r
+                       fHistEta->Fill(vEta);\r
+                       fHistPhi->Fill(vPhi);\r
                        \r
                        // fill charge vector\r
-                       chargeVector[0]->push_back(v_charge);\r
-                       chargeVector[1]->push_back(v_y);\r
-                       chargeVector[2]->push_back(v_eta);\r
-                       chargeVector[3]->push_back(v_phi);\r
-                       chargeVector[4]->push_back(v_p[0]);\r
-                       chargeVector[5]->push_back(v_p[1]);\r
-                       chargeVector[6]->push_back(v_p[2]);\r
-                       chargeVector[7]->push_back(v_pt);\r
-                       chargeVector[8]->push_back(v_E);\r
+                       chargeVector[0]->push_back(vCharge);\r
+                       chargeVector[1]->push_back(vY);\r
+                       chargeVector[2]->push_back(vEta);\r
+                       chargeVector[3]->push_back(vPhi);\r
+                       chargeVector[4]->push_back(vP[0]);\r
+                       chargeVector[5]->push_back(vP[1]);\r
+                       chargeVector[6]->push_back(vP[2]);\r
+                       chargeVector[7]->push_back(vPt);\r
+                       chargeVector[8]->push_back(vE);\r
                        \r
                        \r
                        \r
@@ -637,7 +645,7 @@ void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
                      // calculate balance function for each track in main event\r
                      iMainTrackUsed++; // is needed to do no double counting in Balance Function calculation   \r
                      if(iMainTrackUsed >= (Int_t)chargeVector[0]->size()) break; //do not allow more tracks than in mixed event!\r
-                     fBalance->CalculateBalance(fCentrality,chargeVector,iMainTrackUsed);\r
+                     fBalance->CalculateBalance(fCentrality,chargeVector,iMainTrackUsed,bSign);\r
                      // clean charge vector afterwards\r
                      for(Int_t i = 0; i < 9; i++){                    \r
                        chargeVector[i]->clear();\r
@@ -648,15 +656,16 @@ void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
       //                 }//Vz cut\r
       //               }//Vy cut\r
       //             }//Vx cut\r
-      //           }//proper vertex resolution\r
-      //         }//proper number of contributors\r
-      //       }//vertex object valid\r
-      // }//triggered event \r
+      //           }//proper vertexresolution\r
+      //    }//proper number of contributors\r
+      // }//vertex object valid\r
+      }//triggered event \r
     }//AOD analysis\r
   }\r
 }\r
 \r
 AliMixInputEventHandler *AliAnalysisTaskEventMixingBF::SetupEventsForMixing() {\r
+  //sets the input handlers for event mixing\r
 \r
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
   AliMultiInputEventHandler *inEvHMain = dynamic_cast<AliMultiInputEventHandler *>(mgr->GetInputEventHandler());\r