]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
update to AOD analysis
authorlmilano <lmilano@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Apr 2012 16:09:17 +0000 (16:09 +0000)
committerlmilano <lmilano@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Apr 2012 16:09:17 +0000 (16:09 +0000)
PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraAOD.cxx
PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODHistoManager.cxx
PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODHistoManager.h
PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODTrackCuts.cxx
PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODTrackCuts.h
PWGLF/SPECTRA/PiKaPr/TestAOD/runAODProof.C

index 19c39031fc6e279693792212b57f878e8ecfdd0f..5f1dcfbc767eb30a0cd2f5a9525a75627691e07b 100644 (file)
 #include "AliPID.h"\r
 #include "AliVEvent.h"\r
 #include "AliPIDResponse.h"\r
+#include "AliStack.h"\r
 #include <iostream>\r
 \r
+\r
+\r
+\r
 using namespace AliSpectraNameSpace;\r
 using namespace std;\r
 \r
@@ -50,208 +54,311 @@ ClassImp(AliAnalysisTaskSpectraAOD) // EX1 This stuff tells root to implement th
 //________________________________________________________________________\r
 AliAnalysisTaskSpectraAOD::AliAnalysisTaskSpectraAOD(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0), fIsMC(0), fPIDResponse(0), fNSigmaPID(0), fYCut(0)\r
 {\r
-   // Default constructor\r
-   fNSigmaPID = 3.;\r
-   fYCut = .5;\r
-   DefineInput(0, TChain::Class());\r
-   DefineOutput(1, AliSpectraAODHistoManager::Class());\r
-   DefineOutput(2, AliSpectraAODEventCuts::Class());\r
-   DefineOutput(3, AliSpectraAODTrackCuts::Class());\r
+  // Default constructor\r
+  fNSigmaPID = 3.;\r
+  fYCut = .5;\r
+  DefineInput(0, TChain::Class());\r
+  DefineOutput(1, AliSpectraAODHistoManager::Class());\r
+  DefineOutput(2, AliSpectraAODEventCuts::Class());\r
+  DefineOutput(3, AliSpectraAODTrackCuts::Class());\r
 \r
 }\r
 //________________________________________________________________________\r
 Bool_t AliAnalysisTaskSpectraAOD::CheckYCut(AODParticleSpecies_t species, AliAODTrack* track) const\r
 {\r
-    // check if the rapidity is within the set range\r
-    // note: masses are hardcoded for now. we could look them up in the pdg database, but that would mean accecing it 100k+ times per run ...\r
-    Double_t y;\r
-    if (species == kProton) { y = track->Y(9.38271999999999995e-01); }\r
-    if ( species == kKaon ) { y = track->Y(4.93676999999999977e-01); }\r
-    if ( species == kPion)  { y = track->Y(1.39570000000000000e-01); }\r
-    if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;\r
-    return kTRUE;\r
+  // check if the rapidity is within the set range\r
+  // note: masses are hardcoded for now. we could look them up in the pdg database, but that would mean accecing it 100k+ times per run ...\r
+  Double_t y;\r
+  if (species == kProton) { y = track->Y(9.38271999999999995e-01); }\r
+  if ( species == kKaon ) { y = track->Y(4.93676999999999977e-01); }\r
+  if ( species == kPion)  { y = track->Y(1.39570000000000000e-01); }\r
+  if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;\r
+  return kTRUE;\r
 }\r
 //____________________________________________________________________________\r
 Bool_t AliAnalysisTaskSpectraAOD::CheckYCut(AliAODMCParticle* particle) const\r
 {\r
-    // check if the rapidity is within the set range\r
-    Double_t y = particle->Y();\r
-    if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;\r
-    return kTRUE;\r
+  // check if the rapidity is within the set range\r
+  Double_t y = particle->Y();\r
+  if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;\r
+  return kTRUE;\r
 }\r
 //________________________________________________________________________\r
 void AliAnalysisTaskSpectraAOD::UserCreateOutputObjects()\r
 {\r
   // create output objects\r
-    fHistMan = new AliSpectraAODHistoManager("SpectraHistos");\r
+  fHistMan = new AliSpectraAODHistoManager("SpectraHistos");\r
 \r
-   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
-   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
+  if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
+  if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
 \r
-   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();\r
-   AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());\r
-   fPIDResponse = inputHandler->GetPIDResponse();\r
+  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();\r
+  AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());\r
+  fPIDResponse = inputHandler->GetPIDResponse();\r
 \r
-   PostData(1, fHistMan);\r
-   PostData(2, fEventCuts);\r
-   PostData(3, fTrackCuts);\r
+  PostData(1, fHistMan);\r
+  PostData(2, fEventCuts);\r
+  PostData(3, fTrackCuts);\r
 \r
 }\r
 //________________________________________________________________________\r
 void AliAnalysisTaskSpectraAOD::UserExec(Option_t *)\r
 {\r
   // main event loop\r
-   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);\r
-   if (strcmp(fAOD->ClassName(), "AliAODEvent"))\r
-   {\r
+  fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);\r
+  if (strcmp(fAOD->ClassName(), "AliAODEvent"))\r
+    {\r
       AliFatal("Not processing AODs");\r
-   }\r
-\r
-   if(!fPIDResponse) {\r
-     AliError("Cannot get pid response");\r
-     return;\r
-   }\r
-\r
-   if (!fEventCuts->IsSelected(fAOD)) return;\r
-\r
-   //AliCentrality fAliCentral*;\r
-//   if ((fAOD->GetCentrality())->GetCentralityPercentile("V0M") > 5.) return;\r
-   \r
-   // First do MC to fill up the MC particle array, such that we can use it later\r
-   TClonesArray *arrayMC = 0;\r
-   if (fIsMC)\r
-   {\r
-      arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
-      if (!arrayMC) {\r
-         AliFatal("Error: MC particles branch not found!\n");\r
-      }\r
-      Int_t nMC = arrayMC->GetEntries();\r
-      for (Int_t iMC = 0; iMC < nMC; iMC++)\r
+    }\r
+  \r
+  if(!fPIDResponse) {\r
+    AliError("Cannot get pid response");\r
+    return;\r
+  }\r
+  \r
+  //Selection on QVector, before ANY other selection on the event\r
+  //Spectra MUST be normalized wrt events AFTER the selection on Qvector\r
+  // Can we include this in fEventCuts\r
+  Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;\r
+  Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;\r
+  Int_t multPos = 0;\r
+  Int_t multNeg = 0;\r
+  for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {\r
+    AliAODTrack* aodTrack = fAOD->GetTrack(iT);\r
+    if (!fTrackCuts->IsSelected(aodTrack)) continue;\r
+    if (aodTrack->Eta() >= 0){\r
+      multPos++;\r
+      Qx2EtaPos += TMath::Cos(2*aodTrack->Phi()); \r
+      Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());\r
+    } else {\r
+      multNeg++;\r
+      Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi()); \r
+      Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());\r
+    }\r
+  } \r
+  Double_t qPos = TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);\r
+  Double_t qNeg = TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);\r
+  if((qPos>fTrackCuts->GetQvecMin() && qPos<fTrackCuts->GetQvecMax()) || (qNeg>fTrackCuts->GetQvecMin() && qNeg<fTrackCuts->GetQvecMax())){\r
+    \r
+    //check on centrality distribution\r
+    fHistMan->GetPtHistogram("CentCheck")->Fill(fAOD->GetCentrality()->GetCentralityPercentile("V0M"),fAOD->GetHeader()->GetCentralityP()->GetCentralityPercentileUnchecked("V0M"));\r
+    \r
+    if (!fEventCuts->IsSelected(fAOD)) return;\r
+    \r
+    //fill q distributions vs centrality, after all event selection\r
+    fHistMan->GetqVecHistogram(kHistqVecPos)->Fill(qPos,fAOD->GetCentrality()->GetCentralityPercentile("V0M"));  // qVector distribution\r
+    fHistMan->GetqVecHistogram(kHistqVecNeg)->Fill(qNeg,fAOD->GetCentrality()->GetCentralityPercentile("V0M"));  // qVector distribution\r
+    \r
+    //AliCentrality fAliCentral*;\r
+    //   if ((fAOD->GetCentrality())->GetCentralityPercentile("V0M") > 5.) return;\r
+    \r
+    // First do MC to fill up the MC particle array, such that we can use it later\r
+    TClonesArray *arrayMC = 0;\r
+    if (fIsMC)\r
       {\r
-         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
-         fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt());\r
-         \r
-        if(TMath::Abs(partMC->Eta()) > fTrackCuts->GetEta()) continue;\r
-         // check for true PID + and fill P_t histos \r
-         if (partMC->IsPhysicalPrimary() && CheckYCut(partMC) ) {// only primary vertices and y cut satisfied\r
-            if ( partMC->PdgCode() == 2212) { fHistMan->GetHistogram(kHistPtGenTruePrimaryProtonPlus)->Fill(partMC->Pt()); } \r
-            if ( partMC->PdgCode() == -2212) { fHistMan->GetHistogram(kHistPtGenTruePrimaryProtonMinus)->Fill(partMC->Pt()); } \r
-            if ( partMC->PdgCode() == 321)  { fHistMan->GetHistogram(kHistPtGenTruePrimaryKaonPlus)->Fill(partMC->Pt()); } \r
-            if ( partMC->PdgCode() == -321) { fHistMan->GetHistogram(kHistPtGenTruePrimaryKaonMinus)->Fill(partMC->Pt()); } \r
-            if ( partMC->PdgCode() == 211)  { fHistMan->GetHistogram(kHistPtGenTruePrimaryPionPlus)->Fill(partMC->Pt()); } \r
-            if ( partMC->PdgCode() == -211) { fHistMan->GetHistogram(kHistPtGenTruePrimaryPionMinus)->Fill(partMC->Pt()); } \r
-         }\r
+       arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
+       if (!arrayMC) {\r
+         AliFatal("Error: MC particles branch not found!\n");\r
+       }\r
+       Int_t nMC = arrayMC->GetEntries();\r
+       for (Int_t iMC = 0; iMC < nMC; iMC++)\r
+         {\r
+           AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
+           fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
+           \r
+           if(TMath::Abs(partMC->Eta()) > fTrackCuts->GetEta()) continue;\r
+           // check for true PID + and fill P_t histos \r
+           //if (partMC->IsPhysicalPrimary() && CheckYCut(partMC) ) {// only primary vertices and y cut satisfied\r
+           if (CheckYCut(partMC) ){// only primary vertices and y cut satisfied\r
+             if ( partMC->PdgCode() == 2212) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryProtonPlus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
+             if ( partMC->PdgCode() == -2212) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryProtonMinus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
+             if ( partMC->PdgCode() == 321)  { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryKaonPlus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
+             if ( partMC->PdgCode() == -321) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryKaonMinus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
+             if ( partMC->PdgCode() == 211)  { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryPionPlus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
+             if ( partMC->PdgCode() == -211) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryPionMinus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
+           }\r
+         }\r
       }\r
-   }\r
-   \r
-   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)\r
-   {\r
-\r
-      AliAODTrack* track = fAOD->GetTrack(iTracks);\r
-      if (!fTrackCuts->IsSelected(track)) continue;\r
-\r
-      fHistMan->GetHistogram(kHistPtRec)->Fill(track->Pt());  // PT histo\r
-      fHistMan->GetPIDHistogram(kHistPIDTPC)->Fill(track->P(), track->GetTPCsignal()); // PID histo\r
-      fHistMan->GetPIDHistogram(kHistPIDTPCPt)->Fill(track->Pt(), track->GetTPCsignal());\r
-      \r
-      // MF 22/02/2012\r
-      // Fill DCA vs pt histos\r
-\r
-      // Sigma identify particles (for both MC and regular data):\r
-      // note: as of now, very crude identification. how many sigmas are allowed? what to do with high pt?\r
-      AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(track);\r
-      Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)); \r
-      Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
-      Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)); \r
-      if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton )) { \r
+    \r
+    \r
+    \r
+    //main loop on tracks\r
+    for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)\r
+      {\r
        \r
-          if ((nsigmaTPCkKaon > fNSigmaPID) || (!CheckYCut(kKaon, track) ) ) continue;\r
-          if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaKaonPlus)->Fill(track->Pt()); } \r
-          else { fHistMan->GetHistogram(kHistPtRecSigmaKaonMinus)->Fill(track->Pt()); } \r
-      }\r
-      if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
-          if ( nsigmaTPCkProton > fNSigmaPID || (!CheckYCut(kProton, track) ) )  continue;\r
-          if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaProtonPlus)->Fill(track->Pt()); }\r
-          else { fHistMan->GetHistogram(kHistPtRecSigmaProtonMinus)->Fill(track->Pt()); }\r
-      }\r
-      if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
-          if (nsigmaTPCkPion > fNSigmaPID || (!CheckYCut(kPion, track) ) ) continue;\r
-          if ( track->Charge() > 0 )  { fHistMan->GetHistogram(kHistPtRecSigmaPionPlus)->Fill(track->Pt()); }\r
-          else  { fHistMan->GetHistogram(kHistPtRecSigmaPionMinus)->Fill(track->Pt()); }\r
-      }\r
-      /* MC Part */\r
-      // MF 22/02/2012\r
-      // fill mc DCA vs pt\r
-      if (arrayMC) {\r
-         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
-         if (!partMC) { \r
-             AliError("Cannot get MC particle");\r
-             continue; }\r
-         if (CheckYCut(partMC)) {\r
-         // primaries, true pid\r
-         if ( partMC->PdgCode() == 2212) { fHistMan->GetHistogram(kHistPtRecTrueProtonPlus)->Fill(partMC->Pt()); \r
-            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryProtonPlus)->Fill(partMC->Pt()); }}\r
-         if ( partMC->PdgCode() == -2212) { fHistMan->GetHistogram(kHistPtRecTrueProtonMinus)->Fill(partMC->Pt()); \r
-            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryProtonMinus)->Fill(partMC->Pt()); }}\r
-         if ( partMC->PdgCode() == 321) { fHistMan->GetHistogram(kHistPtRecTrueKaonPlus)->Fill(partMC->Pt()); \r
-            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryKaonPlus)->Fill(partMC->Pt()); }}\r
-         if ( partMC->PdgCode() == -321) { fHistMan->GetHistogram(kHistPtRecTrueKaonMinus)->Fill(partMC->Pt()); \r
-            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryKaonMinus)->Fill(partMC->Pt()); }}\r
-         if ( partMC->PdgCode() == 211) { fHistMan->GetHistogram(kHistPtRecTruePionPlus)->Fill(partMC->Pt()); \r
-            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryPionPlus)->Fill(partMC->Pt()); }}\r
-         if ( partMC->PdgCode() == -211) { fHistMan->GetHistogram(kHistPtRecTruePionMinus)->Fill(partMC->Pt());  \r
-            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryPionMinus)->Fill(partMC->Pt()); }}\r
-         }\r
-\r
-         // primaries, sigma pid \r
-        // DCA vs pt not needed?\r
-         if (partMC->IsPhysicalPrimary()) { \r
-            if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) {\r
-               if ( (nsigmaTPCkKaon > fNSigmaPID ) || (!CheckYCut(kKaon, track) ) ) continue; \r
-               if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryKaonPlus)->Fill(track->Pt()); } \r
-               else { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryKaonMinus)->Fill(track->Pt()); } \r
-            }\r
-            if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
-               if ( (nsigmaTPCkProton > fNSigmaPID ) || (!CheckYCut(kProton, track) ) ) continue;\r
-               if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryProtonPlus)->Fill(track->Pt()); }\r
-               else { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryProtonMinus)->Fill(track->Pt()); }\r
-            }\r
-            if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
-               if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
-               if ( track->Charge() > 0 )  { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryPionPlus)->Fill(track->Pt()); }\r
-               else  { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryPionMinus)->Fill(track->Pt()); }\r
-            }\r
-         }\r
-         //secondaries\r
-        // MF 22/02/2012\r
-        // Distinguish weak decay and material\r
-         else {\r
-            if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) { \r
-               if ( (nsigmaTPCkKaon > fNSigmaPID )  || (!CheckYCut(kKaon, track) ) ) continue;\r
-               if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryKaonPlus)->Fill(track->Pt()); } \r
-               else { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryKaonMinus)->Fill(track->Pt()); } \r
-            }\r
-            if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
-               if ( (nsigmaTPCkProton > fNSigmaPID )  || (!CheckYCut(kProton, track) ) ) continue;\r
-               if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryProtonPlus)->Fill(track->Pt()); }\r
-               else { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryProtonMinus)->Fill(track->Pt()); }\r
-            }\r
-            if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
-               if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
-               if ( track->Charge() > 0 )  { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryPionPlus)->Fill(track->Pt()); }\r
-               else  { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryPionMinus)->Fill(track->Pt()); }\r
-            }\r
-         }\r
-      }\r
-   }\r
+       AliAODTrack* track = fAOD->GetTrack(iTracks);\r
+       if (!fTrackCuts->IsSelected(track)) continue;\r
+       \r
+       //cut on q vectors track-by-track\r
+       if ((qPos>fTrackCuts->GetQvecMin() && qPos<fTrackCuts->GetQvecMax() && track->Eta()<0) || (qNeg>fTrackCuts->GetQvecMin() && qNeg<fTrackCuts->GetQvecMax() && track->Eta()>=0)){\r
+         \r
+         //calculate DCA for AOD track\r
+         Double_t d[2], covd[3];\r
+         AliAODTrack* track_clone=(AliAODTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters\r
+         Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);\r
+         delete track_clone;\r
+         if(!isDCA)d[0]=-999;\r
 \r
-   PostData(1, fHistMan);\r
-   PostData(2, fEventCuts);\r
-   PostData(3, fTrackCuts);\r
+         fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),d[0]);  // PT histo\r
+         //Response\r
+         fHistMan->GetPIDHistogram(kHistPIDTPC)->Fill(track->GetTPCmomentum(), track->GetTPCsignal()*track->Charge()); // PID histo\r
+         fHistMan->GetPIDHistogram(kHistPIDTOF)->Fill(track->P(),track->GetTOFsignal()); // PID histo\r
+         \r
+         AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(track);\r
+         Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
+         Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)); \r
+         Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)); \r
+         Double_t nsigmaTOFkProton=0,nsigmaTOFkKaon=0,nsigmaTOFkPion=0;\r
+         if(track->Pt()>fTrackCuts->GetPtTOFMatching()){\r
+           nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
+           nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)); \r
+           nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)); \r
+           \r
+           //TOF\r
+           fHistMan->GetPtHistogram(kHistNSigProtonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
+           fHistMan->GetPtHistogram(kHistNSigKaonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));\r
+           fHistMan->GetPtHistogram(kHistNSigPionTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));\r
+           fHistMan->GetPtHistogram(kHistNSigProtonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
+           fHistMan->GetPtHistogram(kHistNSigKaonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));\r
+           fHistMan->GetPtHistogram(kHistNSigPionPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));\r
+         }\r
+         \r
+         Double_t nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);\r
+         Double_t nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);\r
+         Double_t nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);\r
+         \r
+         //TPC\r
+         fHistMan->GetPtHistogram(kHistNSigProtonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
+         fHistMan->GetPtHistogram(kHistNSigKaonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));\r
+         fHistMan->GetPtHistogram(kHistNSigPionTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));\r
+         fHistMan->GetPtHistogram(kHistNSigProtonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
+         fHistMan->GetPtHistogram(kHistNSigKaonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));\r
+         fHistMan->GetPtHistogram(kHistNSigPionPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));\r
+         //TPCTOF\r
+         fHistMan->GetPtHistogram(kHistNSigProtonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkProton);\r
+         fHistMan->GetPtHistogram(kHistNSigKaonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkKaon);\r
+         fHistMan->GetPtHistogram(kHistNSigPionTPCTOF)->Fill(track->P(),nsigmaTPCTOFkPion);\r
+        fHistMan->GetPtHistogram(kHistNSigProtonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkProton);\r
+         fHistMan->GetPtHistogram(kHistNSigKaonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkKaon);\r
+         fHistMan->GetPtHistogram(kHistNSigPionPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkPion);\r
+       \r
+         \r
+         if( ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkProton )) { \r
+           if ((nsigmaTPCTOFkKaon > fNSigmaPID) || (!CheckYCut(kKaon, track) ) ) continue;\r
+           if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaKaonPlus)->Fill(track->Pt(),d[0]); } \r
+           else { fHistMan->GetPtHistogram(kHistPtRecSigmaKaonMinus)->Fill(track->Pt(),d[0]); } \r
+         }\r
+         if( ( nsigmaTPCTOFkProton < nsigmaTPCTOFkKaon ) && ( nsigmaTPCTOFkProton < nsigmaTPCTOFkPion ) ) {\r
+           if ( nsigmaTPCTOFkProton > fNSigmaPID || (!CheckYCut(kProton, track) ) )  continue;\r
+           if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaProtonPlus)->Fill(track->Pt(),d[0]); }\r
+           else { fHistMan->GetPtHistogram(kHistPtRecSigmaProtonMinus)->Fill(track->Pt(),d[0]); }\r
+         }\r
+         if( (nsigmaTPCTOFkPion < nsigmaTPCTOFkProton ) && ( nsigmaTPCTOFkPion < nsigmaTPCTOFkKaon ) ) {\r
+           if (nsigmaTPCTOFkPion > fNSigmaPID || (!CheckYCut(kPion, track) ) ) continue;\r
+           if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaPionPlus)->Fill(track->Pt(),d[0]); }\r
+           else  { fHistMan->GetPtHistogram(kHistPtRecSigmaPionMinus)->Fill(track->Pt(),d[0]); }\r
+         }\r
+         /* MC Part */\r
+         // MF 22/02/2012\r
+         // fill mc DCA vs pt\r
+         if (arrayMC) {\r
+           AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
+           if (!partMC) { \r
+             AliError("Cannot get MC particle");\r
+             continue; }\r
+           if (CheckYCut(partMC)) {\r
+             // primaries, true pid\r
+             if ( partMC->PdgCode() == 2212) { fHistMan->GetPtHistogram(kHistPtRecTrueProtonPlus)->Fill(track->Pt(),d[0]); \r
+               if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryProtonPlus)->Fill(track->Pt(),d[0]); }}\r
+             if ( partMC->PdgCode() == -2212) { fHistMan->GetPtHistogram(kHistPtRecTrueProtonMinus)->Fill(track->Pt(),d[0]); \r
+               if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryProtonMinus)->Fill(track->Pt(),d[0]); }}\r
+             if ( partMC->PdgCode() == 321) { fHistMan->GetPtHistogram(kHistPtRecTrueKaonPlus)->Fill(track->Pt(),d[0]); \r
+               if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryKaonPlus)->Fill(track->Pt(),d[0]); }}\r
+             if ( partMC->PdgCode() == -321) { fHistMan->GetPtHistogram(kHistPtRecTrueKaonMinus)->Fill(track->Pt(),d[0]); \r
+               if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryKaonMinus)->Fill(track->Pt(),d[0]); }}\r
+             if ( partMC->PdgCode() == 211) { fHistMan->GetPtHistogram(kHistPtRecTruePionPlus)->Fill(track->Pt(),d[0]); \r
+               if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionPlus)->Fill(track->Pt(),d[0]); }}\r
+             if ( partMC->PdgCode() == -211) { fHistMan->GetPtHistogram(kHistPtRecTruePionMinus)->Fill(track->Pt(),d[0]);  \r
+               if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionMinus)->Fill(track->Pt(),d[0]); }}\r
+             \r
+             // primaries, sigma pid \r
+             if (partMC->IsPhysicalPrimary()) { \r
+               if( ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkProton ) ) {\r
+                 if ( (nsigmaTPCTOFkKaon > fNSigmaPID ) || (!CheckYCut(kKaon, track) ) ) continue; \r
+                 if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryKaonPlus)->Fill(track->Pt(),d[0]); } \r
+                 else { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryKaonMinus)->Fill(track->Pt(),d[0]); } \r
+               }\r
+               if( ( nsigmaTPCTOFkProton < nsigmaTPCTOFkKaon ) && ( nsigmaTPCTOFkProton < nsigmaTPCTOFkPion ) ) {\r
+                 if ( (nsigmaTPCTOFkProton > fNSigmaPID ) || (!CheckYCut(kProton, track) ) ) continue;\r
+                 if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryProtonPlus)->Fill(track->Pt(),d[0]); }\r
+                 else { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryProtonMinus)->Fill(track->Pt(),d[0]); }\r
+               }\r
+               if( (nsigmaTPCTOFkPion < nsigmaTPCTOFkProton ) && ( nsigmaTPCTOFkPion < nsigmaTPCTOFkKaon ) ) {\r
+                 if ( ( nsigmaTPCTOFkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
+                 if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryPionPlus)->Fill(track->Pt(),d[0]); }\r
+                 else  { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryPionMinus)->Fill(track->Pt(),d[0]); }\r
+               }\r
+             }//end if primaries\r
+             // MF 22/02/2012\r
+             // Distinguish weak decay and material\r
+             //done quickly, code can be improved\r
+             else {//to be added in a separate class\r
+               Int_t mfl=-999,uniqueID=-999;\r
+               Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()\r
+               if(indexMoth>=0){//is not fake\r
+                 AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);\r
+                 Float_t codemoth = TMath::Abs(moth->GetPdgCode());\r
+                 mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));\r
+               }\r
+               uniqueID = partMC->GetUniqueID();\r
+               //if(mfl==3 && uniqueID == kPDecay){//strangeness KPDecay not working on AOD, to be understood\r
+               if(mfl==3){//strangeness\r
+                 if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) { \r
+                   if ( (nsigmaTPCkKaon > fNSigmaPID )  || (!CheckYCut(kKaon, track) ) ) continue;\r
+                   if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayKaonPlus)->Fill(track->Pt(),d[0]); } \r
+                   else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayKaonMinus)->Fill(track->Pt(),d[0]); } \r
+                 }\r
+                 if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
+                   if ( (nsigmaTPCkProton > fNSigmaPID )  || (!CheckYCut(kProton, track) ) ) continue;\r
+                   if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayProtonPlus)->Fill(track->Pt(),d[0]); }\r
+                   else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayProtonMinus)->Fill(track->Pt(),d[0]); }\r
+                 }\r
+                 if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
+                   if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
+                   if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayPionPlus)->Fill(track->Pt(),d[0]); }\r
+                   else  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayPionMinus)->Fill(track->Pt(),d[0]); }\r
+                 }\r
+               }//end if strangeness\r
+               else{//material\r
+                 if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) { \r
+                   if ( (nsigmaTPCkKaon > fNSigmaPID )  || (!CheckYCut(kKaon, track) ) ) continue;\r
+                   if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialKaonPlus)->Fill(track->Pt(),d[0]); } \r
+                   else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialKaonMinus)->Fill(track->Pt(),d[0]); } \r
+                 }\r
+                 if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
+                   if ( (nsigmaTPCkProton > fNSigmaPID )  || (!CheckYCut(kProton, track) ) ) continue;\r
+                   if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialProtonPlus)->Fill(track->Pt(),d[0]); }\r
+                   else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialProtonMinus)->Fill(track->Pt(),d[0]); }\r
+                 }\r
+                 if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
+                   if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
+                   if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialPionPlus)->Fill(track->Pt(),d[0]); }\r
+                   else  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialPionMinus)->Fill(track->Pt(),d[0]); }\r
+                 }\r
+               }//end if material\r
+             }//end else (IsPhysPrim)\r
+           }//end if(CheckY)\r
+         }//end if(arrayMC)\r
+       }//end if on q vector (track)\r
+      } // end loop on tracks\r
+  }//end if on q vector (event)\r
+  PostData(1, fHistMan);\r
+  PostData(2, fEventCuts);\r
+  PostData(3, fTrackCuts);\r
 }\r
-  \r
+\r
 //_________________________________________________________________\r
 void   AliAnalysisTaskSpectraAOD::Terminate(Option_t *)\r
 {\r
-   // Terminate\r
+  // Terminate\r
 }\r
index fa7faaeac194df0ef18c574f4a88b42d119af23e..78f37e14a34b811d126aa14cd5dacc069a1b9db3 100644 (file)
@@ -52,27 +52,65 @@ AliSpectraAODHistoManager::AliSpectraAODHistoManager(const char *name): TNamed(n
    TH1::AddDirectory(kFALSE);
    for (Int_t ihist  = 0; ihist < kNHist ; ihist++)
    {
-      if (ihist <= kNPtHist) BookPtHistogram(kHistName[ihist]);  // PT histos
-      if (ihist > kNPtHist) BookPIDHistogram(kHistName[ihist]);  // PID histos
+      if (ihist <= kNPtGenHist) BookPtGenHistogram(kHistName[ihist]);  // PT histos
+      if (ihist > kNPtGenHist && ihist <= kNPtRecHist) BookPtRecHistogram(kHistName[ihist]);  // PT histos
+      if (ihist > kNPtRecHist && ihist <= kNHistPID) BookPIDHistogram(kHistName[ihist]);  // PID histos
+      if (ihist > kNHistPID && ihist <= kNHistNSig) BookNSigHistogram(kHistName[ihist]);  // NSigmaSep histos
+      if (ihist > kNHistNSig) BookqVecHistogram(kHistName[ihist]);  // qDistr histos
    }
+   //adding quickly o plot to check the centrality distr in two different ways
+   TH2F *hist=new TH2F("CentCheck","CentCheck",1000,0,100,1000,0,100);
+   hist->SetXTitle("fAOD->GetCentrality()->GetCentralityPercentile(V0M)");
+   hist->SetYTitle("fAOD->GetHeader()->GetCentralityP()->GetCentralityPercentileUnchecked(V0M)");
+   hist->Sumw2();
+   fOutputList->Add(hist);
+   
    TH1::AddDirectory(oldStatus);
 
 }
+
 //_______________________________________________________
-TH1F* AliSpectraAODHistoManager::BookPtHistogram(const char * name)
+
+TH2F* AliSpectraAODHistoManager::BookPtGenHistogram(const char * name)
 {
    // Return a pt histogram with predefined binning, set the ID and add it to the output list
    AliInfo(Form("Booking pt histogram %s", name));
+   
+   //standard histo
+   const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
+   Int_t nbinsTempl=52;
+   
+   TH2F * hist = new TH2F(name,Form("P_{T} distribution (%s)", name),nbinsTempl,templBins,2,-0.5,1.5);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
+   hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
+   hist->GetYaxis()->SetTitle("DCA xy");
+   hist->SetMarkerStyle(kFullCircle);
+   hist->Sumw2();
+   fOutputList->Add(hist);
+   
+   return hist;
+}
 
-   TH1F * hist = new TH1F(name, Form("P_{T} distribution (%s)", name), 40, 0, 1.2);
+
+//_______________________________________________________
+TH2F* AliSpectraAODHistoManager::BookPtRecHistogram(const char * name)
+{
+   // Return a pt histogram with predefined binning, set the ID and add it to the output list
+   AliInfo(Form("Booking pt histogram %s", name));
+   
+   //standard histo
+   const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
+   Int_t nbinsTempl=52;
+   
+   TH2F * hist = new TH2F(name,Form("reconstructed P_{T} distribution (%s)", name),nbinsTempl,templBins,3000,-3,3);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
    hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
-   hist->GetYaxis()->SetTitle("No. of events");
+   hist->GetYaxis()->SetTitle("DCA xy");
    hist->SetMarkerStyle(kFullCircle);
    hist->Sumw2();
    fOutputList->Add(hist);
 
    return hist;
 }
+
 //_____________________________________________________________________________
 
 TH2F* AliSpectraAODHistoManager::BookPIDHistogram(const char * name)
@@ -80,15 +118,85 @@ TH2F* AliSpectraAODHistoManager::BookPIDHistogram(const char * name)
    // Return a pt histogram with predefined binning, set the ID and add it to the output list
    AliInfo(Form("Booking pt histogram %s", name));
 
-   TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 100, 0, 1.2, 1000, 0, 1000);
-   hist->GetXaxis()->SetTitle("(GeV / c)");
-   hist->GetYaxis()->SetTitle("TPC");
+   TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 100, 0, 1.5, 2000, -1000, 1000);
+   hist->GetXaxis()->SetTitle("(GeV / c)");
+   hist->GetYaxis()->SetTitle("PID signal");
 //  hist->Sumw2();
    fOutputList->Add(hist);
 
    return hist;
 }
 
+//_____________________________________________________________________________
+
+TH2F* AliSpectraAODHistoManager::BookNSigHistogram(const char * name)
+{
+  // Return a pt histogram with predefined binning, set the ID and add it to the output list
+  AliInfo(Form("Booking pt histogram %s", name));
+  
+  TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 200, 0, 1.5, 4000,-40, 40);
+  hist->GetXaxis()->SetTitle("P (GeV / c)");
+  hist->GetYaxis()->SetTitle("TPC");
+  //hist->Sumw2();
+  fOutputList->Add(hist);
+  
+  return hist;
+}
+
+//_____________________________________________________________________________
+
+TH2F* AliSpectraAODHistoManager::BookqVecHistogram(const char * name)
+{
+  // Return a pt histogram with predefined binning, set the ID and add it to the output list
+  AliInfo(Form("Booking q Vector histogram %s", name));
+  
+  TH2F * hist = new TH2F(name, Form("q Vector distribution vs Centrality (%s)", name), 100, 0, 10, 100, 0, 100);
+  hist->GetXaxis()->SetTitle("q vector");
+  hist->GetYaxis()->SetTitle("Centrality (V0)");
+  //  hist->Sumw2();
+  fOutputList->Add(hist);
+  
+  return hist;
+}
+
+
+//_____________________________________________________________________________
+
+TH1F* AliSpectraAODHistoManager::GetPtHistogram1D(const char * name,Double_t minDCA,Double_t maxDCA)
+{
+  //   //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
+  //   //if minDCA=-1 && maxDCA=-1 the projection is done using the full DCA range
+  TH2F *hist=(TH2F*)fOutputList->FindObject(name);
+  TH1F *outhist=0x0;
+  if(minDCA==-1 && maxDCA==-1)outhist=(TH1F*)hist->ProjectionX("",0,-1,"e");
+  else {
+    Int_t firstbin=hist->GetYaxis()->FindBin(minDCA);
+    Int_t lastbin=hist->GetYaxis()->FindBin(maxDCA);
+    Printf("firstbin: %d lastbin: %d",firstbin,lastbin);
+    outhist=(TH1F*)hist->ProjectionX("",firstbin,lastbin,"e");
+  }
+  Printf("Entries outhist: %.0f   Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries());
+  return outhist;
+}
+
+//_____________________________________________________________________________
+
+TH1F* AliSpectraAODHistoManager::GetDCAHistogram1D(const char * name,Double_t minPt,Double_t maxPt)
+{
+  //   //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
+  //   //if minPt=-1 && maxPt=-1 the projection is done using the full DCA range
+  TH2F *hist=(TH2F*)fOutputList->FindObject(name);
+  TH1F *outhist=0x0;
+  if(minPt==-1 && maxPt==-1)outhist=(TH1F*)hist->ProjectionY("",0,-1,"e");
+  else {
+    Int_t firstbin=hist->GetXaxis()->FindBin(minPt);
+    Int_t lastbin=hist->GetXaxis()->FindBin(maxPt);
+    outhist=(TH1F*)hist->ProjectionY("",firstbin,lastbin,"e");
+  }
+  return outhist;
+}
+
+//_____________________________________________________________________________
 
 Long64_t AliSpectraAODHistoManager::Merge(TCollection* list)
 {
index acb6c3c81e33756f26e8473ff3f733e8904503a0..24e62c62698b9eff7f0a73b6cfa95b60f1a5620b 100644 (file)
@@ -27,8 +27,17 @@ namespace AliSpectraNameSpace
      // Add histograms 2D DCA_xy (-3,3) vs pt (0, 3)
      // For Rec data/MC, Rec MC primaries, Rec MC secondaries weak decay, Rec MC secondaries material (x6, each particle hypothesis)
 
+      // 6 Pt Generated True Primary
+      kHistPtGenTruePrimaryProtonPlus=0,          // Pt histo for protons +, generated tracks, true ID, primary Event
+      kHistPtGenTruePrimaryKaonPlus,            // Pt histo for kaons +, generated tracks, true ID, primary Event
+      kHistPtGenTruePrimaryPionPlus,            // Pt histo for pions +, generated tracks, true ID, primary Event
+      kHistPtGenTruePrimaryProtonMinus,         // Pt histo for protons -, generated tracks, true ID, primary Event
+      kHistPtGenTruePrimaryKaonMinus,           // Pt histo for kaons -, generated tracks, true ID, primary Event
+      kHistPtGenTruePrimaryPionMinus,           // Pt histo for pions -, generated tracks, true ID, primary Event
+      kNPtGenHist = kHistPtGenTruePrimaryPionMinus,                    // Number of ptGen-likehistos histos
+      
       // 6 Pt Reconstructed Sigma
-      kHistPtRecSigmaProtonPlus = 0,            // Pt histo for protons +, reconstructed tracks, sigma ID
+      kHistPtRecSigmaProtonPlus,            // Pt histo for protons +, reconstructed tracks, sigma ID
       kHistPtRecSigmaKaonPlus,                  // Pt histo for kaons +, reconsructed tracks, sigma ID
       kHistPtRecSigmaPionPlus,                  // Pt histo for pions +, reconstructed tracks, sigma ID
       kHistPtRecSigmaProtonMinus,               // Pt histo for protons -, reconstructed tracks, sigma ID
@@ -51,23 +60,22 @@ namespace AliSpectraNameSpace
       kHistPtRecSigmaPrimaryKaonMinus,          // Pt histo for kaons -, reconstructed tracks, sigma ID, primary Event
       kHistPtRecSigmaPrimaryPionMinus,          // Pt histo for pions -, reconstructed tracks, sigma ID, primary Event
             
-      // 6 Pt Reconstructed Sigma Secondary
-      kHistPtRecSigmaSecondaryProtonPlus,       // Pt histo for protons +, reconstructed tracks, sigma ID, secondary Event
-      kHistPtRecSigmaSecondaryKaonPlus,         // Pt histo for kaons +, reconsructed tracks, sigma ID, secondary Event
-      kHistPtRecSigmaSecondaryPionPlus,         // Pt histo for pions +, reconstructed tracks, sigma ID, secondary Event
-      kHistPtRecSigmaSecondaryProtonMinus,      // Pt histo for protons -, reconstructed tracks, sigma ID, secondary Event
-      kHistPtRecSigmaSecondaryKaonMinus,        // Pt histo for kaons -, reconstructed tracks, sigma ID, secondary Event
-      kHistPtRecSigmaSecondaryPionMinus,        // Pt histo for pions -, reconstructed tracks, sigma ID, secondary Event
+      // 6 Pt Reconstructed Sigma Secondary Material
+      kHistPtRecSigmaSecondaryMaterialProtonPlus,       // Pt histo for protons +, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryMaterialKaonPlus,         // Pt histo for kaons +, reconsructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryMaterialPionPlus,         // Pt histo for pions +, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryMaterialProtonMinus,      // Pt histo for protons -, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryMaterialKaonMinus,        // Pt histo for kaons -, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryMaterialPionMinus,        // Pt histo for pions -, reconstructed tracks, sigma ID, secondary Event
+
+      // 6 Pt Reconstructed Sigma Secondary WeakDecay
+      kHistPtRecSigmaSecondaryWeakDecayProtonPlus,       // Pt histo for protons +, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryWeakDecayKaonPlus,         // Pt histo for kaons +, reconsructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryWeakDecayPionPlus,         // Pt histo for pions +, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryWeakDecayProtonMinus,      // Pt histo for protons -, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryWeakDecayKaonMinus,        // Pt histo for kaons -, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryWeakDecayPionMinus,        // Pt histo for pions -, reconstructed tracks, sigma ID, secondary Event
 
-      // 6 Pt Generated True Primary
-      kHistPtGenTruePrimaryProtonPlus,          // Pt histo for protons +, generated tracks, true ID, primary Event
-      kHistPtGenTruePrimaryKaonPlus,            // Pt histo for kaons +, generated tracks, true ID, primary Event
-      kHistPtGenTruePrimaryPionPlus,            // Pt histo for pions +, generated tracks, true ID, primary Event
-      kHistPtGenTruePrimaryProtonMinus,         // Pt histo for protons -, generated tracks, true ID, primary Event
-      kHistPtGenTruePrimaryKaonMinus,           // Pt histo for kaons -, generated tracks, true ID, primary Event
-      kHistPtGenTruePrimaryPionMinus,           // Pt histo for pions -, generated tracks, true ID, primary Event
-      kNPtSpecies = kHistPtGenTruePrimaryPionMinus,
-      
       // 6 Pt Reconstructed True Primary
       kHistPtRecTruePrimaryProtonPlus,          // Pt histo for protons +, reconstructed tracks, true ID, primary event
       kHistPtRecTruePrimaryKaonPlus,            // Pt histo for kaons +, reconsructed tracks, true ID, primary event
@@ -79,14 +87,49 @@ namespace AliSpectraNameSpace
       // Rest
       kHistPtRec,                               // Pt histo for all particles, reconstructed tracks
       kHistPtGen,                               // Pt histo for all particles, generated tracks
-      kNPtHist = kHistPtGen,                    // Number of pt-likehistos histos
+      kNPtRecHist = kHistPtGen,                    // Number of ptRec-likehistos histos
+      
       kHistPIDTPC,                              // Particle Identification histo
-      kHistPIDTPCPt,
+      kHistPIDTOF,                              
+      kNHistPID =kHistPIDTOF,                           
+      
+      kHistNSigProtonTPC,                       // NSigma separation plot    
+      kHistNSigKaonTPC,                              
+      kHistNSigPionTPC,                              
+      kHistNSigProtonPtTPC,                              
+      kHistNSigKaonPtTPC,                              
+      kHistNSigPionPtTPC,                              
+      
+      kHistNSigProtonTOF,                              
+      kHistNSigKaonTOF,                              
+      kHistNSigPionTOF,                              
+      kHistNSigProtonPtTOF,                              
+      kHistNSigKaonPtTOF,                              
+      kHistNSigPionPtTOF,                              
+     
+      kHistNSigProtonTPCTOF,                             
+      kHistNSigKaonTPCTOF,                              
+      kHistNSigPionTPCTOF,                              
+      kHistNSigProtonPtTPCTOF,                              
+      kHistNSigKaonPtTPCTOF,                              
+      kHistNSigPionPtTPCTOF,
+      kNHistNSig=kHistNSigPionPtTPCTOF,                              
+      
+      kHistqVecPos,
+      kHistqVecNeg,
       kNHist,                                   // Total number of histos
    };  // Type of events plotted in Pt Histogram
 
    const char * kHistName[] =
    {
+      // 6 Pt Reconstructed Sigma Primary
+      "histPtGenTruePrimaryProtonPlus",         // Pt histo for protons +, generated tracks, sigma ID, primary Event
+      "histPtGenTruePrimaryKaonPlus",           // Pt histo for kaons +, generated tracks, sigma ID, primary Event
+      "histPtGenTruePrimaryPionPlus",           // Pt histo for pions +, generated tracks, sigma ID, primary Event
+      "histPtGenTruePrimaryProtonMinus",          // Pt histo for protons -, generated tracks, sigma ID, primary Event
+      "histPtGenTruePrimaryKaonMinus",            // Pt histo for kaons -, generated tracks, sigma ID, primary Event
+      "histPtGenTruePrimaryPionMinus",            // Pt histo for pions -, generated tracks, sigma ID, primary Event
+      
       // 6 Pt Reconstructed Sigma
       "histPtRecSigmaProtonPlus",               // Pt histo for protons +, reconstructed tracks, sigma ID
       "histPtRecSigmaKaonPlus",                 // Pt histo for kaons +, reconsructed tracks, sigma ID
@@ -112,21 +155,21 @@ namespace AliSpectraNameSpace
       "histPtRecSigmaPrimaryPionMinus",         // Pt histo for pions -, reconstructed tracks, sigma ID, primary Event
       
       // 6 Pt Reconstructed Sigma Seconday
-      "histPtRecSigmaSecondaryProtonPlus",      // Pt histo for protons +, reconstructed tracks, sigma ID, secondary Event
-      "histPtRecSigmaSecondaryKaonPlus",        // Pt histo for kaons +, reconsructed tracks, sigma ID, secondary Event
-      "histPtRecSigmaSecondaryPionPlus",        // Pt histo for pions +, reconstructed tracks, sigma ID, secondary Event
-      "histPtRecSigmaSecondaryProtonMinus",     // Pt histo for protons -, reconstructed tracks, sigma ID, secondary Event
-      "histPtRecSigmaSecondaryKaonMinus",       // Pt histo for kaons -, reconstructed tracks, sigma ID, secondary Event
-      "histPtRecSigmaSecondaryPionMinus",       // Pt histo for pions -, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryMaterialProtonPlus",      // Pt histo for protons +, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryMaterialKaonPlus",        // Pt histo for kaons +, reconsructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryMaterialPionPlus",        // Pt histo for pions +, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryMaterialProtonMinus",     // Pt histo for protons -, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryMaterialKaonMinus",       // Pt histo for kaons -, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryMaterialPionMinus",       // Pt histo for pions -, reconstructed tracks, sigma ID, secondary Event
       
-      // 6 Pt Reconstructed Sigma Primary
-      "histPtGenTruePrimaryProtonPlus",         // Pt histo for protons +, generated tracks, sigma ID, primary Event
-      "histPtGenTruePrimaryKaonPlus",           // Pt histo for kaons +, generated tracks, sigma ID, primary Event
-      "histPtGenTruePrimaryPionPlus",           // Pt histo for pions +, generated tracks, sigma ID, primary Event
-      "histPtGenTruePrimaryProtonMinus",          // Pt histo for protons -, generated tracks, sigma ID, primary Event
-      "histPtGenTruePrimaryKaonMinus",            // Pt histo for kaons -, generated tracks, sigma ID, primary Event
-      "histPtGenTruePrimaryPionMinus",            // Pt histo for pions -, generated tracks, sigma ID, primary Event
-
+      // 6 Pt Reconstructed Sigma Seconday
+      "histPtRecSigmaSecondaryWeakDecayProtonPlus",      // Pt histo for protons +, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryWeakDecayKaonPlus",        // Pt histo for kaons +, reconsructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryWeakDecayPionPlus",        // Pt histo for pions +, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryWeakDecayProtonMinus",     // Pt histo for protons -, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryWeakDecayKaonMinus",       // Pt histo for kaons -, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryWeakDecayPionMinus",       // Pt histo for pions -, reconstructed tracks, sigma ID, secondary Event
+        
       // 6 Pt Reconstructed True
       "histPtRecTruePrimaryProtonPlus",         // Pt histo for protons +, reconstructed tracks, true ID, primary event
       "histPtRecTruePrimaryKaonPlus",           // Pt histo for kaons +, reconsructed tracks, true ID, primary event
@@ -138,16 +181,41 @@ namespace AliSpectraNameSpace
       // Rest
       "histPtRec",                              // Pt histo for all particles, reconstructed tracks
       "histPtGen",                              // Pt histo for all particles, generated tracks
+     
       "histPIDTPC",                             // Particle Identification histo
-      "histPIDTPCPt",
+      "histPIDTOF",                             
+     
+      "histNSigProtonTPC",                      // NSigma Separation plot
+      "histNSigKaonTPC",
+      "histNSigPionTPC",
+      "histNSigProtonPtTPC",
+      "histNSigKaonPtTPC",
+      "histNSigPionPtTPC",
+      
+      "histNSigProtonTOF",
+      "histNSigKaonTOF",
+      "histNSigPionTOF",
+      "histNSigProtonPtTOF",
+      "histNSigKaonPtTOF",
+      "histNSigPionPtTOF",
+      
+      "histNSigProtonTPCTOF",
+      "histNSigKaonTPCTOF",
+      "histNSigPionTPCTOF",
+      "histNSigProtonPtTPCTOF",
+      "histNSigKaonPtTPCTOF",
+      "histNSigPionPtTPCTOF",
+      
+      "histqPos",                             // qVecVsCentrality
+      "histqNeg"
    };
-
+   
    enum AODParticleSpecies_t
    {
-       kProton = 0,
-       kKaon,
-       kPion,
-       kNSpecies = kPion,
+     kProton = 0,
+     kKaon,
+     kPion,
+     kNSpecies = kPion,
    }; // Particle species used in plotting
 
    const char * kParticleSpecies[] =
@@ -178,17 +246,29 @@ public:
    virtual  ~AliSpectraAODHistoManager() {}
 
 
-   TH1F*   BookPtHistogram(const char * name);
+   TH2F*   BookPtGenHistogram(const char * name);
+   TH2F*   BookPtRecHistogram(const char * name);
    TH2F*   BookPIDHistogram(const char * name);
-   TH1*     GetHistogram(UInt_t id)      {      return (TH1*) fOutputList->At(id);   }
-   TH1*     GetPtHistogram(UInt_t id)    {      return (TH1*) fOutputList->At(id);   }
-   TH1*     GetPtHistogram(const char * name)   {      return (TH1*) fOutputList->FindObject(name);   }
-   TH1*     GetPtHistogramByName(UInt_t id)     {      return (TH1*) fOutputList->FindObject(kHistName[id]); }  // Use this if you want to read a file saved with a different histo list   
+   TH2F*   BookNSigHistogram(const char * name);
+   TH2F*   BookqVecHistogram(const char * name);
+   TH1F*   GetPtHistogram1D(const char * name,Double_t minDCA,Double_t maxDCA);
+   TH1F*   GetDCAHistogram1D(const char * name,Double_t minPt,Double_t maxPt);
+   TH2*     GetHistogram(UInt_t id)      {      return (TH2*) fOutputList->At(id);   }
+   TH2*     GetPtHistogram(UInt_t id)    {      return (TH2*) fOutputList->At(id);   }
+   TH2*     GetPtHistogram(const char * name)   {      return (TH2*) fOutputList->FindObject(name);   }
+   TH2*     GetPtHistogramByName(UInt_t id)     {      return (TH2*) fOutputList->FindObject(kHistName[id]); }  // Use this if you want to read a file saved with a different histo list   
    TH2*     GetPIDHistogram(UInt_t id)   {      return (TH2*) fOutputList->At(id);   }
    TH2*     GetPIDHistogram(const char * name)  {      return (TH2*) fOutputList->FindObject(name);   }
    TH2*     GetPIDHistogramByName(UInt_t id)    {      return (TH2*) fOutputList->FindObject(kHistName[id]);  }// Use this if you want to read a file saved with a different histo list
-   TH1F*   GetTH1F(UInt_t id)            {      return (TH1F*) GetPtHistogram(id);   }
-   TH2F*   GetTH2F(UInt_t id)            {      return (TH2F*) GetPIDHistogram(id);   }
+   TH2*     GetNSigHistogram(UInt_t id)   {      return (TH2*) fOutputList->At(id);   }
+   TH2*     GetNSigHistogram(const char * name)  {      return (TH2*) fOutputList->FindObject(name);   }
+   TH2*     GetNSigHistogramByName(UInt_t id)    {      return (TH2*) fOutputList->FindObject(kHistName[id]);  }// Use this if you want to read a file saved with a different histo list
+   TH2*     GetqVecHistogram(UInt_t id)   {      return (TH2*) fOutputList->At(id);   }
+   TH2*     GetqVecHistogram(const char * name)  {      return (TH2*) fOutputList->FindObject(name);   }
+   TH2*     GetqVecHistogramByName(UInt_t id)    {      return (TH2*) fOutputList->FindObject(kHistName[id]);  }// Use this if you want to read a file saved with a different histo list
+  
+   //TH1F*   GetTH1F(UInt_t id)            {      return (TH1F*) GetPtHistogram(id);   }
+   //TH2F*   GetTH2F(UInt_t id)            {      return (TH2F*) GetPIDHistogram(id);   }
 
   TList * GetOutputList() {return fOutputList;}
 
index de04d4d4395287bae4e926cdb6a6c407a651fe92..b431c183c3320e1d3081f14eb6b9157d286a5ddd 100644 (file)
@@ -41,7 +41,7 @@ using namespace std;
 ClassImp(AliSpectraAODTrackCuts)
 
 
-AliSpectraAODTrackCuts::AliSpectraAODTrackCuts(const char *name) : TNamed(name, "AOD Track Cuts"), fIsSelected(0), fTrackBits(0), fEtaCut(0), fDCACut(0), fPCut(0), fPtCut(0), fHistoCuts(0), fTrack(0)
+AliSpectraAODTrackCuts::AliSpectraAODTrackCuts(const char *name) : TNamed(name, "AOD Track Cuts"), fIsSelected(0), fTrackBits(0), fEtaCut(0), fDCACut(0), fPCut(0), fPtCut(0),fQvecCutMin(0),fQvecCutMax(0), fHistoCuts(0), fTrack(0)
 
 {
    // Constructor
@@ -50,6 +50,8 @@ AliSpectraAODTrackCuts::AliSpectraAODTrackCuts(const char *name) : TNamed(name,
    fDCACut = 100000.0; // default value of dca cut ~ no cut
    fPCut = 100000.0; // default value of p cut ~ no cut
    fPtCut = 100000.0; // default value of pt cut ~ no cut 
+   fQvecCutMin = -100000.0; // default value of qvec cut ~ no cut 
+   fQvecCutMax = 100000.0; // default value of qvec cut ~ no cut 
    
 }
 
@@ -63,7 +65,7 @@ Bool_t AliSpectraAODTrackCuts::IsSelected(AliAODTrack * track)
       return kFALSE;
    }
    fTrack = track;
-   fIsSelected = (CheckTrackType() && CheckEtaCut() && CheckDCACut() && CheckPCut() && CheckPtCut());
+   fIsSelected = (CheckTrackType() && CheckEtaCut() && CheckDCACut() && CheckPCut() && CheckPtCut() && CheckTOFMatching());
    return fIsSelected ;
 }
 //_________________________________________________________
@@ -108,6 +110,21 @@ Bool_t AliSpectraAODTrackCuts::CheckPtCut()
     fHistoCuts->Fill(kTrkPt);
     return kFALSE;
 }
+
+//_______________________________________________________
+Bool_t AliSpectraAODTrackCuts::CheckTOFMatching()
+{
+  // check Pt cut
+  //    if ((fTrack->Pt() < fPtCut) && (fTrack->Pt() > 0.3 )) return kTRUE;
+  
+  if (fTrack->Pt() < fPtCutTOFMatching) return kTRUE;
+  else{
+    UInt_t status; 
+    status=fTrack->GetStatus();
+    if((status&AliAODTrack::kTOFout)==0 || (status&AliAODTrack::kTIME)==0) return kFALSE; //tof matching and PID
+    return kTRUE;
+  }
+}
 //_______________________________________________________
 void AliSpectraAODTrackCuts::PrintCuts() const
 {
@@ -118,6 +135,8 @@ void AliSpectraAODTrackCuts::PrintCuts() const
     cout << " > DCA cut\t" << fDCACut << endl;
     cout << " > P cut\t" << fPCut << endl;
     cout << " > Pt cut \t" << fPtCut << endl;
+    cout << " > Q vactor Min \t" << fQvecCutMin << endl;
+    cout << " > Q vactor Max \t" << fQvecCutMax << endl;
 }
 //_______________________________________________________
 void AliSpectraAODTrackCuts::SetTrackType(UInt_t bit)
index b00f92ab9ece16cca2ddcd09675c06efeb334d22..9452adcb7e9c27100cceee6bd8901197ab034574 100644 (file)
@@ -24,7 +24,7 @@ public:
    enum { kTrkBit = 0, kTrkEta, kTrkDCA, kTrkP, kTrkPt, kNTrkCuts};
 
 
-   AliSpectraAODTrackCuts() : TNamed(), fIsSelected(0), fTrackBits(0), fEtaCut(0), fPCut(0), fPtCut(0), fHistoCuts(0), fTrack(0) {}
+   AliSpectraAODTrackCuts() : TNamed(), fIsSelected(0), fTrackBits(0), fEtaCut(0), fPCut(0), fPtCut(0), fPtCutTOFMatching(0), fQvecCutMin(0), fQvecCutMax(0), fHistoCuts(0), fTrack(0) {}
 
    AliSpectraAODTrackCuts(const char *name);
    virtual  ~AliSpectraAODTrackCuts() {} // To be implemented
@@ -37,6 +37,7 @@ public:
    Bool_t CheckDCACut();
    Bool_t CheckPCut();
    Bool_t CheckPtCut();
+   Bool_t CheckTOFMatching();
    void PrintCuts() const;
 
    UInt_t GetTrackType()  const    { return fTrackBits; }
@@ -45,29 +46,38 @@ public:
    void SetDCA(Float_t dca)   { fDCACut = dca; }
    void SetP(Float_t p)       { fPCut = p; }
    void SetPt(Float_t pt)     { fPtCut = pt; }
+   void SetPtTOFMatching(Float_t pt)     { fPtCutTOFMatching = pt; }
+   void SetQvecMin(Float_t qvecmin)     { fQvecCutMin = qvecmin; }
+   void SetQvecMax(Float_t qvecmax)     { fQvecCutMax = qvecmax; }
    Float_t GetEta()       const    { return fEtaCut; }
    Float_t GetDCA()       const    { return fDCACut; }
    Float_t GetP()         const    { return fPCut; }
    Float_t GetPt()        const    { return fPtCut; }
-
-  Long64_t Merge(TCollection* list);
-
-
-private:
-
+   Float_t GetPtTOFMatching()        const    { return fPtCutTOFMatching; }
+   Float_t GetQvecMin()        const    { return fQvecCutMin; }
+   Float_t GetQvecMax()        const    { return fQvecCutMax; }
+    
+   Long64_t Merge(TCollection* list);
+   
+   
+ private:
+   
    Bool_t         fIsSelected;      // True if cuts are selected
    UInt_t         fTrackBits;       // Type of track to be used
    Float_t        fEtaCut;          // Allowed absolute maximum value of Eta
    Float_t        fDCACut;          // Maximum value of DCA
    Float_t        fPCut;            // Maximum value of P
    Float_t        fPtCut;           // Maximum value of Pt
-
+   Float_t        fPtCutTOFMatching;           // TOF Matching
+   Float_t        fQvecCutMin;           // Minimum value of Qvec, done in the analysis task
+   Float_t        fQvecCutMax;           // Minimum value of Qvec, done in the analysis task
+   
    TH1I *         fHistoCuts;       // Cuts statistics
    AliAODTrack *  fTrack;           //! Track pointer
    
    AliSpectraAODTrackCuts(const AliSpectraAODTrackCuts&);
    AliSpectraAODTrackCuts& operator=(const AliSpectraAODTrackCuts&);
-
+   
    ClassDef(AliSpectraAODTrackCuts, 1);
 };
 #endif
index acd4f8c2b455235c6629febc489ba5f4b9c46822..f6af08a139d7647142bfbfedfc2c52bab61c976f 100644 (file)
@@ -1,8 +1,7 @@
-void runAODProof(const char * proofMode = "full")
+void runAODProof(Int_t c=1, const char * proofMode = "full")
 {
 
    gEnv->SetValue("XSec.GSI.DelegProxy", "2");
-   //  TProof::Open("rbertens@alice-caf.cern.ch");
 
    gSystem->Load("libTree.so");
    gSystem->Load("libGeom.so");
@@ -21,17 +20,28 @@ void runAODProof(const char * proofMode = "full")
    handler->SetOverwriteMode();
    handler->SetRunMode(proofMode);
    handler->SetProofReset(0);
-   handler->SetAliROOTVersion("v4-21-29-AN");
+   //handler->SetROOTVersion("v5-33-02a");
+   //handler->SetAliROOTVersion("v5-03-11-AN");
+   handler->SetAliROOTVersion("v5-03-13-AN");
+
    handler->SetProofCluster(Form("%s@alice-caf.cern.ch", gSystem->Getenv("CAFUSER")));
-   //handler->SetProofCluster("rbertens@skaf.saske.sk");
-   //   handler->SetProofDataSet("/alice/data/LHC10h_000138653_p2_AOD049#aodTree");
-   handler->SetProofDataSet("/alice/sim/LHC11a10a_000138653_AOD048#aodTree");
+//   handler->SetProofCluster(Form("%s@skaf.saske.sk",gSystem->Getenv("CAFUSER")));
+
+   // Set handler for Real DATA:
+   if (c == 1)
+   {
+     handler->SetProofDataSet("/alice/data/LHC10h_000138653_p2_AOD049#aodTree");
+ //     handler->SetProofDataSet("/alice/sim/LHC11a10a_000139107_AOD048#aodTree|alice/sim/LHC11a10a_000138653_AOD048#aodTree");
+   }
+   if (c == 2)
+   {
+      handler->SetProofDataSet("/alice/sim/LHC11a10a_000139107_AOD048#aodTree|/alice/sim/LHC11a10a_000138653_AOD048#aodTree|/alice/sim/LHC11a10a_000139110_AOD048#aodTree|/alice/sim/LHC11a10a_000138662_AOD048#aodTree|/alice/sim/LHC11a10a_000138666_AOD048#aodTree|/alice/sim/LHC11a10a_000138795_AOD048#aodTree");      
+   }
    handler->SetNproofWorkersPerSlave(1);
    handler->SetAliRootMode("default");
    handler->SetAdditionalLibs("AliSpectraAODHistoManager.cxx AliSpectraAODHistoManager.h AliSpectraAODEventCuts.cxx AliSpectraAODEventCuts.h AliSpectraAODTrackCuts.cxx AliSpectraAODTrackCuts.h AliAnalysisTaskSpectraAOD.cxx AliAnalysisTaskSpectraAOD.h");
    handler->SetAnalysisSource("AliSpectraAODHistoManager.cxx+ AliSpectraAODEventCuts.cxx+ AliSpectraAODTrackCuts.cxx+ AliAnalysisTaskSpectraAOD.cxx+");
    handler->SetFileForTestMode("filelist.txt"); // list of local files for testing
-
    //  handler->SetAliRootMode("");
    handler->SetClearPackages();
 
@@ -46,39 +56,81 @@ void runAODProof(const char * proofMode = "full")
    gROOT->LoadMacro("AliSpectraAODHistoManager.cxx+g");
    gROOT->LoadMacro("AliAnalysisTaskSpectraAOD.cxx+g");
 
+
    // Add PID task
    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
-   AliAnalysisTask * taskPID = AddTaskPIDResponse();
+   Bool_t isMC = kFALSE;
+   if (c == 2 || c == 3) isMC = kTRUE;   
+   AliAnalysisTask * taskPID = AddTaskPIDResponse(isMC);
    mgr->AddTask(taskPID);
 
-   AliAnalysisTaskSpectraAOD *task = new AliAnalysisTaskSpectraAOD("TaskAODExercise");
-   task->SetIsMC(1);
-   mgr->AddTask(task);
-
-   // Set the cuts
-   AliSpectraAODEventCuts * vcuts = new AliSpectraAODEventCuts("Event Cuts");
-   AliSpectraAODTrackCuts  * tcuts = new AliSpectraAODTrackCuts("Track Cuts");
-   tcuts->SetTrackType(6);
-   tcuts->SetEta(1.);
-   tcuts->SetP(.7);
-   // vcuts->SetCentralityCutMin(0.0) // default
-   // vcuts->SetCentralityCutMax(0.0) // default
-   task->SetEventCuts(vcuts);
-   task->SetTrackCuts(tcuts);
-   task->SetNSigmaForIdentification(3.);
-
-   AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
-   AliAnalysisDataContainer *coutputpt1 = mgr->CreateContainer("chistpt", AliSpectraAODHistoManager::Class(),  AliAnalysisManager::kOutputContainer, "Pt.AOD.1.root");
-   AliAnalysisDataContainer *coutputpt2 = mgr->CreateContainer("cvcutpt", AliSpectraAODEventCuts::Class(),    AliAnalysisManager::kOutputContainer, "Pt.AOD.1.root");
-   AliAnalysisDataContainer *coutputpt3 = mgr->CreateContainer("ctcutpt", AliSpectraAODTrackCuts::Class(),     AliAnalysisManager::kOutputContainer, "Pt.AOD.1.root");
-
-   mgr->ConnectInput(task, 0, cinput);
-   mgr->ConnectOutput(task, 1, coutputpt1);
-   mgr->ConnectOutput(task, 2, coutputpt2);
-   mgr->ConnectOutput(task, 3, coutputpt3);
+   //LOOP OVER SELECTION
+   Double_t CentCutMin[4]={0,0,20,20};
+   Double_t CentCutMax[4]={100,5,40,40};
+   Double_t QvecCutMin[4]={0,0,0,3};
+   Double_t QvecCutMax[4]={100,100,2,100};
+   for(Int_t iCut=0;iCut<4;iCut++){
+     AliAnalysisTaskSpectraAOD *task = new AliAnalysisTaskSpectraAOD("TaskAODExercise");
+     mgr->AddTask(task);
+     //physics selection
+     task->SelectCollisionCandidates();     
+     // Set the cuts
+     AliSpectraAODEventCuts * vcuts = new AliSpectraAODEventCuts("Event Cuts");
+     AliSpectraAODTrackCuts  * tcuts = new AliSpectraAODTrackCuts("Track Cuts");
+     tcuts->SetTrackType(5);
+     tcuts->SetEta(.8);
+     //tcuts->SetDCA(.1);
+     tcuts->SetPt(1.2);
+     tcuts->SetPtTOFMatching(0.6);   
+     tcuts->SetQvecMin(QvecCutMin[iCut]);   
+     tcuts->SetQvecMax(QvecCutMax[iCut]);    
+     vcuts->SetCentralityCutMax(CentCutMax[iCut]);  
+     vcuts->SetCentralityCutMin(CentCutMin[iCut]);
+     task->SetEventCuts(vcuts);
+     task->SetTrackCuts(tcuts);
+     task->SetNSigmaForIdentification(3.); // FIXME
+     task->SetYCut(.5);
+     vcuts->PrintCuts();
+     tcuts->PrintCuts();
+     
+     // check for MC or real data
+     if (c == 2 || c == 3)
+       {
+        task->SetIsMC(kTRUE);
+        AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+        AliAnalysisDataContainer *coutputpt1 = mgr->CreateContainer(Form("chistpt%d",iCut), AliSpectraAODHistoManager::Class(),  AliAnalysisManager::kOutputContainer, 
+                                                                    Form("Pt.AOD.1._MC_Cent%.0fto%.0f_QVec%.1fto%.1f.root",CentCutMin[iCut],CentCutMax[iCut],QvecCutMin[iCut],QvecCutMax[iCut]));
+        AliAnalysisDataContainer *coutputpt2 = mgr->CreateContainer(Form("cvcutpt%d",iCut), AliSpectraAODEventCuts::Class(),    AliAnalysisManager::kOutputContainer, 
+                                                                    Form("Pt.AOD.1._MC_Cent%.0fto%.0f_QVec%.1fto%.1f.root",CentCutMin[iCut],CentCutMax[iCut],QvecCutMin[iCut],QvecCutMax[iCut]));
+        AliAnalysisDataContainer *coutputpt3 = mgr->CreateContainer(Form("ctcutpt%d",iCut), AliSpectraAODTrackCuts::Class(),     AliAnalysisManager::kOutputContainer, 
+                                                                    Form("Pt.AOD.1._MC_Cent%.0fto%.0f_QVec%.1fto%.1f.root",CentCutMin[iCut],CentCutMax[iCut],QvecCutMin[iCut],QvecCutMax[iCut]));
+       }
+     if (c == 1)
+       {
+        AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+        AliAnalysisDataContainer *coutputpt1 = mgr->CreateContainer(Form("chistpt%d",iCut), AliSpectraAODHistoManager::Class(),  AliAnalysisManager::kOutputContainer, 
+                                                                    Form("Pt.AOD.1._data_ptcut_Cent%.0fto%.0f_QVec%.1fto%.1f.root",CentCutMin[iCut],CentCutMax[iCut],QvecCutMin[iCut],QvecCutMax[iCut]));
+        AliAnalysisDataContainer *coutputpt2 = mgr->CreateContainer(Form("cvcutpt%d",iCut), AliSpectraAODEventCuts::Class(),    AliAnalysisManager::kOutputContainer, 
+                                                                    Form("Pt.AOD.1._data_ptcut_Cent%.0fto%.0f_QVec%.1fto%.1f.root",CentCutMin[iCut],CentCutMax[iCut],QvecCutMin[iCut],QvecCutMax[iCut]));
+        AliAnalysisDataContainer *coutputpt3 = mgr->CreateContainer(Form("ctcutpt%d",iCut), AliSpectraAODTrackCuts::Class(),     AliAnalysisManager::kOutputContainer, 
+                                                                    Form("Pt.AOD.1._data_ptcut_Cent%.0fto%.0f_QVec%.1fto%.1f.root",CentCutMin[iCut],CentCutMax[iCut],QvecCutMin[iCut],QvecCutMax[iCut]));
+        
+       }
+     mgr->ConnectInput(task, 0, cinput);
+     mgr->ConnectOutput(task, 1, coutputpt1);
+     mgr->ConnectOutput(task, 2, coutputpt2);
+     mgr->ConnectOutput(task, 3, coutputpt3);
+   }
    mgr->SetDebugLevel(2);
-
+   
    if (!mgr->InitAnalysis()) return;
    mgr->PrintStatus();
-   mgr->StartAnalysis("proof");
+   if (c == 3)
+   {
+      mgr->StartAnalysis("local");
+   }
+   if (c != 3)
+   {
+      mgr->StartAnalysis("proof");
+   }
 }