]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/JetTasks/AliPWG4HighPtSpectra.cxx
New analysis devoted to shower shape studies
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliPWG4HighPtSpectra.cxx
index e0add9661df035e211056ddf53b77cfbaa806717..dab1d9f2084ad74974ba6946b33c408173c6bb7a 100644 (file)
@@ -65,7 +65,7 @@ AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpe
   fCFManagerNeg(0x0),
   fESD(0),
   fTrackCuts(0),
-  fTrigger(0),
+  fTrackCutsTPConly(0),
   fHistList(0),
   fNEventAll(0),
   fNEventSel(0)
@@ -82,7 +82,7 @@ AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
   fCFManagerNeg(0x0),
   fESD(0),
   fTrackCuts(),
-  fTrigger(0),
+  fTrackCutsTPConly(0),
   fHistList(0),
   fNEventAll(0),
   fNEventSel(0)
@@ -90,61 +90,29 @@ AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
   //
   // Constructor. Initialization of Inputs and Outputs
   //
-  AliDebug(2,Form("AliPWG4HighPtSpectra","Calling Constructor"));
+  AliDebug(2,Form("AliPWG4HighPtSpectra Calling Constructor"));
   // Input slot #0 works with a TChain ESD
   DefineInput(0, TChain::Class());
+  // Output slot #0 writes into a TList
   DefineOutput(0,TList::Class());
+  // Output slot #1, #2 writes into a AliCFContainer
   DefineOutput(1,AliCFContainer::Class());
   DefineOutput(2,AliCFContainer::Class());
+  // Output slot #3 writes into a AliESDtrackCuts
+  DefineOutput(3, AliESDtrackCuts::Class());
+  DefineOutput(4, AliESDtrackCuts::Class());
+}
+
+//________________________________________________________________________
+void AliPWG4HighPtSpectra::LocalInit()
+{
+  //
+  // Only called once at beginning
+  //
+  PostData(3,fTrackCuts);
+  PostData(4,fTrackCutsTPConly);
 }
 
-// //___________________________________________________________________________
-// AliPWG4HighPtSpectra& AliPWG4HighPtSpectra::operator=(const AliPWG4HighPtSpectra& c) 
-// {
-//   //
-//   // Assignment operator
-//   //
-//   if (this!=&c) {
-//     AliAnalysisTask::operator=(c) ;
-//     fReadAODData = c.fReadAODData ;
-//     fCFManagerPos  = c.fCFManagerPos;
-//     fCFManagerNeg  = c.fCFManagerNeg;
-//     fHistList = c.fHistList;
-//     fNEventAll = c.fNEventAll;
-//     fNEventSel = c.fNEventSel;
-//   }
-//   return *this;
-// }
-
-// //___________________________________________________________________________
-// AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const AliPWG4HighPtSpectra& c) :
-//   AliAnalysisTask(c),
-//   fReadAODData(c.fReadAODData),
-//   fCFManagerPos(c.fCFManagerPos),
-//   fCFManagerNeg(c.fCFManagerNeg),
-//   fESD(c.fESD),
-//   fTrackCuts(c.fTrackCuts),
-//   fTrigger(c.fTrigger),
-//   fHistList(c.fHistList),
-//   fNEventAll(c.fNEventAll),
-//   fNEventSel(c.fNEventSel)
-// {
-//   //
-//   // Copy Constructor
-//   //
-// }
-
-// //___________________________________________________________________________
-// AliPWG4HighPtSpectra::~AliPWG4HighPtSpectra() {
-//   //
-//   //destructor
-//   //
-//   Info("~AliPWG4HighPtSpectra","Calling Destructor");
-//   if (fCFManagerPos)           delete fCFManagerPos ;
-//   if (fCFManagerNeg)           delete fCFManagerNeg ;
-//   if (fNEventAll)              delete fNEventAll ;
-//   if (fNEventSel)              delete fNEventSel ;
-// }
 //________________________________________________________________________
 void AliPWG4HighPtSpectra::ConnectInputData(Option_t *) 
 {
@@ -188,8 +156,8 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
     return;
   }
 
-  Bool_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
-  if(!isSelected) { //Select collison candidates
+  UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+  if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
     PostData(0,fHistList);
     PostData(1,fCFManagerPos->GetParticleContainer());
@@ -201,7 +169,6 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
   // This handler can return the current MC event
   
   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-  //  AliMCEventHandler* eventHandler = (AliMCEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
   
   AliStack* stack = 0x0;
   AliMCEvent* mcEvent = 0x0;
@@ -212,8 +179,8 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
       AliDebug(2,Form("ERROR: Could not retrieve MC event"));
       PostData(0,fHistList);
       PostData(1,fCFManagerPos->GetParticleContainer());
-    PostData(2,fCFManagerNeg->GetParticleContainer());
-    return;
+      PostData(2,fCFManagerNeg->GetParticleContainer());
+      return;
     }
     
     AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));
@@ -226,11 +193,18 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
   // Need vertex cut
-  if (vtx->GetNContributors() < 2) {
-    PostData(0,fHistList);
-    PostData(1,fCFManagerPos->GetParticleContainer());
-    PostData(2,fCFManagerNeg->GetParticleContainer());
-    return;
+  TString vtxName(vtx->GetName());
+  if(vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
+    // SPD vertex
+    vtx = fESD->GetPrimaryVertexSPD();
+    if(vtx->GetNContributors()<2) {
+      vtx = 0x0;
+      // Post output data
+      PostData(0,fHistList);
+      PostData(1,fCFManagerPos->GetParticleContainer());
+      PostData(2,fCFManagerNeg->GetParticleContainer());
+      return;
+    }
   }
   
   double primVtx[3];
@@ -264,9 +238,12 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
   fNEventSel->Fill(0.);
   
   
-  Double_t containerInputRec[5] ;
-  Double_t containerInputTPConly[5];
-  Double_t containerInputMC[5];
+  Double_t containerInputRec[3]       = {0.,0.,0.};
+  Double_t containerInputTPConly[3]   = {0.,0.,0.};
+  Double_t containerInputMC[3]        = {0.,0.,0.};
+  Double_t containerInputRecMC[3]     = {0.,0.,0.};
+  Double_t containerInputTPConlyMC[3] = {0.,0.,0.};
+
   //Now go to rec level
   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++) 
     {   
@@ -276,60 +253,72 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
       AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)track->GetTPCInnerParam();
       if(!track || !trackTPC) continue;
 
-      Float_t dca2D, dcaZ;
-      track->GetImpactParameters(dca2D,dcaZ);
-      Float_t dca2DTPC, dcaZTPC;
-      track->GetImpactParametersTPC(dca2DTPC,dcaZTPC); 
-      Float_t chi2PerClusterTPC = -1.;
-      Float_t nClustersTPC = track->GetTPCNcls();//track->GetTPCclusters(0);
-      if(nClustersTPC>0.) chi2PerClusterTPC = track->GetTPCchi2()/(2.*nClustersTPC-5.);
-      Float_t chi2PerClusterTPCIter1 = -1.;
-      Float_t nClustersTPCIter1 = track->GetTPCNclsIter1();   
-      if(nClustersTPCIter1>0.) chi2PerClusterTPCIter1 = track->GetTPCchi2Iter1()/(2.*nClustersTPCIter1-5.);
-
       //fill the container
       containerInputRec[0] = track->Pt();
       containerInputRec[1] = track->Phi();
       containerInputRec[2] = track->Eta();
-      containerInputRec[3] = dca2D;
-      containerInputRec[4] = chi2PerClusterTPC;
-
+    
+      //Store TPC Inner Params for TPConly tracks
       containerInputTPConly[0] = trackTPC->Pt();
       containerInputTPConly[1] = trackTPC->Phi();
       containerInputTPConly[2] = trackTPC->Eta();
-      containerInputTPConly[3] = dca2DTPC/10.; //Divide by 10 in order to store in same containter. Should be corrected back when looking at output.
-      containerInputTPConly[4] = chi2PerClusterTPCIter1;//TPC;
 
-      if (fTrackCuts->AcceptTrack(track)) {
-       if(track->GetSign()>0.) {
-         fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
-         fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
-       }
-       if(track->GetSign()<0.) {
-         fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
-         fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
+      AliESDtrack* trackTPCESD = fTrackCutsTPConly->GetTPCOnlyTrack(fESD, iTrack);
+      if(trackTPCESD) {
+       if (fTrackCutsTPConly->AcceptTrack(trackTPCESD)) {
+         if(trackTPC->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
+         if(trackTPC->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
+
+         //Only fill the MC containers if MC information is available
+         if(eventHandler) {
+           Int_t label = TMath::Abs(track->GetLabel());
+           TParticle *particle = stack->Particle(label) ;
+           if(!particle) continue;
+           
+           containerInputTPConlyMC[0] = particle->Pt();      
+           containerInputTPConlyMC[1] = particle->Phi();      
+           containerInputTPConlyMC[2] = particle->Eta();  
+           
+           //Container with primaries
+           if(stack->IsPhysicalPrimary(label)) {
+             if(particle->GetPDG()->Charge()>0.) {
+               fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedTPCOnlyMC);
+             }
+             if(particle->GetPDG()->Charge()<0.) {
+               fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedTPCOnlyMC);
+             }
+           }
+         }
+         
        }
-       
+      }
+
+      if (fTrackCuts->AcceptTrack(track)) {
+       if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
+       if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
+
        
-       //Only fill the secondary particle container if MC information is available
+       //Only fill the MC containers if MC information is available
        if(eventHandler) {
          Int_t label = TMath::Abs(track->GetLabel());
          TParticle *particle = stack->Particle(label) ;
          if(!particle) continue;
 
-         containerInputMC[0] = particle->Pt();      
-         containerInputMC[1] = particle->Phi();      
-         containerInputMC[2] = particle->Eta();  
-         containerInputMC[3] = 0.0;      
-         containerInputMC[4] = 0.0;  
+         containerInputRecMC[0] = particle->Pt();      
+         containerInputRecMC[1] = particle->Phi();      
+         containerInputRecMC[2] = particle->Eta();  
 
-         if(particle->GetPDG()->Charge()>0.) {
-           fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
-         }
-         if(particle->GetPDG()->Charge()<0.) {
-           fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
+         //Container with primaries
+         if(stack->IsPhysicalPrimary(label)) {
+           if(particle->GetPDG()->Charge()>0.) {
+             fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
+           }
+           if(particle->GetPDG()->Charge()<0.) {
+             fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
+           }
          }
 
+         //Container with secondaries
          if (!stack->IsPhysicalPrimary(label) ) {
            if(particle->GetPDG()->Charge()>0.) {
              fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
@@ -340,34 +329,26 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
          }
        }
        
-      }
-      
-    }
+      }//trackCuts
+
+      delete trackTPCESD;
+    }//track loop
   
 
   //Fill MC containters if particles are findable
   if(eventHandler) {
-    for(int iPart = 1; iPart<(mcEvent->GetNumberOfTracks()); iPart++)//stack->GetNprimary();
+    for(int iPart = 1; iPart<(mcEvent->GetNumberOfPrimaries()); iPart++)//stack->GetNprimary();
       {
        AliMCParticle *mcPart  = (AliMCParticle*)mcEvent->GetTrack(iPart);
        if(!mcPart) continue;
-       
        //fill the container
        containerInputMC[0] = mcPart->Pt();
        containerInputMC[1] = mcPart->Phi();      
        containerInputMC[2] = mcPart->Eta();  
-       containerInputMC[3] = 0.0;
-       containerInputMC[4] = 0.0;
-
-       int counter;
-       Float_t trackLengthTPC = 0.;
-       if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(3,mcPart)) {
-         trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
-         if(trackLengthTPC > (float)fTrackCuts->GetMinNClusterTPC()) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
-       }
-       if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(3,mcPart)) {
-         trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
-         if(trackLengthTPC > (float)fTrackCuts->GetMinNClusterTPC()) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
+
+       if(stack->IsPhysicalPrimary(iPart)) {
+         if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
+         if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
        }
       }
   }
@@ -385,6 +366,7 @@ void AliPWG4HighPtSpectra::Terminate(Option_t*)
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
+
 }
 
 //___________________________________________________________________________
@@ -392,7 +374,7 @@ void AliPWG4HighPtSpectra::CreateOutputObjects() {
   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
   //TO BE SET BEFORE THE EXECUTION OF THE TASK
   //
-  AliDebug(2,Form("CreateOutputObjects","CreateOutputObjects of task %s", GetName()));
+  AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
 
   Bool_t oldStatus = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kFALSE);