]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/AliAnalysisTaskJetCorePP.cxx
AliAODEvent::GetHeader() returns AliVHeader
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCorePP.cxx
index 4d23c729ac785ecdb07d872142ec98ddbdcdf924..678856c00ebe8e65d63a7e0dad59f45f8c3dbb6f 100644 (file)
@@ -62,6 +62,7 @@
 #include "AliAODJet.h"
 #include "AliVVertex.h"
 #include "AliAnalysisTaskJetCorePP.h"
+#include "AliHeader.h" //KF//
 
 using std::cout;
 using std::endl;
@@ -136,6 +137,8 @@ fhDphiTriggerJet(0x0),
 fhDphiTriggerJetAccept(0x0),
 fhCentrality(0x0),
 fhCentralityAccept(0x0),
+fhNofMultipleTriggers(0x0),
+fhDeltaRMultTriggers(0x0),
 //fHJetPtRaw(0x0),
 //fHLeadingJetPtRaw(0x0), 
 //fHDphiVsJetPtAll(0x0), 
@@ -164,6 +167,8 @@ fhEntriesToMedian(0x0),
 fhEntriesToMedianGen(0x0),
 fhCellAreaToMedian(0x0),
 fhCellAreaToMedianGen(0x0),
+fhNofMultipleTriggersGen(0x0),
+fhDeltaRMultTriggersGen(0x0),
 fIsChargedMC(0),
 fIsKine(0),
 fIsFullMC(0),
@@ -265,6 +270,8 @@ fhDphiTriggerJet(0x0),
 fhDphiTriggerJetAccept(0x0),
 fhCentrality(0x0),
 fhCentralityAccept(0x0),
+fhNofMultipleTriggers(0x0),
+fhDeltaRMultTriggers(0x0),
 //fHJetPtRaw(0x0),
 //fHLeadingJetPtRaw(0x0), 
 //fHDphiVsJetPtAll(0x0), 
@@ -293,6 +300,8 @@ fhEntriesToMedian(0x0),
 fhEntriesToMedianGen(0x0),
 fhCellAreaToMedian(0x0),
 fhCellAreaToMedianGen(0x0),
+fhNofMultipleTriggersGen(0x0),
+fhDeltaRMultTriggersGen(0x0),
 fIsChargedMC(0),
 fIsKine(0),
 fIsFullMC(0),
@@ -330,7 +339,7 @@ fDoubleBinning(kFALSE)
    TString dummy(name);
    if(dummy.Contains("KINE")){
       DefineInput(1, TClonesArray::Class());
-      //XXXX//DefineInput(2, TClonesArray::Class());
+      DefineInput(2, TClonesArray::Class());
    }
 }
 
@@ -400,6 +409,8 @@ fhDphiTriggerJet(a.fhDphiTriggerJet),
 fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
 fhCentrality(a.fhCentrality),
 fhCentralityAccept(a.fhCentralityAccept),
+fhNofMultipleTriggers(a.fhNofMultipleTriggers),
+fhDeltaRMultTriggers(a.fhDeltaRMultTriggers),
 //fHJetPtRaw(a.fHJetPtRaw),
 //fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
 //fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
@@ -428,6 +439,8 @@ fhEntriesToMedian(a.fhEntriesToMedian),
 fhEntriesToMedianGen(a.fhEntriesToMedianGen),
 fhCellAreaToMedian(a.fhCellAreaToMedian),
 fhCellAreaToMedianGen(a.fhCellAreaToMedianGen),
+fhNofMultipleTriggersGen(a.fhNofMultipleTriggersGen),
+fhDeltaRMultTriggersGen(a.fhDeltaRMultTriggersGen),
 fIsChargedMC(a.fIsChargedMC),
 fIsKine(a.fIsKine),
 fIsFullMC(a.fIsFullMC),
@@ -493,49 +506,52 @@ Bool_t AliAnalysisTaskJetCorePP::Notify()
    //and number of trials from pyxsec.root
    //inspired by AliAnalysisTaskJetSpectrum2::Notify()
    if(!(fIsChargedMC || fIsKine)) return kFALSE; 
+   Float_t xsection = 0;
+   Float_t trials  = 1;
+   fAvgTrials = 1;
 
-   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
-   if(!fESD){
-      if(fDebug>1) AliError("ESD not available");
-      fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
-   } 
+   if(fIsChargedMC){ 
+      fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+      if(!fESD){
+         if(fDebug>1) AliError("ESD not available");
+         fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
+      } 
  
-   fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
+      fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 
 
-   if(fNonStdFile.Length()!=0){
-      // case that we have an AOD extension we can fetch the jets from the extended output
-      AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
-      fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
-      if(!fAODExtension){
-         if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
-      } 
-   }
+      if(fNonStdFile.Length()!=0){
+         // case that we have an AOD extension we can fetch the jets from the extended output
+         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+         fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
+         if(!fAODExtension){
+            if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
+         
+      }
  
-   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
-   Float_t xsection = 0;
-   Float_t ftrials  = 1;
+      TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 
-   fAvgTrials = 1;
-   if(tree){
-      TFile *curfile = tree->GetCurrentFile();
-      if(!curfile) {
-         Error("Notify","No current file");
-         return kFALSE;
-      }
-      if(!fh1Xsec || !fh1Trials){
-         Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
-         return kFALSE;
-      }
-      AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
-      fh1Xsec->Fill("<#sigma>",xsection);
-      // construct a poor man average trials
-      Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
-      if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries;
-      fh1Trials->Fill("#sum{ntrials}",ftrials);
-   }  
+      if(tree){
+         TFile *curfile = tree->GetCurrentFile();
+         if(!curfile) {
+            Error("Notify","No current file");
+            return kFALSE;
+         }
+         if(!fh1Xsec || !fh1Trials){
+            Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
+            return kFALSE;
+         }
+         AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
+         fh1Xsec->Fill("<#sigma>",xsection);
+         // construct a poor man average trials
+         Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
+         if(trials>=nEntries && nEntries>0.) fAvgTrials = trials/nEntries;
+         fh1Trials->Fill("#sum{ntrials}",trials);
+      }  
+
+      if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
+   }
 
-   if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
 
    return kTRUE;
 }
@@ -571,10 +587,6 @@ void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
 
    fRandom = new TRandom3(0);
 
-   //if(fIsKine){ //?????????????????????????????????????????
-   //}
-
-
    if(fIsChargedMC || fIsKine){   //full MC or pure kine
       fListJetsGen   = new TList(); //generator level charged antikt jets
       fListJetsBgGen = new TList(); //generator level jets to be removed from bg 
@@ -711,6 +723,9 @@ void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
    fhEntriesToMedian = new TH1D("fhEntriesToMedian","fhEntriesToMedian",30,0,30);
    fhCellAreaToMedian =  new TH1D("fhCellAreaToMedian", "fhCellAreaToMedian", 75,0,1.5);
 
+   fhNofMultipleTriggers = new TH1D("fhNofMultipleTriggers","fhNofMultipleTriggers",100,0,100);
+   fhDeltaRMultTriggers = new  TH1D("fhDeltaRMultTriggers","fhDeltaRMultTriggers", 100,0,4);  
+
    if(!fIsKine){
       fOutputList->Add(fhJetPhi);
       fOutputList->Add(fhTriggerPhi);
@@ -728,6 +743,8 @@ void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
       fOutputList->Add(fhCentralityAccept);
       fOutputList->Add(fhEntriesToMedian);
       fOutputList->Add(fhCellAreaToMedian);
+      fOutputList->Add(fhNofMultipleTriggers);
+      fOutputList->Add(fhDeltaRMultTriggers);
    }
    // raw spectra of INCLUSIVE jets  
    //Centrality, pTjet, A
@@ -855,6 +872,13 @@ void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
       fhCellAreaToMedianGen  = (TH1D*) fhCellAreaToMedian->Clone("fhCellAreaToMedianGen");
       fhCellAreaToMedianGen->SetTitle(Form("%s Gen MC", fhCellAreaToMedian->GetTitle())); 
       fOutputList->Add(fhCellAreaToMedianGen);
+
+      fhNofMultipleTriggersGen = (TH1D*) fhNofMultipleTriggers->Clone("fhNofMultipleTriggersGen"); 
+      fOutputList->Add(fhNofMultipleTriggersGen);
+
+      fhDeltaRMultTriggersGen  = (TH1D*) fhDeltaRMultTriggers->Clone("fhDeltaRMultTriggersGen");
+      fOutputList->Add(fhDeltaRMultTriggersGen);
+
    }
    //-------------------------------------
    //     pythia histograms
@@ -906,6 +930,9 @@ void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
 
 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
 {
+   //User Exec
+   
+
    //Event loop
    Double_t eventW  = 1.0;
    Double_t ptHard  = 0.0;
@@ -932,26 +959,43 @@ void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
          return;
       }
    }else{  //Kine
-      //XXXX// if(!strlen(fJetBranchNameBgKine.Data())){
-      //XXXX//   AliError("Name of jet bg branch for kine not set.");
-      //XXXX//   return;
-      //XXXX// }
+       if(!strlen(fJetBranchNameBgKine.Data())){
+         AliError("Name of jet bg branch for kine not set.");
+         return;
+       }
 
       Init(); 
       if(fMcHandler){
          fMcEvent = fMcHandler->MCEvent(); 
       }else{
-         if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Handler() fMcHandler=NULL\n");
+         if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcHandler=NULL\n");
          PostData(1, fOutputList);
          return;
       } 
       if(!fMcEvent){
-         if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Exec()   fMcEvent=NULL \n");
+         if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcEvent=NULL \n");
          PostData(1, fOutputList);
          return;
       }
+      Float_t xsection = 0;
+      Float_t trials  = 0;
+
+      AliGenPythiaEventHeader *genPH =
+         dynamic_cast<AliGenPythiaEventHeader*> (fMcEvent->GenEventHeader()); 
+      if(genPH){
+         xsection = genPH->GetXsection();
+         trials   = genPH->Trials();
+         ptHard   = genPH->GetPtHard();
+      }
+      fh1Xsec->Fill("<#sigma>",xsection);
+      fh1Trials->Fill("#sum{ntrials}",trials);
+      fh1PtHard->Fill(ptHard,eventW);
+      fh1PtHardNoW->Fill(ptHard,1);
+      fh1PtHardTrials->Fill(ptHard,trials);
    }
 
+
    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
    if(!fESD){
       if(fDebug>1) AliError("ESD not available");
@@ -1058,7 +1102,7 @@ void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
          cent = fESD->GetCentrality();
          if(cent) centValue = cent->GetCentralityPercentile("V0M");
       }else{
-         centValue = aod->GetHeader()->GetCentrality();
+         centValue = ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
       }   
       if(fDebug) printf("centrality: %f\n", centValue);
       //Input events
@@ -1076,7 +1120,10 @@ void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
  
    //-----------------select disjunct event subsamples ----------------
    if(!fIsKine){ //reconstructed data
-      Int_t eventnum  = aod->GetHeader()->GetEventNumberESDFile();
+      AliAODHeader * header = dynamic_cast<AliAODHeader*>(aod->GetHeader());
+      if(!header) AliFatal("Not a standard AOD");
+
+      Int_t eventnum  = header->GetEventNumberESDFile();
       Int_t lastdigit = eventnum % 10;
       if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
          fHistEvtSelection->Fill(5);
@@ -1114,14 +1161,13 @@ void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
    //==============  analyze generator level MC  ================ 
    TList particleListGen; //list of tracks in MC
 
-
    if(fIsChargedMC || fIsKine){
 
       if(fIsKine){  //pure kine
 
          //================= generated charged antikt jets ================
          ReadTClonesArray(fJetBranchNameKine.Data(),   fListJetsGen); 
-         //XXXX//ReadTClonesArray(fJetBranchNameBgKine.Data(), fListJetsBgGen); 
+         ReadTClonesArray(fJetBranchNameBgKine.Data(), fListJetsBgGen); 
       }else{  
          //================= generated charged antikt jets ================
          ReadTClonesArray(fJetBranchNameChargMC.Data(),   fListJetsGen); 
@@ -1228,7 +1274,6 @@ void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
       EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
 
       //Estimate rho from cell median minus jets
-      if(!fIsKine)  //XXXX//
       EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
 
       //============  Generator trigger+jet ==================
@@ -1236,6 +1281,25 @@ void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
          if(ntriggersMC>0){ //there is at least one trigger 
             Int_t rnd     = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
             indexTriggGen = triggersMC[rnd];
+
+            fhNofMultipleTriggersGen->Fill(ntriggersMC-1);
+            if(ntriggersMC>1){
+               Double_t deltaPhi, deltaEta, deltaR;
+               for(Int_t ia=0; ia<ntriggersMC-1; ia++){
+                  AliVParticle* tGenI = (AliVParticle*) particleListGen.At(ia);  
+                  if(!tGenI) continue;
+                  for(Int_t ib=ia+1; ib<ntriggersMC; ib++){
+                     AliVParticle* tGenII = (AliVParticle*) particleListGen.At(ib);  
+                     if(!tGenII) continue;
+                        
+                        deltaPhi = RelativePhi(tGenI->Phi(),tGenII->Phi());
+                        deltaEta = tGenI->Eta()-tGenII->Eta(); 
+                        deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
+                        fhDeltaRMultTriggersGen->Fill(deltaR);
+                  }
+               }
+            } 
+
          }else{
             indexTriggGen = -1; //trigger not found
          }
@@ -1723,6 +1787,7 @@ Int_t  AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
    Double_t ptmax = -10;
    Int_t    triggers[200];
    Int_t    ntriggers = 0; //index in triggers array
 
    for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
       AliAODTrack *tr = aodevt->GetTrack(it);
@@ -1755,6 +1820,23 @@ Int_t  AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
    if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
       Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
       index     = triggers[rnd];
+
+      fhNofMultipleTriggers->Fill(ntriggers-1);
+      if(ntriggers>1){
+         Double_t deltaPhi, deltaEta, deltaR;
+         for(Int_t ia=0; ia<ntriggers-1; ia++){
+            AliVParticle* tGeni = (AliVParticle*) list->At(ia);  
+            if(!tGeni) continue;
+            for(Int_t ib=ia+1; ib<ntriggers; ib++){
+               AliVParticle* tGenii = (AliVParticle*) list->At(ib);  
+               if(!tGenii) continue;
+               deltaPhi = RelativePhi(tGeni->Phi(),tGenii->Phi());
+               deltaEta = tGeni->Eta()-tGenii->Eta(); 
+               deltaR   = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
+               fhDeltaRMultTriggers->Fill(deltaR);
+            }
+         }
+      }
    }
 
    return index;
@@ -2061,9 +2143,9 @@ void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
       if(fJetBranchNameKine.Length()>0 && bname.CompareTo(fJetBranchNameKine.Data())==0){
          array = dynamic_cast<TClonesArray*>(GetInputData(1)); //connect exchange container slot 1
       }
-      //XXXX// if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
-      //XXXX//   array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
-      //XXXX// }
+      if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
+        array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
+      }
    }else{ //take input from AOD
       if(fAODOut&&!array){
          array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));