]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/AliPWG4HighPtSpectra.cxx
ghost area to 0.005 by default
[u/mrichter/AliRoot.git] / PWGJE / AliPWG4HighPtSpectra.cxx
index 631be43cf315568651f0c495d85edcc0dc84e1e7..5d440599280679e0c1cad3b56394a1d06a945435 100644 (file)
@@ -54,6 +54,7 @@
 #include "AliMCEventHandler.h"
 #include "AliCFContainer.h"
 #include "AliGenPythiaEventHeader.h"
+#include "AliGenHijingEventHeader.h"
 #include "AliGenCocktailEventHeader.h"
 //#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
 
@@ -71,11 +72,13 @@ AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpe
   fMC(0x0),
   fStack(0x0),
   fVtx(0x0),
+  fTriggerMask(AliVEvent::kMB),
   fIsPbPb(0),
   fCentClass(10),
   fTrackType(0),
   fTrackCuts(0x0),
   fTrackCutsReject(0x0),
+  fbSelectHIJING(kFALSE),
   fSigmaConstrainedMax(100.),
   fAvgTrials(1),
   fHistList(0),
@@ -104,11 +107,13 @@ AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
   fMC(0x0),
   fStack(0x0),
   fVtx(0x0),
+  fTriggerMask(AliVEvent::kMB),
   fIsPbPb(0),
   fCentClass(10),
   fTrackType(0),
   fTrackCuts(0x0),
   fTrackCutsReject(0x0),
+  fbSelectHIJING(kFALSE),
   fSigmaConstrainedMax(100.),
   fAvgTrials(1),
   fHistList(0),
@@ -203,7 +208,7 @@ Bool_t AliPWG4HighPtSpectra::SelectEvent() {
 
   //Trigger
   UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
-  if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
+  if(fTriggerMask != AliVEvent::kAny && !(isSelected&fTriggerMask)) { //Select collison candidates
     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
     fNEventReject->Fill("Trigger",1);
     selectEvent = kFALSE;
@@ -307,7 +312,6 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
   fNEventAll->Fill(0.);
 
   if(!SelectEvent()) {
-    fNEventReject->Fill("NTracks<2",1);
     // Post output data
     PostData(0,fHistList);
     PostData(1,fCFManagerPos->GetParticleContainer());
@@ -347,117 +351,179 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
   //Now go to rec level
   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++) 
     {   
+   
+
       //Get track for analysis
       AliESDtrack *track = 0x0;
       AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
-      if(!esdtrack) continue;
+      if(!esdtrack)
+       continue;
+      
 
-      if(fTrackType==1) {
-       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
-       if(!track) continue;
+      if(fTrackType==4) {
+       if (!(fTrackCuts->AcceptTrack(esdtrack)))
+         continue;
       }
-      else if(fTrackType==2) {
-       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
-       if(!track) continue;
 
+      if(fTrackType==1)
+       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
+      else if(fTrackType==2 || fTrackType==4) {
+       track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(fESD),esdtrack->GetID());
+       if(!track)
+         continue;
+       
        AliExternalTrackParam exParam;
        Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
        if( !relate ) {
-         delete track;
+         if(track) delete track;
          continue;
        }
        track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
       }
+      else if(fTrackType==5 || fTrackType==6) {
+       if(fTrackCuts->AcceptTrack(esdtrack)) {
+         continue;
+       }
+       else {
+         if( !(fTrackCutsReject->AcceptTrack(esdtrack)) && fTrackCuts->AcceptTrack(esdtrack) ) {
+
+           if(fTrackType==5) {
+             //use TPConly constrained track
+             track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
+             if(!track)
+               continue;
+             
+             AliExternalTrackParam exParam;
+             Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
+             if( !relate ) {
+               if(track) delete track;
+               continue;
+             }
+             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
+           }
+           else if(fTrackType==6) {
+             //use global constrained track
+             track = new AliESDtrack(*esdtrack);
+             track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
+
+           }
+         }
+       }
+      }
       else if(fTrackType==7) {
        //use global constrained track
        track = new AliESDtrack(*esdtrack);
-       //      track = esdtrack;
-       //      track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
-
       }
       else
        track = esdtrack;
     
-      if(!track) continue;
-      if(fTrackType==2) {
+      if(!track)
+       continue;
+      
+
+      if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
        //Cut on chi2 of constrained fit
-       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
-         delete track;
+       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
+         if(track) delete track;
          continue;
        }
       }
+      if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
+       if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
+         if(track) delete track;
+       }
+       continue;
+      }
 
-      if (fTrackCuts->AcceptTrack(track)) {
-
-       if(fTrackType==7) {
-         if(fTrackCutsReject ) {
-           if(fTrackCutsReject->AcceptTrack(track) ) {
-             if(track) delete track;
-             continue;
-           }
+      if(fTrackType==7) {
+       if(fTrackCutsReject ) {
+         if(fTrackCutsReject->AcceptTrack(track) ) {
+           if(track) delete track;
+           continue;
          }
-         
-         if(esdtrack->GetConstrainedParam()) 
-           track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
        }
+      
+       if(esdtrack->GetConstrainedParam()) 
+         track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
+      }
 
-       //fill the container
-       containerInputRec[0] = track->Pt();
-       containerInputRec[1] = track->Phi();
-       containerInputRec[2] = track->Eta();
-       containerInputRec[3] = track->GetTPCNclsIter1();
+      if(!track) {
+       if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
+         if(track) delete track;
+       }
+       continue;
+      }
 
-       if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
-       if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
+      //fill the container
+      containerInputRec[0] = track->Pt();
+      containerInputRec[1] = track->Phi();
+      containerInputRec[2] = track->Eta();
+      containerInputRec[3] = track->GetTPCNclsIter1();
+
+      if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
+      if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 
        
-       //Only fill the MC containers if MC information is available
-       if(fMC) {
-         Int_t label = TMath::Abs(track->GetLabel());
-         TParticle *particle = fStack->Particle(label) ;
-         if(!particle) {
+      //Only fill the MC containers if MC information is available
+      if(fMC) {
+       Int_t label = TMath::Abs(track->GetLabel());
+       if(label>fStack->GetNtrack()) {
+         if(fTrackType==1 || fTrackType==2 || fTrackType==7)
+           delete track;
+         continue;
+       }
+       //Only select particles generated by HIJING if requested
+       if(fbSelectHIJING) {
+         if(!IsHIJINGParticle(label)) {
            if(fTrackType==1 || fTrackType==2 || fTrackType==7)
              delete track;
            continue;
          }
-         containerInputRecMC[0] = particle->Pt();      
-         containerInputRecMC[1] = particle->Phi();      
-         containerInputRecMC[2] = particle->Eta();  
-         containerInputRecMC[3] = track->GetTPCNclsIter1();
-
-         //Container with primaries
-         if(fStack->IsPhysicalPrimary(label)) {
-           if(particle->GetPDG()->Charge()>0.) {
-             fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
-           }
-           if(particle->GetPDG()->Charge()<0.) {
-             fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
-           }
+       }
+       TParticle *particle = fStack->Particle(label) ;
+       if(!particle) {
+         if(fTrackType==1 || fTrackType==2 || fTrackType==7)
+           delete track;
+         continue;
+       }
+       containerInputRecMC[0] = particle->Pt();      
+       containerInputRecMC[1] = particle->Phi();      
+       containerInputRecMC[2] = particle->Eta();  
+       containerInputRecMC[3] = track->GetTPCNclsIter1();
+
+       //Container with primaries
+       if(fStack->IsPhysicalPrimary(label)) {
+         if(particle->GetPDG()->Charge()>0.) {
+           fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
+         }
+         if(particle->GetPDG()->Charge()<0.) {
+           fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
+         }
 
-           //Fill pT resolution plots for primaries
-           fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
+         //Fill pT resolution plots for primaries
+         fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
 
-         }
+       }
 
-         //Container with secondaries
-         if (!fStack->IsPhysicalPrimary(label) ) {
-           if(particle->GetPDG()->Charge()>0.) {
-             fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
-           }
-           if(particle->GetPDG()->Charge()<0.) {
-             fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
-           }
-           //Fill pT resolution plots for primaries
-           fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
+       //Container with secondaries
+       if (!fStack->IsPhysicalPrimary(label) ) {
+         if(particle->GetPDG()->Charge()>0.) {
+           fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
+         }
+         if(particle->GetPDG()->Charge()<0.) {
+           fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
          }
+         //Fill pT resolution plots for primaries
+         fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
        }
+      }
        
-      }//trackCuts global tracks
 
-      if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
+      if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
        if(track) delete track;
       }
+
+
     }//track loop
   
 
@@ -467,6 +533,12 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
       AliMCParticle *mcPart  = (AliMCParticle*)fMC->GetTrack(iPart);
       if(!mcPart) continue;
 
+      //Only select particles generated by HIJING if requested
+      if(fbSelectHIJING) {
+       if(!IsHIJINGParticle(iPart))
+         continue;
+      }
+
       Int_t pdg = mcPart->PdgCode();
       
       // select charged pions, protons, kaons , electrons, muons
@@ -521,7 +593,7 @@ Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fX
     // next trial fetch the histgram file
     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
     if(!fxsec){
-       // not a severe condition but inciate that we have no information
+      // not a severe condition but indicates that we have no information
       return kFALSE;
     }
     else{
@@ -618,6 +690,36 @@ AliGenPythiaEventHeader*  AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent
 
 }
 
+//________________________________________________________________________
+AliGenHijingEventHeader*  AliPWG4HighPtSpectra::GetHijingEventHeader(AliMCEvent *mcEvent){
+  
+  if(!mcEvent)return 0;
+  AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
+  AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
+  if(!hijingGenHeader){
+    // cocktail ??
+    AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
+    
+    if (!genCocktailHeader) {
+      //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
+      //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
+      return 0;
+    }
+    TList* headerList = genCocktailHeader->GetHeaders();
+    for (Int_t i=0; i<headerList->GetEntries(); i++) {
+      hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
+      if (hijingGenHeader)
+        break;
+    }
+    if(!hijingGenHeader){
+      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Hijing event header not found");
+      return 0;
+    }
+  }
+  return hijingGenHeader;
+
+}
+
 
 //___________________________________________________________________________
 void AliPWG4HighPtSpectra::Terminate(Option_t*)
@@ -706,4 +808,31 @@ void AliPWG4HighPtSpectra::CreateOutputObjects() {
 
 }
 
+//________________________________________________________________________
+Bool_t AliPWG4HighPtSpectra::IsHIJINGParticle(Int_t label) {
+  //
+  // Return kTRUE in case particle is from HIJING event
+  //
+
+  AliGenHijingEventHeader*  hijingHeader = GetHijingEventHeader(fMC);
+
+  Int_t nproduced = hijingHeader->NProduced();
+
+  TParticle * mom = fStack->Particle(label);
+  Int_t iMom = label;
+  Int_t iParent = mom->GetFirstMother();
+  while(iParent!=-1){
+    iMom = iParent;
+    mom = fStack->Particle(iMom);
+    iParent = mom->GetFirstMother();
+  }
+
+  if(iMom<nproduced) {
+    return kTRUE;
+  } else {
+    return kFALSE;
+  }
+
+}
+
 #endif