fEvent(0x0),
fEventESD(0x0),
fEventAOD(0x0),
- fMCStack(0x0),
fRunNumber(-999),
fInternalRunNumber(0),
fPHOSGeo(0),
TArrayI centNMixed(nbins, nMixed);
SetCentralityBinning(centEdges, centNMixed);
+ for(int mod=1; mod <= kNMod; ++mod)
+ fModuleEnabled[mod-1] = kTRUE;
+
for(Int_t i=0;i<kNCenBins;i++){
for(Int_t j=0;j<2; j++)
for(Int_t k=0; k<2; k++) {
fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
}
+ this->MakeMCHistograms();
+
// Setup photon lists
Int_t kapacity = kNVtxZBins * GetNumberOfCentralityBins() * fNEMRPBins;
fCaloPhotonsPHOSLists = new TObjArray(kapacity);
PostData(1, fOutputContainer);
}
+void AliAnalysisTaskPi0Flow::MakeMCHistograms()
+{
+ //empty, MC extensions should override
+}
+
+void AliAnalysisTaskPi0Flow::DoMC()
+{
+ //empty, MC extensions should override
+}
+
+
//________________________________________________________________________
void AliAnalysisTaskPi0Flow::UserExec(Option_t *)
{
fEvent = GetEvent();
fEventESD = dynamic_cast<AliESDEvent*> (fEvent);
fEventAOD = dynamic_cast<AliAODEvent*> (fEvent);
- fMCStack = GetMCStack();
LogProgress(0);
ConsiderPi0sMix();
LogProgress(10);
- // Step 11: Update lists
+ // Step 11: MC
+ this->DoMC();
+
+ // Step 12: Update lists
UpdateLists();
LogProgress(11);
fCentNMixed = nMixed;
}
+//_____________________________________________________________________________
+void AliAnalysisTaskPi0Flow::SetEnablePHOSModule(int module, Bool_t enable)
+{
+ if( module < 1 || kNMod < module )
+ AliFatal(Form("PHOS Module must be between 1 and %i", kNMod));
+ else
+ fModuleEnabled[module-1] = enable;
+}
//________________________________________________________________________
Int_t mod = relId[0] ;
Int_t cellX = relId[2];
Int_t cellZ = relId[3] ;
+ if ( ! fModuleEnabled[mod-1] )
+ continue;
if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
continue ; // reject if not.
fCaloPhotonsPHOS->Add(new AliCaloPhoton(lorentzMomentum.X(),lorentzMomentum.Py(),lorentzMomentum.Z(),lorentzMomentum.E()) );
AliCaloPhoton * ph = (AliCaloPhoton*) fCaloPhotonsPHOS->At( fCaloPhotonsPHOS->GetLast() );
+ ph->SetCluster(clu);
ph->SetModule(mod) ;
lorentzMomentum*= ecore/lorentzMomentum.E() ;
//Evaluate CoreDispersion
Double_t m02=0.,m20=0. ;
EvalCoreLambdas(clu, cells, m02, m20) ;
- ph->SetDisp2Bit(TestCoreLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
+ ph->SetDisp2Bit(TestCoreLambda(clu->E(),m20,m02)) ; //Correct order m20,m02
+// ph->SetDisp2Bit(TestCoreLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
if(ph->IsDispOK()){
FillHistogram(Form("hCluDispM%d",mod),cellX,cellZ,1.);
}
AliError(Form("can not find histogram (of instance TH2) <%s> ",key)) ;
}
+void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z, Double_t w) const{
+ //Fills 1D histograms with key
+ TObject * obj = fOutputContainer->FindObject(key);
+
+ TH3 * th3 = dynamic_cast<TH3*> (obj);
+ if(th3) {
+ th3->Fill(x, y, z, w) ;
+ return;
+ }
+
+ AliError(Form("can not find histogram (of instance TH3) <%s> ",key)) ;
+}
+
+
//_____________________________________________________________________________
AliVEvent* AliAnalysisTaskPi0Flow::GetEvent()
{
//___________________________________________________________________________
-AliStack* AliAnalysisTaskPi0Flow::GetMCStack()
-{
- fMCStack = 0;
- AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
- if(eventHandler){
- AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
- if( mcEventHandler)
- fMCStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
- }
- return fMCStack;
-}
+// AliStack* AliAnalysisTaskPi0Flow::GetMCStack()
+// {
+// fMCStack = 0;
+// AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
+// if(eventHandler){
+// AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
+// if( mcEventHandler)
+// fMCStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
+// }
+// return fMCStack;
+// }
//___________________________________________________________________________
Int_t AliAnalysisTaskPi0Flow::GetCentralityBin(Float_t centralityV0M)
{
Double_t averageRP;
if(fHaveTPCRP)
- averageRP = (fRPV0A+fRPV0C+fRP) /3.;
+ averageRP = fRP ;// If possible, it is better to have EP bin from TPC
+ // to have similar events for miximng (including jets etc) (fRPV0A+fRPV0C+fRP) /3.;
else
averageRP = (fRPV0A+fRPV0C) /2.;
default : return 199;
}
}
+ if( kLHC13 == fPeriod ) {
+ switch(run){
+ case 195344 : return 1;
+ case 195346 : return 2;
+ case 195351 : return 3;
+ case 195389 : return 4;
+ case 195390 : return 5;
+ case 195391 : return 6;
+ case 195478 : return 7;
+ case 195479 : return 8;
+ case 195480 : return 9;
+ case 195481 : return 10;
+ case 195482 : return 11;
+ case 195483 : return 12;
+ case 195529 : return 13;
+ case 195531 : return 14;
+ case 195532 : return 15;
+ case 195566 : return 16;
+ case 195567 : return 17;
+ case 195568 : return 18;
+ case 195592 : return 19;
+ case 195593 : return 20;
+ case 195596 : return 21;
+ case 195633 : return 22;
+ case 195635 : return 23;
+ case 195644 : return 24;
+ case 195673 : return 25;
+ case 195675 : return 26;
+ case 195676 : return 27;
+ case 195677 : return 28;
+ case 195681 : return 29;
+ case 195682 : return 30;
+ case 195720 : return 31;
+ case 195721 : return 32;
+ case 195722 : return 33;
+ case 195724 : return 34;
+ case 195725 : return 34;
+ case 195726 : return 35;
+ case 195727 : return 36;
+ case 195760 : return 37;
+ case 195761 : return 38;
+ case 195765 : return 39;
+ case 195767 : return 40;
+ case 195783 : return 41;
+ case 195787 : return 42;
+ case 195826 : return 43;
+ case 195827 : return 44;
+ case 195829 : return 45;
+ case 195830 : return 46;
+ case 195831 : return 47;
+ case 195867 : return 48;
+ case 195869 : return 49;
+ case 195871 : return 50;
+ case 195872 : return 51;
+ case 195873 : return 52;
+ case 195935 : return 53;
+ case 195949 : return 54;
+ case 195950 : return 55;
+ case 195954 : return 56;
+ case 195955 : return 57;
+ case 195958 : return 58;
+ case 195989 : return 59;
+ case 195994 : return 60;
+ case 195998 : return 61;
+ case 196000 : return 62;
+ case 196006 : return 63;
+ case 196085 : return 64;
+ case 196089 : return 65;
+ case 196090 : return 66;
+ case 196091 : return 67;
+ case 196099 : return 68;
+ case 196105 : return 69;
+ case 196107 : return 70;
+ case 196185 : return 71;
+ case 196187 : return 72;
+ case 196194 : return 73;
+ case 196197 : return 74;
+ case 196199 : return 75;
+ case 196200 : return 76;
+ case 196201 : return 77;
+ case 196203 : return 78;
+ case 196208 : return 79;
+ case 196214 : return 80;
+ case 196308 : return 81;
+ case 196309 : return 82;
+ case 196310 : return 83;
+ case 196311 : return 84;
+ case 196433 : return 85;
+ case 196474 : return 86;
+ case 196475 : return 87;
+ case 196477 : return 88;
+ case 196528 : return 89;
+ case 196533 : return 90;
+ case 196535 : return 91;
+ case 196563 : return 92;
+ case 196564 : return 93;
+ case 196566 : return 94;
+ case 196568 : return 95;
+ case 196601 : return 96;
+ case 196605 : return 97;
+ case 196608 : return 98;
+ case 196646 : return 99;
+ case 196648 : return 100;
+ case 196701 : return 101;
+ case 196702 : return 102;
+ case 196703 : return 103;
+ case 196706 : return 104;
+ case 196714 : return 105;
+ case 196720 : return 106;
+ case 196721 : return 107;
+ case 196722 : return 108;
+ case 196772 : return 109;
+ case 196773 : return 110;
+ case 196774 : return 111;
+ case 196869 : return 112;
+ case 196870 : return 113;
+ case 196874 : return 114;
+ case 196876 : return 115;
+ case 196965 : return 116;
+ case 196967 : return 117;
+ case 196972 : return 118;
+ case 196973 : return 119;
+ case 196974 : return 120;
+ case 197003 : return 121;
+ case 197011 : return 122;
+ case 197012 : return 123;
+ case 197015 : return 124;
+ case 197027 : return 125;
+ case 197031 : return 126;
+ case 197089 : return 127;
+ case 197090 : return 128;
+ case 197091 : return 129;
+ case 197092 : return 130;
+ case 197094 : return 131;
+ case 197098 : return 132;
+ case 197099 : return 133;
+ case 197138 : return 134;
+ case 197139 : return 135;
+ case 197142 : return 136;
+ case 197143 : return 137;
+ case 197144 : return 138;
+ case 197145 : return 139;
+ case 197146 : return 140;
+ case 197147 : return 140;
+ case 197148 : return 141;
+ case 197149 : return 142;
+ case 197150 : return 143;
+ case 197152 : return 144;
+ case 197153 : return 145;
+ case 197184 : return 146;
+ case 197189 : return 147;
+ case 197247 : return 148;
+ case 197248 : return 149;
+ case 197254 : return 150;
+ case 197255 : return 151;
+ case 197256 : return 152;
+ case 197258 : return 153;
+ case 197260 : return 154;
+ case 197296 : return 155;
+ case 197297 : return 156;
+ case 197298 : return 157;
+ case 197299 : return 158;
+ case 197300 : return 159;
+ case 197302 : return 160;
+ case 197341 : return 161;
+ case 197342 : return 162;
+ case 197348 : return 163;
+ case 197349 : return 164;
+ case 197351 : return 165;
+ case 197386 : return 166;
+ case 197387 : return 167;
+ case 197388 : return 168;
+ default : return 199;
+ }
+ }
if((fPeriod == kUndefinedPeriod) && (fDebug >= 1) ) {
AliWarning("Period not defined");
}
const Int_t mulDigit=clu->GetNCells() ;
Double_t xc[mulDigit] ;
Double_t zc[mulDigit] ;
- Double_t ei[mulDigit] ;
+ Double_t wi[mulDigit] ;
Double_t x = 0 ;
Double_t z = 0 ;
const Double_t logWeight=4.5 ;
fPHOSGeo->RelPosInModule(relid, xi, zi);
xc[iDigit]=xi ;
zc[iDigit]=zi ;
- ei[iDigit] = elist[iDigit]*cells->GetCellAmplitude(absId) ;
- if (clu->E()>0 && ei[iDigit]>0) {
- Float_t w = TMath::Max( 0., logWeight + TMath::Log( ei[iDigit] / clu->E() ) ) ;
+ Double_t ei = elist[iDigit]*cells->GetCellAmplitude(absId) ;
+ wi[iDigit]=0. ;
+ if (clu->E()>0 && ei>0) {
+ wi[iDigit] = TMath::Max( 0., logWeight + TMath::Log( ei / clu->E() ) ) ;
+ Double_t w=wi[iDigit];
x += xc[iDigit] * w ;
z += zc[iDigit] * w ;
wtot += w ;
Double_t xCut = 0. ;
Double_t zCut = 0. ;
for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
- if (clu->E()>0 && ei[iDigit]>0.) {
- Double_t w = TMath::Max( 0., logWeight + TMath::Log( ei[iDigit] / clu->E() ) ) ;
+ Double_t w=wi[iDigit];
+ if (w>0.) {
Double_t xi= xc[iDigit] ;
Double_t zi= zc[iDigit] ;
if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){