--- /dev/null
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TList.h"
+#include "TCanvas.h"
+
+//#include <iostream.h>
+//#include <TVector3.h>
+//#include <TGeoHMatrix.h>
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+//#include "AliESDCaloCluster.h"
+//#include "AliESDCaloCells.h"
+//#include "AliESDHeader.h"
+#include "AliESDVZERO.h"
+#include "AliESDTZERO.h"
+#include "AliESDZDC.h"
+
+#include "AliAnalysisTaskPt.h"
+//#include "AliEMCALGeoUtils.h"
+
+// example of an analysis task creating a p_t spectrum
+// Authors: Panos Cristakoglou, Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing
+// Reviewed: A.Gheata (19/02/10)
+
+ClassImp(AliAnalysisTaskPt)
+
+//________________________________________________________________________
+AliAnalysisTaskPt::AliAnalysisTaskPt(const char *name)
+ : AliAnalysisTaskSE(name), fESD(0), fOutputList(0), fMyTr(0), jev(0), iev(0)
+{
+ // Constructor
+
+ // Define input and output slots here
+ // Input slot #0 works with a TChain
+ DefineInput(0, TChain::Class());
+ // Output slot #0 id reserved by the base class for AOD
+ // Output slot #1 writes into a TH1 container
+ DefineOutput(1, TList::Class());
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPt::UserCreateOutputObjects()
+{
+ // Create histograms
+ // Called once
+ // fEMCALGeo = new AliEMCALGeoUtils("EMCAL_FIRSTYEAR");
+
+ // float pi=acos(-1.0);
+ fOutputList = new TList();
+/*
+ fHstPt = new TH1F("fHstPt", "Pt distribution", 100, 0.0, 20.0);
+ fEMCe = new TH1F("fEMCe", "E distribution", 100, 0.0, 10.0);
+ fEMCt = new TH1F("fEMCt", "T distribution", 100, 0.0, 1000.0);
+ fEMCn = new TH1F("fEMCn", "N distribution", 100, -0.5, 99.5);
+ fEMCm = new TH2F("fEMCm", "phi-zps map", 50, 0, pi, 50, -300.0, 300.0);
+ fCelle = new TH1F("fCelle", "Cell E distribution", 100, 0.0, 10.0);
+ fCellf = new TH1F("fCellf", "Cell E fraction", 100, -0.1, 1.1);
+ fCellt = new TH1F("fCellt", "Cell T distribution", 100, 0.0, 1000.0);
+ fCellc = new TH2F("fCellc", "T-ID map", 5000, -200.0, 4800.0, 50, 0.0, 1000.0);
+ fCelld = new TH2F("fCelld", "E-ID map", 5000, -200.0, 4800.0, 50, 0.0, 10.0);
+ fCellm = new TH2F("fCellm", "phi-eta map", 50, 0, pi, 50, -300.0, 300.0);
+
+ fOutputList->Add(fHstPt);
+ fOutputList->Add(fEMCe);
+ fOutputList->Add(fEMCt);
+ fOutputList->Add(fEMCn);
+ fOutputList->Add(fEMCm);
+ fOutputList->Add(fCelle);
+ fOutputList->Add(fCellf);
+ fOutputList->Add(fCellt);
+ fOutputList->Add(fCellc);
+ fOutputList->Add(fCelld);
+ fOutputList->Add(fCellm);
+ iii=0;
+*/
+ jev=0;
+ iev=0;
+ fMyTr = new TTree("fMyTr","beam tree");
+ fOutputList->Add(fMyTr);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPt::UserExec(Option_t *)
+{
+ // Main loop
+ // Called for each event
+
+ fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+ if (!fESD) {
+ printf("ERROR: fESD not available\n");
+ return;
+ }
+
+//AliESDHeader *header = fESD->GetHeader();
+//fEMCALGeo->SetMisalMatrixes(header->GetEMCALMisalMatrix());
+//
+//for(Int_t mod=0; mod<(fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
+// if(fESD){
+// const TGeoHMatrix* m=fESD->GetEMCALMatrix(mod) ;
+// fEMCALGeo->SetMisalMatrix(m, mod) ;
+// }
+//}
+//
+
+ int ntr, etr;
+ float z1n, z2n, z1p, z2p, zvt;
+ float z1s, z2s, w1s, w2s, t1s, t2s, v1s, v2s;
+ float z1nt[10], z2nt[10], z1pt[10], z2pt[10];
+ float t0amp[24], t0tim[24], v0mul[64];
+ ntr = (int)(fESD->GetNumberOfTracks());
+ if (ntr>0 && iev==0) {
+ fMyTr->Branch("jev",&jev,"jev/I");
+ fMyTr->Branch("iev",&iev,"iev/I");
+ fMyTr->Branch("ntr",&ntr,"ntr/I");
+ fMyTr->Branch("etr",&etr,"etr/I");
+ fMyTr->Branch("z1n",&z1n,"z1n/F");
+ fMyTr->Branch("z2n",&z2n,"z2n/F");
+ fMyTr->Branch("z1p",&z1p,"z1p/F");
+ fMyTr->Branch("z2p",&z2p,"z2p/F");
+ fMyTr->Branch("zvt",&zvt,"zvt/F");
+ fMyTr->Branch("z1s",&z1s,"z1s/F");
+ fMyTr->Branch("z2s",&z2s,"z2s/F");
+ fMyTr->Branch("w1s",&w1s,"w1s/F");
+ fMyTr->Branch("w2s",&w2s,"w2s/F");
+ fMyTr->Branch("t1s",&t1s,"t1s/F");
+ fMyTr->Branch("t2s",&t2s,"t2s/F");
+ fMyTr->Branch("v1s",&v1s,"v1s/F");
+ fMyTr->Branch("v2s",&v2s,"v2s/F");
+ fMyTr->Branch("z1nt",z1nt,"z1nt[10]/F");
+ fMyTr->Branch("z2nt",z2nt,"z2nt[10]/F");
+ fMyTr->Branch("z1pt",z1pt,"z1pt[10]/F");
+ fMyTr->Branch("z2pt",z2pt,"z2pt[10]/F");
+ fMyTr->Branch("t0amp",t0amp,"t0amp[24]/F");
+ fMyTr->Branch("t0tim",t0tim,"t0tim[24]/F");
+ fMyTr->Branch("v0mul",v0mul,"v0mul[64]/F");
+ printf("my tree branches are created\n");
+ }
+ const AliESDZDC* zdcData = fESD->GetESDZDC();
+ const AliESDTZERO* tzrData = fESD->GetESDTZERO();
+ const AliESDVZERO* vzrData = fESD->GetVZEROData();
+//const AliESDVZERO* vzrData = (const_cast<AliESDEvent*>(fESD))->GetVZEROData();
+ const AliESDVertex* vtxData = fESD->GetPrimaryVertex();
+
+ etr = (int)(fESD->GetNumberOfCaloClusters());
+ zvt = (float)(vtxData->GetZ());
+ z1n = (float)(zdcData->GetZDCN1Energy());
+ z2n = (float)(zdcData->GetZDCN2Energy());
+ z1p = (float)(zdcData->GetZDCP1Energy());
+ z2p = (float)(zdcData->GetZDCP2Energy());
+ const Double_t *z1ntTmp,*z2ntTmp,*z1ptTmp,*z2ptTmp;
+ const Double_t *z1ntTmq,*z2ntTmq,*z1ptTmq,*z2ptTmq;
+ const Double32_t *t0timTmp,*t0ampTmp;
+ z1ntTmp = zdcData->GetZN1TowerEnergyLR();
+ z2ntTmp = zdcData->GetZN2TowerEnergyLR();
+ z1ptTmp = zdcData->GetZP1TowerEnergyLR();
+ z2ptTmp = zdcData->GetZP2TowerEnergyLR();
+ z1ntTmq = zdcData->GetZN1TowerEnergy();
+ z2ntTmq = zdcData->GetZN2TowerEnergy();
+ z1ptTmq = zdcData->GetZP1TowerEnergy();
+ z2ptTmq = zdcData->GetZP2TowerEnergy();
+ t0timTmp = tzrData->GetT0time();
+ t0ampTmp = tzrData->GetT0amplitude();
+ for (int i=0; i<5; i++) {
+ z1nt[i]=(float)(z1ntTmp[i]); z1nt[i+5]=(float)(z1ntTmq[i]);
+ z2nt[i]=(float)(z2ntTmp[i]); z2nt[i+5]=(float)(z2ntTmq[i]);
+ z1pt[i]=(float)(z1ptTmp[i]); z1pt[i+5]=(float)(z1ptTmq[i]);
+ z2pt[i]=(float)(z2ptTmp[i]); z2pt[i+5]=(float)(z2ptTmq[i]);
+ }
+ for (int i=0; i<24; i++) {
+ t0tim[i] = (float)(t0timTmp[i]);
+ t0amp[i] = (float)(t0ampTmp[i]);
+ }
+ for (int i=0; i<64; i++) {
+ v0mul[i] = (float)(vzrData->GetMultiplicity(i));
+ }
+ z1s=z1nt[1]+z1nt[2]+z1nt[3]+z1nt[4];
+ z2s=z2nt[1]+z2nt[2]+z2nt[3]+z2nt[4];
+ w1s=z1pt[1]+z1pt[2]+z1pt[3]+z1pt[4];
+ w2s=z2pt[1]+z2pt[2]+z2pt[3]+z2pt[4];
+ t1s=0; t2s=0; v1s=0; v2s=0;
+ for (Int_t n=0; n<24; n++) {
+ if (t0tim[n]>100 && t0amp[n]>0) {
+ if (n<12) t1s+=t0amp[n];
+ else t2s+=t0amp[n];
+ }
+ }
+ for (Int_t n=0; n<64; n++) {
+ if (v0mul[n]>0) {
+ if (n<32) v1s+=v0mul[n];
+ else v2s+=v0mul[n];
+ }
+ }
+ jev++;
+ if (ntr>0) {
+ fMyTr->Fill();
+ iev++;
+ printf(": %d %d %d %d %f : %f %f %f %f :\n",jev,iev,ntr,etr,zvt,z1s+z2s,w1s+w2s,t1s+t2s,v1s+v2s);
+ }
+
+/*
+ if(Ntrack==0 || nCluster==0) {
+ printf("no trk or clus %d %d\n",Ntrack,nCluster);
+ return;
+ }
+
+ Double_t bfield = fESD->GetMagneticField();
+ printf("bfield %f\n",bfield);
+
+ // Track loop to fill a pT spectrum
+ Double_t posI[5],posO[5],posE[5];
+ Double_t minR=300.0;
+ Double_t maxR=500.0;
+ Double_t stpR=10.0;
+ Float_t arr[5][22][1000];
+ int k=0;
+ for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
+ AliESDtrack* track = fESD->GetTrack(iTracks);
+ if (!track) {
+ printf("ERROR: Could not receive track %d\n", iTracks);
+ continue;
+ }
+ float pt=track->Pt();
+ fHstPt->Fill(pt);
+ float px=track->Px();
+ float py=track->Py();
+ float ph=atan2(py,px);
+ track->GetInnerXYZ(posI);
+ posI[3] = sqrt(posI[0]*posI[0]+posI[1]*posI[1]);
+ posI[4] = atan2(posI[1],posI[0]);
+ track->GetOuterXYZ(posO);
+ posO[3] = sqrt(posO[0]*posO[0]+posO[1]*posO[1]);
+ posO[4] = atan2(posO[1],posO[0]);
+ if (pt>1.0 && ph>1.0 && ph<2.5 &&
+ posI[3]>80.0 && posI[3]<90.0 &&
+ posO[3]>280.0 && posO[3]<300.0) {
+ for (int i=0; i<5; i++) arr[i][0][k]=posI[i];
+ for (int i=0; i<5; i++) arr[i][1][k]=posO[i];
+ int j=0;
+ for (Double_t rad=minR; rad<maxR; rad+=stpR) {
+ track->GetOuterParam()->GetXYZAt(rad,bfield,posE);
+ posE[3] = sqrt(posE[0]*posE[0]+posE[1]*posE[1]);
+ posE[4] = atan2(posE[1],posE[0]);
+ for (int i=0; i<5; i++) arr[i][2+j][k]=posE[i];
+ j++;
+ }
+ if (k<1000) k++;
+ else printf("error max trk 1000\n");
+ }
+ } //track loop
+ int nTrkEMC=k;
+
+ // EMC loop
+ float arr0[10000][7],arr1[10000][7];
+ float scl=1.0E+09;
+// TVector3 CaloCellPos;
+ AliESDCaloCells &cells= *(fESD->GetEMCALCells());
+ int totalCell = cells.GetNumberOfCells();
+ int totCell = 0;
+ int totClus = 0;
+ for (Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++) {
+ AliESDCaloCluster * clust = fESD->GetCaloCluster(iCluster);
+ if (clust->IsEMCAL())
+ {
+ Int_t nCells = clust->GetNCells();
+ UShort_t * index = clust->GetCellsAbsId() ;
+ Double_t * fraction = clust->GetCellsAmplitudeFraction() ;
+ Float_t pos[3]={0,0,0};
+ clust->GetPosition(pos);
+ float rad=sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
+ float phi=atan2(pos[1],pos[0]);
+ float zps=pos[2];
+ float tof = scl*clust->GetTOF();
+ float ene = clust->E();
+ fEMCe->Fill(ene);
+ fEMCt->Fill(tof);
+ fEMCn->Fill((float)nCells);
+ fEMCm->Fill(phi,zps);
+ arr0[totClus][0]=ene;
+ arr0[totClus][1]=tof;
+ arr0[totClus][2]=rad;
+ arr0[totClus][3]=phi;
+ arr0[totClus][4]=zps;
+ arr0[totClus][5]=totCell;
+ arr0[totClus][6]=totCell+nCells;
+ if (totClus<10000) totClus++;
+ else printf("error max cls 10000\n");
+ for(Int_t i = 0; i < nCells ; i++)
+ {
+ Int_t absId = index[i]; // or clus->GetCellNumber(i) ;
+ Double_t cpos[3]={0,0,0};
+//
+ if (fEMCALGeo) {
+// fEMCALGeo->GetGlobal(absId,CaloCellPos);
+ fEMCALGeo->GetGlobal(absId,cpos);
+ } else {
+ printf("Error, fEMCALGeo=0\n");
+ }
+// Double_t ccphi=CaloCellPos.Phi();
+// Double_t cceta=CaloCellPos.Eta();
+//
+ Double_t ccrad=sqrt(cpos[0]*cpos[0]+cpos[1]*cpos[1]);
+ Double_t ccphi=atan2(cpos[1],cpos[0]);
+ Double_t cczps=cpos[2];
+
+ Double_t ampFract = fraction[i];
+ Float_t amp = cells.GetCellAmplitude(absId) ;
+ Float_t time = scl*cells.GetCellTime(absId);
+ fCelle->Fill(amp);
+ fCellf->Fill(ampFract);
+ fCellt->Fill(time);
+ fCellc->Fill((float)absId,time);
+ fCelld->Fill((float)absId,amp);
+ fCellm->Fill(ccphi,cczps);
+ arr1[totCell][0]=i;
+ arr1[totCell][1]=absId;
+ arr1[totCell][2]=amp;
+ arr1[totCell][3]=time;
+ arr1[totCell][4]=ccrad;
+ arr1[totCell][5]=ccphi;
+ arr1[totCell][6]=cczps;
+ if (totCell<10000) totCell++;
+ else printf("error max cel 10000\n");
+ } // cell loop
+
+ } // end of cluster loop
+
+// Printf("+++++ # EMC cluster E = %f", mesE);
+
+ }
+ if(nTrkEMC==0 || totClus==0) {
+ printf("no emc trk or emc clus %d %d\n",nTrkEMC,totClus);
+ return;
+ }
+ FILE *ofs;
+ char opt[4];
+ if (iii==0) sprintf(opt,"w");
+ else sprintf(opt,"a");
+ if ((ofs = fopen("mydata.dat",opt)) == NULL) {
+ printf("error in fopen\n");
+ }
+ iii++;
+ printf("%d %s : %d %d : %d %d : %d %d\n", iii, opt,
+ Ntrack, nTrkEMC, nCluster, totClus, totalCell, totCell);
+//fprintf(ofs,"%d %d %d %d %d %d\n",
+//Ntrack, nTrkEMC, nCluster, totClus, totalCell, totCell);
+ for (int itrk=0; itrk<nTrkEMC; itrk++) {
+ for (int ihit=0; ihit<22; ihit++) {
+// fprintf(ofs,"%f %f %f\n",
+// arr[3][ihit][itrk],arr[4][ihit][itrk],arr[2][ihit][itrk]);
+ }
+ }
+ for (int icls=0; icls<totClus; icls++) {
+// for (int ihit=0; ihit<5; ihit++) fprintf(ofs,"%f ",arr0[icls][ihit]);
+// for (int ihit=5; ihit<7; ihit++) fprintf(ofs,"%d ",(int)arr0[icls][ihit]);
+// fprintf(ofs,"\n");
+ }
+ for (int icel=0; icel<totCell; icel++) {
+// for (int ihit=0; ihit<2; ihit++) fprintf(ofs,"%d ",(int)arr1[icel][ihit]);
+// for (int ihit=2; ihit<7; ihit++) fprintf(ofs,"%f ",arr1[icel][ihit]);
+// fprintf(ofs,"\n");
+ }
+ fclose(ofs);
+*/
+ // Post output data.
+ PostData(1, fOutputList);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPt::Terminate(Option_t *)
+{
+ // Draw result to the screen
+ // Called once at the end of the query
+
+/*
+ fOutputList = dynamic_cast<TList*> (GetOutputData(1));
+ if (!fOutputList) {
+ printf("ERROR: Output list not available\n");
+ return;
+ }
+
+ fHstPt = dynamic_cast<TH1F*> (fOutputList->At(0));
+ if (!fHstPt) {
+ printf("ERROR: fHstPt not available\n");
+ return;
+ }
+
+ TCanvas *c1 = new TCanvas("AliAnalysisTaskPt","Pt",10,10,510,510);
+ c1->cd(1)->SetLogy();
+ fHstPt->DrawCopy("E");
+*/
+
+}
--- /dev/null
+#ifndef AliAnalysisTaskPt_cxx\r
+#define AliAnalysisTaskPt_cxx\r
+\r
+// example of an analysis task creating a p_t spectrum\r
+// Authors: Panos Cristakoglou, Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing\r
+\r
+class TH1F;\r
+class TH2F;\r
+class AliESDEvent;\r
+// class AliEMCALGeometry;\r
+// class AliEMCALGeoUtils;\r
+\r
+// #include "AliEMCALGeometry.h"\r
+// #include "AliEMCALGeoUtils.h"\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class AliAnalysisTaskPt : public AliAnalysisTaskSE {\r
+ public:\r
+ AliAnalysisTaskPt() : AliAnalysisTaskSE(), fESD(0),\r
+ fOutputList(0), fMyTr(0), jev(0), iev(0) {} \r
+\r
+ AliAnalysisTaskPt(const char *name);\r
+ virtual ~AliAnalysisTaskPt() {}\r
+ \r
+ virtual void UserCreateOutputObjects();\r
+ virtual void UserExec(Option_t *option);\r
+ virtual void Terminate(Option_t *);\r
+ \r
+ private:\r
+ AliESDEvent *fESD; //! ESD object\r
+ TList *fOutputList; //! Output list\r
+ TTree *fMyTr;\r
+ int jev, iev;\r
+/*\r
+ TH1F *fHstPt;\r
+ TH1F *fEMCe;\r
+ TH1F *fEMCt;\r
+ TH1F *fEMCn;\r
+ TH2F *fEMCm;\r
+ TH1F *fCelle;\r
+ TH1F *fCellf;\r
+ TH1F *fCellt;\r
+ TH2F *fCellc;\r
+ TH2F *fCelld;\r
+ TH2F *fCellm;\r
+//AliEMCALGeometry *fEMCALGeo;\r
+ AliEMCALGeoUtils *fEMCALGeo;\r
+*/\r
+ \r
+ AliAnalysisTaskPt(const AliAnalysisTaskPt&); // not implemented\r
+ AliAnalysisTaskPt& operator=(const AliAnalysisTaskPt&); // not implemented\r
+ \r
+ ClassDef(AliAnalysisTaskPt, 1); // example of analysis\r
+};\r
+\r
+#endif\r
--- /dev/null
+#define AnaEveMyTree_cxx
+#include "AnaEveMyTree.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+
+void AnaEveMyTree::Loop()
+{
+// Root > .L AnaEveMyTree.C
+// Root > AnaEveMyTree t
+// Root > t.Loop(); // Loop on all entries
+//
+// T0C -3.3 ~ -2.9
+// T0A 4.5 ~ 5.0
+//
+// V0C -1.7 ~ -2.2 ~ -2.7 ~ -3.2 ~ -3.7
+// V0A 2.8 ~ 3.4 ~ 3.9 ~ 4.5 ~ 5.1
+//
+ float pi=acos(-1.0);
+ float t0phi[24];
+ for(int it=0; it<12; it++) {
+ t0phi[it]=(it+1.0)*2.0*pi/12.0;
+ t0phi[it+12]=it*2.0*pi/12.0;
+ }
+ for(int it=0; it<24; it++) {
+ float phi=t0phi[it];
+ t0phi[it]=atan2(sin(phi),cos(phi));
+ }
+ float v0phi[64];
+ int is=0;
+ for(int iRing = 0; iRing<4 ; iRing++){
+ for(int iSec = 0; iSec<8 ; iSec++){
+ v0phi[is] = 22.5 + 45. * iSec;
+ is++;
+ }
+ }
+ for(int iRing = 0; iRing<4 ; iRing++){
+ for(int iSec = 0; iSec<8 ; iSec++){
+ v0phi[is] = 22.5 + 45. * iSec;
+ is++;
+ }
+ }
+ for (int i=0; i<64; i++) {
+ float phi=v0phi[i]*pi/180.0;
+ v0phi[i]=atan2(sin(phi),cos(phi));
+ cout << v0phi[i] << " ";
+ if (i%8==7) cout << endl;
+ }
+ float z0gai[8],t0gai[24],v0gai[64];
+ float mean[10][10][4][27][2];
+ float widt[10][10][4][27][2];
+ float four[10][10][4][27][2][8];
+ float f0,f1,f2,f3,f4,f5,f6,f7;
+ ifstream ifs;
+ ifs.open("AnaEveMyTree.cal");
+ for (int i=0; i<1; i++) {
+ ifs>>f0>>f1>>f2>>f3>>f4>>f5>>f6>>f7;
+ z0gai[i*8+0]=f0; z0gai[i*8+1]=f1;
+ z0gai[i*8+2]=f2; z0gai[i*8+3]=f3;
+ z0gai[i*8+4]=f4; z0gai[i*8+5]=f5;
+ z0gai[i*8+6]=f6; z0gai[i*8+7]=f7;
+ cout << f0 << " " << f1 << " " << f2 << " " << f3 << " "
+ << f4 << " " << f5 << " " << f6 << " " << f7 << endl;
+ }
+ for (int i=0; i<3; i++) {
+ ifs>>f0>>f1>>f2>>f3>>f4>>f5>>f6>>f7;
+ t0gai[i*8+0]=f0; t0gai[i*8+1]=f1;
+ t0gai[i*8+2]=f2; t0gai[i*8+3]=f3;
+ t0gai[i*8+4]=f4; t0gai[i*8+5]=f5;
+ t0gai[i*8+6]=f6; t0gai[i*8+7]=f7;
+ cout << f0 << " " << f1 << " " << f2 << " " << f3 << " "
+ << f4 << " " << f5 << " " << f6 << " " << f7 << endl;
+ }
+ for (int i=0; i<8; i++) {
+ ifs>>f0>>f1>>f2>>f3>>f4>>f5>>f6>>f7;
+ v0gai[i*8+0]=f0; v0gai[i*8+1]=f1;
+ v0gai[i*8+2]=f2; v0gai[i*8+3]=f3;
+ v0gai[i*8+4]=f4; v0gai[i*8+5]=f5;
+ v0gai[i*8+6]=f6; v0gai[i*8+7]=f7;
+ cout << f0 << " " << f1 << " " << f2 << " " << f3 << " "
+ << f4 << " " << f5 << " " << f6 << " " << f7 << endl;
+ }
+ for (int ic=0; ic<10; ic++) {
+ for (int iz=0; iz<10; iz++) {
+ for (int ih=0; ih<4; ih++) {
+ for (int id=0; id<27; id++) {
+ for (int ib=0; ib<2; ib++) {
+ mean[ic][iz][ih][id][ib]=0.0;
+ widt[ic][iz][ih][id][ib]=1.0;
+ }
+ ifs >> f0 >> f1 >> f2 >> f3;
+ if (f1==0.0 || f1<0.0) f1=1.0;
+ if (f3==0.0 || f3<0.0) f3=1.0;
+ mean[ic][iz][ih][id][0]=f0;
+ widt[ic][iz][ih][id][0]=f1;
+ mean[ic][iz][ih][id][1]=f2;
+ widt[ic][iz][ih][id][1]=f3;
+ if (ic==4 && iz==4 && ih==0) {
+ cout << f0 << " " << f1 << " " << f2 << " " << f3 << endl;
+ }
+ for (int ib=0; ib<2; ib++) {
+ for (int io=0; io<8; io++) {
+ four[ic][iz][ih][id][ib][io]=0.0;
+ }
+ ifs >> f0 >> f1 >> f2 >> f3 >> f4 >> f5 >> f6 >> f7;
+ four[ic][iz][ih][id][ib][0]=f0;
+ four[ic][iz][ih][id][ib][1]=f1;
+ four[ic][iz][ih][id][ib][2]=f2;
+ four[ic][iz][ih][id][ib][3]=f3;
+ four[ic][iz][ih][id][ib][4]=f4;
+ four[ic][iz][ih][id][ib][5]=f5;
+ four[ic][iz][ih][id][ib][6]=f6;
+ four[ic][iz][ih][id][ib][7]=f7;
+ if (ic==4 && iz==4 && ih==0) {
+ cout << f0 << " " << f1 << " " << f2 << " " << f3 << " "
+ << f4 << " " << f5 << " " << f6 << " " << f7 << endl;
+ }
+ }
+ }
+ }
+ }
+ }
+ ifs.close();
+ float av1=0;
+ float no1=0;
+ for (int i=0; i<8; i++) {
+ if (z0gai[i]>3) {
+ av1+=z0gai[i];
+ no1+=1.0;
+ } else {
+ z0gai[i]=-1.0;
+ }
+ }
+ av1/=no1;
+ float av2=0;
+ float no2=0;
+ for (int i=0; i<24; i++) {
+ if (t0gai[i]>3) {
+ av2+=t0gai[i];
+ no2+=1.0;
+ } else {
+ t0gai[i]=-1.0;
+ }
+ }
+ av2/=no2;
+ float av3=0;
+ float no3=0;
+ for (int i=0; i<64; i++) {
+ if (v0gai[i]>3) {
+ av3+=v0gai[i];
+ no3+=1.0;
+ } else {
+ v0gai[i]=-1.0;
+ }
+ }
+ av3/=no3;
+ cout << av1 << " " << av2 << " " << av3 << endl;
+ TFile *tf = new TFile("AnaEveMyTree.root","recreate");
+ TH2D *cnzv = new TH2D("cnzv","",50,0.0,11000.0,50,-15.0,15.0);
+ TProfile *det0g = new TProfile("det0g","", 16,0.0, 16.0 ,-600.0,6000.0);
+ TH2D *det0m = new TH2D ("det0m","", 16,0.0, 16.0,50,-300.0,3000.0);
+ TProfile *det1g = new TProfile("det1g","", 48,0.0, 48.0 , -20.0, 200.0);
+ TH2D *det1m = new TH2D ("det1m","", 48,0.0, 48.0,50, -10.0, 100.0);
+ TProfile *det2g = new TProfile("det2g","",128,0.0,128.0 ,-100.0,1000.0);
+ TH2D *det2m = new TH2D ("det2m","",128,0.0,128.0,50, -50.0, 500.0);
+ char name[80];
+ TH2D *zdcxy[2][5];
+ for (int i=0; i<2; i++) {
+ for (int j=0; j<5; j++) {
+ sprintf(name,"zdcxy%d%d",i,j);
+ zdcxy[i][j] = new TH2D(name,"",50,-2,2,50,-2,2);
+ }
+ }
+ float rng0=11000.0;
+ float rng1=av3*64.0*4.5;
+ float rng2=av2*12.0*7.0;
+ float rng3=av1* 8.0*2.2;
+ TH2D *cor00 = new TH2D("corr00","",50,0.0,rng0,50,0.0,rng1);
+ TH2D *cor01 = new TH2D("corr01","",50,0.0,rng0,50,0.0,rng2);
+ TH2D *cor02 = new TH2D("corr02","",50,0.0,rng0,50,0.0,rng3);
+ TH2D *cor03 = new TH2D("corr03","",50,0.0,rng1,50,0.0,rng2);
+ TH2D *cor04 = new TH2D("corr04","",50,0.0,rng1,50,0.0,rng3);
+ TH2D *cor05 = new TH2D("corr05","",50,0.0,rng2,50,0.0,rng3);
+ TH2D *cor10 = new TH2D("corr10","",50,0.0,rng0,50,0.0,rng1);
+ TH2D *cor11 = new TH2D("corr11","",50,0.0,rng0,50,0.0,rng2);
+ TH2D *cor12 = new TH2D("corr12","",50,0.0,rng0,50,0.0,rng3);
+ TH2D *cor13 = new TH2D("corr13","",50,0.0,rng1,50,0.0,rng2);
+ TH2D *cor14 = new TH2D("corr14","",50,0.0,rng1,50,0.0,rng3);
+ TH2D *cor15 = new TH2D("corr15","",50,0.0,rng2,50,0.0,rng3);
+ TProfile *ave[10][10][4][27];
+ TProfile *flt[10][10][4][27];
+ for (int ic=0; ic<10; ic++) {
+ for (int iz=0; iz<10; iz++) {
+ for (int ih=0; ih<4; ih++) {
+ for (int id=0; id<27; id++) {
+ sprintf(name,"ave%d%d%d%d",ic,iz,ih,id);
+ ave[ic][iz][ih][id] = new TProfile(name,name,4,-0.5,3.5,-1.0E+06,1.0E+06,"S");
+ sprintf(name,"flt%d%d%d%d",ic,iz,ih,id);
+ flt[ic][iz][ih][id] = new TProfile(name,name,32,-0.5,31.5,-1.0E+06,1.0E+06);
+ }
+ }
+ }
+ }
+ TProfile *res[4][22][2];
+ TH2D *cor[4][22][5];
+ for (int ih=0; ih<4; ih++) {
+ for (int id=0; id<22; id++) {
+ for (int ia=0; ia<2; ia++) {
+ sprintf(name,"res%d%d%d",ih,id,ia);
+ res[ih][id][ia] = new TProfile(name,name,10,-0.5,9.5,-1.1,1.1);
+ }
+ for (int ia=0; ia<5; ia++) {
+ sprintf(name,"cor%d%d%d",ih,id,ia);
+ cor[ih][id][ia] = new TH2D(name,name,50,-pi,pi,50,-pi,pi);
+ }
+ }
+ }
+ TH1D *t0d[24][5],*t0e[24][5];
+ for (int i=0; i<24; i++) {
+ for (int j=0; j<5; j++) {
+ if (j<2) {
+ sprintf(name,"t0d%d_%d",i,j);
+ t0d[i][j] = new TH1D(name,"",24,-pi,pi);
+ sprintf(name,"t0e%d_%d",i,j);
+ t0e[i][j] = new TH1D(name,"",24,-pi,pi);
+ } else {
+ sprintf(name,"t0d%d_%d",i,j);
+ t0d[i][j] = new TH1D(name,"",24,-pi/2.,pi/2.);
+ sprintf(name,"t0e%d_%d",i,j);
+ t0e[i][j] = new TH1D(name,"",24,-pi/2.,pi/2.);
+ }
+ }
+ }
+ TH2D *tvcr0[10][4],*tvcr1[10][4];
+ for (int i=0; i<10; i++) {
+ for (int j=0; j<4; j++) {
+ sprintf(name,"tvcr0_%d_%d",i,j);
+ tvcr0[i][j] = new TH2D(name,"",48,0,48,24,-pi,pi);
+ sprintf(name,"tvcr1_%d_%d",i,j);
+ tvcr1[i][j] = new TH2D(name,"",128,0,128,24,-pi,pi);
+ }
+ }
+ int neve[10][10];
+ float buft[10][10][3][24];
+ float bufv[10][10][3][64];
+ for (int i=0; i<10; i++) {
+ for (int j=0; j<10; j++) {
+ neve[i][i]=0;
+ }
+ }
+ if (fChain == 0) return;
+ Long64_t nentries = fChain->GetEntriesFast();
+ cout << nentries << endl;
+ Long64_t nbytes = 0, nb = 0;
+ int ievent==0;
+ for (Long64_t jentry=0; jentry<nentries;jentry++) {
+ Long64_t ientry = LoadTree(jentry);
+ if (ientry < 0) break;
+ nb = fChain->GetEntry(jentry); nbytes += nb;
+ if (t1s > 1 && t2s > 1
+ && v1s > 10 && v2s > 10
+ && z1s > 100 && z2s >100
+ && zvt > -10 && zvt < 10
+ && v2s < 0.8*ntr+1000
+ && v2s > 0.8*ntr-1000
+ && t1s > 0.02*ntr) {
+ cnzv->Fill(ntr,zvt);
+ int icen=(int)(ntr/900.0);
+ if (icen<0) icen=0;
+ if (icen>9) icen=9;
+ int izvt=(int)((zvt+10.0)/2.0);
+ if (izvt<0) izvt=0;
+ if (izvt>9) izvt=9;
+
+ float sumxy[4][27][4];
+ for (int ih=0; ih<4; ih++) {
+ for (int id=0; id<27; id++) {
+ for (int iv=0; iv<4; iv++) {
+ sumxy[ih][id][iv]=0;
+ }
+ }
+ }
+ float z0co[8];
+ float z0sm=0;
+ float z0x[8] = {-1, 1, -1, 1, 1, -1, 1, -1};
+ float z0y[8] = {-1, -1, 1, 1, -1, -1, 1, 1};
+ for (int i=0; i<8; i++) {
+ float gain=z0gai[i];
+ if (gain<0) gain=0;
+ else gain=av1/gain;
+ if (i<4) { // for C-side
+ det0g->Fill(i+0.5,z1nt[i+1]);
+ det0m->Fill(i+0.5,z1nt[i+1]);
+ det0g->Fill(i+8.5,z1nt[i+1]*gain);
+ det0m->Fill(i+8.5,z1nt[i+1]*gain);
+ z0co[i]=z1nt[i+1]*gain;
+ z0sm+=z1nt[i+1]*gain;
+ } else { // for A-side
+ det0g->Fill(i+0.5,z2nt[i-3]);
+ det0m->Fill(i+0.5,z2nt[i-3]);
+ det0g->Fill(i+8.5,z2nt[i-3]*gain);
+ det0m->Fill(i+8.5,z2nt[i-3]*gain);
+ z0co[i]=z2nt[i-3]*gain;
+ z0sm+=z2nt[i-3]*gain;
+ }
+ sumxy[0][i/4][0]+=z0co[i]*z0x[i];
+ sumxy[0][i/4][1]+=z0co[i]*z0y[i];
+ sumxy[0][i/4][2]+=z0co[i];
+ float sgn=-1.0*(1-(int)(i/4)); // sign flip for C-side
+ sumxy[0][2][0]+=z0co[i]*z0x[i]*sgn;
+ sumxy[0][2][1]+=z0co[i]*z0y[i]*sgn;
+ sumxy[0][2][2]+=z0co[i];
+ }
+ for (int i=0; i<3; i++) {
+ sumxy[0][i][0]/=sumxy[0][i][2];
+ sumxy[0][i][1]/=sumxy[0][i][2];
+ }
+ zdcxy[0][icen/2]->Fill(sumxy[0][0][0],sumxy[0][1][0]);
+ zdcxy[1][icen/2]->Fill(sumxy[0][0][1],sumxy[0][1][1]);
+ float t0co[24];
+ float t0sm=0;
+ for (int i=0; i<24; i++) {
+ float gain=t0gai[i];
+ if (gain<0) gain=0;
+ else gain=av2/gain;
+ if (t0tim[i]>100 && t0amp[i]>0) {
+ det1g->Fill(i+0.5,t0amp[i]);
+ det1m->Fill(i+0.5,t0amp[i]);
+ det1g->Fill(i+24.5,t0amp[i]*gain);
+ det1m->Fill(i+24.5,t0amp[i]*gain);
+ t0co[i]=t0amp[i]*gain;
+ t0sm+=t0amp[i]*gain;
+ } else {
+ t0co[i]=0.0;
+ }
+ }
+ float v0co[64];
+ float v0sm=0;
+ for (int i=0; i<64; i++) {
+ float gain=v0gai[i];
+ if (gain<0) gain=0;
+ else gain=av3/gain;
+ if (v0mul[i]>0) {
+ det2g->Fill(i+0.5,v0mul[i]);
+ det2m->Fill(i+0.5,v0mul[i]);
+ det2g->Fill(i+64.5,v0mul[i]*gain);
+ det2m->Fill(i+64.5,v0mul[i]*gain);
+ v0co[i]=v0mul[i]*gain;
+ v0sm+=v0mul[i]*gain;
+ } else {
+ v0co[i]=0.0;
+ }
+ }
+ for (int i=0; i<14; i++) {
+ int id=(int)(i/12);
+ for (int ih=0; ih<4; ih++) {
+ float sgn=1;
+ if (ih%2==0) sgn=-1.0*(1-id); // sign flip for C-side
+ sumxy[ih][3+id][0]+=t0co[i]*cos((ih+1.0)*t0phi[i]);
+ sumxy[ih][3+id][1]+=t0co[i]*sin((ih+1.0)*t0phi[i]);
+ sumxy[ih][3+id][2]+=t0co[i];
+ sumxy[ih][5][0]+=t0co[i]*cos((ih+1.0)*t0phi[i]);
+ sumxy[ih][5][1]+=t0co[i]*sin((ih+1.0)*t0phi[i]);
+ sumxy[ih][5][2]+=t0co[i];
+ sumxy[ih][6][0]+=t0co[i]*cos((ih+1.0)*t0phi[i])*sgn;
+ sumxy[ih][6][1]+=t0co[i]*sin((ih+1.0)*t0phi[i])*sgn;
+ sumxy[ih][6][2]+=t0co[i];
+ }
+ }
+ for (int i=0; i<64; i++) {
+ int id=(int)(i/32);
+ int jd=(int)(i/8);
+ int kd=jd%4;
+ for (int ih=0; ih<4; ih++) {
+ float sgn=1;
+ if (ih%2==0) sgn=-1.0*(1-id); // sign flip for C-side
+ sumxy[ih][7+id][0]+=v0co[i]*cos((ih+1.0)*v0phi[i]);
+ sumxy[ih][7+id][1]+=v0co[i]*sin((ih+1.0)*v0phi[i]);
+ sumxy[ih][7+id][2]+=v0co[i];
+ sumxy[ih][9][0]+=v0co[i]*cos((ih+1.0)*v0phi[i]);
+ sumxy[ih][9][1]+=v0co[i]*sin((ih+1.0)*v0phi[i]);
+ sumxy[ih][9][2]+=v0co[i];
+ sumxy[ih][10][0]+=v0co[i]*cos((ih+1.0)*v0phi[i])*sgn;
+ sumxy[ih][10][1]+=v0co[i]*sin((ih+1.0)*v0phi[i])*sgn;
+ sumxy[ih][10][2]+=v0co[i];
+ sumxy[ih][11+4*kd+id][0]+=v0co[i]*cos((ih+1.0)*v0phi[i]);
+ sumxy[ih][11+4*kd+id][1]+=v0co[i]*sin((ih+1.0)*v0phi[i]);
+ sumxy[ih][11+4*kd+id][2]+=v0co[i];
+ sumxy[ih][13+4*kd][0]+=v0co[i]*cos((ih+1.0)*v0phi[i]);
+ sumxy[ih][13+4*kd][1]+=v0co[i]*sin((ih+1.0)*v0phi[i]);
+ sumxy[ih][13+4*kd][2]+=v0co[i];
+ sumxy[ih][14+4*kd][0]+=v0co[i]*cos((ih+1.0)*v0phi[i])*sgn;
+ sumxy[ih][14+4*kd][1]+=v0co[i]*sin((ih+1.0)*v0phi[i])*sgn;
+ sumxy[ih][14+4*kd][2]+=v0co[i];
+ }
+ }
+ if (icen>0 && icen<7) {
+ for (int i=0; i<24; i++) {
+ float amp1=t0co[i]/av2;
+ float phi1=t0phi[i];
+ for (int j=0; j<64; j++) {
+ int id=(int)(j/32);
+ float amp2=v0co[j]/av3;
+ float phi2=v0phi[j];
+ float dphi=phi2-phi1;
+ dphi=atan2(sin(dphi),cos(dphi));
+ t0d[i][id]->Fill(phi2,amp1*amp2);
+ t0e[i][id]->Fill(dphi,amp1*amp2);
+ }
+ }
+ }
+ int ieve=neve[icen][izvt];
+ for (int i=0; i<24; i++) buft[icen][izvt][ieve][i]=t0co[i]/av2;
+ for (int i=0; i<64; i++) bufv[icen][izvt][ieve][i]=v0co[i]/av3;
+ neve[icen][izvt]++;
+ if (neve[icen][izvt]==3) {
+ neve[icen][izvt]=0;
+ for (int ie=0; ie<3; ie++) {
+ for (int je=0; je<3; je++) {
+ int mix=1;
+ if (ie==je) mix=0;
+ for (int ip=0; ip<24; ip++) {
+ int iq=(int)(ip/12);
+ float phi1=t0phi[ip];
+ float amp1=buft[icen][izvt][ie][ip];
+ for (int jp=0; jp<64; jp++) {
+ int jq=(int)(jp/32);
+ float phi2=v0phi[jp];
+ float amp2=bufv[icen][izvt][je][jp];
+ float dphi=phi2-phi1;
+ dphi=atan2(sin(dphi),cos(dphi));
+ tvcr0[icen][iq*2+jq]->Fill(mix*24+ip+0.5,dphi,amp1*amp2);
+ tvcr1[icen][iq*2+jq]->Fill(mix*64+jp+0.5,dphi,amp1*amp2);
+ }
+ }
+ }
+ }
+ }
+
+// r.p. flattening // -----------------------------------
+ int calFlag=1;
+ for (int ih=0; ih<4; ih++) {
+ for (int id=0; id<27; id++) {
+ for (int ib=0; ib<2; ib++) {
+ if (calFlag>0) ave[icen][izvt][ih][id]->Fill(ib+0.0,sumxy[ih][id][ib]);
+ float sxy=sumxy[ih][id][ib];
+ float mxy=mean[icen][izvt][ih][id][ib];
+ float wxy=widt[icen][izvt][ih][id][ib];
+ sumxy[ih][id][ib]=(sxy-mxy)/wxy;
+ if (calFlag>0) ave[icen][izvt][ih][id]->Fill(ib+2.0,sumxy[ih][id][ib]);
+ }
+ if (sumxy[ih][id][2]>0.0) {
+ sumxy[ih][id][3]=atan2(sumxy[ih][id][1],sumxy[ih][id][0])/(ih+1.0);
+ float psi=sumxy[ih][id][3]*(ih+1.0);
+ float dp=0.0;
+ for (int io=0; io<8; io++) {
+ float cc=cos((io+1.0)*psi);
+ float ss=sin((io+1.0)*psi);
+ if (calFlag>0) flt[icen][izvt][ih][id]->Fill(io+0.0,cc);
+ if (calFlag>0) flt[icen][izvt][ih][id]->Fill(io+8.0,ss);
+ float aa=four[icen][izvt][ih][id][0][io]; // mean cos
+ float bb=four[icen][izvt][ih][id][1][io]; // mean sin
+ dp+=(aa*ss-bb*cc)*2.0/(io+1.0);
+ }
+ psi+=dp;
+ for (int io=0; io<8; io++) {
+ float cc=cos((io+1.0)*psi);
+ float ss=sin((io+1.0)*psi);
+ if (calFlag>0) flt[icen][izvt][ih][id]->Fill(io+16.0,cc);
+ if (calFlag>0) flt[icen][izvt][ih][id]->Fill(io+24.0,ss);
+ }
+ sumxy[ih][id][3]=atan2(sin(psi),cos(psi))/(ih+1.0);
+ } else {
+ sumxy[ih][id][3]=-9999.9;
+ }
+ }
+ }
+
+ if (icen>0 && icen<7) {
+ for (int i=0; i<24; i++) {
+ float amp1=t0co[i]/av2;
+ float phi1=t0phi[i];
+ for (int id=7; id<10; id++) {
+ float phi2=sumxy[1][id][3];
+ if (phi2>-9000) {
+ float dphi=phi2-phi1;
+ dphi=atan2(sin(2.0*dphi),cos(2.0*dphi))/2.0;
+ t0d[i][id-5]->Fill(phi2,amp1);
+ t0e[i][id-5]->Fill(dphi,amp1);
+ }
+ }
+ }
+ }
+
+// r.p. correlation // -----------------------------------
+ for (int ih=0; ih<4; ih++) {
+ for (int id=0; id<22; id++) {
+ int id1=0;
+ int id2=0;
+ int ih1=ih;
+ int ih2=ih;
+ if (id==0) {id1=0; id2=1; ih1=0; ih2=0;} // z0c-z0c (1st)
+ if (id==1) {id1=3; id2=4;} // t0c-t0a
+ if (id==2) {id1=7; id2=8;} // v0c-v0a
+ if (id==3) {id1=2; id2=5; ih1=0;} // zd1-t0ac 6 for t0ac flip
+ if (id==4) {id1=2; id2=9; ih1=0;} // zd1-v0ac 10 for v0ac flip
+ if (id==5) {id1=3; id2=7;} // t0c-v0c
+ if (id==6) {id1=3; id2=8;} // t0c-v0a
+ if (id==7) {id1=4; id2=7;} // t0a-v0c
+ if (id==8) {id1=4; id2=8;} // t0a-v0a
+ if (id==9) {id1=5; id2=9;} // t0ac-v0ac
+ if (id==10) {id1=11; id2=12;} // v0c0-v0a0
+ if (id==11) {id1=15; id2=16;} // v0c1-v0a1
+ if (id==12) {id1=19; id2=20;} // v0c2-v0a2
+ if (id==13) {id1=23; id2=24;} // v0c3-v0a3
+ if (id==14) {id1=3; id2=12;} // t0c-v0a0
+ if (id==15) {id1=3; id2=16;} // t0c-v0a1
+ if (id==16) {id1=3; id2=20;} // t0c-v0a2
+ if (id==17) {id1=3; id2=24;} // t0c-v0a3
+ if (id==18) {id1=4; id2=11;} // t0a-v0c0
+ if (id==19) {id1=4; id2=15;} // t0a-v0c1
+ if (id==20) {id1=4; id2=19;} // t0a-v0c2
+ if (id==21) {id1=4; id2=23;} // t0a-v0c3
+ if (sumxy[ih1][id1][3]>-9000 && sumxy[ih2][id2][3]>-9000) {
+ float dph=sumxy[ih1][id1][3]-sumxy[ih2][id2][3];
+ res[ih][id][0]->Fill((float)icen,cos((ih+1.0)*dph));
+ res[ih][id][1]->Fill((float)icen,sin((ih+1.0)*dph));
+ float phi1=sumxy[ih1][id1][3];
+ float phi2=sumxy[ih2][id2][3];
+ if (ih>ih1) phi1*=(ih+1.0);
+ else phi1*=(ih1+1.0);
+ if (ih>ih2) phi2*=(ih+1.0);
+ else phi2*=(ih2+1.0);
+ phi1=atan2(sin(phi1),cos(phi1));
+ phi2=atan2(sin(phi2),cos(phi2));
+ cor[ih][id][icen/2]->Fill(phi1,phi2);
+ }
+ }
+ }
+
+ corr00->Fill(1.0*ntr,v1s+v2s);
+ corr01->Fill(1.0*ntr,t1s+t2s);
+ corr02->Fill(1.0*ntr,z1s+z2s);
+ corr03->Fill(v1s+v2s,t1s+t2s);
+ corr04->Fill(v1s+v2s,z1s+z2s);
+ corr05->Fill(t1s+t2s,z1s+z2s);
+ corr10->Fill(1.0*ntr,v0sm);
+ corr11->Fill(1.0*ntr,t0sm);
+ corr12->Fill(1.0*ntr,z0sm);
+ corr13->Fill(v0sm,t0sm);
+ corr14->Fill(v0sm,z0sm);
+ corr15->Fill(t0sm,z0sm);
+ if (ievent%1000==0) cout << "eve " << jentry << " " << ievent << " : "
+ << jev << " " << iev << " : " << ntr << " : "
+ << z1s+z2s << " " << t1s+t2s << " " << v1s+v2s << " : "
+ << z0sm << " " << t0sm << " " << v0sm << endl;
+ ievent++;
+ }
+ }
+ ofstream ofs;
+ ofs.open("AnaEveMyTree.cal");
+ for (int i=0; i<8; i++) {
+ ofs << det0g->GetBinContent(i+1) << " ";
+ if (i%8==7) ofs << endl;
+ cout << det0g->GetBinContent(i+1) << " ";
+ if (i%8==7) cout << endl;
+ }
+ for (int i=0; i<24; i++) {
+ ofs << det1g->GetBinContent(i+1) << " ";
+ if (i%8==7) ofs << endl;
+ cout << det1g->GetBinContent(i+1) << " ";
+ if (i%8==7) cout << endl;
+ }
+ for (int i=0; i<64; i++) {
+ ofs << det2g->GetBinContent(i+1) << " ";
+ if (i%8==7) ofs << endl;
+ cout << det2g->GetBinContent(i+1) << " ";
+ if (i%8==7) cout << endl;
+ }
+ for (int ic=0; ic<10; ic++) {
+ for (int iz=0; iz<10; iz++) {
+ for (int ih=0; ih<4; ih++) {
+ for (int id=0; id<27; id++) {
+ for (int ib=0; ib<2; ib++) {
+ ofs << ave[ic][iz][ih][id]->GetBinContent(ib+1) << " ";
+ ofs << ave[ic][iz][ih][id]->GetBinError (ib+1) << " ";
+ if (ic==4 && iz==4 && ih==0) {
+ cout << ave[ic][iz][ih][id]->GetBinContent(ib+1) << " ";
+ cout << ave[ic][iz][ih][id]->GetBinError (ib+1) << " ";
+ }
+ }
+ ofs << endl;
+ if (ic==4 && iz==4 && ih==0) cout << endl;
+ for (int ib=0; ib<2; ib++) {
+ for (int io=0; io<8; io++) {
+ ofs << flt[ic][iz][ih][id]->GetBinContent(ib*8+io+1) << " ";
+ if (ic==4 && iz==4 && ih==0) {
+ cout << flt[ic][iz][ih][id]->GetBinContent(ib*8+io+1) << " ";
+ }
+ }
+ ofs << endl;
+ if (ic==4 && iz==4 && ih==0) cout << endl;
+ }
+ }
+ }
+ }
+ }
+ ofs.close();
+ tf->Write();
+ tf->Close();
+}
--- /dev/null
+#! /bin/csh -f
+
+source root.setup
+@ ii = 0
+# foreach run (137549)
+# foreach file (1 4 5 6 7 8 9 10 11 12 14 15 16 18 19 23 24 25 26 27 28 29 30 31 32 33 35 37 38 39 40 42 43 44 45 46 47 49 51 52 53 57 58 60 62 64 65 66 67 68 69 70 71 72 73 74 75)
+foreach run (137431)
+foreach file (1 2 9 11 27 28 39 40 41 99)
+foreach nn (1 2 3 4 5)
+ @ ii = $ii + 1
+ echo "++++++++++++++++++++++++++++++++++++++++++++++"
+ echo "+++++" num $ii run $run file $file step $nn "+++++"
+ echo "++++++++++++++++++++++++++++++++++++++++++++++"
+ date
+ echo "++++++++++++++++++++++++++++++++++++++++++++++"
+
+root.exe -b > AnaEveMyTree_${nn}.log <<EOF
+ .L AnaEveMyTree.C
+ AnaEveMyTree a(${run},${file})
+ a.Loop()
+ .q
+EOF
+
+end
+mv AnaEveMyTree.root root${run}/AnaEveMyTree_${file}.root
+mv AnaEveMyTree.cal cal${run}/AnaEveMyTree_${file}.cal
+mv AnaEveMyTree_1.log log${run}/AnaEveMyTree_${file}_1.log
+mv AnaEveMyTree_2.log log${run}/AnaEveMyTree_${file}_2.log
+mv AnaEveMyTree_3.log log${run}/AnaEveMyTree_${file}_3.log
+mv AnaEveMyTree_4.log log${run}/AnaEveMyTree_${file}_4.log
+mv AnaEveMyTree_5.log log${run}/AnaEveMyTree_${file}_5.log
+end
+end
+
+echo "++++++++++++++++++++++++++++++++++++++++++++++"
+date
+echo "++++++++++++++++++++++++++++++++++++++++++++++"
+echo finished
+echo "++++++++++++++++++++++++++++++++++++++++++++++"
--- /dev/null
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Fri Nov 26 11:48:52 2010 by ROOT version 5.26/00
+// from TTree fMyTr/beam tree
+// found on file: Memory Directory
+//////////////////////////////////////////////////////////
+
+#ifndef AnaEveMyTree_h
+#define AnaEveMyTree_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+
+class AnaEveMyTree {
+public :
+ TTree *fChain; //!pointer to the analyzed TTree or TChain
+ Int_t fCurrent; //!current Tree number in a TChain
+
+ // Declaration of leaf types
+ Int_t jev;
+ Int_t iev;
+ Int_t ntr;
+ Int_t etr;
+ Float_t z1n;
+ Float_t z2n;
+ Float_t z1p;
+ Float_t z2p;
+ Float_t zvt;
+ Float_t z1s;
+ Float_t z2s;
+ Float_t w1s;
+ Float_t w2s;
+ Float_t t1s;
+ Float_t t2s;
+ Float_t v1s;
+ Float_t v2s;
+ Float_t z1nt[10];
+ Float_t z2nt[10];
+ Float_t z1pt[10];
+ Float_t z2pt[10];
+ Float_t t0amp[24];
+ Float_t t0tim[24];
+ Float_t v0mul[64];
+
+ // List of branches
+ TBranch *b_jev; //!
+ TBranch *b_iev; //!
+ TBranch *b_ntr; //!
+ TBranch *b_etr; //!
+ TBranch *b_z1n; //!
+ TBranch *b_z2n; //!
+ TBranch *b_z1p; //!
+ TBranch *b_z2p; //!
+ TBranch *b_zvt; //!
+ TBranch *b_z1s; //!
+ TBranch *b_z2s; //!
+ TBranch *b_w1s; //!
+ TBranch *b_w2s; //!
+ TBranch *b_t1s; //!
+ TBranch *b_t2s; //!
+ TBranch *b_v1s; //!
+ TBranch *b_v2s; //!
+ TBranch *b_z1nt; //!
+ TBranch *b_z2nt; //!
+ TBranch *b_z1pt; //!
+ TBranch *b_z2pt; //!
+ TBranch *b_t0amp; //!
+ TBranch *b_t0tim; //!
+ TBranch *b_v0mul; //!
+
+ AnaEveMyTree(int irun=137431, int ifile=2, TTree *tree=0);
+ virtual ~AnaEveMyTree();
+ virtual Int_t Cut(Long64_t entry);
+ virtual Int_t GetEntry(Long64_t entry);
+ virtual Long64_t LoadTree(Long64_t entry);
+ virtual void Init(TTree *tree);
+ virtual void Loop();
+ virtual Bool_t Notify();
+ virtual void Show(Long64_t entry = -1);
+};
+
+#endif
+
+#ifdef AnaEveMyTree_cxx
+AnaEveMyTree::AnaEveMyTree(int irun, int ifile, TTree *tree)
+{
+// if parameter tree is not specified (or zero), connect the ifile
+// used to generate this class and read the Tree.
+ char name[80];
+ sprintf(name,"%06d/EveMyTree%03d.root",irun,ifile);
+ cout << "open file :" << name << endl;
+ if (tree == 0) {
+ TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject(name);
+ if (!f) {
+ f = new TFile(name);
+ }
+ tree = (TTree*)gROOT->FindObject("chist")->FindObject("fMyTr");
+ }
+ Init(tree);
+}
+
+AnaEveMyTree::~AnaEveMyTree()
+{
+ if (!fChain) return;
+ delete fChain->GetCurrentFile();
+}
+
+Int_t AnaEveMyTree::GetEntry(Long64_t entry)
+{
+// Read contents of entry.
+ if (!fChain) return 0;
+ return fChain->GetEntry(entry);
+}
+Long64_t AnaEveMyTree::LoadTree(Long64_t entry)
+{
+// Set the environment to read one entry
+ if (!fChain) return -5;
+ Long64_t centry = fChain->LoadTree(entry);
+ if (centry < 0) return centry;
+ if (!fChain->InheritsFrom(TChain::Class())) return centry;
+ TChain *chain = (TChain*)fChain;
+ if (chain->GetTreeNumber() != fCurrent) {
+ fCurrent = chain->GetTreeNumber();
+ Notify();
+ }
+ return centry;
+}
+
+void AnaEveMyTree::Init(TTree *tree)
+{
+ // The Init() function is called when the selector needs to initialize
+ // a new tree or chain. Typically here the branch addresses and branch
+ // pointers of the tree will be set.
+ // It is normally not necessary to make changes to the generated
+ // code, but the routine can be extended by the user if needed.
+ // Init() will be called many times when running on PROOF
+ // (once per file to be processed).
+
+ // Set branch addresses and branch pointers
+ if (!tree) return;
+ fChain = tree;
+ fCurrent = -1;
+ fChain->SetMakeClass(1);
+
+ fChain->SetBranchAddress("jev", &jev, &b_jev);
+ fChain->SetBranchAddress("iev", &iev, &b_iev);
+ fChain->SetBranchAddress("ntr", &ntr, &b_ntr);
+ fChain->SetBranchAddress("etr", &etr, &b_etr);
+ fChain->SetBranchAddress("z1n", &z1n, &b_z1n);
+ fChain->SetBranchAddress("z2n", &z2n, &b_z2n);
+ fChain->SetBranchAddress("z1p", &z1p, &b_z1p);
+ fChain->SetBranchAddress("z2p", &z2p, &b_z2p);
+ fChain->SetBranchAddress("zvt", &zvt, &b_zvt);
+ fChain->SetBranchAddress("z1s", &z1s, &b_z1s);
+ fChain->SetBranchAddress("z2s", &z2s, &b_z2s);
+ fChain->SetBranchAddress("w1s", &w1s, &b_w1s);
+ fChain->SetBranchAddress("w2s", &w2s, &b_w2s);
+ fChain->SetBranchAddress("t1s", &t1s, &b_t1s);
+ fChain->SetBranchAddress("t2s", &t2s, &b_t2s);
+ fChain->SetBranchAddress("v1s", &v1s, &b_v1s);
+ fChain->SetBranchAddress("v2s", &v2s, &b_v2s);
+ fChain->SetBranchAddress("z1nt", z1nt, &b_z1nt);
+ fChain->SetBranchAddress("z2nt", z2nt, &b_z2nt);
+ fChain->SetBranchAddress("z1pt", z1pt, &b_z1pt);
+ fChain->SetBranchAddress("z2pt", z2pt, &b_z2pt);
+ fChain->SetBranchAddress("t0amp", t0amp, &b_t0amp);
+ fChain->SetBranchAddress("t0tim", t0tim, &b_t0tim);
+ fChain->SetBranchAddress("v0mul", v0mul, &b_v0mul);
+ Notify();
+}
+
+Bool_t AnaEveMyTree::Notify()
+{
+ // The Notify() function is called when a new file is opened. This
+ // can be either for a new TTree in a TChain or when when a new TTree
+ // is started when using PROOF. It is normally not necessary to make changes
+ // to the generated code, but the routine can be extended by the
+ // user if needed. The return value is currently not used.
+
+ return kTRUE;
+}
+
+void AnaEveMyTree::Show(Long64_t entry)
+{
+// Print contents of entry.
+// If entry is not specified, print current entry
+ if (!fChain) return;
+ fChain->Show(entry);
+}
+Int_t AnaEveMyTree::Cut(Long64_t entry)
+{
+// This function may be called from Loop.
+// returns 1 if entry is accepted.
+// returns -1 otherwise.
+ return 1;
+}
+#endif // #ifdef AnaEveMyTree_cxx
--- /dev/null
+// 115393 115401 115413 115414 115514
+// 117048 117050 117052 117053 117054
+// 117059 117060 117063 117065 117077
+// 117082 117086 117092
+AliAnalysisGrid* CreateAlienHandler()
+{
+// Check if user has a valid token, otherwise make one. This has limitations.
+// One can always follow the standard procedure of calling alien-token-init then
+// source /tmp/gclient_env_$UID in the current shell.
+// if (!AliAnalysisGrid::CreateToken()) return NULL;
+ AliAnalysisAlien *plugin = new AliAnalysisAlien();
+// Set the run mode (can be "full", "test", "offline", "submit" or "terminate")
+// plugin->SetRunMode("test");
+ plugin->SetRunMode("full");
+// plugin->SetRunMode("terminate");
+// Set versions of used packages
+ plugin->SetAPIVersion("V1.1x");
+// plugin->SetROOTVersion("v5-26-00b");
+ plugin->SetROOTVersion("v5-27-06b");
+// plugin->SetAliROOTVersion("v4-19-19-AN");
+// plugin->SetAliROOTVersion("v4-19-Rev-25");
+// plugin->SetAliROOTVersion("v4-21-05-AN");
+ plugin->SetAliROOTVersion("v4-21-06-AN");
+// Declare input data to be processed.
+// Method 1: Create automatically XML collections using alien 'find' command.
+// Define production directory LFN
+// plugin->SetGridDataDir("/alice/sim/LHC10a6");
+ // On real reconstructed data:
+ plugin->SetGridDataDir("/alice/data/2010/LHC10h");
+// Set data search pattern
+ // plugin->SetDataPattern("*ESDs.root");
+// Data pattern for reconstructed data
+// plugin->SetDataPattern("*ESDs/pass1_plusplusplus/*ESDs.root");
+ plugin->SetDataPattern("*ESDs/pass1_4plus/*ESDs.root");
+ plugin->SetRunPrefix("000"); // real data
+// ...then add run numbers to be considered
+// plugin->AddRunNumber(125020);
+// plugin->AddRunNumber(137161); // 137161 137162 real data
+ plugin->AddRunNumber(137609); // 137161 *137431 *137549 *137595 *137609 137638 137639 137693
+// plugin->SetOutputSingleFolder("output");
+// plugin->SetOutputToRunNo();
+// Method 2: Declare existing data files (raw collections, xml collections, root file)
+// If no path mentioned data is supposed to be in the work directory (see SetGridWorkingDir())
+// XML collections added via this method can be combined with the first method if
+// the content is compatible (using or not tags)
+// plugin->AddDataFile("tag.xml");
+// plugin->AddDataFile("/alice/data/2008/LHC08c/000057657/raw/Run57657.Merged.RAW.tag.root");
+// Define alien work directory where all files will be copied. Relative to alien $HOME.
+ plugin->SetGridWorkingDir("work_137609");
+// Declare alien output directory. Relative to working directory.
+ plugin->SetGridOutputDir("output"); // In this case will be $HOME/work/output
+// Declare the analysis source files names separated by blancs. To be compiled runtime
+// using ACLiC on the worker nodes.
+ plugin->SetAnalysisSource("AliAnalysisTaskPt.cxx");
+// Declare all libraries (other than the default ones for the framework. These will be
+// loaded by the generated analysis macro. Add all extra files (task .cxx/.h) here.
+ plugin->SetAdditionalLibs("AliAnalysisTaskPt.h AliAnalysisTaskPt.cxx");
+// Declare the output file names separated by blancs.
+// (can be like: file.root or file.root@ALICE::Niham::File)
+// plugin->SetDefaultOutputs(kTRUE);
+ plugin->SetDefaultOutputs(kFALSE);
+ plugin->SetOutputFiles("EveMyTree.root");
+// Optionally define the files to be archived.
+// plugin->SetOutputArchive("log_archive.zip:stdout,stderr@ALICE::NIHAM::File root_archive.zip:*.root@ALICE::NIHAM::File");
+// plugin->SetOutputArchive("log_archive.zip:mydata.dat");
+// Optionally set a name for the generated analysis macro (default MyAnalysis.C)
+ plugin->SetAnalysisMacro("TaskPt.C");
+// Optionally set maximum number of input files/subjob (default 100, put 0 to ignore)
+ plugin->SetSplitMaxInputFileNumber(100);
+// Optionally modify the executable name (default analysis.sh)
+ plugin->SetExecutableCommand("aliroot -b -q");
+ plugin->SetExecutable("TaskPt.sh");
+// Optionally set number of failed jobs that will trigger killing waiting sub-jobs.
+// plugin->SetMaxInitFailed(5);
+// Optionally resubmit threshold.
+// plugin->SetMasterResubmitThreshold(90);
+// Optionally set time to live (default 30000 sec)
+ plugin->SetTTL(30000);
+// Optionally set input format (default xml-single)
+ plugin->SetInputFormat("xml-single");
+// Optionally modify the name of the generated JDL (default analysis.jdl)
+ plugin->SetJDLName("TaskPt.jdl");
+// Optionally modify job price (default 1)
+ plugin->SetPrice(1);
+// Optionally modify split mode (default 'se')
+ plugin->SetSplitMode("se");
+ return plugin;
+}
--- /dev/null
+const char *anatype = "ESD";
+
+void TaskPt()
+{
+// Analysis using ESD data
+// Automatically generated analysis steering macro executed in grid subjobs
+
+ TStopwatch timer;
+ timer.Start();
+
+// Set temporary merging directory to current one
+ gSystem->Setenv("TMPDIR", gSystem->pwd());
+
+// include path
+ gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+
+// Load analysis framework libraries
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+ gSystem->Load("libCORRFW");
+
+// Add aditional AliRoot libraries
+
+// analysis source to be compiled at runtime (if any)
+ gROOT->ProcessLine(".L AliAnalysisTaskPt.cxx+g");
+
+// connect to AliEn and make the chain
+ if (!TGrid::Connect("alien://")) return;
+// read the analysis manager from file
+ TFile *file = TFile::Open("TaskPt.root");
+ if (!file) return;
+ TIter nextkey(file->GetListOfKeys());
+ AliAnalysisManager *mgr = 0;
+ TKey *key;
+ while ((key=(TKey*)nextkey())) {
+ if (!strcmp(key->GetClassName(), "AliAnalysisManager"))
+ mgr = (AliAnalysisManager*)file->Get(key->GetName());
+ };
+ if (!mgr) {
+ ::Error("TaskPt", "No analysis manager found in file TaskPt.root");
+ return;
+ }
+
+ mgr->PrintStatus();
+ AliLog::SetGlobalLogLevel(AliLog::kError);
+ TChain *chain = CreateChain("wn.xml", anatype);
+
+ mgr->StartAnalysis("localfile", chain);
+ timer.Stop();
+ timer.Print();
+}
+
+//________________________________________________________________________________
+TChain* CreateChain(const char *xmlfile, const char *type="ESD")
+{
+// Create a chain using url's from xml file
+ TString filename;
+ Int_t run = 0;
+ TString treename = type;
+ treename.ToLower();
+ treename += "Tree";
+ printf("***************************************\n");
+ printf(" Getting chain of trees %s\n", treename.Data());
+ printf("***************************************\n");
+ TAlienCollection *coll = TAlienCollection::Open(xmlfile);
+ if (!coll) {
+ ::Error("CreateChain", "Cannot create an AliEn collection from %s", xmlfile);
+ return NULL;
+ }
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ TChain *chain = new TChain(treename);
+ coll->Reset();
+ while (coll->Next()) {
+ filename = coll->GetTURL();
+ if (mgr) {
+ Int_t nrun = AliAnalysisManager::GetRunFromAlienPath(filename);
+ if (nrun && nrun != run) {
+ printf("### Run number detected from chain: %d\n", nrun);
+ mgr->SetRunFromPath(nrun);
+ run = nrun;
+ }
+ }
+ chain->Add(filename);
+ }
+ if (!chain->GetNtrees()) {
+ ::Error("CreateChain", "No tree found from collection %s", xmlfile);
+ return NULL;
+ }
+ return chain;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff