Int_t ntracks=event->GetNumberOfTracks();
UInt_t specie = event->GetEventSpecie(); // skip laser events
if (specie==AliRecoParam::kCalib) return;
-
+ Int_t ntracksFriend = esdFriend->GetNumberOfTracks();
for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
//if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
// const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
AliESDfriendTrack* friendTrack0=NULL;
- if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack0 = esdFriend->GetTrack(itrack0);} //this guy can be NULL
+ if (esdFriend &&!esdFriend->TestSkipBit()){
+ if (itrack0<ntracksFriend){
+ friendTrack0 = esdFriend->GetTrack(itrack0);
+ } //this guy can be NULL
+ }
for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
AliESDtrack *track1 = event->GetTrack(itrack1);
AliESDfriendTrack* friendTrack1=NULL;
- if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack1 = esdFriend->GetTrack(itrack1);} //this guy can be NULL
+ if (esdFriend &&!esdFriend->TestSkipBit()){
+ if (itrack1<ntracksFriend){
+ friendTrack1 = esdFriend->GetTrack(itrack1);
+ } //this guy can be NULL
+ }
+
//
AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
AliESDfriendTrack* friendTrack1=NULL;
if (esdFriend) {
if (!esdFriend->TestSkipBit()){
- friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL
- friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL
+ Int_t ntracksFriend = esdFriend->GetNumberOfTracks();
+ if (v0->GetIndex(0)<ntracksFriend){
+ friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL
+ }
+ if (v0->GetIndex(1)<ntracksFriend){
+ friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL
+ }
}
}
if (track0->GetSign()<0) {
AliESDtrack *track = esdEvent->GetTrack(iTrack);
if(!track) continue;
AliESDfriendTrack* friendTrack=NULL;
- if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
+ if (esdFriend && !esdFriend->TestSkipBit()) {
+ Int_t ntracksFriend = esdFriend->GetNumberOfTracks();
+ if (iTrack<ntracksFriend){
+ friendTrack = esdFriend->GetTrack(iTrack);
+ } //this guy can be NULL
+ }
+
if(track->Charge()==0) continue;
if(!esdTrackCuts->AcceptTrack(track)) continue;
if(!accCuts->AcceptTrack(track)) continue;
local ftype=${3}
local outFile=${4}
shift 4
- if [ ! -f ${inFile} ] ; then
- echo ${inFile} not found!
- return 1
- fi
+# if [ ! -f ${inFile} ] ; then
+# echo ${inFile} not found!
+# return 1
+# fi
if [ -f ${outFile} ] ; then
>${outFile}
fi
#include "TGraphErrors.h"
#include "TCut.h"
#include "TCanvas.h"
+#include "TLinearFitter.h"
TObjArray arrayFit(3);
-Int_t kMarkers[10]={25,24,20,21};
-Int_t kColors[10]={1,2,4,3};
-const char * names[10] = {"LTM","LThisto","LTMhisto1"};
-const char * distNames[10] = {"Gauss","Gauss+flat","Gauss(m)+0.2*Gauss(m+5#sigma)","Gauss(m)+0.3*Gauss(m+3#sigma)"};
+Int_t kMarkers[10]={25,24,20,21,22};
+Int_t kColors[10]={1,2,4,3,6};
+const char * names[10] = {"LTM","LThisto","LTMhistoPar","Gaus fit", "Gaus in range"};
+const char * distNames[10] = {"Gauss","Gauss+flat","Gauss(m)+0.2*Gauss(m+5#sigma)","Gauss(m)+0.4*Gauss(m+5#sigma)"};
void TestLTM(Int_t nevents=10000){
TVectorD sigmaLTMHisto(6);
TVectorD meanLTMHisto1(6);
TVectorD sigmaLTMHisto1(6);
+ TVectorD meanLTMHistoPar(6);
+ TVectorD sigmaLTMHistoPar(6);
+ TVectorD meanLTMHistoGausLim(6);
+ TVectorD sigmaLTMHistoGausLim(6);
TVectorD paramLTM(10);
//
for (Int_t iltm=0; iltm<6; iltm++) vecLTM[iltm]=0.99999*(0.5+iltm/10.);
}
}
if (distType==3) {
- if (gRandom->Rndm()>0.3) { // gauss + second gaus 4 sigma away
+ if (gRandom->Rndm()>0.4) { // gauss + second gaus 5 sigma away
value=gRandom->Gaus(mean,sigma);
}else{
value=gRandom->Gaus(mean+5*sigma,sigma);
TStatToolkit::LTMHisto(&histo, paramLTM, vecLTM[iltm]);
meanLTMHisto[iltm]=paramLTM[1];
sigmaLTMHisto[iltm]=paramLTM[2];
+ //
+ // graph fit
+ //
+ Int_t binC=histo.GetXaxis()->FindBin(paramLTM[1]);
+ Int_t bin0=TMath::Max(histo.GetXaxis()->FindBin(paramLTM[1]-5*paramLTM[2]),1);
+ Int_t bin1=TMath::Min(histo.GetXaxis()->FindBin(paramLTM[1]+5*paramLTM[2]),histo.GetXaxis()->GetNbins()) ;
+ if ( (bin0-bin1) <3){
+ bin0=TMath::Max(binC-2,1);
+ bin1=TMath::Min(binC+2, histo.GetXaxis()->GetNbins());
+ }
+ //
+ /* test input:
+ TH1F histo("aaa","aaa",100,-10,10);for (Int_t i=0; i<10000; i++) histo.Fill(gRandom->Gaus(2,1));
+ */
+
+ TLinearFitter fitter(3,"hyp2");
+ for (Int_t ibin=bin0; ibin<bin1; ibin++){
+ Double_t x=histo.GetXaxis()->GetBinCenter(ibin);
+ Double_t y=histo.GetBinContent(ibin);
+ Double_t sy=histo.GetBinError(ibin);
+ Double_t xxx[3]={x,x*x};
+ if (y>3.*sy){
+ fitter.AddPoint(xxx,TMath::Log(y),sy/y);
+ }
+ }
+ if (fitter.GetNpoints()>3){
+ fitter.Eval();
+ // f(x)=[0]+[1]*x+[2]*x*x;
+ // f(x)=s*(x0-x0)^2+s
+ Double_t parSigma=TMath::Sqrt(2*TMath::Abs(fitter.GetParameter(2)));
+ Double_t parMean=-fitter.GetParameter(1)/(2.*fitter.GetParameter(2));
+ histo.Fit(&fg,"QN","QN", paramLTM[1]-5*paramLTM[2], paramLTM[1]+5*paramLTM[2]);
+ meanLTMHistoPar[iltm]=parMean;
+ sigmaLTMHistoPar[iltm]=parSigma;
+ meanLTMHistoGausLim[iltm]=fg.GetParameter(1);
+ sigmaLTMHistoGausLim[iltm]=fg.GetParameter(2);
+ }else{
+ meanLTMHistoPar[iltm]=paramLTM[1];
+ sigmaLTMHistoPar[iltm]=paramLTM[2];
+ meanLTMHistoGausLim[iltm]=meanG;
+ sigmaLTMHistoGausLim[iltm]=rmsG;
+ }
+
}
(*pcstream)<<"ltm"<<
//Input
"meanLTM.="<<&meanLTM<<
"meanLTMHisto.="<<&meanLTMHisto<<
"meanLTMHisto1.="<<&meanLTMHisto1<<
+ "meanLTMHistoPar.="<<&meanLTMHistoPar<<
+ "meanLTMHistoGausLim.="<<&meanLTMHistoGausLim<<
//
"sigmaLTM.="<<&sigmaLTM<<
"sigmaLTMHisto.="<<&sigmaLTMHisto<<
- "sigmaLTMHisto1.="<<&sigmaLTMHisto1<<
+ "sigmaLTMHistoPar.="<<&sigmaLTMHistoPar<<
+ "sigmaLTMHistoGausLim.="<<&sigmaLTMHistoGausLim<<
"\n";
}
delete pcstream;
TH1 * hisResol[10]={0};
TGraphErrors * grSigma[10]={0};
TCanvas *canvasLTM= new TCanvas("canvasLTM","canvasLTM",800,700);
- canvasLTM->Divide(2,2);
+ canvasLTM->Divide(2,2,0,0);
for (Int_t itype=0; itype<4; itype++){
canvasLTM->cd(itype+1);
TCut cutType=TString::Format("distType==0%d",itype).Data();
hisLTMTrunc[0]= (TH2*)(tree->GetHistogram()->Clone());
tree->Draw("sqrt(npoints-2)*(meanLTMHisto.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHisto(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
hisLTMTrunc[1]= (TH2*)tree->GetHistogram()->Clone();
- tree->Draw("sqrt(npoints-2)*(meanLTMHisto1.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
+ tree->Draw("sqrt(npoints-2)*(meanLTMHistoPar.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
hisLTMTrunc[2]= (TH2*)tree->GetHistogram()->Clone();
- TLegend * legend = new TLegend(0.5,0.7,0.9,0.9,distNames[itype]);
+ tree->Draw("sqrt(npoints-2)*(meanG-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
+ hisLTMTrunc[3]= (TH2*)tree->GetHistogram()->Clone();
+ tree->Draw("sqrt(npoints-2)*(meanLTMHistoGausLim.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
+ hisLTMTrunc[4]= (TH2*)tree->GetHistogram()->Clone();
+ TLegend * legend = new TLegend(0.5,0.7,0.88,0.88,distNames[itype]);
legend->SetBorderSize(0);
- for (Int_t ihis=0; ihis<3; ihis++){
+ for (Int_t ihis=0; ihis<5; ihis++){
// MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor);
- grSigma[ihis]=TStatToolkit::MakeStat1D( hisLTMTrunc[ihis],0,0.99,5,kMarkers[ihis],kColors[ihis]);
+ grSigma[ihis]=TStatToolkit::MakeStat1D( hisLTMTrunc[ihis],0,0.99,1,kMarkers[ihis],kColors[ihis]);
+ grSigma[ihis]->SetMaximum(5.0);
+ grSigma[ihis]->SetMinimum(0.5);
grSigma[ihis]->GetXaxis()->SetTitle("LTM fraction");
- grSigma[ihis]->GetYaxis()->SetTitle("#sigma_{meas}/sigma_{gauss}");
+ grSigma[ihis]->GetYaxis()->SetTitle("#sigma_{meas}/#sigma_{Theor.}");
if (ihis==0)grSigma[ihis]->Draw("alp");
if (ihis>0)grSigma[ihis]->Draw("lp");
legend->AddEntry(grSigma[ihis],names[ihis],"p");
Bool_t DiffObjects(const char *cdbFile1, const char *cdbFile2) const;
void ExtractBaseFolder(TString& url); // remove everything but the url from OCDB path
+
protected:
static TString fgkCondUri; // URI of the Conditions data base folder
TString cdbPath="";
TObjString *ostr;
AliCDBEntry *cdbEntry=0;
+ TGrid *myGrid = NULL;
UInt_t hash;
TMessage * file;
Int_t size;
if(!ostr) ostr = (TObjString*)cdbMap0->GetValue("default");
cdbPath = ostr->GetString();
if(cdbPath.Contains("local://"))cdbPath=cdbPath(8,cdbPath.Length()).Data();
+ if(!myGrid && cdbPath.Contains("alien://")){ //check if connection to alien is initialized
+ myGrid = TGrid::Connect("alien://"); //Oddly this will return also a pointer if connection fails
+ if(myGrid->GetPort()==0){ //if connection fails port 0 is saved, using this to check for successful connection
+ cerr << "Cannot connect to grid!" << endl;
+ continue;
+ }
+ }
try {
cdbEntry = (AliCDBEntry*) man->Get(*CDBId,kTRUE);
}catch(const exception &e){
// DumpOCDBFile("$ALICE_ROOT/OCDB/ITS/Align/Data/Run0_999999999_v0_s0.root", "ITS_Align_Data_Run0_999999999_v0_s0.dump")
//
if (finput==0) return ;
+ if (TString(finput).Contains("alien://") && gGrid==0){
+ TGrid *myGrid = TGrid::Connect("alien://"); //Oddly this will return also a pointer if connection fails
+ if(myGrid->GetPort()==0){ //if connection fails port 0 is saved, using this to check for successful connection
+ cerr << "Cannot connect to grid!" << endl;
+ return;
+ }
+ }
TFile *falignITS = TFile::Open(finput);
AliCDBEntry *entry = (AliCDBEntry*)falignITS->Get("AliCDBEntry");
if (!entry) return;
if (!fIp) return 0;
return fTPCdEdxInfo->GetdEdxInfo(fIp, regionID, calibID, qID, valueID);
}
+
+
+Double_t AliESDtrack::GetdEdxInfoTRD(Int_t method, Double_t p0, Double_t p1, Double_t p2){
+ //
+ // Methods
+ // mean values:
+ // 0.) linear
+ // 1.) logarithmic
+ // 2.) 1/sqrt
+ // 3.) power()
+ // time COG:
+ // 4.) linear
+ // 5.) logarithmic
+ // 6.) square
+ Int_t nSlicesPerLayer=GetNumberOfTRDslices();
+ Int_t nSlicesAll=GetNumberOfTRDslices()*kTRDnPlanes;
+
+ if (method<=3){
+ Double_t sumAmp=0;
+ Int_t sumW=0;
+ for (Int_t ibin=0; ibin<nSlicesAll; ibin++){
+ if (fTRDslices[ibin]<=0) continue;
+ sumW++;
+ if (method==0) sumAmp+=fTRDslices[ibin];
+ if (method==1) sumAmp+=TMath::Log(TMath::Abs(fTRDslices[ibin])+p0);
+ if (method==2) sumAmp+=1/TMath::Sqrt(TMath::Abs(fTRDslices[ibin])+p0);
+ if (method==3) sumAmp+=TMath::Power(TMath::Abs(fTRDslices[ibin])+p0,p1);
+ }
+ if (sumW==0) return 0;
+ Double_t dEdx=sumAmp/sumW;
+ if (method==1) dEdx= TMath::Exp(dEdx);
+ if (method==2) dEdx= 1/(dEdx*dEdx);
+ if (method==3) dEdx= TMath::Power(dEdx,1/p1);
+ return dEdx;
+ }
+ if (method>3){
+ Double_t sumWT=0;
+ Double_t sumW=0;
+ for (Int_t ibin=0; ibin<nSlicesAll; ibin++){
+ if (fTRDslices[ibin]<=0) continue;
+ Double_t time=(ibin%nSlicesPerLayer);
+ Double_t weight=fTRDslices[ibin];
+ if (method==5) weight=TMath::Log((weight+p0)/p0);
+ if (method==6) weight=TMath::Power(weight+p0,p1);
+ sumWT+=time*weight;
+ sumW+=weight;
+ }
+ if (sumW<=0) return 0;
+ Double_t meanTime=sumWT/sumW;
+ return meanTime;
+ }
+ return 0;
+}
}
void SetTPCdEdxInfo(AliTPCdEdxInfo * dEdxInfo);
Double_t GetdEdxInfo(Int_t regionID, Int_t calibID, Int_t qID,Int_t valueID);
+ Double_t GetdEdxInfoTRD(Int_t method, Double_t p0, Double_t p1, Double_t p2);
AliTPCdEdxInfo * GetTPCdEdxInfo() const {return fTPCdEdxInfo;}
Double_t GetTPCsignal() const {return fTPCsignal;}
AliCodeTimerStart(Form("creating summable digits for %s", det->GetName()));
det->Hits2SDigits();
AliCodeTimerStop(Form("creating summable digits for %s", det->GetName()));
- AliSysInfo::AddStamp(Form("Digit_%s_%d",det->GetName(),eventNr), 0,1, eventNr);
+ AliSysInfo::AddStamp(Form("SDigit_%s_%d",det->GetName(),eventNr), 0,1, eventNr);
}
}
const char* excludeDetectors)
{
// run the digitization and produce digits from sdigits
-
AliCodeTimerAuto("",0)
// initialize CDB storage, run number, set CDB lock
}
detArr.AddLast(digitizer);
AliInfo(Form("Created digitizer from SDigits -> Digits for %s", det->GetName()));
+
}
//
if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
digInp.InitEvent(); //this must be after call of Connect Input tress.
if (outRl) outRl->SetEventNumber(eventsCreated-1);
static_cast<AliStream*>(digInp.GetInputStream(0))->ImportgAlice(); // use gAlice of the first input stream
- for (int id=0;id<ndigs;id++) ((AliDigitizer*)detArr[id])->Digitize("");
+ for (int id=0;id<ndigs;id++) {
+ ((AliDigitizer*)detArr[id])->Digitize("");
+ AliSysInfo::AddStamp(Form("Digit_%s_%d",detArr[id]->GetName(),eventsCreated), 0,2, eventsCreated);
+ }
digInp.FinishEvent();
};
digInp.FinishGlobal();
- //
+ //
return kTRUE;
}
//5.) and now see results in file testredirector.root
}
-void TTreeSRedirector::UnitTest(){
+void TTreeSRedirector::UnitTest(Int_t testEntries){
//
//
//
- UnitTestSparse(0.5);
- UnitTestSparse(0.1);
- UnitTestSparse(0.01);
+ UnitTestSparse(0.5,testEntries);
+ UnitTestSparse(0.1,testEntries);
+ UnitTestSparse(0.01,testEntries);
}
-void TTreeSRedirector::UnitTestSparse(Double_t scale){
+void TTreeSRedirector::UnitTestSparse(Double_t scale, Int_t testEntries){
//
// Unit test for the TTreeSRedirector
// 1.) Test TTreeRedirector
if (scale<=0) scale=1;
if (scale>1) scale=1;
TTreeSRedirector *pcstream = new TTreeSRedirector("testpcstreamSparse.root","recreate");
- for (Int_t ientry=0; ientry<20000; ientry++){
- TVectorD vecRandom(50);
- TVectorD vecZerro(50); // zerro vector
- for (Int_t j=0; j<50; j++) vecRandom[j]=ientry+gRandom->Rndm();
+ for (Int_t ientry=0; ientry<testEntries; ientry++){
+ TVectorD vecRandom(200);
+ TVectorD vecZerro(200); // zerro vector
+ for (Int_t j=0; j<200; j++) vecRandom[j]=j+ientry+0.1*gRandom->Rndm();
Bool_t isSelected= (gRandom->Rndm()<scale);
TVectorD *pvecFull = &vecRandom;
TVectorD *pvecSparse = isSelected ? &vecRandom:0;
printf("#UnitTest:\tTTreeSRedirector::TestSparse(%f)\tRatioSkip0\t%f\n",scale,ratio0);
printf("#UnitTest:\tTTreeSRedirector::TestSparse(%f)\tRatioZerro\t%f\n",scale,ratio1);
// b.) Integrity
- Int_t outlyersSparseSkip=treeSparseSkip->Draw("1","(vec.fElements-ientry-0.5)>0.5","goff");
- Int_t outlyersSparseSkip0=treeSparseSkip0->Draw("1","(vec.fElements-ientry-0.5)>0.5","goff");
- printf("#UnitTest:\tTTreeSRedirector::TestSparse(%f)\tOutlyersSkip\t%d\n",scale,outlyersSparseSkip);
- printf("#UnitTest:\tTTreeSRedirector::TestSparse(%f)\tOutlyersSkip0\t%d\n",scale,outlyersSparseSkip);
-
- Bool_t isOK=(ratio<1.2)&&(outlyersSparseSkip0);
+ Int_t outlyersSparseSkip=treeSparseSkip->Draw("1","(vec.fElements-ientry-Iteration$-0.5)>0.5","goff");
+ Int_t outlyersSparseSkip0=treeSparseSkip0->Draw("1","(vec.fElements-ientry-Iteration$-0.5)>0.5","goff");
+ printf("#UnitTest:\tTTreeSRedirector::TestSparse(%f)\tOutlyersSkip\t%d\n",scale,outlyersSparseSkip!=0);
+ printf("#UnitTest:\tTTreeSRedirector::TestSparse(%f)\tOutlyersSkip0\t%d\n",scale,outlyersSparseSkip0!=0);
+ // c.) Number of entries
+ //
+ Int_t entries=treeFull->GetEntries();
+ Int_t entries0=treeSparseSkip0->GetEntries();
+ Bool_t isOKStat =(entries==entries0);
+ printf("#UnitTest:\tTTreeSRedirector::TestSparse(%f)\tEntries\t%d\n",scale,isOKStat);
+ //
+ // d.)Reading test
+ TVectorD *pvecRead = 0;
+ treeSparseSkip0->SetBranchAddress("vec.",&pvecRead);
+ Bool_t readOK=kTRUE;
+ for (Int_t ientry=0; ientry<testEntries; ientry++){
+ if (!pvecRead) continue;
+ if (pvecRead->GetNrows()==0) continue;
+ if (TMath::Abs((*pvecRead)[0]-ientry)>0.5) readOK=kFALSE;
+ }
+ printf("#UnitTest:\tTTreeSRedirector::TestSparse(%f)\tReadOK\t%d\n",scale,readOK);
+ //
+ // e.)Global test
+ Bool_t isOK=(outlyersSparseSkip0==0)&&isOKStat&&readOK;
printf("#UnitTest:\tTTreeSRedirector::TestSparse(%f)\tisOk\t%d\n",scale,isOK);
+
}
TTreeSRedirector::TTreeSRedirector(const char *fname,const char * option) :
if (element->fClass==0) {
element->fClass=pClass;
}else{
- if (element->fClass!=pClass){
+ if (element->fClass!=pClass && pClass!=0){
fStatus++;
return 1; //mismatched data element
}
// Build the Tree
//
//if (fTree && fTree->GetEntries()>0) return;
- if (!fTree) fTree = new TTree(GetName(),GetName());
+ Int_t entriesFilled=0;
+ if (!fTree) {
+ fTree = new TTree(GetName(),GetName());
+ }else{
+ entriesFilled=fTree->GetEntries();
+ }
Int_t entries = fElements->GetEntriesFast();
if (!fBranches) fBranches = new TObjArray(entries);
if (element->fClass){
if (element->fClass->GetBaseClass("TClonesArray")){
TBranch * br = fTree->Branch(bname1,element->fClass->GetName(),&(element->fPointer));
+ if (entriesFilled!=0) {
+ br->SetAddress(0);
+ for (Int_t ientry=0; ientry<entriesFilled;ientry++) br->Fill();
+ br->SetAddress(&(element->fPointer));
+ }
fBranches->AddAt(br,i);
}else
{
TBranch * br = fTree->Branch(bname1,element->fClass->GetName(),&(element->fPointer));
+ if (entriesFilled!=0) {
+ br->SetAddress(0);
+ for (Int_t ientry=0; ientry<entriesFilled;ientry++) br->Fill();
+ br->SetAddress(&(element->fPointer));
+ }
fBranches->AddAt(br,i);
}
}
char bname2[1000];
snprintf(bname2,1000,"B%d/%c",i,element->GetType());
TBranch * br = fTree->Branch(bname1,element->fPointer,bname2);
+ if (entriesFilled!=0) {
+ br->SetAddress(0);
+ for (Int_t ientry=0; ientry<entriesFilled;ientry++) br->Fill();
+ br->SetAddress(element->fPointer);
+ }
+
fBranches->AddAt(br,i);
}
}
void Close();
static void Test();
static void Test2();
- static void UnitTestSparse(Double_t scale);
- static void UnitTest();
+ static void UnitTestSparse(Double_t scale, Int_t testEntries);
+ static void UnitTest(Int_t testEntries=5000);
void StoreObject(TObject* object);
TFile * GetFile() {return fDirectory->GetFile();}
TDirectory * GetDirectory() {return fDirectory;}
segRowUp = 48;
} else if ( wireSegmentID == 3 || wireSegmentID == 7 ) {
segRowDown = 48;
- segRowUp = 64;
+ segRowUp = 63;
} else if ( wireSegmentID == 8 ) {
segRowDown = 64;
segRowUp = 75;
// Analyze gain as a function of multiplicity and produce calibration graphs
// padRegion -- 0: short, 1: medium, 2: long, 3: absolute calibration of full track
//
+ Int_t kMarkers[10]={25,24,20,21,22};
+ Int_t kColors[10]={1,2,4,3,6};
if (!fGainMult) return kFALSE;
if (!(fGainMult->GetHistTopology())) return kFALSE;
+ const Double_t kMinStat=100;
//
// "dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength"
TObjArray arrMax;
histQtot->SetName("fGainMult_GetHistPadEqual_00");
}
//
- histQmax->FitSlicesY(0,0,-1,0,"QNR",&arrMax);
- TH1D * corrMax = (TH1D*)arrMax.At(1)->Clone();
- histQtot->FitSlicesY(0,0,-1,0,"QNR",&arrTot);
- TH1D * corrTot = (TH1D*)arrTot.At(1)->Clone();
- AliInfo(Form("hisQtot.GetEntries()=%d",histQtot->GetEntries()));
- AliInfo(Form("hisQmax.GetEntries()=%d",histQmax->GetEntries()));
- if (histQmax->GetMean(2)<=0 || histQtot->GetMean(2)) {
+
+ if (histQmax->GetEntries()<=kMinStat || histQtot->GetEntries()<=kMinStat) {
+ AliError(Form("hisQtot.GetEntries()=%f",histQtot->GetEntries()));
+ AliError(Form("hisQmax.GetEntries()=%f",histQmax->GetEntries()));
return kFALSE;
}
- corrMax->Scale(1./histQmax->GetMean(2));
- corrTot->Scale(1./histQtot->GetMean(2));
+
+ TGraphErrors * graphMax = TStatToolkit::MakeStat1D( histQmax,0,0.8,4,kMarkers[padRegion],kColors[padRegion]);
+ TGraphErrors * graphTot = TStatToolkit::MakeStat1D( histQtot,0,0.8,4,kMarkers[padRegion],kColors[padRegion]);
//
const char* names[4]={"SHORT","MEDIUM","LONG","ABSOLUTE"};
//
- TGraphErrors * graphMax = new TGraphErrors(corrMax);
- TGraphErrors * graphTot = new TGraphErrors(corrTot);
Double_t meanMax = TMath::Mean(graphMax->GetN(), graphMax->GetY());
Double_t meanTot = TMath::Mean(graphTot->GetN(), graphTot->GetY());
+ if (meanMax<=0 || meanTot<=0){
+ AliError(Form("meanMax=%f",meanMax));
+ AliError(Form("meanTot=%f",meanTot));
+ return kFALSE;
+ }
//
for (Int_t ipoint=0; ipoint<graphMax->GetN(); ipoint++) {graphMax->GetY()[ipoint]/=meanMax;}
for (Int_t ipoint=0; ipoint<graphTot->GetN(); ipoint++) {graphTot->GetY()[ipoint]/=meanTot;}
-
//
graphMax->SetNameTitle(Form("TGRAPHERRORS_QMAX_DIPANGLE_%s_BEAM_ALL",names[padRegion]),
Form("TGRAPHERRORS_QMAX_DIPANGLE_%s_BEAM_ALL",names[padRegion]));
TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
}
}
+ printf("\ndEdx correction matrix used in GetQnormCorr\n");
+ fQNormCorr->Print();
+
}
}
}
-void AliTPCClusterParam::SetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType, Float_t val){
+void AliTPCClusterParam::SetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType, Float_t val, Int_t mode){
//
// ipad - pad type
// itype - 0- qtot 1-qmax
// rows
// itype*3+ipad - itype=0 qtot itype=1 qmax ipad=0
//
+ if (mode==1) {
+ // mode introduced in 20.07.2014 - git describe ~ vAN-20140703-48-g3449a97 - to keep back compatibility with o
+ (*fQNormCorr)(itype*3+ipad, corrType) *= val; // set value
+ return;
+ }
if (itype<2) (*fQNormCorr)(itype*3+ipad, corrType) *= val; // multiplicative correction
if (itype>=2) (*fQNormCorr)(itype*3+ipad, corrType)+= val; // additive correction
}
void FitResol(TTree * tree);
void FitRMS(TTree * tree);
void SetQnorm(Int_t ipad, Int_t itype, const TVectorD *const norm);
- void SetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType, Float_t val);
+ void SetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType, Float_t val, Int_t mode=1);
Double_t GetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType) const;
TMatrixD *GetQnormCorrMatrix(){return fQNormCorr;};
void ResetQnormCorr();
fUseMultiplicityCorrectionDedx(kTRUE), // use Dedx multiplicity correction
fUseAlignmentTime(kTRUE), // use time dependent alignment correction
fUseIonTailCorrection(0), // no ion tail correction for now
- fCrosstalkCorrection(0), // crosstalk correction factor (fro each signal substracted by (mean signal in wite patch)xfCrosstalkCorrection) - Effect important only after removing oc capacitors in 2012
- //
+ fCrosstalkCorrection(0), // crosstalk correction factor (from each signal substracted by (mean signal in wite patch)xfCrosstalkCorrection) - Effect important only after removing oc capacitors in 2012
+ fCrosstalkCorrectionMissingCharge(1), // crosstalk correction factor - missing charge (from each signal substracted by (mean signal in wite patch)xfCrosstalkCorrection) - Effect important only after removing oc capacitors in 2012
+ //
fUseTotCharge(kTRUE), // switch use total or max charge
fMinFraction(0.01), // truncated mean - lower threshold
fMaxFaction(0.7), // truncated mean - upper threshold
void SetUseTOFCorrection(Bool_t flag) {fUseTOFCorrection = flag;}
void SetUseIonTailCorrection(Int_t flag) {fUseIonTailCorrection = flag;}
void SetCrosstalkCorrection(Float_t crosstalkCorrection) {fCrosstalkCorrection= crosstalkCorrection; }
- //
+ void SetCrosstalkCorrectionMissingCharge(Float_t crosstalkCorrection) {fCrosstalkCorrectionMissingCharge= crosstalkCorrection; }
+ //
Int_t GetUseFieldCorrection() const {return fUseFieldCorrection;}
Int_t GetUseComposedCorrection() const {return fUseComposedCorrection;}
Int_t GetUseRPHICorrection() const {return fUseRPHICorrection;}
Bool_t GetUseTOFCorrection() {return fUseTOFCorrection;}
Int_t GetUseIonTailCorrection() const {return fUseIonTailCorrection;}
Double_t GetCrosstalkCorrection() const {return fCrosstalkCorrection;}
+ Double_t GetCrosstalkCorrectionMissingCharge() const {return fCrosstalkCorrectionMissingCharge;}
Bool_t GetUseMultiplicityCorrectionDedx() const {return fUseMultiplicityCorrectionDedx;}
Int_t GetGainCorrectionHVandPTMode() const { return fGainCorrectionHVandPTMode;}
Bool_t fUseAlignmentTime; // use time dependent alignment correction
Int_t fUseIonTailCorrection; // use ion tail correction
Double_t fCrosstalkCorrection; // crosstalk correction factor (fro each signal substracted by (mean signal in wite patch)xfCrosstalkCorrection) - Effect important only after removing oc capacitors in 2012
- //
+ Double_t fCrosstalkCorrectionMissingCharge; // crosstalk correction factor - missing charge factor (from each signal substracted by (mean signal in wite patch)xfCrosstalkCorrection) - Effect important only after removing capacitors in 2012
+ //
// dEdx switches
//
Bool_t fUseTotCharge; // switch use total or max charge
// to be switched off for pass 0 reconstruction
// Use static function, other option will be to use
// additional specific storage ?
- ClassDef(AliTPCRecoParam, 19)
+ ClassDef(AliTPCRecoParam, 21)
};
#include <TFile.h>
#include <TObjArray.h>
#include <TTree.h>
+#include <TMatrixD.h>
#include <TGraphErrors.h>
#include <TTimeStamp.h>
#include "AliLog.h"
fkParam(0),
fDebugStreamer(0),
fUseHLTClusters(4),
+ fCrossTalkSignalArray(0),
fSeedsPool(0),
fFreeSeedsID(500),
fNFreeSeeds(0),
fYMax[irow]=0;
fPadLength[irow]=0;
}
-
}
//_____________________________________________________________________
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"Update"<<
"cl.="<<c<<
+ "event="<<event<<
"track.="<<¶m<<
"\n";
}
fSeeds(0),
fIteration(0),
fkParam(0),
- fDebugStreamer(0),
- fUseHLTClusters(4),
- fSeedsPool(0),
+ fDebugStreamer(0),
+ fUseHLTClusters(4),
+ fCrossTalkSignalArray(0),
+ fSeedsPool(0),
fFreeSeedsID(500),
fNFreeSeeds(0),
fLastSeedID(-1)
}
//
fSeedsPool = new TClonesArray("AliTPCseed",1000);
+
+ // crosstalk array and matrix initialization
+ Int_t nROCs = 72;
+ Int_t nTimeBinsAll = par->GetMaxTBin();
+ Int_t nWireSegments = 11;
+ fCrossTalkSignalArray = new TObjArray(nROCs); // for 36 sectors
+ fCrossTalkSignalArray->SetOwner(kTRUE);
+ for (Int_t isector=0; isector<nROCs; isector++){
+ TMatrixD * crossTalkSignal = new TMatrixD(nWireSegments,nTimeBinsAll);
+ for (Int_t imatrix = 0; imatrix<11; imatrix++)
+ for (Int_t jmatrix = 0; jmatrix<nTimeBinsAll; jmatrix++){
+ (*crossTalkSignal)[imatrix][jmatrix]=0.;
+ }
+ fCrossTalkSignalArray->AddAt(crossTalkSignal,isector);
+ }
+
}
//________________________________________________________________________
AliTPCtracker::AliTPCtracker(const AliTPCtracker &t):
fSeeds(0),
fIteration(0),
fkParam(0),
- fDebugStreamer(0),
- fUseHLTClusters(4),
- fSeedsPool(0),
+ fDebugStreamer(0),
+ fUseHLTClusters(4),
+ fCrossTalkSignalArray(0),
+ fSeedsPool(0),
fFreeSeedsID(500),
fNFreeSeeds(0),
fLastSeedID(-1)
fSeeds->Clear();
delete fSeeds;
}
+ if (fCrossTalkSignalArray) delete fCrossTalkSignalArray;
if (fDebugStreamer) delete fDebugStreamer;
if (fSeedsPool) delete fSeedsPool;
}
//
// TTree * tree = fClustersArray.GetTree();
AliInfo("LoadClusters()\n");
+ // reset crosstalk matrix
+ //
+ for (Int_t isector=0; isector<72; isector++){ //set all ellemts of crosstalk matrix to 0
+ TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
+ if (crossTalkMatrix)(*crossTalkMatrix)*=0;
+ }
TTree * tree = fInput;
TBranch * br = tree->GetBranch("Segment");
// Conversion of pad, row coordinates in local tracking coords.
// Could be skipped here; is already done in clusterfinder
+
Int_t j=Int_t(tree->GetEntries());
for (Int_t i=0; i<j; i++) {
br->GetEntry(i);
//
Int_t sec,row;
fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
+
+ // wire segmentID and nPadsPerSegment to be used for Xtalk calculation
+ Int_t wireSegmentID = fkParam->GetWireSegment(sec,row);
+ Float_t nPadsPerSegment = (Float_t)(fkParam->GetNPadsPerSegment(wireSegmentID));
+ TMatrixD &crossTalkSignal = *((TMatrixD*)fCrossTalkSignalArray->At(sec));
+ Int_t nCols=crossTalkSignal.GetNcols();
+ Double_t missingChargeFactor= AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrectionMissingCharge();
for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
- Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
+ AliTPCclusterMI *clXtalk= static_cast<AliTPCclusterMI*>(clrow->GetArray()->At(icl));
+ Transform((AliTPCclusterMI*)(clXtalk));
+
+ Int_t timeBinXtalk = clXtalk->GetTimeBin();
+ Double_t qTotXtalk = 0.;
+ Double_t rmsTime2 = clXtalk->GetSigmaZ2()/(fkParam->GetZWidth()*fkParam->GetZWidth());
+ Double_t rmsPad2 = clXtalk->GetSigmaY2()/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec));
+
+ Double_t norm= 2.*TMath::Exp(-1.0/(2.*rmsTime2))+2.*TMath::Exp(-4.0/(2.*rmsTime2))+1.;
+ for (Int_t itb=timeBinXtalk-2, idelta=-2; itb<=timeBinXtalk+2; itb++,idelta++) {
+ if (itb<0 || itb>=nCols) continue;
+ Double_t missingCharge=0;
+ Double_t trf= TMath::Exp(-idelta*idelta/(2.*rmsTime2));
+ if (missingChargeFactor>0) for (Int_t dpad=-2; dpad<=2; dpad++){
+ Double_t qPad = clXtalk->GetMax()*TMath::Exp(-dpad*dpad/(2*rmsPad2));
+ if (qPad*trf<fkParam->GetZeroSup()){
+ missingCharge+=qPad*missingChargeFactor;
+ }
+ }
+ qTotXtalk = (clXtalk->GetQ()+missingCharge)*trf/norm;
+ crossTalkSignal[wireSegmentID][itb]+= qTotXtalk/nPadsPerSegment;
+ }
}
//
AliTPCtrackerRow * tpcrow=0;
clrow->Clear("C");
LoadOuterSectors();
LoadInnerSectors();
+
+ cout << " =================================================================================================================================== " << endl;
+ cout << " AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() << endl;
+ cout << " AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() << endl;
+ cout << " =================================================================================================================================== " << endl;
+
if (AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection()) ApplyTailCancellation();
+ if (AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()!=0.) ApplyXtalkCorrection();
return 0;
}
}
}
+void AliTPCtracker::ApplyXtalkCorrection(){
+ //
+ // ApplyXtalk correction
+ // Loop over all clusters
+ // add to each cluster signal corresponding to common Xtalk mode for given time bin at given wire segment
+ // cluster loop
+ for (Int_t isector=0; isector<36; isector++){ //loop tracking sectors
+ for (Int_t iside=0; iside<2; iside++){ // loop over sides A/C
+ AliTPCtrackerSector §or= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
+ Int_t nrows = sector.GetNRows();
+ for (Int_t row = 0;row<nrows;row++){ // loop over rows
+ AliTPCtrackerRow& tpcrow = sector[row]; // row object
+ Int_t ncl = tpcrow.GetN1(); // number of clusters in the row
+ if (iside>0) ncl=tpcrow.GetN2();
+ Int_t xSector=0; // sector number in the TPC convention 0-72
+ if (isector<18){ //if IROC
+ xSector=isector+(iside>0)*18;
+ }else{
+ xSector=isector+18; // isector -18 +36
+ if (iside>0) xSector+=18;
+ }
+ TMatrixD &crossTalkMatrix= *((TMatrixD*)fCrossTalkSignalArray->At(xSector));
+ Int_t wireSegmentID = fkParam->GetWireSegment(xSector,row);
+ for (Int_t i=0;i<ncl;i++) {
+ AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
+ Int_t iTimeBin=TMath::Nint(cluster->GetTimeBin());
+ Double_t xTalk= crossTalkMatrix[wireSegmentID][iTimeBin];
+ cluster->SetMax(cluster->GetMax()+xTalk);
+ const Double_t kDummy=4;
+ Double_t sumxTalk=xTalk*kDummy; // should be calculated via time response function
+ cluster->SetQ(cluster->GetQ()+sumxTalk);
+
+
+ if (AliTPCReconstructor::StreamLevel()==1) {
+ TTreeSRedirector &cstream = *fDebugStreamer;
+ if (gRandom->Rndm() > 0.){
+ cstream<<"Xtalk"<<
+ "isector=" << isector << // sector [0,36]
+ "iside=" << iside << // side A or C
+ "row=" << row << // padrow
+ "i=" << i << // index of the cluster
+ "xSector=" << xSector << // sector [0,72]
+ "wireSegmentID=" << wireSegmentID << // anode wire segment id [0,10]
+ "iTimeBin=" << iTimeBin << // timebin of the corrected cluster
+ "xTalk=" << xTalk << // Xtalk contribution added to Qmax
+ "sumxTalk=" << sumxTalk << // Xtalk contribution added to Qtot (roughly 3*Xtalk)
+ "cluster.=" << cluster << // corrected cluster object
+ "\n";
+ }
+ }// dump the results to the debug streamer if in debug mode
+
+
+
+
+
+ }
+ }
+ }
+ }
+}
+
void AliTPCtracker::ApplyTailCancellation(){
//
// Correct the cluster charge for the ion tail effect
if (!cl0) continue;
Int_t nclPad=0;
- for (Int_t icl1=0; icl1<ncl;icl1++){ // second loop over clusters
-
+ for (Int_t icl1=0; icl1<ncl;icl1++){ // second loop over clusters
AliTPCclusterMI *cl1= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl1]));
if (!cl1) continue;
if (TMath::Abs(cl0->GetPad()-cl1->GetPad())>4) continue; // no contribution if far away in pad direction
//
// Function in order to calculate the amount of the correction to be added for a given cluster, return values are ionTailTaoltal and ionTailMax
+ // Parameters:
+ // cl0 - cluster to be modified
+ // cl1 - source cluster ion tail of this cluster will be added to the cl0 (accroding time and pad response function)
//
-
const Double_t kMinPRF = 0.5; // minimal PRF width
ionTailTotal = 0.; // correction value to be added to Qtot of cl0
ionTailMax = 0.; // correction value to be added to Qmax of cl0
const Int_t deltaTimebin = TMath::Nint(TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin()))+12; //distance between pads of cl1 and cl0 increased by 12 bins
Double_t rmsPad1 = (cl1->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl1->GetSigmaY2())/padWidth);
Double_t rmsPad0 = (cl0->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl0->GetSigmaY2())/padWidth);
-
+
+
Double_t sumAmp1=0.;
for (Int_t idelta =-2; idelta<=2;idelta++){
+
#ifndef ALITPCTRACKER_H
#define ALITPCTRACKER_H
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
//-------------------------------------------------------
#include <TArrayI.h>
+#include <TMatrixD.h>
#include "AliTracker.h"
#include "AliTPCreco.h"
#include "AliTPCclusterMI.h"
virtual void FillClusterArray(TObjArray* array) const;
void Transform(AliTPCclusterMI * cluster);
void ApplyTailCancellation();
+ void ApplyXtalkCorrection();
void GetTailValue(Float_t ampfactor,Double_t &ionTailMax,Double_t &ionTailTotal,TGraphErrors **graphRes,Float_t *indexAmpGraphs,AliTPCclusterMI *cl0,AliTPCclusterMI *cl1);
//
void FillESD(const TObjArray* arr);
TTreeSRedirector *fDebugStreamer; //!debug streamer
Int_t fUseHLTClusters; // use HLT clusters instead of offline clusters
//
+
+ TObjArray * fCrossTalkSignalArray; // for 36 sectors
TClonesArray* fSeedsPool; //! pool of seeds
TArrayI fFreeSeedsID; //! array of ID's of freed seeds
Int_t fNFreeSeeds; //! number of seeds freed in the pool
// ctor which should be used
//
AliDebug(2,"(AliDigitizationInput* digInput) was processed");
- if (AliTPCReconstructor::StreamLevel()>0) fDebugStreamer = new TTreeSRedirector("TPCDigitDebug.root");
+ if (AliTPCReconstructor::StreamLevel()>0) fDebugStreamer = new TTreeSRedirector("TPCDigitDebug.root","recreate");
}
//------------------------------------------------------------------------
void AliTPCDigitizer::Digitize(Option_t* option)
{
- DigitizeFast(option);
- //DigitizeWithTailAndCrossTalk(option);
+ //DigitizeFast(option);
+ DigitizeWithTailAndCrossTalk(option);
}
//------------------------------------------------------------------------
// Wire segmentationn is obtatined from the
// AliTPCParam::GetWireSegment(Int_t sector, Int_t row); // to be implemented
// AliTPCParam::GetNPadsPerSegment(Int_t segmentID); // to be implemented
+ // 3.) Substract form the signal contribution from the previous tracks - Ion tail in case speified in the AliTPCRecoParam
+ // AliTPCRecoParam::GetUseIonTailCorrection()
//
// Ion tail simulation:
// 1.) Needs signal from pad+-1, taking signal from history
// merge input tree's with summable digits
// output stored in TreeTPCD
- //
-
+ //
+ AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
+ AliTPCRecoParam *recoParam = calib->GetRecoParam(0);
+ AliDebug(1, Form(" recoParam->GetCrosstalkCorrection() = %f", recoParam->GetCrosstalkCorrection()));
+ AliDebug(1,Form(" recoParam->GetUseIonTailCorrection() = %f ", recoParam->GetUseIonTailCorrection()));
Int_t nROCs = 72;
char s[100];
char ss[100];
TObject *rocFactorOROC = ionTailArr->FindObject("factorOROC");
Float_t factorIROC = (atof(rocFactorIROC->GetTitle()));
Float_t factorOROC = (atof(rocFactorOROC->GetTitle()));
+ Int_t nIonTailBins =0;
TObjArray timeResFunc(nROCs);
for (Int_t isec = 0;isec<nROCs;isec++){ //loop overs sectors
// Array of TGraphErrors for a given sector
TGraphErrors ** graphRes = new TGraphErrors *[20];
- Float_t * indexAmpGraphs = new Float_t[20];
+ Float_t * trfIndexArr = new Float_t[20];
for (Int_t icache=0; icache<20; icache++)
{
graphRes[icache] = NULL;
- indexAmpGraphs[icache] = 0;
+ trfIndexArr[icache] = 0;
}
- if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(isec,graphRes,indexAmpGraphs)) continue;
+ if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(isec,graphRes,trfIndexArr)) continue;
+
// fill all TGraphErrors of trfs (time response function) of a given sector to a TObjArray
TObjArray *timeResArr = new TObjArray(20); timeResArr -> SetOwner(kTRUE);
for (Int_t ires = 0;ires<20;ires++) timeResArr->AddAt(graphRes[ires],ires);
timeResFunc.AddAt(timeResArr,isec); // Fill all trfs into a single TObjArray
- delete timeResArr;
+ nIonTailBins = graphRes[3]->GetN();
}
//
TObjArray crossTalkSignalArray(nROCs); // for 36 sectors
TVectorD * qTotSectorOld = new TVectorD(nROCs);
- TVectorD * qTotSectorNew = new TVectorD(nROCs);
Float_t qTotTPC = 0.;
+ Float_t qTotPerSector = 0.;
Int_t nTimeBinsAll = 1100;
Int_t nWireSegments=11;
// 1.a) crorstalk matrix initialization
continue;
}
// Calculate number of pads in a anode wire segment for normalization
- Int_t anodeSegmentID = param->GetWireSegment(sector,padRow);
- Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(anodeSegmentID));
+ Int_t wireSegmentID = param->GetWireSegment(sector,padRow);
+ Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(wireSegmentID));
// structure with mean signal per pad to be filled for each timebin in first loop (11 anodeWireSegment and 1100 timebin)
TMatrixD &crossTalkSignal = *((TMatrixD*)crossTalkSignalArray.At(sector));
AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
- digrow->SetID(globalRowID);
+ // digrow->SetID(globalRowID);
Int_t nTimeBins = 0;
Int_t nPads = 0;
Bool_t digitize = kFALSE;
if (GetRegionOfInterest() && !digitize) break;
}
if (!digitize) continue;
- digrow->Allocate(nTimeBins,nPads);
+ //digrow->Allocate(nTimeBins,nPads);
Float_t q = 0;
Int_t labptr = 0;
Int_t nElems = nTimeBins*nPads; // element is a unit of a given row's pad-timebin space
q += *(pdig[i]);
pdig[i]++;
}
+ if (q<=0) continue;
Int_t padNumber = elem/nTimeBins;
Int_t timeBin = elem%nTimeBins;
- q/=16.; //conversion factor
Float_t gain = gainROC->GetValue(padRow,padNumber); // get gain for given - pad-row pad
q*= gain;
-
- crossTalkSignal[anodeSegmentID][timeBin]+= q/nPadsPerSegment; // Qtot per segment for a given timebin
+ crossTalkSignal[wireSegmentID][timeBin]+= q/nPadsPerSegment; // Qtot per segment for a given timebin
qTotSectorOld -> GetMatrixArray()[sector] += q; // Qtot for each sector
qTotTPC += q; // Qtot for whole TPC
} // end of q loop
-
- cout << " sector = " << sector << " row = " << padRow << " SegmentID = " << anodeSegmentID ;
- cout << " nPadsPerSegment = " << nPadsPerSegment << " QtotSector = " << qTotSectorOld -> GetMatrixArray()[sector] << endl;
-
} // end of global row loop
cerr<<"AliTPC warning: invalid segment ID ! "<<globalRowID<<endl;
continue;
}
- Int_t anodeSegmentID = param->GetWireSegment(sector,padRow);
+ TObjArray *arrTRF = (TObjArray*)timeResFunc.At(sector);
+ TGraphErrors *graphTRF = (TGraphErrors*)arrTRF->At(1);
+ Int_t wireSegmentID = param->GetWireSegment(sector,padRow);
+ Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(wireSegmentID));
const Float_t ampfactor = (sector<36)?factorIROC:factorOROC; // factor for the iontail which is ROC type dependent
AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector); // noise per given sector
digrow->Allocate(nTimeBins,nPads);
digrow->AllocateTrack(3);
+ Int_t localPad = 0;
Float_t q = 0.;
- Float_t qX = 0.;
- Float_t qIon = 0.;
- Float_t qold = 0.;
+ Float_t qXtalk = 0.;
+ Float_t qIonTail = 0.;
+ Float_t qOrig = 0.;
Int_t label[1000]; //stack for 300 events
Int_t labptr = 0;
Int_t nElems = nTimeBins*nPads; // element is a unit of a given row's pad-timebin space
Short_t *pdig1= digrow->GetDigits();
Int_t *ptr1= digrow->GetTracks() ;
// loop over elements i.e pad-timebin space of a row
- for (Int_t elem=0;elem<nElems; elem++)
- {
+ for (Int_t elem=0;elem<nElems; elem++) {
q=0;
labptr=0;
// looop over digits
- for (Int_t i=0;i<nInputs; i++) if (active[i])
- {
- // q += digarr[i]->GetDigitFast(rows,col);
+ for (Int_t i=0;i<nInputs; i++) if (active[i]){
q += *(pdig[i]);
-
- for (Int_t tr=0;tr<3;tr++)
- {
- // Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
+ for (Int_t tr=0;tr<3;tr++) {
Int_t lab = ptr[i][tr*nElems];
- if ( (lab > 1) && *(pdig[i])>zerosup)
- {
+ if ( (lab > 1) && *(pdig[i])>zerosup) {
label[labptr]=lab+masks[i];
labptr++;
}
}
Int_t padNumber = elem/nTimeBins;
Int_t timeBin = elem%nTimeBins;
-
- q/=16.; //conversion factor
+ localPad = padNumber-nPads/2;
+
Float_t gain = gainROC->GetValue(padRow,padNumber); // get gain for given - pad-row pad
//if (gain<0.5){
//printf("problem\n");
//}
q*= gain;
- qold = q;
- // Crosstalk correction
- qX = (*(TMatrixD*)crossTalkSignalArray.At(sector))[anodeSegmentID][timeBin];
-
- // Ion tail correction:
- // elem=padNumber*nTimeBins+timeBin;
- // lowerElem=elem-nIonTailBins;
- // if (lowerElem<0) lowerElem=0;
- // if (lowerElem in previospad) lowerElem = padNumber*nTimeBins;
- //
- // for (Int_t celem=elem-1; celem>lowerElem; celem--){
- // Int_t deltaT=elem-celem;
- //
- // }
- //
+ qOrig = q;
Float_t noisePad = noiseROC->GetValue(padRow,padNumber);
- // Float_t noise = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
- Float_t noise = pTPC->GetNoise();
- q+=noise*noisePad;
+ Float_t noise = pTPC->GetNoise()*noisePad;
+ if ( (q/16.+noise)> zerosup){
+ // Crosstalk correction
+ qXtalk = (*(TMatrixD*)crossTalkSignalArray.At(sector))[wireSegmentID][timeBin];
+ qTotPerSector = qTotSectorOld -> GetMatrixArray()[sector];
+
+ // Ion tail correction: being elem=padNumber*nTimeBins+timeBin;
+ Int_t lowerElem=elem-nIonTailBins;
+ Int_t zeroElem =(elem/nTimeBins)*nTimeBins;
+ if (zeroElem<0) zeroElem=0;
+ if (lowerElem<zeroElem) lowerElem=zeroElem;
+ //
+ qIonTail=0;
+ if (q>0 && recoParam->GetUseIonTailCorrection()){
+ for (Int_t i=0;i<nInputs; i++) if (active[i]){
+ Short_t *pdigC= digarr[i]->GetDigits();
+ for (Int_t celem=elem-1; celem>lowerElem; celem--){
+ //for Mesut - her we substact the ion tail
+ Double_t qCElem=pdigC[celem];
+ if (graphTRF->GetY()[elem-celem]<0)qIonTail+=qCElem*graphTRF->GetY()[elem-celem];
+ }
+ }
+ }
+ }
+ //
+ q -= qXtalk*recoParam->GetCrosstalkCorrection();
+ q+=qIonTail;
+ q/=16.; //conversion factor
+ q+=noise;
q=TMath::Nint(q); // round to the nearest integer
// fill info for degugging
- if (AliTPCReconstructor::StreamLevel()==1) {
+ if (AliTPCReconstructor::StreamLevel()==1 && qOrig > zerosup ) {
TTreeSRedirector &cstream = *fDebugStreamer;
cstream <<"ionTailXtalk"<<
- "sec=" << sector << //
- "row=" << globalRowID <<
- "pad=" << padNumber <<
- "tb=" << timeBin <<
- // "qsecOld.=" << qTotSectorOld << // vector of total charge in sector => number total charge in given sector
- //"qsecNew.=" << qTotSectorNew << //
- "qtpc=" << qTotTPC << // acumulated charge without crosstalk and ion tail in full TPC
- "qold=" << qold << // charge in given pad-row,pad,time-bin
- "qX=" << qX << // crosstal contribtion at given position
- "qIon=" << qIon << // ion tail cotribution from past signal
- "q=" << q << // q=qold-qX-qIon - to check sign of the effects
- // "mult=" << sec <<
+ "sector="<< sector<<
+ "globalRowID="<<globalRowID<<
+ "padRow="<< padRow<< //pad row
+ "wireSegmentID="<< wireSegmentID<< //wire segment 0-11, 0-3 in IROC 4-10 in OROC
+ "localPad="<<localPad<< // pad number -npads/2 .. npads/2
+ "padNumber="<<padNumber<< // pad number 0..npads
+ "timeBin="<< timeBin<< // time bin
+ "nPadsPerSegment="<<nPadsPerSegment<<// number of pads per wire segment
+ "qTotPerSector="<<qTotPerSector<< // total charge in sector
+ //
+ "qTotTPC="<<qTotTPC<< // acumulated charge without crosstalk and ion tail in full TPC
+ "qOrig="<< qOrig<< // charge in given pad-row,pad,time-bin
+ "q="<<q<< // q=qOrig-qXtalk-qIonTail - to check sign of the effects
+ "qXtalk="<<qXtalk<< // crosstal contribtion at given position
+ "qIonTail="<<qIonTail<< // ion tail cotribution from past signal
"\n";
}// dump the results to the debug streamer if in debug mode
pdig1++;
ptr1++;
}
+
+ if (AliTPCReconstructor::StreamLevel()==1 && qOrig > zerosup ) {
+ cout << " sector = " << sector << " row = " << padRow << " localPad = " << localPad << " SegmentID = " << wireSegmentID << " nPadsPerSegment = " << nPadsPerSegment << endl;
+ cout << " qXtalk = " << qXtalk ;
+ cout << " qOrig = " << qOrig ;
+ cout << " q = " << q ;
+ cout << " qsec = " << qTotPerSector ;
+ cout << " qTotTPC "<< qTotTPC << endl;
+ }
//
// glitch filter
//
#include "TH1.h"
#include "TClonesArray.h"
#include "TTreeStream.h"
+#include "TGrid.h"
class AliTPCclusterFast: public TObject {
public:
//
//
if (!fCl) fCl = new TClonesArray("AliTPCclusterFast",160);
+ //
+ // 0.) Init data structure
+ //
for (Int_t i=0;i<fN;i++){
AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
if (!cluster) cluster = new ((*fCl)[i]) AliTPCclusterFast;
cluster->Init();
}
-
+ //
+ // 1.) Create hits - with crosstalk diffusion
+ //
for (Int_t i=0;i<fN;i++){
Double_t tY = i*fAngleY;
Double_t tZ = i*fAngleZ;
cluster->SetParam(fMNprim,fDiff, fDiffLong, posY,posZ,fAngleY,fAngleZ);
//
cluster->GenerElectrons(cluster, clusterm, clusterp);
+ }
+ //
+ // 2.) make digitization
+ //
+ for (Int_t i=0;i<fN;i++){
+ AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
cluster->Digitize();
}
+
}
Double_t AliTPCtrackFast::CookdEdxNtot(Double_t f0,Float_t f1){
fPosY[i]=0;
fPosZ[i]=0;
fGain[i]=0;
+ fSec[i]=0;
}
}
// Generate number of secondary electrons
// copy of procedure implemented in geant
//
- const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2, EEND1=1E-6;
+ const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2; // EEND1=1E-6;
const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
const Double_t W=20.77E-9;
Float_t RAN = gRandom->Rndm();
//
const Int_t knMax=1000;
cl0->fNprim = gRandom->Poisson(cl0->fMNprim); //number of primary electrons
- cl0->fNtot=0; //total number of electrons
- cl0->fQtot=0; //total number of electrons after gain multiplification
+ // cl0->fNtot=0; //total number of electrons
+ // cl0->fQtot=0; //total number of electrons after gain multiplification
//
Double_t sumQ=0;
Double_t sumYQ=0;
Double_t sumZQ=0;
Double_t sumY2Q=0;
Double_t sumZ2Q=0;
- for (Int_t i=0;i<knMax;i++){
- cl0->fSec[i]=0;
- }
+ // for (Int_t i=0;i<knMax;i++){
+ // cl0->fSec[i]=0;
+ //}
for (Int_t iprim=0; iprim<cl0->fNprim;iprim++){
Float_t dN = cl0->GetNsec();
cl0->fSec[iprim]=dN;
cl->fQtot+=gg;
cl->fNtot++;
//
- // cl->sumQ+=gg;
-// cl->sumYQ+=gg*y;
-// cl->sumY2Q+=gg*y*y;
-// cl->sumZQ+=gg*z;
-// cl->sumZ2Q+=gg*z*z;
+ // cl->sumQ+=gg;
+ // cl->sumYQ+=gg*y;
+ // cl->sumY2Q+=gg*y*y;
+ // cl->sumZQ+=gg*z;
+ // cl->sumZ2Q+=gg*z*z;
if (cl->fNtot>=knMax) continue;
}
if (cl0->fNtot>=knMax) break;
+++ /dev/null
-Scripts to test reconstruction/calibration chain
-marian.ivanov@cern.ch
-
-$ALICE_ROOT/TPC/scripts/calibPassX
-Content:
-README.PassX
-runTrainPassX.sh
-runPassXJob.sh
-runPassX.sh
-recPass0GSI.C
-
-
-
-runTrainPassX.sh - test of the PassX reconstruction/calibration chain
- - filtering of the log files
- - merging of results
- - It is pseudo code
-
-submitPass0.sh - submit reconstruction/calibration job
-
-runPassXJob.sh - wrapper around the script $ALICE_ROOT/TPC/scripts/runPassX.sh
- to be run
- - intermediate results on the local machine
-
-runPassX.sh - copy from $ALICE_ROOT/ANALYSIS/macros/runPassX.sh
- - Modification needed:
- Provide the run number and number of events as argument
-
-recPass0GSI.C - reconstruction script
+++ /dev/null
-#
-# recursive merging
-#
-maxMerge=$1
-queue="$2"
-mask=$3
-output=$4
-reject=$5
-#
-counter=0;
-counter2=0;
-wdir=`pwd`
-rm -rf merge*
-mkdir merge$counter2
-cd merge$counter
-for a in `cat ../calib.list`; do
- let counter=counter+1;
- echo $counter $counter2
- echo $a >>calib.list
- if [ $counter -gt $maxMerge ] ; then
- echo bsub -q $queue -oo outMerge.log $ALICE_ROOT/ANALYSIS/CalibMacros/MergeCalibration/mergeCustom.C\(\"$output\",\"$mask\",\"$5\"\);
- cat calib.list
- bsub -q $queue -oo outMerge.log aliroot -b -q $ALICE_ROOT/ANALYSIS/CalibMacros/MergeCalibration/mergeCustom.C\(\"calib.list\",\"$output\",\"$mask\",\"$5\"\)
- let counter2=counter2+1;
- let counter=0;
- cd $wdir
- mkdir merge$counter2
- cd merge$counter2
- if [ -e calib.list ]; then
- rm calib.list
- fi;
- fi;
-done;
-
-bsub -q $queue -oo outMerge.log aliroot -b -q $ALICE_ROOT/ANALYSIS/CalibMacros/MergeCalibration/mergeCustom.C\(\"calib.list\",\"$output\",\"$mask\",\"$5\"\)
-
+++ /dev/null
-//
-// rec.C to be used for pass0
-//
-
-void rec(const char *filename="raw.root",Int_t nevents=-1)
-{
- // Load some system libs for Grid and monitoring
- // Set the CDB storage location
- AliCDBManager * man = AliCDBManager::Instance();
- //man->SetDefaultStorage("raw://");
- man->SetDefaultStorage("local:///lustre/alice/alien/alice/data/2010/OCDB/");
- // Reconstruction settings
- AliReconstruction rec;
-
- // Set protection against too many events in a chunk (should not happen)
- if (nevents>0) rec.SetEventRange(0,nevents);
-
- // Switch off HLT until the problem with schema evolution resolved
- //rec.SetRunReconstruction("ALL-HLT");
- //
- // QA options
- //
- AliQAManager *qam = AliQAManager::QAManager(AliQAv1::kRECMODE) ;
- rec.SetRunQA(":");
- rec.SetRunGlobalQA(kFALSE);
-
- // AliReconstruction settings
- rec.SetWriteESDfriend(kTRUE);
- rec.SetWriteAlignmentData();
- rec.SetInput(filename);
- rec.SetUseTrackingErrorsForAlignment("ITS");
- rec.SetRunReconstruction("ITS TPC TRD TOF");
- rec.SetFillESD("ITS TPC TRD TOF");
-
- // switch off cleanESD
- rec.SetCleanESD(kFALSE);
-
- AliLog::Flush();
- rec.Run();
-
-}
-
-
+++ /dev/null
-#!/bin/bash
-
-# Script to run:
-# 1. reconstruction
-# 2. calibration and friend track filtering
-# 3. tag creation
-#
-# Files assumed to be in working directory:
-# rec.C - reconstruction macro
-# runCalibTrain.C - calibration/filtering macro
-# Arguments:
-# 1 - raw data file name
-# 2 - number of events to be processed
-# 3 - run number
-
-# example:
-# runPassX.sh raw.root 50 104892
-
-#ALIEN setting
-#entries=1000
-# $1 = raw input filename
-#runnum=`echo $1 | cut -d "/" -f 6`
-
-#Local setting : setting variables
-
-entries=$2
-runnum=$3
-source $HOME/alienSetup.sh
-
-echo File to be processed $1
-echo Number of events to be processed $entries
-echo Run mumber $runnum
-
-echo ALICE_ROOT = $ALICE_ROOT
-echo AliROOT = $AliROOT
-cp $ALICE_ROOT/.rootrc ~/.rootrc
-cp $ALICE_ROOT/.rootrc $HOME
-#cat $HOME/.rootrc
-export GRID_TOKEN=OK
-
-echo ">>>>>>>>> PATH is..."
-echo $PATH
-echo ">>>>>>>>> LD_LIBRARY_PATH is..."
-echo $LD_LIBRARY_PATH
-echo ">>>>>>>>> rec.C is..."
-cat rec.C
-echo
-
-
-echo
-echo ">>>>>>> Running AliRoot to reconstruct $1. Run number is $runnum..."
-echo
-if [ -e AliESDs.root ]; then
- echo AliESDs.root exist
- ls -al AliESD*
-else
- echo aliroot -l -b -q rec.C\(\"$1\",$2\) 2>&1 | tee rec.log
- aliroot -l -b -q rec.C\(\"$1\",$2\) 2>&1 | tee rec.log
- echo aliroot -l -b -q tag.C\(\) 2>&1 | tee tag.log
- aliroot -l -b -q tag.C\(\) 2>&1 | tee tag.log
-fi
-
-echo
-echo ">>>>>>> Running AliRoot to make calibration..."
-echo
-echo aliroot -l -b -q runCalibTrain.C\($runnum\) 2>&1 | tee calib.log
-aliroot -l -b -q runCalibTrain.C\($runnum\) 2>&1 | tee calib.log
-
-echo
-echo ">>>>>>> Running AliRoot to generate Tags..."
-echo
-echo aliroot -l -b -q tag.C\(\) 2>&1 | tee tag.log
-#aliroot -l -b -q tag.C\(\) 2>&1 | tee tag.log
+++ /dev/null
-#
-# run PassX recosntruction/calibration jobs
-#
-# Needed in order to run the jobs locally
-#
-# 1 - raw data file name
-# 2 - number of events to be processed
-# 3 - run number
-
-
-echo Input file $1
-echo Run number $3
-echo Events to reconstruct/calibrate $2
-echo HOSTNAME $HOSTNAME
-
-outdir=`pwd`
-tmpdir=/tmp/$USER/`echo $outdir | sed s_/_x_g`
-mkdirhier $tmpdir
-echo tmpdir = $tmpdir
-#
-cp $outdir/*.C $tmpdir/
-cp $outdir/*.sh $tmpdir/
-cp $outdir/AliESDs.root $tmpdir/
-cp $outdir/AliESDfriends.root $tmpdir/
-
-cd $tmpdir
-#
-#
-#
-echo tmpdir = `pwd`
-echo outdir = $outdir
-ls
-
-cp $1 .
-
-if [ -e root_archive.zip ]; then
- echo unzip reconstructed file
- unzip root_archive.zip
- rm -f AliESDfriends_v1.root
- rm root_archive.zip
-fi;
-
-echo runPassX.sh $1 $2 $3
-runPassX.sh $1 $2 $3
-
-ls -alrt
-rm *ebug.root
-cp -r $tmpdir/* $outdir/
-rm -rf $tmpdir
-
-
+++ /dev/null
-#
-# Sequence of steps to test Pass0 and PassX reconstruction/calibration which run on GRID
-# by default
-#
-# Semi automatic test performed on the batch farm
-# Important features:
-# 1. Parsing of the log files
-# 2. Parsing stack traces
-
-# author: marian.ivanov@cern.ch
-
-
-#
-# 0. copy a template setup
-# To be modified if necessary
-# "standard scripts
-cp $ALICE_ROOT/ANALYSIS/macros/runCalibTrain.C .
-# scipts to run on batch farm
-cp $ALICE_ROOT/TPC/scripts/calibPassX/recPass0GSI.C rec.C
-cp $ALICE_ROOT/TPC/scripts/calibPassX/runPassX.sh .
-cp $ALICE_ROOT/TPC/scripts/calibPassX/submitPass0.sh .
-cp $ALICE_ROOT/TPC/scripts/calibPassX/runPassXJob.sh .
-
-cp ../lists/run.list .
-cp ../lists/raw.list .
-
-#
-# 1. Make workspace - directory structure with run and raw lists
-#
-$ALICE_ROOT/TPC/scripts/makeWorkspace.sh run.list
-#
-# 1.a clean the workspace all
-#
-find `pwd` | grep AliESD | xargs rm
-find `pwd` | grep out | xargs rm
-find `pwd` | grep log | xargs rm
-#
-# 1.b clean workscpace - rm root files if not tags found
-#
-find `pwd` | grep AliESDs.root > esd.txt
-for efile in `cat esd.txt`; do
- dname=`dirname $efile`
- status=`ls $dname/*tag.root`;
- echo CHECK $efile $status
- if [ -z $status ]; then
- echo NON OK rm -r $dname/*.root
- rm -r $dname/*.root
- else
- echo IS OK $status
- fi;
-done
-
-#
-# 2. Run reconstruction/calibration
-#
-bgroup=/recPass0/`pwd | xargs basename`
-bgadd $bgroup
-nEvents=500000
-submitPass0.sh run.list "alice-t3_8h -g $bgroup" $nEvents | tee submit.log
-
-
-#
-# 3. Run merging - run level
-#
-for arun in `cat run.list`; do
- cd $arun
- find `pwd`/ | grep AliESDfriends_v1.root > calib.list
- echo bsub -q alice-t3_8h -o outMerge.log aliroot -b -q $ALICE_ROOT/ANALYSIS/macros/mergeCalibObjects.C
- bsub -q alice-t3_8h -o outMerge.log aliroot -b -q $ALICE_ROOT/ANALYSIS/macros/mergeCalibObjects.C
- cd ../
-done;
-
-#
-# 4. Merge all
-#
-ls `pwd`/*/CalibObjects.root > calib.list
-aliroot -b -q $ALICE_ROOT/ANALYSIS/macros/mergeCalibObjects.C
-
-#
-# 5. filter reconstruction logs, make statistic of problems
-#
-mkdirhier logs/reco
-cd logs/reco
-find `pwd`/../../ | grep rec.log > errRec.log
-$ALICE_ROOT/TPC/scripts/filterRecLog.sh
-cd ../..
-
-#
-# 6. filter calibration logs, make statistic of problems
-#
-mkdirhier logs/calib
-cd logs/calib
-find `pwd`/../../ | grep calib.log > errRec.log
-$ALICE_ROOT/TPC/scripts/filterRecLog.sh
-cd ../..
-
-
-#
-# 7. filter debug streamer
-#
-mkdir debug
-cd debug
-find `pwd`/../*/ | grep V6 | grep .root > debugall.txt
-cat debugall.txt | grep calibTimeDebug > timeitstpc.txt
-
-
-
-
+++ /dev/null
-#!/bin/bash
-#
-# Submit jobs for reconstruction/calibration
-# Test - Equivalent of pass0/passX on the grid
-#
-
-# Parameters:
-# 1 - queue name
-# 2 - run number
-# 3 - number of events to be reconstructed
-# 4 - run.list
-# Example: submitPass0.sh run.list "alice-t3_8h" 1000
-
-runList=$1
-qName=$2
-nEvents=$3
-workDir=`pwd`
-for arun in `cat $runList`; do
- echo
- echo Run number $arun
- echo
- #
- cd $workDir/$arun
- cp $workDir/*.C .
- cp $workDir/*.sh .
- cat ../../lists/raw.list | grep $arun >raw.list
- #
- #
- for afile in `cat raw.list | grep -v tag`; do
- echo $afile;
- cdir0=`basename $afile | sed s_.root__`
- mkdir $cdir0
- cd $cdir0
- cp $workDir/$arun/*.C .
- cp $workDir/$arun/*.sh .
- echo
- echo
- pwd
- echo bsub -q $qName -o out${cdir0}.log runPassXJob.sh $afile $nEvents $arun
- bsub -q $qName -o out${cdir0}.log runPassXJob.sh $afile $nEvents $arun
- cd ../
- cd $workDir/$arun
- done;
- cd $workDir
-done
-
-
-
+++ /dev/null
-#
-# parameters:
-# 1 - basedir
-# Example:
-# source /usr/local/grid/AliRoot/HEAD0108/TPC/scripts/tpcPass0Env.sh `pwd`
-export balice=/u/miranov/.balice
-export aliensetup=$HOME/alienSetup.sh
-export PASS0_DIR=/usr/local/grid/AliRoot/HEAD0108
-source $balice
-#source $aliensetup >aliensetup.log
-#
-# Test setup
-#
-export workdir=$1
-if [ ! -n length ]; then
- echo \############################
- echo Directory was not specified. Exiting
- echo \############################
- return;
-fi;
-if [ ! -r $workdir/lists/esd.list ] ; then
- echo \############################
- echo File esd list does not exist. Exiting
- echo \############################
- return;
-fi;
-if [ ! -r $workdir/lists/run.list ] ; then
- echo \############################
- echo File run list does not exist. Exiting
- echo \############################
- return;
-fi;
-
-#
-# Make directories
-#
-cd $workdir
-chgrp -R alice $workdir
-chmod -R g+rwx $workdir
-chmod -R o+rx $workdir
-mkdirhier $workdir/calibNoDrift
-mkdirhier $workdir/calibNoRefit
-mkdirhier $workdir/calibQA
-#
-#modify ConfigOCDB.C
-#
-# copy predefined Config files
-#
-cp $PASS0_DIR/TPC/macros/CalibrateTPC.C calibNoDrift/CalibrateTPC.C
-cat $PASS0_DIR/TPC/macros/CalibrateTPC.C | grep -v AddCalibCalib\(task\) > calibNoRefit/CalibrateTPC.C
-cp $PASS0_DIR/TPC/macros/CalibrateTPC.C calibQA/CalibrateTPC.C
-cp $PASS0_DIR/TPC/macros/ConfigOCDBNoDrift.C calibNoDrift/ConfigOCDB.C
-cp $PASS0_DIR/TPC/macros/ConfigOCDBNoRefit.C calibNoRefit/ConfigOCDB.C
-cp $PASS0_DIR/TPC/macros/ConfigOCDBQA.C calibQA/ConfigOCDB.C
-cp lists/*.list calibNoDrift/
-cp lists/*.list calibNoRefit/
-cp lists/*.list calibQA/
-ln -sf $balice calibNoDrift/balice.sh
-ln -sf $balice calibNoRefit/balice.sh
-ln -sf $balice calibQA/balice.sh
-ln -sf $aliensetup calibNoDrift/alienSetup.sh
-ln -sf $aliensetup calibNoRefit/alienSetup.sh
-ln -sf $aliensetup calibQA/alienSetup.sh
-# make workspaces
-#
-cd $workdir/calibNoDrift
-$PASS0_DIR/TPC/scripts/makeWorkspace.sh run.list
-$PASS0_DIR/TPC/scripts/submitCalib.sh run.list alice-t3 20
-cd $workdir/calibNoRefit
-$PASS0_DIR/TPC/scripts/makeWorkspace.sh run.list
-$PASS0_DIR/TPC/scripts/submitCalib.sh run.list alice-t3 20
-cd $workdir/calibQA
-$PASS0_DIR/TPC/scripts/makeWorkspace.sh run.list
-$PASS0_DIR/TPC/scripts/submitCalib.sh run.list alice-t3 20
-cd $workdir/
-#
-#
-#
Double_t ratioPointsHighPt = highPt->GetBranch("friendTrack.fCalibContainer")->GetZipBytes()/Double_t(0.00001+highPt->GetZipBytes());
printf("#UnitTest:\tAliAnalysisTaskFiltered\tRatioPointsV0\t%f\n",ratioPointsV0);
printf("#UnitTest:\tAliAnalysisTaskFiltered\tRatioPointsHighPt\t%f\n",ratioPointsHighPt);
+ //
+ // a.) Check track correspondence
+ //
+ Int_t entries= highPt->Draw("(friendTrack.fTPCOut.fP[3]-esdTrack.fIp.fP[3])/sqrt(friendTrack.fTPCOut.fC[9]+esdTrack.fIp.fC[9])","friendTrack.fTPCOut.fP[3]!=0","");
+ // here we should check if the tracks
+ Double_t pulls=TMath::RMS(entries, highPt->GetV1());
+ printf("#UnitTest:\tAliAnalysisTaskFiltered\tFriendPull\t%2.4f\n",pulls);
+ printf("#UnitTest:\tAliAnalysisTaskFiltered\tFriendOK\t%d\n",pulls<10);
}
#get path to input list
inputListfiles=$TestData_pPb
#get scale number for tracks
- filterT=${2-100}
+ filterT=${2-20}
#get scale number for V0s
- filterV=${3-10}
-#get scale number of riends
- filterFriend=${4--10}
+ filterV=${3-2}
+#get scale number of friends
+ filterFriend=${4--1}
#get OCDB path (optional)
OCDBpath=${5-"\"$OCDBPath_pPb\""}
#get max number of files
--- /dev/null
+add_test("testdEdx" runtest.sh)
\ No newline at end of file
--- /dev/null
+//Configuration of simulation
+
+enum PprRun_t
+{
+ test50, testdEdx,
+ kParam_8000, kParam_4000, kParam_2000,
+ kHijing_cent1, kHijing_cent2,
+ kHijing_per1, kHijing_per2, kHijing_per3, kHijing_per4, kHijing_per5,
+ kHijing_jj25, kHijing_jj50, kHijing_jj75, kHijing_jj100, kHijing_jj200,
+ kHijing_gj25, kHijing_gj50, kHijing_gj75, kHijing_gj100, kHijing_gj200,
+ kHijing_pA, kPythia6,
+ kPythia6Jets20_24, kPythia6Jets24_29, kPythia6Jets29_35,
+ kPythia6Jets35_42, kPythia6Jets42_50, kPythia6Jets50_60,
+ kPythia6Jets60_72, kPythia6Jets72_86, kPythia6Jets86_104,
+ kPythia6Jets104_125, kPythia6Jets125_150, kPythia6Jets150_180,
+ kD0PbPb5500, kCharmSemiElPbPb5500, kBeautySemiElPbPb5500,
+ kCocktailTRD, kPyJJ, kPyGJ,
+ kMuonCocktailCent1, kMuonCocktailPer1, kMuonCocktailPer4,
+ kMuonCocktailCent1HighPt, kMuonCocktailPer1HighPt, kMuonCocktailPer4HighPt,
+ kMuonCocktailCent1Single, kMuonCocktailPer1Single, kMuonCocktailPer4Single,
+ kFlow_2_2000, kFlow_10_2000, kFlow_6_2000, kFlow_6_5000,
+ kHIJINGplus, kRunMax
+};
+
+const char* pprRunName[] = {
+ "test50", "testdEdx",
+ "kParam_8000", "kParam_4000", "kParam_2000",
+ "kHijing_cent1", "kHijing_cent2",
+ "kHijing_per1", "kHijing_per2", "kHijing_per3", "kHijing_per4",
+ "kHijing_per5",
+ "kHijing_jj25", "kHijing_jj50", "kHijing_jj75", "kHijing_jj100",
+ "kHijing_jj200",
+ "kHijing_gj25", "kHijing_gj50", "kHijing_gj75", "kHijing_gj100",
+ "kHijing_gj200", "kHijing_pA", "kPythia6",
+ "kPythia6Jets20_24", "kPythia6Jets24_29", "kPythia6Jets29_35",
+ "kPythia6Jets35_42", "kPythia6Jets42_50", "kPythia6Jets50_60",
+ "kPythia6Jets60_72", "kPythia6Jets72_86", "kPythia6Jets86_104",
+ "kPythia6Jets104_125", "kPythia6Jets125_150", "kPythia6Jets150_180",
+ "kD0PbPb5500", "kCharmSemiElPbPb5500", "kBeautySemiElPbPb5500",
+ "kCocktailTRD", "kPyJJ", "kPyGJ",
+ "kMuonCocktailCent1", "kMuonCocktailPer1", "kMuonCocktailPer4",
+ "kMuonCocktailCent1HighPt", "kMuonCocktailPer1HighPt", "kMuonCocktailPer4HighPt",
+ "kMuonCocktailCent1Single", "kMuonCocktailPer1Single", "kMuonCocktailPer4Single",
+ "kFlow_2_2000", "kFlow_10_2000", "kFlow_6_2000", "kFlow_6_5000", "kHIJINGplus"
+};
+
+enum PprRad_t
+{
+ kGluonRadiation, kNoGluonRadiation
+};
+
+enum PprTrigConf_t
+{
+ kDefaultPPTrig, kDefaultPbPbTrig
+};
+
+const char * pprTrigConfName[] = {
+ "p-p","Pb-Pb"
+};
+
+// This part for configuration
+
+//static PprRun_t srun = kPythia6;
+static PprRun_t srun = testdEdx;
+static PprRad_t srad = kGluonRadiation;
+static AliMagF::BMap_t smag = AliMagF::k5kG;
+static Int_t sseed = 0; //Set 0 to use the current time
+static PprTrigConf_t strig = kDefaultPPTrig; // default pp trigger configuration
+
+// Comment line
+static TString comment;
+
+// Functions
+Float_t EtaToTheta(Float_t arg);
+AliGenerator* GeneratorFactory(PprRun_t srun);
+AliGenHijing* HijingStandard();
+AliGenGeVSim* GeVSimStandard(Float_t, Float_t);
+void ProcessEnvironmentVars();
+
+void Config()
+{
+ // ThetaRange is (0., 180.). It was (0.28,179.72) 7/12/00 09:00
+ // Theta range given through pseudorapidity limits 22/6/2001
+
+ // Get settings from environment variables
+ ProcessEnvironmentVars();
+
+ // Set Random Number seed
+ gRandom->SetSeed(sseed);
+ cout<<"Seed for random number generation= "<<gRandom->GetSeed()<<endl;
+
+
+ // libraries required by geant321 and Pythia: loaded in sim.C
+
+ new TGeant3TGeo("C++ Interface to Geant3");
+
+ // Output every 100 tracks
+
+ TVirtualMC * vmc = TVirtualMC::GetMC();
+
+ ((TGeant3*)vmc)->SetSWIT(4,100);
+
+ AliRunLoader* rl=0x0;
+
+ AliLog::Message(AliLog::kInfo, "Creating Run Loader", "", "", "Config()"," ConfigPPR.C", __LINE__);
+
+ rl = AliRunLoader::Open("galice.root",
+ AliConfig::GetDefaultEventFolderName(),
+ "recreate");
+ if (rl == 0x0)
+ {
+ gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
+ return;
+ }
+ rl->SetCompressionLevel(2);
+ rl->SetNumberOfEventsPerFile(100);
+ gAlice->SetRunLoader(rl);
+
+ // Set the trigger configuration
+ AliSimulation::Instance()->SetTriggerConfig(pprTrigConfName[strig]);
+ cout<<"Trigger configuration is set to "<<pprTrigConfName[strig]<<endl;
+
+ //
+ // Set External decayer
+ AliDecayer *decayer = new AliDecayerPythia();
+
+
+ switch (srun) {
+ case kD0PbPb5500:
+ decayer->SetForceDecay(kHadronicD);
+ break;
+ case kCharmSemiElPbPb5500:
+ decayer->SetForceDecay(kSemiElectronic);
+ break;
+ case kBeautySemiElPbPb5500:
+ decayer->SetForceDecay(kSemiElectronic);
+ break;
+ default:
+ decayer->SetForceDecay(kAll);
+ break;
+ }
+ decayer->Init();
+ vmc->SetExternalDecayer(decayer);
+ //
+ //
+ //=======================================================================
+ //
+ //=======================================================================
+ // ************* STEERING parameters FOR ALICE SIMULATION **************
+ // --- Specify event type to be tracked through the ALICE setup
+ // --- All positions are in cm, angles in degrees, and P and E in GeV
+
+ vmc->SetProcess("DCAY",1);
+ vmc->SetProcess("PAIR",1);
+ vmc->SetProcess("COMP",1);
+ vmc->SetProcess("PHOT",1);
+ vmc->SetProcess("PFIS",0);
+ vmc->SetProcess("DRAY",0);
+ vmc->SetProcess("ANNI",1);
+ vmc->SetProcess("BREM",1);
+ vmc->SetProcess("MUNU",1);
+ vmc->SetProcess("CKOV",1);
+ vmc->SetProcess("HADR",1);
+ vmc->SetProcess("LOSS",2);
+ vmc->SetProcess("MULS",1);
+ vmc->SetProcess("RAYL",1);
+
+ Float_t cut = 1.e-3; // 1MeV cut by default
+ Float_t tofmax = 1.e10;
+
+ vmc->SetCut("CUTGAM", cut);
+ vmc->SetCut("CUTELE", cut);
+ vmc->SetCut("CUTNEU", cut);
+ vmc->SetCut("CUTHAD", cut);
+ vmc->SetCut("CUTMUO", cut);
+ vmc->SetCut("BCUTE", cut);
+ vmc->SetCut("BCUTM", cut);
+ vmc->SetCut("DCUTE", cut);
+ vmc->SetCut("DCUTM", cut);
+ vmc->SetCut("PPCUTM", cut);
+ vmc->SetCut("TOFMAX", tofmax);
+
+ // Generator Configuration
+ AliGenerator* gener = GeneratorFactory(srun);
+ gener->SetOrigin(0, 0, 0); // vertex position
+ gener->SetSigma(0, 0, 5.3); // Sigma in (X,Y,Z) (cm) on IP position
+ gener->SetCutVertexZ(1.); // Truncate at 1 sigma
+ gener->SetVertexSmear(kPerEvent);
+ gener->SetTrackingFlag(1);
+ gener->Init();
+
+ if (smag == AliMagF::k2kG) {
+ comment = comment.Append(" | L3 field 0.2 T");
+ } else if (smag == AliMagF::k5kG) {
+ comment = comment.Append(" | L3 field 0.5 T");
+ }
+
+
+ if (srad == kGluonRadiation)
+ {
+ comment = comment.Append(" | Gluon Radiation On");
+
+ } else {
+ comment = comment.Append(" | Gluon Radiation Off");
+ }
+
+ printf("\n \n Comment: %s \n \n", comment.Data());
+
+
+// Field
+ TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., smag));
+
+ rl->CdGAFile();
+//
+ Int_t iABSO = 1;
+ Int_t iDIPO = 1;
+ Int_t iFMD = 1;
+ Int_t iFRAME = 1;
+ Int_t iHALL = 1;
+ Int_t iITS = 1;
+ Int_t iMAG = 1;
+ Int_t iMUON = 1;
+ Int_t iPHOS = 1;
+ Int_t iPIPE = 1;
+ Int_t iPMD = 1;
+ Int_t iHMPID = 1;
+ Int_t iSHIL = 1;
+ Int_t iT0 = 1;
+ Int_t iTOF = 1;
+ Int_t iTPC = 1;
+ Int_t iTRD = 1;
+ Int_t iZDC = 1;
+ Int_t iEMCAL = 1;
+ Int_t iVZERO = 1;
+ Int_t iACORDE = 1;
+ Int_t iAD = 0;
+
+ //=================== Alice BODY parameters =============================
+ AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
+
+
+ if (iMAG)
+ {
+ //=================== MAG parameters ============================
+ // --- Start with Magnet since detector layouts may be depending ---
+ // --- on the selected Magnet dimensions ---
+ AliMAG *MAG = new AliMAG("MAG", "Magnet");
+ }
+
+
+ if (iABSO)
+ {
+ //=================== ABSO parameters ============================
+ AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");
+ }
+
+ if (iDIPO)
+ {
+ //=================== DIPO parameters ============================
+
+ AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");
+ }
+
+ if (iHALL)
+ {
+ //=================== HALL parameters ============================
+
+ AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");
+ }
+
+
+ if (iFRAME)
+ {
+ //=================== FRAME parameters ============================
+
+ AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
+ FRAME->SetHoles(1);
+ }
+
+ if (iSHIL)
+ {
+ //=================== SHIL parameters ============================
+
+ AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");
+ }
+
+
+ if (iPIPE)
+ {
+ //=================== PIPE parameters ============================
+
+ AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
+ }
+
+ if (iITS)
+ {
+ //=================== ITS parameters ============================
+
+ AliITS *ITS = new AliITSv11("ITS","ITS v11");
+ }
+
+ if (iTPC)
+ {
+ //============================ TPC parameters =====================
+ AliTPC *TPC = new AliTPCv2("TPC", "Default");
+ }
+
+
+ if (iTOF) {
+ //=================== TOF parameters ============================
+ AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
+ }
+
+
+ if (iHMPID)
+ {
+ //=================== HMPID parameters ===========================
+ AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");
+
+ }
+
+
+ if (iZDC)
+ {
+ //=================== ZDC parameters ============================
+
+ AliZDC *ZDC = new AliZDCv4("ZDC", "normal ZDC");
+ }
+
+ if (iTRD)
+ {
+ //=================== TRD parameters ============================
+
+ AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
+ }
+
+ if (iFMD)
+ {
+ //=================== FMD parameters ============================
+ AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
+ }
+
+ if (iMUON)
+ {
+ //=================== MUON parameters ===========================
+ // New MUONv1 version (geometry defined via builders)
+ AliMUON *MUON = new AliMUONv1("MUON", "default");
+ }
+ //=================== PHOS parameters ===========================
+
+ if (iPHOS)
+ {
+ AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");
+ }
+
+
+ if (iPMD)
+ {
+ //=================== PMD parameters ============================
+ AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
+ }
+
+ if (iT0)
+ {
+ //=================== T0 parameters ============================
+ AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
+ }
+
+ if (iEMCAL)
+ {
+ //=================== EMCAL parameters ============================
+ AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETEV1");
+ }
+
+ if (iACORDE)
+ {
+ //=================== ACORDE parameters ============================
+ AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
+ }
+
+ if (iVZERO)
+ {
+ //=================== VZERO parameters ============================
+ AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
+ }
+
+ if (iAD)
+ {
+ //=================== AD parameters ============================
+ AliAD *AD = new AliADv1("AD", "normal AD test");
+ AD->SetADAToInstalled(kTRUE);
+ AD->SetADCToInstalled(kTRUE);
+ }
+
+
+}
+
+Float_t EtaToTheta(Float_t arg){
+ return (180./TMath::Pi())*2.*atan(exp(-arg));
+}
+
+
+
+AliGenerator* GeneratorFactory(PprRun_t srun) {
+ Int_t isw = 3;
+ if (srad == kNoGluonRadiation) isw = 0;
+
+
+ AliGenerator * gGener = 0x0;
+ switch (srun) {
+ case test50:
+ {
+ comment = comment.Append(":HIJINGparam test 50 particles");
+ AliGenHIJINGpara *gener = new AliGenHIJINGpara(50);
+ gener->SetMomentumRange(0, 999999.);
+ gener->SetPhiRange(0., 360.);
+ // Set pseudorapidity range from -8 to 8.
+ Float_t thmin = EtaToTheta(8); // theta min. <---> eta max
+ Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min
+ gener->SetThetaRange(thmin,thmax);
+ gGener=gener;
+ }
+ break;
+ case testdEdx:
+ {
+ // generator for the dEdx simulation study
+ gSystem->Getenv("TestdEdxNTracks");
+ Int_t ntracks =50;
+ if (gSystem->Getenv("TestdEdxNTracks")) ntracks= atoi(gSystem->Getenv("TestdEdxNTracks"));
+ comment = comment.Append(":HIJINGparam test N particles");
+ AliGenHIJINGpara *gener = new AliGenHIJINGpara(ntracks);
+ gener->SetMomentumRange(0, 999999.);
+ gener->SetPhiRange(0., 360.);
+ // Set pseudorapidity range from -8 to 8.
+ Float_t thmin = EtaToTheta(2); // theta min. <---> eta max
+ Float_t thmax = EtaToTheta(-2); // theta max. <---> eta min
+ gener->SetThetaRange(thmin,thmax);
+ gGener=gener;
+ }
+ break;
+
+ case kParam_8000:
+ {
+ comment = comment.Append(":HIJINGparam N=8000");
+ AliGenHIJINGpara *gener = new AliGenHIJINGpara(86030);
+ gener->SetMomentumRange(0, 999999.);
+ gener->SetPhiRange(0., 360.);
+ // Set pseudorapidity range from -8 to 8.
+ Float_t thmin = EtaToTheta(8); // theta min. <---> eta max
+ Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min
+ gener->SetThetaRange(thmin,thmax);
+ gGener=gener;
+ }
+ break;
+ case kParam_4000:
+ {
+ comment = comment.Append("HIJINGparam N=4000");
+ AliGenHIJINGpara *gener = new AliGenHIJINGpara(43015);
+ gener->SetMomentumRange(0, 999999.);
+ gener->SetPhiRange(0., 360.);
+ // Set pseudorapidity range from -8 to 8.
+ Float_t thmin = EtaToTheta(8); // theta min. <---> eta max
+ Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min
+ gener->SetThetaRange(thmin,thmax);
+ gGener=gener;
+ }
+ break;
+ case kParam_2000:
+ {
+ comment = comment.Append("HIJINGparam N=2000");
+ AliGenHIJINGpara *gener = new AliGenHIJINGpara(21507);
+ gener->SetMomentumRange(0, 999999.);
+ gener->SetPhiRange(0., 360.);
+ // Set pseudorapidity range from -8 to 8.
+ Float_t thmin = EtaToTheta(8); // theta min. <---> eta max
+ Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min
+ gener->SetThetaRange(thmin,thmax);
+ gGener=gener;
+ }
+ break;
+//
+// Hijing Central
+//
+ case kHijing_cent1:
+ {
+ comment = comment.Append("HIJING cent1");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ gGener=gener;
+ }
+ break;
+ case kHijing_cent2:
+ {
+ comment = comment.Append("HIJING cent2");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 2.);
+ gGener=gener;
+ }
+ break;
+//
+// Hijing Peripheral
+//
+ case kHijing_per1:
+ {
+ comment = comment.Append("HIJING per1");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(5., 8.6);
+ gGener=gener;
+ }
+ break;
+ case kHijing_per2:
+ {
+ comment = comment.Append("HIJING per2");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(8.6, 11.2);
+ gGener=gener;
+ }
+ break;
+ case kHijing_per3:
+ {
+ comment = comment.Append("HIJING per3");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(11.2, 13.2);
+ gGener=gener;
+ }
+ break;
+ case kHijing_per4:
+ {
+ comment = comment.Append("HIJING per4");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(13.2, 15.);
+ gGener=gener;
+ }
+ break;
+ case kHijing_per5:
+ {
+ comment = comment.Append("HIJING per5");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(15., 100.);
+ gGener=gener;
+ }
+ break;
+//
+// Jet-Jet
+//
+ case kHijing_jj25:
+ {
+ comment = comment.Append("HIJING Jet 25 GeV");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ // trigger
+ gener->SetTrigger(1);
+ gener->SetPtJet(25.);
+ gener->SetRadiation(isw);
+ gener->SetSimpleJets(!isw);
+ gener->SetJetEtaRange(-0.3,0.3);
+ gener->SetJetPhiRange(75., 165.);
+ gGener=gener;
+ }
+ break;
+
+ case kHijing_jj50:
+ {
+ comment = comment.Append("HIJING Jet 50 GeV");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ // trigger
+ gener->SetTrigger(1);
+ gener->SetPtJet(50.);
+ gener->SetRadiation(isw);
+ gener->SetSimpleJets(!isw);
+ gener->SetJetEtaRange(-0.3,0.3);
+ gener->SetJetPhiRange(75., 165.);
+ gGener=gener;
+ }
+ break;
+
+ case kHijing_jj75:
+ {
+ comment = comment.Append("HIJING Jet 75 GeV");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ // trigger
+ gener->SetTrigger(1);
+ gener->SetPtJet(75.);
+ gener->SetRadiation(isw);
+ gener->SetSimpleJets(!isw);
+ gener->SetJetEtaRange(-0.3,0.3);
+ gener->SetJetPhiRange(75., 165.);
+ gGener=gener;
+ }
+ break;
+
+ case kHijing_jj100:
+ {
+ comment = comment.Append("HIJING Jet 100 GeV");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ // trigger
+ gener->SetTrigger(1);
+ gener->SetPtJet(100.);
+ gener->SetRadiation(isw);
+ gener->SetSimpleJets(!isw);
+ gener->SetJetEtaRange(-0.3,0.3);
+ gener->SetJetPhiRange(75., 165.);
+ gGener=gener;
+ }
+ break;
+
+ case kHijing_jj200:
+ {
+ comment = comment.Append("HIJING Jet 200 GeV");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ // trigger
+ gener->SetTrigger(1);
+ gener->SetPtJet(200.);
+ gener->SetRadiation(isw);
+ gener->SetSimpleJets(!isw);
+ gener->SetJetEtaRange(-0.3,0.3);
+ gener->SetJetPhiRange(75., 165.);
+ gGener=gener;
+ }
+ break;
+//
+// Gamma-Jet
+//
+ case kHijing_gj25:
+ {
+ comment = comment.Append("HIJING Gamma 25 GeV");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ // trigger
+ gener->SetTrigger(2);
+ gener->SetPtJet(25.);
+ gener->SetRadiation(isw);
+ gener->SetSimpleJets(!isw);
+ gener->SetJetEtaRange(-0.12, 0.12);
+ gener->SetJetPhiRange(220., 320.);
+ gGener=gener;
+ }
+ break;
+
+ case kHijing_gj50:
+ {
+ comment = comment.Append("HIJING Gamma 50 GeV");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ // trigger
+ gener->SetTrigger(2);
+ gener->SetPtJet(50.);
+ gener->SetRadiation(isw);
+ gener->SetSimpleJets(!isw);
+ gener->SetJetEtaRange(-0.12, 0.12);
+ gener->SetJetPhiRange(220., 320.);
+ gGener=gener;
+ }
+ break;
+
+ case kHijing_gj75:
+ {
+ comment = comment.Append("HIJING Gamma 75 GeV");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ // trigger
+ gener->SetTrigger(2);
+ gener->SetPtJet(75.);
+ gener->SetRadiation(isw);
+ gener->SetSimpleJets(!isw);
+ gener->SetJetEtaRange(-0.12, 0.12);
+ gener->SetJetPhiRange(220., 320.);
+ gGener=gener;
+ }
+ break;
+
+ case kHijing_gj100:
+ {
+ comment = comment.Append("HIJING Gamma 100 GeV");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ // trigger
+ gener->SetTrigger(2);
+ gener->SetPtJet(100.);
+ gener->SetRadiation(isw);
+ gener->SetSimpleJets(!isw);
+ gener->SetJetEtaRange(-0.12, 0.12);
+ gener->SetJetPhiRange(220., 320.);
+ gGener=gener;
+ }
+ break;
+
+ case kHijing_gj200:
+ {
+ comment = comment.Append("HIJING Gamma 200 GeV");
+ AliGenHijing *gener = HijingStandard();
+// impact parameter range
+ gener->SetImpactParameterRange(0., 5.);
+ // trigger
+ gener->SetTrigger(2);
+ gener->SetPtJet(200.);
+ gener->SetRadiation(isw);
+ gener->SetSimpleJets(!isw);
+ gener->SetJetEtaRange(-0.12, 0.12);
+ gener->SetJetPhiRange(220., 320.);
+ gGener=gener;
+ }
+ break;
+ case kHijing_pA:
+ {
+ comment = comment.Append("HIJING pA");
+
+ AliGenCocktail *gener = new AliGenCocktail();
+
+ AliGenHijing *hijing = new AliGenHijing(-1);
+// centre of mass energy
+ hijing->SetEnergyCMS(TMath::Sqrt(82./208.) * 14000.);
+// impact parameter range
+ hijing->SetImpactParameterRange(0., 15.);
+// reference frame
+ hijing->SetReferenceFrame("CMS");
+ hijing->SetBoostLHC(1);
+// projectile
+ hijing->SetProjectile("P", 1, 1);
+ hijing->SetTarget ("A", 208, 82);
+// tell hijing to keep the full parent child chain
+ hijing->KeepFullEvent();
+// enable jet quenching
+ hijing->SetJetQuenching(0);
+// enable shadowing
+ hijing->SetShadowing(1);
+// Don't track spectators
+ hijing->SetSpectators(0);
+// kinematic selection
+ hijing->SetSelectAll(0);
+//
+ AliGenSlowNucleons* gray = new AliGenSlowNucleons(1);
+ AliSlowNucleonModel* model = new AliSlowNucleonModelExp();
+ gray->SetSlowNucleonModel(model);
+ gray->SetDebug(1);
+ gener->AddGenerator(hijing,"Hijing pPb", 1);
+ gener->AddGenerator(gray, "Gray Particles",1);
+ gGener=gener;
+ }
+ break;
+ case kPythia6:
+ {
+ comment = comment.Append(":Pythia p-p @ 14 TeV");
+ AliGenPythia *gener = new AliGenPythia(-1);
+ gener->SetMomentumRange(0,999999);
+ gener->SetThetaRange(0., 180.);
+ gener->SetYRange(-12,12);
+ gener->SetPtRange(0,1000);
+ gener->SetProcess(kPyMb);
+ gener->SetEnergyCMS(14000.);
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets20_24:
+ {
+ comment = comment.Append(":Pythia jets 20-24 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(20., 24.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets24_29:
+ {
+ comment = comment.Append(":Pythia jets 24-29 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(24., 29.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets29_35:
+ {
+ comment = comment.Append(":Pythia jets 29-35 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(29., 35.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets35_42:
+ {
+ comment = comment.Append(":Pythia jets 35-42 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(35., 42.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets42_50:
+ {
+ comment = comment.Append(":Pythia jets 42-50 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(42., 50.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets50_60:
+ {
+ comment = comment.Append(":Pythia jets 50-60 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(50., 60.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets60_72:
+ {
+ comment = comment.Append(":Pythia jets 60-72 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(60., 72.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets72_86:
+ {
+ comment = comment.Append(":Pythia jets 72-86 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(72., 86.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets86_104:
+ {
+ comment = comment.Append(":Pythia jets 86-104 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(86., 104.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets104_125:
+ {
+ comment = comment.Append(":Pythia jets 105-125 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(104., 125.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets125_150:
+ {
+ comment = comment.Append(":Pythia jets 125-150 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(125., 150.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gGener=gener;
+ }
+ break;
+ case kPythia6Jets150_180:
+ {
+ comment = comment.Append(":Pythia jets 150-180 GeV @ 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);// Centre of mass energy
+ gener->SetProcess(kPyJets);// Process type
+ gener->SetJetEtaRange(-0.5, 0.5);// Final state kinematic cuts
+ gener->SetJetPhiRange(0., 360.);
+ gener->SetJetEtRange(10., 1000.);
+ gener->SetGluonRadiation(1,1);
+ // gener->SetPtKick(0.);
+ // Structure function
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(150., 180.);// Pt transfer of the hard scattering
+ gener->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
+ gener->SetForceDecay(kAll);// Decay type (semielectronic, etc.)
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kD0PbPb5500:
+ {
+ comment = comment.Append(" D0 in Pb-Pb at 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(10);
+ gener->SetProcess(kPyD0PbPbMNR);
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(2.1,-1.0);
+ gener->SetEnergyCMS(5500.);
+ gener->SetNuclei(208,208);
+ gener->SetForceDecay(kHadronicD);
+ gener->SetYRange(-2,2);
+ gener->SetFeedDownHigherFamily(kFALSE);
+ gener->SetStackFillOpt(AliGenPythia::kParentSelection);
+ gener->SetCountMode(AliGenPythia::kCountParents);
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kCharmSemiElPbPb5500:
+ {
+ comment = comment.Append(" Charm in Pb-Pb at 5.5 TeV");
+ AliGenPythia * gener = new AliGenPythia(10);
+ gener->SetProcess(kPyCharmPbPbMNR);
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(2.1,-1.0);
+ gener->SetEnergyCMS(5500.);
+ gener->SetNuclei(208,208);
+ gener->SetForceDecay(kSemiElectronic);
+ gener->SetYRange(-2,2);
+ gener->SetFeedDownHigherFamily(kFALSE);
+ gener->SetCountMode(AliGenPythia::kCountParents);
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kBeautySemiElPbPb5500:
+ {
+ comment = comment.Append(" Beauty in Pb-Pb at 5.5 TeV");
+ AliGenPythia *gener = new AliGenPythia(10);
+ gener->SetProcess(kPyBeautyPbPbMNR);
+ gener->SetStrucFunc(kCTEQ4L);
+ gener->SetPtHard(2.75,-1.0);
+ gener->SetEnergyCMS(5500.);
+ gener->SetNuclei(208,208);
+ gener->SetForceDecay(kSemiElectronic);
+ gener->SetYRange(-2,2);
+ gener->SetFeedDownHigherFamily(kFALSE);
+ gener->SetCountMode(AliGenPythia::kCountParents);
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kCocktailTRD:
+ {
+ comment = comment.Append(" Cocktail for TRD at 5.5 TeV");
+ AliGenCocktail *gener = new AliGenCocktail();
+
+ AliGenParam *phi = new AliGenParam(10,
+ new AliGenMUONlib(),
+ AliGenMUONlib::kPhi,
+ "Vogt PbPb");
+
+ phi->SetPtRange(0, 100);
+ phi->SetYRange(-1., +1.);
+ phi->SetForceDecay(kDiElectron);
+
+ AliGenParam *omega = new AliGenParam(10,
+ new AliGenMUONlib(),
+ AliGenMUONlib::kOmega,
+ "Vogt PbPb");
+
+ omega->SetPtRange(0, 100);
+ omega->SetYRange(-1., +1.);
+ omega->SetForceDecay(kDiElectron);
+
+ AliGenParam *jpsi = new AliGenParam(10,
+ new AliGenMUONlib(),
+ AliGenMUONlib::kJpsiFamily,
+ "Vogt PbPb");
+
+ jpsi->SetPtRange(0, 100);
+ jpsi->SetYRange(-1., +1.);
+ jpsi->SetForceDecay(kDiElectron);
+
+ AliGenParam *ups = new AliGenParam(10,
+ new AliGenMUONlib(),
+ AliGenMUONlib::kUpsilonFamily,
+ "Vogt PbPb");
+ ups->SetPtRange(0, 100);
+ ups->SetYRange(-1., +1.);
+ ups->SetForceDecay(kDiElectron);
+
+ AliGenParam *charm = new AliGenParam(10,
+ new AliGenMUONlib(),
+ AliGenMUONlib::kCharm,
+ "central");
+ charm->SetPtRange(0, 100);
+ charm->SetYRange(-1.5, +1.5);
+ charm->SetForceDecay(kSemiElectronic);
+
+
+ AliGenParam *beauty = new AliGenParam(10,
+ new AliGenMUONlib(),
+ AliGenMUONlib::kBeauty,
+ "central");
+ beauty->SetPtRange(0, 100);
+ beauty->SetYRange(-1.5, +1.5);
+ beauty->SetForceDecay(kSemiElectronic);
+
+ AliGenParam *beautyJ = new AliGenParam(10,
+ new AliGenMUONlib(),
+ AliGenMUONlib::kBeauty,
+ "central");
+ beautyJ->SetPtRange(0, 100);
+ beautyJ->SetYRange(-1.5, +1.5);
+ beautyJ->SetForceDecay(kBJpsiDiElectron);
+
+ gener->AddGenerator(phi,"Phi",1);
+ gener->AddGenerator(omega,"Omega",1);
+ gener->AddGenerator(jpsi,"J/psi",1);
+ gener->AddGenerator(ups,"Upsilon",1);
+ gener->AddGenerator(charm,"Charm",1);
+ gener->AddGenerator(beauty,"Beauty",1);
+ gener->AddGenerator(beautyJ,"J/Psi from Beauty",1);
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kPyJJ:
+ {
+ comment = comment.Append(" Jet-jet at 5.5 TeV");
+ AliGenPythia *gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);
+ gener->SetProcess(kPyJets);
+ Double_t ptHardMin=10.0, ptHardMax=-1.0;
+ gener->SetPtHard(ptHardMin,ptHardMax);
+ gener->SetYHard(-0.7,0.7);
+ gener->SetJetEtaRange(-0.2,0.2);
+ gener->SetEventListRange(0,1);
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kPyGJ:
+ {
+ comment = comment.Append(" Gamma-jet at 5.5 TeV");
+ AliGenPythia *gener = new AliGenPythia(-1);
+ gener->SetEnergyCMS(5500.);
+ gener->SetProcess(kPyDirectGamma);
+ Double_t ptHardMin=10.0, ptHardMax=-1.0;
+ gener->SetPtHard(ptHardMin,ptHardMax);
+ gener->SetYHard(-1.0,1.0);
+ gener->SetGammaEtaRange(-0.13,0.13);
+ gener->SetGammaPhiRange(210.,330.);
+ gener->SetEventListRange(0,1);
+ gener->SetProjectile("p", 1, 1) ;
+ gener->SetTarget("p", 1, 1) ;
+ gGener=gener;
+ }
+ break;
+ case kMuonCocktailCent1:
+ {
+ comment = comment.Append(" Muon Cocktail Cent1");
+ AliGenMUONCocktail * gener = new AliGenMUONCocktail();
+ gener->SetPtRange(0.4,100.); // Transverse momentum range
+ gener->SetPhiRange(0.,360.); // Azimuthal angle range
+ gener->SetYRange(-4.0,-2.4);
+ gener->SetMuonPtCut(0.8);
+ gener->SetMuonThetaCut(171.,178.);
+ gener->SetMuonMultiplicity(2);
+ gener->SetImpactParameterRange(0.,5.); //Centrality class Cent1 for PDC04
+ gGener=gener;
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ }
+ break;
+ case kMuonCocktailPer1:
+ {
+ comment = comment.Append(" Muon Cocktail Per1");
+ AliGenMUONCocktail * gener = new AliGenMUONCocktail();
+ gener->SetPtRange(0.0,100.); // Transverse momentum range
+ gener->SetPhiRange(0.,360.); // Azimuthal angle range
+ gener->SetYRange(-4.0,-2.4);
+ gener->SetMuonPtCut(0.8);
+ gener->SetMuonThetaCut(171.,178.);
+ gener->SetMuonMultiplicity(2);
+ gener->SetImpactParameterRange(5.,8.6);//Centrality class Per1 for PDC04
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kMuonCocktailPer4:
+ {
+ comment = comment.Append(" Muon Cocktail Per4");
+ AliGenMUONCocktail * gener = new AliGenMUONCocktail();
+ gener->SetPtRange(0.0,100.); // Transverse momentum range
+ gener->SetPhiRange(0.,360.); // Azimuthal angle range
+ gener->SetYRange(-4.0,-2.4);
+ gener->SetMuonPtCut(0.8);
+ gener->SetMuonThetaCut(171.,178.);
+ gener->SetMuonMultiplicity(2);
+ gener->SetImpactParameterRange(13.2,15.0);//Centrality class Per4 for PDC04
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kMuonCocktailCent1HighPt:
+ {
+ comment = comment.Append(" Muon Cocktail HighPt Cent1");
+ AliGenMUONCocktail * gener = new AliGenMUONCocktail();
+ gener->SetPtRange(0.0,100.); // Transverse momentum range
+ gener->SetPhiRange(0.,360.); // Azimuthal angle range
+ gener->SetYRange(-4.0,-2.4);
+ gener->SetMuonPtCut(2.5);
+ gener->SetMuonThetaCut(171.,178.);
+ gener->SetMuonMultiplicity(2);
+ gener->SetImpactParameterRange(0.,5.); //Centrality class Cent1 for PDC04
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kMuonCocktailPer1HighPt :
+ {
+ comment = comment.Append(" Muon Cocktail HighPt Per1");
+ AliGenMUONCocktail * gener = new AliGenMUONCocktail();
+ gener->SetPtRange(0.0,100.); // Transverse momentum range
+ gener->SetPhiRange(0.,360.); // Azimuthal angle range
+ gener->SetYRange(-4.0,-2.4);
+ gener->SetMuonPtCut(2.5);
+ gener->SetMuonThetaCut(171.,178.);
+ gener->SetMuonMultiplicity(2);
+ gener->SetImpactParameterRange(5.,8.6);//Centrality class Per1 for PDC04
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kMuonCocktailPer4HighPt:
+ {
+ comment = comment.Append(" Muon Cocktail HighPt Per4");
+ AliGenMUONCocktail * gener = new AliGenMUONCocktail();
+ gener->SetPtRange(0.0,100.); // Transverse momentum range
+ gener->SetPhiRange(0.,360.); // Azimuthal angle range
+ gener->SetYRange(-4.0,-2.4);
+ gener->SetMuonPtCut(2.5);
+ gener->SetMuonThetaCut(171.,178.);
+ gener->SetMuonMultiplicity(2);
+ gener->SetImpactParameterRange(13.2,15.0);//Centrality class Per4 for PDC04
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kMuonCocktailCent1Single:
+ {
+ comment = comment.Append(" Muon Cocktail Single Cent1");
+ AliGenMUONCocktail * gener = new AliGenMUONCocktail();
+ gener->SetPtRange(0.0,100.); // Transverse momentum range
+ gener->SetPhiRange(0.,360.); // Azimuthal angle range
+ gener->SetYRange(-4.0,-2.4);
+ gener->SetMuonPtCut(0.8);
+ gener->SetMuonThetaCut(171.,178.);
+ gener->SetMuonMultiplicity(1);
+ gener->SetImpactParameterRange(0.,5.); //Centrality class Cent1 for PDC04
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kMuonCocktailPer1Single :
+ {
+ comment = comment.Append(" Muon Cocktail Single Per1");
+ AliGenMUONCocktail * gener = new AliGenMUONCocktail();
+ gener->SetPtRange(0.0,100.); // Transverse momentum range
+ gener->SetPhiRange(0.,360.); // Azimuthal angle range
+ gener->SetYRange(-4.0,-2.4);
+ gener->SetMuonPtCut(0.8);
+ gener->SetMuonThetaCut(171.,178.);
+ gener->SetMuonMultiplicity(1);
+ gener->SetImpactParameterRange(5.,8.6);//Centrality class Per1 for PDC04
+ gener->SetNumberOfParticipants(229.3);//Centrality class Per1 for PDC04
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kMuonCocktailPer4Single:
+ {
+ comment = comment.Append(" Muon Cocktail Single Per4");
+ AliGenMUONCocktail * gener = new AliGenMUONCocktail();
+ gener->SetPtRange(0.0,100.); // Transverse momentum range
+ gener->SetPhiRange(0.,360.); // Azimuthal angle range
+ gener->SetYRange(-4.0,-2.4);
+ gener->SetMuonPtCut(0.8);
+ gener->SetMuonThetaCut(171.,178.);
+ gener->SetMuonMultiplicity(1);
+ gener->SetImpactParameterRange(13.2,15.0);//Centrality class Per4 for PDC04
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener=gener;
+ }
+ break;
+ case kFlow_2_2000:
+ {
+ comment = comment.Append(" Flow with dN/deta = 2000, vn = 2%");
+ gGener = GeVSimStandard(2000., 2.);
+ }
+ break;
+
+ case kFlow_10_2000:
+ {
+ comment = comment.Append(" Flow with dN/deta = 2000, vn = 10%");
+ gGener = GeVSimStandard(2000., 10.);
+ }
+ break;
+
+ case kFlow_6_2000:
+ {
+ comment = comment.Append(" Flow with dN/deta = 2000, vn = 6%");
+ gGener = GeVSimStandard(2000., 6.);
+ }
+ break;
+
+ case kFlow_6_5000:
+ {
+ comment = comment.Append(" Flow with dN/deta = 5000, vn = 6%");
+ gGener = GeVSimStandard(5000., 6.);
+ }
+ break;
+ case kHIJINGplus:
+ {
+ //
+ // The cocktail
+ AliGenCocktail *gener = new AliGenCocktail();
+
+ //
+ // Charm production by Pythia
+ AliGenPythia * genpyc = new AliGenPythia(230);
+ genpyc->SetProcess(kPyCharmPbPbMNR);
+ genpyc->SetStrucFunc(kCTEQ4L);
+ genpyc->SetPtHard(2.1,-1.0);
+ genpyc->SetEnergyCMS(5500.);
+ genpyc->SetNuclei(208,208);
+ genpyc->SetYRange(-999,999);
+ genpyc->SetForceDecay(kAll);
+ genpyc->SetFeedDownHigherFamily(kFALSE);
+ genpyc->SetCountMode(AliGenPythia::kCountParents);
+ //
+ // Beauty production by Pythia
+ AliGenPythia * genpyb = new AliGenPythia(9);
+ genpyb->SetProcess(kPyBeautyPbPbMNR);
+ genpyb->SetStrucFunc(kCTEQ4L);
+ genpyb->SetPtHard(2.75,-1.0);
+ genpyb->SetEnergyCMS(5500.);
+ genpyb->SetNuclei(208,208);
+ genpyb->SetYRange(-999,999);
+ genpyb->SetForceDecay(kAll);
+ genpyb->SetFeedDownHigherFamily(kFALSE);
+ genpyb->SetCountMode(AliGenPythia::kCountParents);
+ //
+ // Hyperons
+ //
+ AliGenSTRANGElib *lib = new AliGenSTRANGElib();
+ Int_t particle;
+ // Xi
+ particle = kXiMinus;
+ AliGenParam *genXi = new AliGenParam(16,lib,particle);
+ genXi->SetPtRange(0., 12.);
+ genXi->SetYRange(-1.1, 1.1);
+ genXi->SetForceDecay(kNoDecay);
+
+ //
+ // Omega
+ particle = kOmegaMinus;
+ AliGenParam *genOmega = new AliGenParam(10,lib,particle);
+ genOmega->SetPtRange(0, 12.);
+ genOmega->SetYRange(-1.1, 1.1);
+ genOmega->SetForceDecay(kNoDecay);
+
+ //
+ // Central Hijing
+ AliGenHijing *genHi = HijingStandard();
+ genHi->SwitchOffHeavyQuarks(kTRUE);
+ genHi->SetImpactParameterRange(0.,5.);
+ //
+ // Add everything to the cocktail and shake ...
+ gener->AddGenerator(genHi, "Hijing cent1", 1);
+ gener->AddGenerator(genpyc, "Extra charm", 1);
+ gener->AddGenerator(genpyb, "Extra beauty", 1);
+ gener->AddGenerator(genXi, "Xi" , 1);
+ gener->AddGenerator(genOmega, "Omega", 1);
+ gener->SetProjectile("A", 208, 82) ;
+ gener->SetTarget("A", 208, 82) ;
+ gGener = gener;
+ }
+ break;
+ default: break;
+ }
+
+ return gGener;
+}
+
+AliGenHijing* HijingStandard()
+{
+ AliGenHijing *gener = new AliGenHijing(-1);
+ // centre of mass energy
+ gener->SetEnergyCMS(5500.);
+ // reference frame
+ gener->SetReferenceFrame("CMS");
+ // projectile
+ gener->SetProjectile("A", 208, 82);
+ gener->SetTarget ("A", 208, 82);
+ // tell hijing to keep the full parent child chain
+ gener->KeepFullEvent();
+ // enable jet quenching
+ gener->SetJetQuenching(1);
+ // enable shadowing
+ gener->SetShadowing(1);
+ // neutral pion and heavy particle decays switched off
+ gener->SetDecaysOff(1);
+ // Don't track spectators
+ gener->SetSpectators(0);
+ // kinematic selection
+ gener->SetSelectAll(0);
+ return gener;
+}
+
+AliGenGeVSim* GeVSimStandard(Float_t mult, Float_t vn)
+{
+ AliGenGeVSim* gener = new AliGenGeVSim(0);
+ //
+ // Mult is the number of charged particles in |eta| < 0.5
+ // Vn is in (%)
+ //
+ // Sigma of the Gaussian dN/deta
+ Float_t sigma_eta = 2.75;
+ //
+ // Maximum eta
+ Float_t etamax = 7.00;
+ //
+ //
+ // Scale from multiplicity in |eta| < 0.5 to |eta| < |etamax|
+ Float_t mm = mult * (TMath::Erf(etamax/sigma_eta/sqrt(2.)) / TMath::Erf(0.5/sigma_eta/sqrt(2.)));
+ //
+ // Scale from charged to total multiplicity
+ //
+ mm *= 1.587;
+ //
+ // Vn
+ vn /= 100.;
+ //
+ // Define particles
+ //
+ //
+ // 78% Pions (26% pi+, 26% pi-, 26% p0) T = 250 MeV
+ AliGeVSimParticle *pp = new AliGeVSimParticle(kPiPlus, 1, 0.26 * mm, 0.25, sigma_eta) ;
+ AliGeVSimParticle *pm = new AliGeVSimParticle(kPiMinus, 1, 0.26 * mm, 0.25, sigma_eta) ;
+ AliGeVSimParticle *p0 = new AliGeVSimParticle(kPi0, 1, 0.26 * mm, 0.25, sigma_eta) ;
+ //
+ // 12% Kaons (3% K0short, 3% K0long, 3% K+, 3% K-) T = 300 MeV
+ AliGeVSimParticle *ks = new AliGeVSimParticle(kK0Short, 1, 0.03 * mm, 0.30, sigma_eta) ;
+ AliGeVSimParticle *kl = new AliGeVSimParticle(kK0Long, 1, 0.03 * mm, 0.30, sigma_eta) ;
+ AliGeVSimParticle *kp = new AliGeVSimParticle(kKPlus, 1, 0.03 * mm, 0.30, sigma_eta) ;
+ AliGeVSimParticle *km = new AliGeVSimParticle(kKMinus, 1, 0.03 * mm, 0.30, sigma_eta) ;
+ //
+ // 10% Protons / Neutrons (5% Protons, 5% Neutrons) T = 250 MeV
+ AliGeVSimParticle *pr = new AliGeVSimParticle(kProton, 1, 0.05 * mm, 0.25, sigma_eta) ;
+ AliGeVSimParticle *ne = new AliGeVSimParticle(kNeutron, 1, 0.05 * mm, 0.25, sigma_eta) ;
+ //
+ // Set Elliptic Flow properties
+
+ Float_t pTsaturation = 2. ;
+
+ pp->SetEllipticParam(vn,pTsaturation,0.) ;
+ pm->SetEllipticParam(vn,pTsaturation,0.) ;
+ p0->SetEllipticParam(vn,pTsaturation,0.) ;
+ pr->SetEllipticParam(vn,pTsaturation,0.) ;
+ ne->SetEllipticParam(vn,pTsaturation,0.) ;
+ ks->SetEllipticParam(vn,pTsaturation,0.) ;
+ kl->SetEllipticParam(vn,pTsaturation,0.) ;
+ kp->SetEllipticParam(vn,pTsaturation,0.) ;
+ km->SetEllipticParam(vn,pTsaturation,0.) ;
+ //
+ // Set Direct Flow properties
+ pp->SetDirectedParam(vn,1.0,0.) ;
+ pm->SetDirectedParam(vn,1.0,0.) ;
+ p0->SetDirectedParam(vn,1.0,0.) ;
+ pr->SetDirectedParam(vn,1.0,0.) ;
+ ne->SetDirectedParam(vn,1.0,0.) ;
+ ks->SetDirectedParam(vn,1.0,0.) ;
+ kl->SetDirectedParam(vn,1.0,0.) ;
+ kp->SetDirectedParam(vn,1.0,0.) ;
+ km->SetDirectedParam(vn,1.0,0.) ;
+ //
+ // Add particles to the list
+ gener->AddParticleType(pp) ;
+ gener->AddParticleType(pm) ;
+ gener->AddParticleType(p0) ;
+ gener->AddParticleType(pr) ;
+ gener->AddParticleType(ne) ;
+ gener->AddParticleType(ks) ;
+ gener->AddParticleType(kl) ;
+ gener->AddParticleType(kp) ;
+ gener->AddParticleType(km) ;
+ //
+ // Random Ev.Plane ----------------------------------
+ TF1 *rpa = new TF1("gevsimPsiRndm","1", 0, 360);
+ // --------------------------------------------------
+ gener->SetPtRange(0., 9.) ; // Use a resonable range! (used for bin size in numerical integration)
+ gener->SetPhiRange(0, 360);
+ //
+ // Set pseudorapidity range
+ Float_t thmin = EtaToTheta(+etamax);
+ Float_t thmax = EtaToTheta(-etamax);
+ gener->SetThetaRange(thmin,thmax);
+ return gener;
+}
+
+
+
+void ProcessEnvironmentVars()
+{
+ // Run type
+ if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
+ for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
+ if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) {
+ srun = (PprRun_t)iRun;
+ cout<<"Run type set to "<<pprRunName[iRun]<<endl;
+ }
+ }
+ }
+
+ // Random Number seed
+ if (gSystem->Getenv("CONFIG_SEED")) {
+ sseed = atoi(gSystem->Getenv("CONFIG_SEED"));
+ }
+}
--- /dev/null
+void DrawSysWatchTime(){
+ TString prefix="/hera/alice/marsland/MAF/MAFbenchmark/mcmaker/workdir/test_6_375000_25_20140713_20/";
+ TString listSyswatch = gSystem->GetFromPipe(Form("ls %s/*/mult_10000/event_1/simwatch.log",prefix.Data()));
+ TObjArray * arraySyswatch= listSyswatch.Tokenize("\n");
+ Int_t nfiles= arraySyswatch->GetEntries();
+ TTree * treeSyswatch[] = new TTree*[nfiles];
+
+
+ for (Int_t ifile=0; ifile<nfiles; ifile++){
+ treeSyswatch[ifile] = AliSysInfo::MakeTree(arraySyswatch->At(ifile)->GetName());
+ treeSyswatch[0]->AddFriend( treeSyswatch[ifile],Form("T%d",ifile));
+ }
+
+ treeSyswatch[0]->Draw("deltaT:sname","deltaT>1","");
+ for (Int_t ifile=1; ifile<nfiles; ifile++){
+ treeSyswatch[0]->SetMarkerStyle(25);
+ treeSyswatch[0]->SetMarkerColor(1+ifile);
+ treeSyswatch[0]->Draw(Form("T%d.stampSec-T%d.stampOldSec:sname",ifile,ifile),"deltaT>1","same");
+ }
+}
+
+
--- /dev/null
+const char * recoStorage="local:///cvmfs/alice.gsi.de/alice/data/2010/OCDB";
+Int_t run=0;
+
+
+void ModifyRecoParam(TObjArray* recoArray, Bool_t useIonTail, Double_t crossTalkCorrection){
+ //
+ // Modify reco param - and store it in the OCDB in local directory
+ //
+ AliCDBManager * man = AliCDBManager::Instance();
+ for (Int_t i=0; i<4; i++){
+ AliTPCRecoParam* p = ( AliTPCRecoParam*)recoArray->At(i);
+ p->SetUseIonTailCorrection(useIonTail);
+ p->SetCrosstalkCorrection(crossTalkCorrection);
+ }
+ TString localStorage = "local://"+gSystem->GetFromPipe("pwd")+"/OCDBrec";
+ AliCDBStorage*pocdbStorage = AliCDBManager::Instance()->GetStorage(localStorage.Data());
+ AliCDBMetaData *metaData= new AliCDBMetaData();
+ AliCDBId* id1=new AliCDBId("TPC/Calib/RecoParam/", man->GetRun(), man->GetRun());
+ pocdbStorage->Put(recoArray, (*id1), metaData);
+ AliCDBManager::Instance()->SetSpecificStorage("TPC/Calib/RecoParam/",localStorage.Data());
+}
+
+
+void rec(Bool_t useIonTail, Double_t crossTalkCorrection) {
+ //
+ // run reconstruction
+ // Parameters:
+ // useIonTail - switch for using ion tail correction - OCDB entry will be overritten in
+ //
+ // stard reco setting
+ //
+ AliCDBManager * man = AliCDBManager::Instance();
+ man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ man->SetSpecificStorage("TPC/Calib/RecoParam/",recoStorage);
+ man->SetRun(run);
+ AliReconstruction reco;
+ reco.SetWriteESDfriend();
+ reco.SetWriteAlignmentData();
+ reco.SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ reco.SetSpecificStorage("GRP/GRP/Data",
+ Form("local://%s",gSystem->pwd()));
+ reco.SetRunPlaneEff(kTRUE);
+ reco.SetRunQA(":");
+ reco.SetRunGlobalQA(kFALSE);
+ reco.ResetCheckRecoCDBvsSimuCDB();
+
+ //
+ //Switch Iontail in RecoParam. Creation of a new OCDB entry
+ //
+ AliCDBEntry* e = man->Get("TPC/Calib/RecoParam/",run); // get default
+ // modify content
+ TObjArray* recoArray = (TObjArray*)e->GetObject();
+ ModifyRecoParam(recoArray, useIonTail, crossTalkCorrection);
+ //
+ //
+ // Run reconstruction
+ //
+ TStopwatch timer;
+ timer.Start();
+ reco.Run();
+ timer.Stop();
+ timer.Print();
+ //
+ // Print the OCDB setup which we used
+ //
+ AliCDBEntry* ocdbEntry = man->Get("TPC/Calib/RecoParam/",run);
+ TObjArray* recoArray = (TObjArray*)ocdbEntry->GetObject();
+ for (Int_t i=0; i<4; i++){
+ AliTPCRecoParam* recoParam = ( AliTPCRecoParam*)recoArray->At(i);
+ printf("ipar=%d\t%d\t%f\n",i,recoParam->GetUseIonTailCorrection(), recoParam->GetCrosstalkCorrection());
+ }
+}
--- /dev/null
+//
+// parameter to take from config file
+//
+const char * recoStorage="local:///cvmfs/alice.gsi.de/alice/data/2010/OCDB";
+Int_t run=0;
+
+
+void ModifyRecoParam(TObjArray* recoArray, Bool_t useIonTail, Double_t crossTalkCorrection){
+ //
+ // Modify reco param - and store it in the OCDB in local directory
+ //
+ AliCDBManager * man = AliCDBManager::Instance();
+ for (Int_t i=0; i<4; i++){
+ AliTPCRecoParam* p = ( AliTPCRecoParam*)recoArray->At(i);
+ p->SetUseIonTailCorrection(useIonTail);
+ p->SetCrosstalkCorrection(crossTalkCorrection);
+ }
+ TString localStorage = "local://"+gSystem->GetFromPipe("pwd")+"/OCDBsim";
+ AliCDBStorage*pocdbStorage = AliCDBManager::Instance()->GetStorage(localStorage.Data());
+ AliCDBMetaData *metaData= new AliCDBMetaData();
+ AliCDBId* id1=new AliCDBId("TPC/Calib/RecoParam/", man->GetRun(), AliCDBRunRange::Infinity());
+ pocdbStorage->Put(recoArray, (*id1), metaData);
+ AliCDBManager::Instance()->SetSpecificStorage("TPC/Calib/RecoParam/",localStorage.Data());
+}
+
+
+void sim(Int_t nev, Bool_t useIonTail, Double_t crossTalkCorrection) {
+ gSystem->Load("liblhapdf");
+ gSystem->Load("libEGPythia6");
+ gSystem->Load("libpythia6");
+ gSystem->Load("libAliPythia6");
+ gSystem->Load("libhijing");
+ gSystem->Load("libTHijing");
+ gSystem->Load("libgeant321");
+
+ if (nev<0){
+ AliCDBManager * man = AliCDBManager::Instance();
+ man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ man->SetSpecificStorage("TPC/Calib/RecoParam/",recoStorage);
+ man->SetRun(run);
+ AliCDBEntry* e = man->Get("TPC/Calib/RecoParam/",run); // get default
+ // modify content
+ TObjArray* recoArray = (TObjArray*)e->GetObject();
+ ModifyRecoParam(recoArray, useIonTail, crossTalkCorrection);
+ return;
+ }
+
+ if (gSystem->Getenv("EVENT")) nev = atoi(gSystem->Getenv("EVENT")) ;
+
+ AliSimulation simulator;
+ simulator.SetMakeSDigits("ITS TPC TRD TOF PHOS HMPID EMCAL MUON FMD ZDC PMD T0");
+ //simulator.SetMakeDigitsFromHits( "ITS TPC");
+ simulator.SetWriteRawData("ALL","raw.root",kTRUE);
+
+ simulator.SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ simulator.SetSpecificStorage("GRP/GRP/Data", Form("local://%s",gSystem->pwd()));
+ TString localStorage = "local://"+gSystem->GetFromPipe("pwd")+"/OCDBsim";
+ simulator.SetSpecificStorage("TPC/Calib/RecoParam/",localStorage.Data());
+
+ simulator.SetRunQA(":") ;
+
+ TStopwatch timer;
+ timer.Start();
+ simulator.Run(nev);
+ timer.Stop();
+ timer.Print();
+ //
+ // Print the OCDB setup which we used
+ //
+ AliCDBManager * man = AliCDBManager::Instance();
+ AliCDBEntry* ocdbEntry = man->Get("TPC/Calib/RecoParam/",run);
+ TObjArray* recoArray = (TObjArray*)ocdbEntry->GetObject();
+ for (Int_t i=0; i<4; i++){
+ AliTPCRecoParam* recoParam = ( AliTPCRecoParam*)recoArray->At(i);
+ printf("ipar=%d\t%d\t%f\n",i,recoParam->GetUseIonTailCorrection(), recoParam->GetCrosstalkCorrection());
+ }
+
+
+
+}
--- /dev/null
+#!/bin/sh
+
+###############################################################################
+#
+# Aim is to submit simulation/reconstruction jobs for multiplicity dependence check
+# For local running one can use the instructions in runSim() function.
+# submitMultiplicityScan() --> main function to submit jobs for each event with a given multiplicity
+# runSim() --> runs the macros for simulation and reconstruction
+# To see how to run these functions look below
+#
+###############################################################################
+
+# global variables to be set
+DATE=`date +%Y_%m_%d_%H`
+nEvents=1 #fixed to 1 to make an event by event job submission
+
+###############################################################################
+submitIonTailXTalkScan()
+{
+
+ #
+ # cals the submitMultiplicityScan fuction for
+ # --> IonTail ON/OFF, XTALK ON/OFF, in reconstruction and MC
+ ###############################################################################
+ ## example --> for 5 multiplicity bins, each most central event having 15000 track multiplicity.
+ ## For each multiplicity bin total statistic is 75000
+ if [ 1 -eq 0 ]; then
+ cd $ALICE_ROOT/test/testdEdx
+ source submitSimJobs.sh
+ submitIonTailXTalkScan 8 1000000 50 /hera/alice/marsland/software/bin/set_private_TPCdev.sh > out.log 2>&1 &
+ fi
+ ###############################################################################
+
+ # inputs
+ nMultBins=$1
+ maxNTracks=$2
+ nEventsCentral=$3
+ alirootSourceFile=$4
+
+ if [ $# -ne 4 ]; then
+ echo "4 parameters needed --> multiplicity bins, maximum number of tracks, maximum number of central events, source file for root"
+ return 1
+ fi
+
+ baseDir=$(pwd)
+ testDir=$baseDir/test_$nMultBins\_$maxNTracks\_$nEventsCentral\_$DATE
+ DirCheckCreate $testDir
+ cd $testDir
+
+ for ((iontail=0; iontail<2 ; iontail=iontail+1))
+ do
+ for ((xtalk=0; xtalk<2 ; xtalk=xtalk+1))
+ do
+
+ MultiplicityScan $nMultBins $maxNTracks $nEventsCentral $alirootSourceFile $iontail $xtalk
+
+ done
+ done
+
+}
+###############################################################################
+runSim(){
+
+ #
+ # Function which runs rec.C and sim.C for given multiplicity and event number.
+ # (if submitMultiplicityScan function is called, this parameter is always fixed to 1, i.e event by event)
+ # Input parameters are
+ # 1) total track multiplicity
+ # 2) number of events to be processed for given total track multiplicity (if submitMultiplicityScan() is called, it is 1)
+ # 3) script to source a specific aliroot
+ #
+ ###############################################################################
+ ## example --> 2 events with total multiplicity of 1000 tracks
+ if [ 1 -eq 0 ]; then
+ cd $ALICE_ROOT/test/testdEdx
+ source submitSimJobs.sh
+ runSim 100 1 0 0 /hera/alice/marsland/software/bin/set_private_TPCdev.sh > out.log 2>&1 &
+ fi
+ ###############################################################################
+
+ export TestdEdxNTracks=$1
+ nevents=$2
+ ionTail=$3
+ xTalk=$4
+ alirootSourceFile=$5
+
+ if [ $# -ne 4 ]; then
+ echo "5 parameters needed --> multiplicity, number of events, source file for aliroot, iontail switch (0 or 1), xTalk switch (0 or 1)"
+ return 1
+ fi
+
+ echo " Running dEdx digitzer test job"
+ echo " NEvents = $nevents"
+ echo " NTracks per event $TestdEdxNTracks"
+
+
+ ## main body of the simulation part
+ rm -rf *.root *.dat *.log fort* hlt hough raw* recraw/*.root recraw/*.log GRP*
+ printf "\n ======================================================================\n\n"
+ echo Running: aliroot -b -q sim.C\($nevents,$ionTail,$xTalk\)
+ aliroot -b -q sim.C\(-1,$ionTail,$xTalk\) 2>&1 | tee sim.log # make a specific OCDB for simulation
+ aliroot -b -q sim.C\($nevents,$ionTail,$xTalk\) 2>&1 | tee sim.log
+ mv syswatch.log simwatch.log
+ printf "\n ======================================================================\n\n"
+ echo Running: aliroot -b -q rec.C\($ionTail\,$xTalk\)
+ aliroot -b -q rec.C\($ionTail\,$xTalk\) 2>&1 | tee rec.log
+ mv syswatch.log recwatch.log
+ ## OCDB entries to be dumped in human readable format
+ source $ALICE_ROOT/PWGPP/CalibMacros/AliOCDBtoolkit.sh
+ printf "\n ======================================================================\n\n"
+ echo Running: ocdbMakeTable AliESDs.root "ESD" OCDBrec.list
+ ocdbMakeTable AliESDs.root "ESD" OCDBrec.list
+ printf "\n ======================================================================\n\n"
+ echo Running: ocdbMakeTable galice.root MC OCDBsim.list
+ ocdbMakeTable galice.root MC OCDBsim.list
+ ocdbFileName=$(cat OCDBrec.list | grep "TPC/Calib/RecoParam" | gawk '{print $2"/"$3}' )
+ printf "\n ======================================================================\n\n"
+ echo Running: dumpObject $ocdbFileName "object" "XML" RecoParam
+ dumpObject $ocdbFileName "object" "XML" RecoParam
+
+ return 1;
+
+}
+###############################################################################
+MultiplicityScan(){
+ #
+ # Here we submit the jobs for the simulation//reconstruction for one setting of IonTail and XTalk configuration
+ # Parameters:
+ # 1. multiplicity bins to be investigated (default 5)
+ # 2. max multiplicity for whole processing (default 75000 tracks --> 5 central PbPb event )
+ # 3. number of central events to be used (default 5)
+ # 4. file to source aliroot
+ # (2)/(3) should be a reasonable multiplicity estimate (e.g. 15000 tracks which is 1 central PbPb event)
+ # Jobs will be submitted per event
+ #
+ # For each setting new directory will be created - indicating muiltiplicity
+ # dir<ntracks>/dir<eventNr>
+ #
+ ###############################################################################
+ ## example --> for 5 multiplicity bins, each most central event having 15000 track multiplicity.
+ ## For each multiplicity bin total statistic is 75000
+ if [ 1 -eq 0 ]; then
+ cd $ALICE_ROOT/test/testdEdx
+ source submitSimJobs.sh
+ MultiplicityScan 8 1000000 50 /hera/alice/marsland/software/bin/set_private_TPCdev.sh 0 1 > out.log 2>&1 &
+ fi
+ ###############################################################################
+
+ # inputs
+ nMultBins=$1
+ maxNTracks=$2
+ nEventsCentral=$3
+ alirootSourceFile=$4
+ ionTail=$5
+ xTalk=$6
+
+ if [ $# -ne 6 ]; then
+ echo "6 parameters needed --> multiplicity bins, maximum number of tracks, maximum number of central events, source file for root, iontail switch (0 or 1), xTalk switch (0 or 1)"
+ return 1
+ fi
+
+ baseDir=$(pwd)
+ testDir=$baseDir/IonTail_XTalk_$ionTail\_$xTalk
+ DirCheckCreate $testDir
+ cd $testDir
+
+ # create multiplicity bins
+ multPerCentralEvent=$(echo $maxNTracks/$nEventsCentral | bc)
+ echo "multiplicity per most central event is $multPerCentralEvent"
+ for ((i=0; i<$nMultBins ; i=i+1))
+ do
+
+ multSteps=$(echo $multPerCentralEvent/$nMultBins | bc)
+ multBin=$(echo $multPerCentralEvent - $multSteps*$i | bc)
+ multBinDir=$testDir/mult_$multBin
+ DirCheckCreate $multBinDir
+ cd $multBinDir
+ echo $multBinDir
+
+ nEventsPerMultBin=$(echo $maxNTracks/$multBin | bc)
+ echo $nEventsPerMultBin
+ for ((j=1; j<$(echo $nEventsPerMultBin+1 | bc) ; j=j+1))
+ do
+
+ eventDir=$multBinDir/event_$j
+ DirCheckCreate $eventDir
+ cd $eventDir
+ cp -r $ALICE_ROOT/test/testdEdx/GRP $ALICE_ROOT/test/testdEdx/OCDB* $ALICE_ROOT/test/testdEdx/*.* .
+ #cp $ALICE_ROOT/test/testdEdx/*.* .
+
+ qsub -V -cwd -l h_rt=24:0:0,h_rss=6G -P alice -b y -r y -o outSim.log -e errSim.log $baseDir/submitSimJobs.sh runSim $multBin $nEvents $4 $5 $6
+
+ done
+ done
+
+ cd $baseDir
+}
+###############################################################################
+DirCheckCreate()
+{
+
+ #
+ # check if the directory exist. If so, delete it
+ #
+
+ dirName=$1
+ if [ -d "$dirName" ]; then
+ echo " $dirName already exist delete it and create a new one "
+ rm -rf $dirName
+ fi
+ mkdir $dirName
+
+}
+###############################################################################
+main()
+{
+ eval "$@"
+}
+
+main $@
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff