]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/AliAnalysisTaskJetCore.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
index d6ef8f9a8ac81df1cfd13d9c4ab9e2006bbb31bd..3562be28467ac0c6e663f2ab41c69a88f062e037 100644 (file)
@@ -1,4 +1,3 @@
-
 // ******************************************
 // This task computes several jet observables like 
 // the fraction of energy in inner and outer coronnas,
@@ -32,7 +31,7 @@
 #include "TH3F.h"
 #include "THnSparse.h"
 #include "TCanvas.h"
-
+#include "TRandom3.h"
 #include "AliLog.h"
 
 #include "AliAnalysisTask.h"
@@ -82,12 +81,20 @@ fCentMin(0.),
 fCentMax(100.),
 fNInputTracksMin(0),
 fNInputTracksMax(-1),
+fRequireITSRefit(0),
+fApplySharedClusterCut(0),
 fAngStructCloseTracks(0),
 fCheckMethods(0),
 fDoEventMixing(0), 
 fFlagPhiBkg(0),
 fFlagEtaBkg(0),
 fFlagJetHadron(0),
+fFrac(0.8),
+fTTLowRef(11),
+fTTUpRef(13),
+fTTLowSig(15),
+fTTUpSig(19),
+fHardest(0),
 fFlagRandom(0),
 fFlagOnlyRecoil(0),
 fFlagOnlyHardest(1),
@@ -143,6 +150,8 @@ fh2AngStructpt1C60(0x0),
 fh2AngStructpt2C60(0x0),
 fh2AngStructpt3C60(0x0),
 fh2AngStructpt4C60(0x0),
+fh1TrigRef(0x0),
+fh1TrigSig(0x0),
 fh2Ntriggers(0x0),
 fh2Ntriggers2C10(0x0),
 fh2Ntriggers2C20(0x0), 
@@ -154,8 +163,8 @@ fh2RPTC10(0x0),
 fh2RPTC20(0x0), 
 fHJetSpec(0x0),
 fhTTPt(0x0),
-fHJetPhiCorr(0x0)
-
+fHJetPhiCorr(0x0),
+fRandom(0x0)
  
 {
    // default Constructor
@@ -203,12 +212,20 @@ fCentMin(0.),
 fCentMax(100.),
 fNInputTracksMin(0),
 fNInputTracksMax(-1),
+fRequireITSRefit(0),
+fApplySharedClusterCut(0),
 fAngStructCloseTracks(0),
 fCheckMethods(0),
 fDoEventMixing(0),
 fFlagPhiBkg(0),
 fFlagEtaBkg(0),
 fFlagJetHadron(0),
+fFrac(0.8),
+fTTLowRef(11),
+fTTUpRef(13),
+fTTLowSig(15),
+fTTUpSig(19),
+fHardest(0),
 fFlagRandom(0),
 fFlagOnlyRecoil(0),
 fFlagOnlyHardest(1),
@@ -264,6 +281,8 @@ fh2AngStructpt1C60(0x0),
 fh2AngStructpt2C60(0x0),
 fh2AngStructpt3C60(0x0),
 fh2AngStructpt4C60(0x0),    
+fh1TrigRef(0x0),
+fh1TrigSig(0x0),
 fh2Ntriggers(0x0),
 fh2Ntriggers2C10(0x0),
 fh2Ntriggers2C20(0x0),
@@ -275,7 +294,8 @@ fh2RPTC10(0x0),
 fh2RPTC20(0x0), 
 fHJetSpec(0x0),
 fhTTPt(0x0),
-fHJetPhiCorr(0x0)
+fHJetPhiCorr(0x0),
+fRandom(0x0)
 
  {
    // Constructor
@@ -332,6 +352,12 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
    TH1::AddDirectory(kFALSE);
 
 
+        // set seed
+   fRandom = new TRandom3(0);
+  
+
+
+
    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
@@ -408,12 +434,13 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
     fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
 
    
-
-    fh2Ntriggers=new TH2F("# of triggers","",10,0.,100.,50,0.,50.);
+    fh1TrigRef=new TH1D("Trig Ref","",10,0.,10);
+    fh1TrigSig=new TH1D("Trig Sig","",10,0.,10);  
+    fh2Ntriggers=new TH2F("# of triggers","",100,0.,100.,50,0.,50.);
     fh2Ntriggers2C10=new TH2F("# of triggers2C10","",50,0.,50.,50,0.,50.);
     fh2Ntriggers2C20=new TH2F("# of triggers2C20","",50,0.,50.,50,0.,50.);
-    fh3JetDensity=new TH3F("Jet density vs mutliplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
-    fh3JetDensityA4=new TH3F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
+    fh3JetDensity=new TH3F("Jet density vs mutliplicity A>0.07","",100,0.,4000.,100,0.,5.,25,0.,50.);
+    fh3JetDensityA4=new TH3F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.,25,0.,50.);
     fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,3.5,100,0.,100.);
     fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,3.5,100,0.,100.); 
     fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,3.5,50,0.,50.); 
@@ -469,7 +496,8 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
 
 
 
+        fOutputList->Add(fh1TrigRef);
+        fOutputList->Add(fh1TrigSig); 
        fOutputList->Add(fh2Ntriggers);
         fOutputList->Add(fh2Ntriggers2C10);
         fOutputList->Add(fh2Ntriggers2C20); 
@@ -481,10 +509,27 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
         fOutputList->Add(fh2RPTC20);
 
         const Int_t dimSpec = 5;
-       const Int_t nBinsSpec[dimSpec]     = {10,100, 140, 50, fNRPBins};
+       const Int_t nBinsSpec[dimSpec]     = {100,6, 140, 50, fNRPBins};
        const Double_t lowBinSpec[dimSpec] = {0,0,-80, 0, 0};
-       const Double_t hiBinSpec[dimSpec]  = {100,1, 200, 50, fNRPBins};
+       const Double_t hiBinSpec[dimSpec]  = {100,1, 200, 50,  static_cast<Double_t>(fNRPBins)};
        fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum",dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
+
+             //change binning in jet area
+     Double_t *xPt6 = new Double_t[7];
+     xPt6[0] = 0.;
+     xPt6[1]=0.07;
+     xPt6[2]=0.2;
+     xPt6[3]=0.4;
+     xPt6[4]=0.6;
+     xPt6[5]=0.8; 
+     xPt6[6]=1;
+    fHJetSpec->SetBinEdges(1,xPt6);
+    delete [] xPt6;
+
+
+
+
+
        fOutputList->Add(fHJetSpec);  
 
 
@@ -557,7 +602,17 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    // -- event selection --
    fHistEvtSelection->Fill(1); // number of events before event selection
 
-   // physics selection
+
+       Bool_t selected=kTRUE;
+       selected = AliAnalysisHelperJetTasks::Selected();
+       if(!selected){
+      // no selection by the service task, we continue
+       PostData(1,fOutputList);
+       return;}
+    
+
+
+  // physics selection: this is now redundant, all should appear as accepted after service task selection
    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
         std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
@@ -568,6 +623,8 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
       return;
    }
 
+
+      
    // vertex selection
    if(!aod){
      if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
@@ -608,7 +665,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    if(fIsPbPb){
    if(fESD) {cent = fESD->GetCentrality();
      if(cent) centValue = cent->GetCentralityPercentile("V0M");}
-   else     centValue=aod->GetHeader()->GetCentrality();
+   else     centValue=((AliVAODHeader*)aod->GetHeader())->GetCentrality();
    
    if(fDebug) printf("centrality: %f\n", centValue);
       if (centValue < fCentMin || centValue > fCentMax){
@@ -642,9 +699,10 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
 
    if(fFlagRandom==0){
      if(externalBackground)rho = externalBackground->GetBackground(0);}
-   if(fFlagRandom==1){
+   if(fFlagRandom==2){
       if(externalBackground)rho = externalBackground->GetBackground(2);}
-
+   if(fFlagRandom==3){
+      if(externalBackground)rho = externalBackground->GetBackground(3);}
    // fetch jets
    TClonesArray *aodJets[2];
    aodJets[0]=0;
@@ -659,14 +717,29 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data()));  } 
 
 
-   //Double_t ptsub[aodJets[0]->GetEntriesFast()];
-   //Int_t inord[aodJets[0]->GetEntriesFast()];
-   //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
-   //  ptsub[n]=0;
-   //  inord[n]=0;}   
 
+   Int_t nT=0;
    TList ParticleList;
-   Int_t nT = GetListOfTracks(&ParticleList);
+   Double_t minT=0;
+   Double_t maxT=0;
+   Int_t number=0;
+   Double_t dice=fRandom->Uniform(0,1);
+   if(dice>fFrac){ minT=fTTLowRef;
+                   maxT=fTTUpRef;}
+   if(dice<=fFrac){minT=fTTLowSig;
+                   maxT=fTTUpSig;} 
+
+
+
+   if(fHardest==1 || fHardest==2) nT = GetListOfTracks(&ParticleList);
+   if(fHardest==0) nT=SelectTrigger(&ParticleList,minT,maxT,number);
+   if(nT<0){  
+   PostData(1, fOutputList);
+   return;}   
+
+      if(dice>fFrac) fh1TrigRef->Fill(number);
+      if(dice<=fFrac)fh1TrigSig->Fill(number);
+
      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
       fListJets[iJetType]->Clear();
       if (!aodJets[iJetType]) continue;
@@ -684,7 +757,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    Double_t phibig=0.;
    Double_t etasmall=0;
    Double_t ptsmall=0;
-   Double_t areasmall=0;
+   //   Double_t areasmall=0;
    Double_t phismall=0.;
          
   
@@ -692,28 +765,40 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    Int_t iCount=0; 
    Int_t trigJet=-1;
    Int_t trigBBTrack=-1;
-   Int_t trigInTrack=-1;
-   fRPAngle = aod->GetHeader()->GetEventplane();     
+          // Int_t trigInTrack=-1;
+   fRPAngle = ((AliVAODHeader*)aod->GetHeader())->GetEventplane();     
 
+   if(fHardest==0 || fHardest==1){
    AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);     
    if(!partback){  
    PostData(1, fOutputList);
    return;}
 
+    if(fSemigoodCorrect){
+   Double_t disthole=RelativePhi(partback->Phi(),fHolePos);
+   if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-0.6){ 
+   PostData(1, fOutputList);
+   return;}} 
+
+
+    }
+
+
+   for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
+     if(fHardest==0||fHardest==1){if(tt!=nT) continue;}
+     AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);     
+     if(!partback) continue;
+     if(partback->Pt()<8) continue;
 
-   //for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
-   //if(fFlagOnlyHardest!=0){if(tt!=nT) continue;}
-   //AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);     
    Double_t accep=2.*TMath::Pi()*1.8;
    Int_t injet4=0;
    Int_t injet=0; 
+
    if(fSemigoodCorrect){
    Double_t disthole=RelativePhi(partback->Phi(),fHolePos);
-   if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-0.6){ 
-   PostData(1, fOutputList);
-   return;}
+   if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-0.6) continue;}
 
-   }
+   
 
    fh2Ntriggers->Fill(centValue,partback->Pt());
    Double_t phiBinT = RelativePhi(partback->Phi(),fRPAngle);
@@ -748,8 +833,10 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
 
 
            if(fFlagJetHadron==0){
-           if(fFlagPhiBkg!=0) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
-           if(fFlagPhiBkg==0) if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;}
+           if(fFlagPhiBkg==1) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
+           if(fFlagPhiBkg==0) if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
+           if(fFlagPhiBkg==2) if(TMath::Abs(dphi)<TMath::Pi()-0.7) continue;
+           if(fFlagPhiBkg==3) if(TMath::Abs(dphi)<TMath::Pi()-0.5) continue;}
  
            if(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
 
@@ -765,7 +852,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
                    if(phitt<0)phitt+=TMath::Pi()*2.; 
                    Int_t phiBintt = GetPhiBin(phitt-fRPAngle);
 
-                  Double_t fillspec[] = {centValue,jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt(),phiBintt};
+                  Double_t fillspec[] = {centValue,jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt(), static_cast<Double_t>(phiBintt)};
                  fHJetSpec->Fill(fillspec);
            
           
@@ -804,7 +891,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
                       if(iCount==0){                        
                         trigJet=i;
                         trigBBTrack=nT;
-                        trigInTrack=ippt;
+                        //                      trigInTrack=ippt;
                         iCount=iCount+1;} 
                       
    
@@ -814,7 +901,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
                           etasmall  = jetsmall->Eta();
                           phismall = jetsmall->Phi();
                           ptsmall   = jetsmall->Pt();
-                          areasmall = jetsmall->EffectiveAreaCharged();
+                          //                      areasmall = jetsmall->EffectiveAreaCharged();
                           Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
                           tmpDeltaR=TMath::Sqrt(tmpDeltaR);
                           //Fraction in the jet core  
@@ -1042,7 +1129,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
 
     
 
-
+   }
 
 
 
@@ -1079,14 +1166,20 @@ Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
 
     
      for(int it = 0;it < aod->GetNumberOfTracks();++it){
-      AliAODTrack *tr = aod->GetTrack(it);
+      AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
+      if(!tr) AliFatal("Not a standard AOD");
       Bool_t bGood = false;
       if(fFilterType == 0)bGood = true;
       else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
       else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();    
-     if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
+      if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
+      if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
       if(bGood==false) continue;
-      if(TMath::Abs(tr->Eta())>0.9)continue;
+      if (fApplySharedClusterCut) {
+           Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
+           if (frac > 0.4) continue;
+      }
+     if(TMath::Abs(tr->Eta())>0.9)continue;
       if(tr->Pt()<0.15)continue;
       list->Add(tr);
       iCount++;
@@ -1125,6 +1218,71 @@ Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
  
 }
 
+
+
+Int_t  AliAnalysisTaskJetCore::SelectTrigger(TList *list,Double_t minT,Double_t maxT,Int_t &number){
+     Int_t iCount = 0;
+     AliAODEvent *aod = 0;
+     if(!fESD)aod = fAODIn;
+     else aod = fAODOut;   
+     if(!aod)return 0;
+     Int_t index=-1;
+     Int_t triggers[100];
+     for(Int_t cr=0;cr<100;cr++){triggers[cr]=-1;}
+     Int_t im=0;
+     for(int it = 0;it < aod->GetNumberOfTracks();++it){
+      AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
+      if(!tr) AliFatal("Not a standard AOD");
+      Bool_t bGood = false;
+      if(fFilterType == 0)bGood = true;
+      else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
+      else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
+      if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
+       if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
+      if(bGood==false) continue;
+      if (fApplySharedClusterCut) {
+           Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
+           if (frac > 0.4) continue;
+      }
+      if(TMath::Abs(tr->Eta())>0.9)continue;
+      if(tr->Pt()<0.15)continue;
+      list->Add(tr);
+      iCount++;
+      
+      if(tr->Pt()>=minT && tr->Pt()<maxT){
+       triggers[im]=iCount-1;
+        im=im+1;}
+
+     }
+      number=im;
+      Int_t rd=0;
+      if(im==0) rd=0;
+      if(im>0) rd=fRandom->Integer(im);
+      index=triggers[rd];
+
+     
+  
+   
+  
+      return index;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
  
     AliAODEvent *aod = 0;
@@ -1133,10 +1291,11 @@ Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
     Int_t index=-1;
     Double_t ptmax=-10;
     Double_t dphi=0;
-    Double_t dif=0;
+    // Double_t dif=0;
     Int_t iCount=0;
     for(int it = 0;it < aod->GetNumberOfTracks();++it){
-      AliAODTrack *tr = aod->GetTrack(it);
+      AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
+      if(!tr) AliFatal("Not a standard AOD");
       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
       if(TMath::Abs(tr->Eta())>0.9)continue;
       if(tr->Pt()<0.15)continue;
@@ -1145,7 +1304,8 @@ Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
       if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
       index=iCount-1;
-      dif=dphi;  }}
+      //  dif=dphi; 
+      }}
   
       return index;
 
@@ -1167,7 +1327,8 @@ Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
      else aod = fAODOut;   
   
       for(int it = 0;it < aod->GetNumberOfTracks();++it){
-      AliAODTrack *tr = aod->GetTrack(it);
+      AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
+      if(!tr) AliFatal("Not a standard AOD");
       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
       if(TMath::Abs(tr->Eta())>0.9)continue;
       if(tr->Pt()<0.15)continue;