+++ /dev/null
-// Config file test for MUON spectormeter
-// Remember to define the directory and option
-// gAlice->SetConfigFunction("Config('$HOME','box');");
-
-void
-Config (char directory[100] = "", char option[6] = "box")
-{
- //
- // Config file for MUON test
- // Gines MARTINEZ, Subatech, mai 2003, august 2003
- //
-
- //=====================================================================
- // Libraries required by geant321
- gSystem->Load ("libgeant321.so");
- new TGeant3 ("C++ Interface to Geant3");
- //=======================================================================
- // Create the output file
- Text_t filename[100];
- sprintf (filename, "%sgalice.root", directory);
- cout << ">>> Output file is " << filename << endl;
- cout << ">>> Config_MUON_test.C: Creating Run Loader ..." << endl;
- AliRunLoader *rl = 0x0;
- rl =
- AliRunLoader::Open (filename, AliConfig::GetDefaultEventFolderName (),
- "recreate");
- if (rl == 0x0)
- {
- gAlice->Fatal ("Config_MUON_test.C",
- "Can not instatiate the Run Loader");
- return;
- }
- rl->SetCompressionLevel (2);
- rl->SetNumberOfEventsPerFile (100);
- gAlice->SetRunLoader (rl);
-
-
- //=======================================================================
- // Set External decayer
- TVirtualMCDecayer *decayer = new AliDecayerPythia ();
- decayer->SetForceDecay (kAll);
- decayer->Init ();
- gMC->SetExternalDecayer (decayer);
-
- //
- //=======================================================================
- // ******* GEANT STEERING parameters FOR ALICE SIMULATION *******
- gMC->SetProcess ("DCAY", 1);
- gMC->SetProcess ("PAIR", 1);
- gMC->SetProcess ("COMP", 1);
- gMC->SetProcess ("PHOT", 1);
- gMC->SetProcess ("PFIS", 0);
- gMC->SetProcess ("DRAY", 0);
- gMC->SetProcess ("ANNI", 1);
- gMC->SetProcess ("BREM", 1);
- gMC->SetProcess ("MUNU", 1);
- gMC->SetProcess ("CKOV", 1);
- gMC->SetProcess ("HADR", 1);
- gMC->SetProcess ("LOSS", 2);
- gMC->SetProcess ("MULS", 1);
- gMC->SetProcess ("RAYL", 1);
-
- Float_t cut = 1.e-3; // 1MeV cut by default
- Float_t tofmax = 1.e10;
-
- gMC->SetCut ("CUTGAM", cut);
- gMC->SetCut ("CUTELE", cut);
- gMC->SetCut ("CUTNEU", cut);
- gMC->SetCut ("CUTHAD", cut);
- gMC->SetCut ("CUTMUO", cut);
- gMC->SetCut ("BCUTE", cut);
- gMC->SetCut ("BCUTM", cut);
- gMC->SetCut ("DCUTE", cut);
- gMC->SetCut ("DCUTM", cut);
- gMC->SetCut ("PPCUTM", cut);
- gMC->SetCut ("TOFMAX", tofmax);
- //
- //=======================================================================
- // ************* STEERING parameters FOR ALICE SIMULATION **************
- // Chamber positions
- // From AliMUONConstants class we get :
- // Position Z (along beam) of the chambers (in cm)
- // (from AliMUONConstants class):
- // 533.5, 546.5, 678.5, 693.5, 964.0, 986.0, 1251.5, 1278.5,
- // 1416.5, 1443.5, 1610, 1625., 1710., 1725.
- // Internal Radius (in cm)
- // 36.4, 46.2, 66.0, 80., 80., 100., 100.
- // External Radius (in cm)
- // 183., 245., 395., 560., 563., 850., 900.
- //=======================================================================
- if (!strcmp (option, "box"))
- {
- AliGenBox *gener = new AliGenBox (1);
- gener->SetMomentumRange (20., 20.1);
- gener->SetPhiRange (0.,180.);
- gener->SetThetaRange (171.000, 178.001);
- gener->SetPart (13); // Muons
- gener->SetOrigin (0., 0., 0.); //vertex position
- gener->SetSigma (0.0, 0.0, 0.0); //Sigma in (X,Y,Z) (cm) on IP position
- }
- if (!strcmp (option, "gun"))
- {
- //*********************************************
- // Example for Fixed Particle Gun *
- //*********************************************
- AliGenFixed *gener = new AliGenFixed (ntracks);
- gener->SetMomentum (10);
- gener->SetPhiRange (0.);
- gener->SetThetaRange (0.);
- gener->SetOrigin (30, 30, 1200); //vertex position
- gener->SetPart (13); //GEANT particle type 13 is muons
- }
- if (!strcmp (option, "scan"))
- {
- AliGenScan *gener = new AliGenScan (-1);
- gener->SetMomentumRange (10, 10);
- gener->SetPhiRange (0, 0);
- gener->SetThetaRange (-180, -180);
- //vertex position
- //gener->SetSigma(1,1,0); //Sigma in (X,Y,Z) (cm) on IP position
- gener->SetPart (kRootino);
- gener->SetRange (100, -300., 300., 100, -300., 300., 1, 2000, 2000);
- }
- if (!strcmp (option, "param"))
- {
- //*******************************************************
- // Example for J/psi or Upsilon Production from Parameterisation *
- //*******************************************************
- AliGenParam *gener = new AliGenParam (1, AliGenMUONlib::kUpsilon);
- gener->SetMomentumRange (0, 999);
- gener->SetPtRange (0, 100.);
- gener->SetPhiRange (0., 360.);
- gener->SetCutOnChild (1);
- gener->SetChildPhiRange (0., 360.);
- gener->SetChildThetaRange (171.0, 178.0);
- gener->SetOrigin (0, 0, 0); //vertex position gener->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position
- gener->SetForceDecay (kDiMuon);
- gener->SetTrackingFlag (1);
- gener->Init ();
- }
- //=============================================================
- // Field (L3 0.4 T)
- AliMagFMaps *field =
- new AliMagFMaps ("Maps", "Maps", 1, 1., 10., AliMagFMaps::k4kG);
- gAlice->SetField (field);
-
- //=================== Alice BODY parameters =============================
- AliBODY *BODY = new AliBODY ("BODY", "Alice envelop");
- //=================== ABSO parameters ============================
- AliABSO *ABSO = new AliABSOv0 ("ABSO", "Muon Absorber");
- //=================== DIPO parameters ============================
- AliDIPO *DIPO = new AliDIPOv2 ("DIPO", "Dipole version 2");
- //================== HALL parameters ============================
- AliHALL *HALL = new AliHALL ("HALL", "Alice Hall");
- //=================== PIPE parameters ============================
- AliPIPE *PIPE = new AliPIPEv0 ("PIPE", "Beam Pipe");
- //=================== SHIL parameters ============================
- AliSHIL *SHIL = new AliSHILv2 ("SHIL", "Shielding Version 2");
- //=================== MUON Subsystem ===========================
- cout << ">>> Config_MUON_test.C: Creating AliMUONv1 ..." << endl;
-
- // Old MUONv1 version (renamed to MUONv3)
- //AliMUONv3 *MUON = new AliMUONv3("MUON","default");
-
- // New MUONv1 version (geometry defined via builders)
- AliMUON *MUON = new AliMUONv1 ("MUON", "default");
-
- //MUON->SetAlign(true);
- // If align = true, the detection elements transformations
- // are taken from the input files and not from the code
-
- //MUON->SetDebug(2);
- // To check setting of transformations from input files
- // set align = true, debug level = 2 and run with scan generator
-
- //MUON->AddGeometryBuilder(new AliMUONSt1GeometryBuilder(MUON));
- MUON->AddGeometryBuilder (new AliMUONSt1GeometryBuilderV2 (MUON));
- MUON->AddGeometryBuilder (new AliMUONSt2GeometryBuilder (MUON));
- MUON->AddGeometryBuilder (new AliMUONSlatGeometryBuilder (MUON));
- MUON->AddGeometryBuilder (new AliMUONTriggerGeometryBuilder (MUON));
-}
-
-Float_t
-EtaToTheta (Float_t arg)
-{
- return (180. / TMath::Pi ()) * 2. * atan (exp (-arg));
-}
+++ /dev/null
-#define DIFFY 0.5 //5000 micron difference in Y
-#define DIFFX 1.0 //10000 micron difference in X
-#define X 56 //No. of Columns in bending plane
-#define Y 128 //No. of Rows in nonbending plane
-#define NUM 35 //Max Padhits in the column or rows
-#define C 200 //Max padhits in one ddl
-#define D 500 //Total Number of padhits in one ddl
-
-void dhlt_v2(Int_t evNumber1=0,Int_t evNumber2=0,Int_t lutNum=13,Int_t dcCut=50)
-{
- Int_t event,dsp,hit,iqpad,header,iPx,iPy,pl,count1=0,count2=0;
- Double_t MeanTime=0.;
- Int_t Data[D];
-
- /********************************************************
- * beginning of file reading - decide which ones and open
- * them.
- *******************************************************/
-
- //Using Lookup Table
-
- if(lutNum==13||lutNum==14||lutNum==15||lutNum==16)
- {
- Int_t lineNum,xvalue[48448],yvalue[48448],plane[48448];
- Float_t xproper[48448],yproper[48448];
- char fnamelut[30];
- sprintf(fnamelut,"lut%d.dat",lutNum);
- FILE *ft = fopen(fnamelut,"r");
- for(Int_t n=0;n<48448;n++)
- {
- fscanf(ft,"%d",&lineNum);
- fscanf(ft,"%d",&xvalue[lineNum]);
- fscanf(ft,"%d",&yvalue[lineNum]);
- fscanf(ft,"%f",&xproper[lineNum]);
- fscanf(ft,"%f",&yproper[lineNum]);
- fscanf(ft,"%d",&plane[lineNum]);
- }
- fclose(ft);
- }
- else if(lutNum==17||lutNum==18||lutNum==19||lutNum==20)
- {
- Int_t lineNum,xvalue[59136],yvalue[59136],plane[59136];
- Float_t xproper[59136],yproper[59136];
- char fnamelut[30];
- sprintf(fnamelut,"lut%d.dat",lutNum);
- FILE *ft1 = fopen(fnamelut,"r");
- for(Int_t n=0;n<59136;n++)
- {
- fscanf(ft1,"%d",&lineNum);
- fscanf(ft1,"%d",&xvalue[lineNum]);
- fscanf(ft1,"%d",&yvalue[lineNum]);
- fscanf(ft1,"%f",&xproper[lineNum]);
- fscanf(ft1,"%f",&yproper[lineNum]);
- fscanf(ft1,"%d",&plane[lineNum]);
- }
- fclose(ft1);
- }
- // else
- // continue;
-
- //Declarations of arrays
-
- Float_t xdatab[X][NUM];
- Float_t ydatab[X][NUM];
- Float_t xdatan[Y][NUM];
- Float_t ydatan[Y][NUM];
-
- Int_t n3=0,n4=0;
- //Dymanic allocation of the arrays
-
- Int_t *counterx;
- counterx=(Int_t *) malloc(X*sizeof(Int_t));
- if(counterx==NULL)exit(1);
- Int_t *countery;
- countery=(Int_t *) malloc(Y*sizeof(Int_t));
- if(countery==NULL)exit(1);
- Int_t **datab;
- datab=(Int_t **) malloc(X*sizeof(Int_t *));
-
- for(n3=0;n3<X;n3++)
- {
- datab[n3]=(Int_t *) malloc(NUM*sizeof(Int_t));
- if(datab[n3]==NULL)exit(1);
- }
- Int_t **datan;
- datan=(Int_t **) malloc(Y*sizeof(Int_t *));
- for(n4=0;n4<Y;n4++)
- {
- datan[n4]=(Int_t *) malloc(NUM*sizeof(Int_t));
- if(datan[n4]==NULL)exit(1);
- }
-
- //declaration of cg arrays
-
- // C is max padhits per ddl.
-
- Float_t ycgnb[C],xcgnb[C];
- Int_t column[C],row[C],chargeb[C],chargenb[C],s1,s2;
- Float_t ycgb[C],xcgb[C];
-
- //Unformatted Binary Input Files and Output Files created
-
- char fnameubd[30];
- sprintf(fnameubd,"ddl%d.dat",lutNum);
- FILE *fp = fopen(fnameubd,"r+");
- char fnamehitb[30];
- sprintf(fnamehitb,"hit%db.dat",lutNum);
- FILE *fr = fopen(fnamehitb,"w+");
- char fnamehitn[30];
- sprintf(fnamehitn,"hit%dn.dat",lutNum);
- FILE *fs = fopen(fnamehitn,"w+");
-
- for(int nev=evNumber1;nev<=evNumber2;nev++) // Start Event Loop
- {
-
- // Initiallisation of the arrays
-
- for(Int_t m=0;m<X;m++)
- { // loop over columns in bending plane
- for(Int_t n=0;n<NUM;n++)
- { // loop over padhits in the column
- datab[m][n]=0;
- xdatab[m][n]=0;
- ydatab[m][n]=0;
- }
- }
- for(Int_t m1=0;m1<Y;m1++)
- { // loop over rows in the nonbending plane
- for(Int_t n1=0;n1<NUM;n1++)
- { // loop over padhits in the bending plane
- datan[m1][n1]=0;
- xdatan[m1][n1]=0;
- ydatan[m1][n1]=0;
- }
- }
-
- // setting counters to zero.
- // x counter and y counter are contained in the LUT's
- for(Int_t i1=0;i1<X;i1++)
- counterx[i1]=0;
-
- for(Int_t j1=0;j1<Y;j1++)
- countery[j1]=0;
- // Data has the size of the total number of padhits in a DDL.
- for(Int_t h=0;h<D;h++)
- Data[h]=0;
-
- // initialising binary data c.o.g
- for(Int_t tsn=0;tsn<C;tsn++)
- {
- ycgnb[tsn]=0;
- xcgnb[tsn]=0;
- row[tsn]=0;
- chargenb[tsn]=0;
- }
-
- for(Int_t tsb=0;tsb<C;tsb++)
- {
- ycgb[tsb]=0;
- xcgb[tsb]=0;
- column[tsb]=0;
- chargeb[tsb]=0;
- }
-
- // Time Starts for the algorithm
-
- TStopwatch timer;
- timer.Start();
-
- // Reading the unformatted data stream from ddl(DDL#).dat generated by
- // rawddl(STATION#).C from digitized data from AliRoot.
-
- // this should probably be changed to case...
- Int_t initial=0,final=0;
- if(lutNum==13)
- {initial=501;final=512;}
- if(lutNum==14)
- {initial=521;final=532;}
- if(lutNum==15)
- {initial=541;final=552;}
- if(lutNum==16)
- {initial=561;final=572;}
- if(lutNum==17)
- {initial=601;final=613;}
- if(lutNum==18)
- {initial=621;final=633;}
- if(lutNum==19)
- {initial=641;final=653;}
- if(lutNum==20)
- {initial=661;final=673;}
-
- for(Int_t p=initial;p<final;p++)
- {//dsp loop start
- fread(&dsp,2,1,fp);
- fread(&hit,2,1,fp); // size_t fread(void *ptr, size_t size, size_t nmemb, FILE *stream);
-
- //reading unformatted binary data stream
-
- for(Int_t g=0;g<hit;g++)
- { // DDL data
- fread(&Data[g],4,1,fp);
- iqpad = (((Data[g]<<17)>>17)&4095); // charge on pad first 12 bits
- if(iqpad>dcCut) // threshold
- {
- header = (((Data[g])>>12)&65535); // next 17 bits
- iPx = xvalue[header]; // ix (logical position)
- iPy = yvalue[header]; // ix (logical position)
- pl=plane[header];
- if(iqpad>=4096)
- iqpad = 4095;
- iPx=abs(iPx);
- if(pl==0)
- {
- count1 = counterx[iPx];
- datab[iPx][count1]=4096*iPy+iqpad; // new data structure: column-wise
- // in order to optimize the search
- // for central pads
- xdatab[iPx][count1]=xproper[header];
- ydatab[iPx][count1]=yproper[header];
- counterx[iPx]+=1;
- }
- else
- {
- count2 = countery[iPy];
- datan[iPy][count2]=4096*iPx+iqpad; // new data structure: column-wise
- // in order to optimize the search
- // for central pads
- xdatan[iPy][count2]=xproper[header];
- ydatan[iPy][count2]=yproper[header];
- countery[iPy]+=1;
- }
- }// DcCut condition
- else continue;
- }//dsp counter loop
- }//dsp loop
-
- s1=0;s2=0;
-
- // Reconstruction of the bending output array-streams
-
- for(Int_t b1=0;b1<X;b1++)
- { //Loop on Columns
- if(counterx[b1]>0){ // Condition on Empty Columns
- Int_t iq=0;
- for(Int_t a1=0;a1<counterx[b1];a1++)
- { // A pad is obtained
- if(datab[b1][a1]<18388608)
- {
- iq = (((*(*(datab+b1)+a1)<<10)>>10)&4095);
- Int_t iyc1 = ((*(*(datab+b1)+a1)>>12)&1023);
-
- // Checking Neibouring Pads of the first pad
- Int_t iq1=0,iq2=0,byp1=0,byn1=0,byp2=0,byn2=0;
- for(Int_t c1=0;c1<counterx[b1];c1++)
- {
- if(((*(*(datab+b1)+c1)>>12)&1023)==iyc1+1)
- {
- iq1 = (((*(*(datab+b1)+c1)<<10)>>10)&4095);
- byp1=c1;
- }
- if(((*(*(datab+b1)+c1)>>12)&1023)==iyc1-1)
- {
- iq2 = (((*(*(datab+b1)+c1)<<10)>>10)&4095);
- byn1=c1;
- }
- if(((*(*(datab+b1)+c1)>>12)&1023)==iyc1+2)
- {
- byp2=c1;
- }
- if(((*(*(datab+b1)+c1)>>12)&1023)==iyc1-2)
- {
- byn2=c1;
- }
- } // Neibouring pads
-
- // Reconstruction of 3pad hits
- if(iq1>0 && iq2>0 && iq>=iq1 && iq>=iq2)
- {
- ycgb[s1] = funcy(ydatab[b1][a1],iq,iq1,iq2);
- xcgb[s1] = xdatab[b1][a1];
- column[s1]=b1;
- chargeb[s1]=iq;
- s1+=1;
- datab[b1][a1]+=18388608;
- datab[b1][byp1]+=18388608;
- datab[b1][byn1]+=18388608;
- if(byp2>0)
- datab[b1][byp2]+=18388608;
- if(byn2>0)
- datab[b1][byn2]+=18388608;
- }// end 3pad hits
-
- // Reconstruction of 2pad hits
- else if(iq2==0 && iq1>0 && iq>=iq1)
- {
- ycgb[s1] = funcy(ydatab[b1][a1],iq,iq1,iq2);
- xcgb[s1] = xdatab[b1][a1];
- column[s1]=b1;
- chargeb[s1]=iq;
- s1+=1;
- datab[b1][a1]+=18388608;
- datab[b1][byp1]+=18388608;
- if(byp2>0)
- datab[b1][byp2]+=18388608;
- }//end 2pad hits1
- else if(iq1==0 && iq2>0 && iq>=iq2)
- {
- ycgb[s1] = funcy(ydatab[b1][a1],iq,iq1,iq2);
- xcgb[s1] = xdatab[b1][a1];
- column[s1]=b1;
- chargeb[s1]=iq;
- s1+=1;
- datab[b1][a1]+=18388608;
- datab[b1][byn1]+=18388608;
- if(byn2>0)
- datab[b1][byn2]+=18388608;
- }//end 2pad hits2
- else continue;
- }//dccut condition
- else continue;
- }//loop 2
- }//end condition1
- }// end loop on columns
- //Marging two consecutive columns in the bending plane
-
- for(Int_t z1b=0;z1b<s1;z1b++)
- {
- for(Int_t z2b=z1b+1;z2b<s1;z2b++)
- {
- if((column[z2b]==column[z1b]+1)&&(abs(ycgb[z1b]-ycgb[z2b])<=DIFFY))
- {
- if(chargeb[z1b]>chargeb[z2b])
- {
- ycgb[z1b]=(ycgb[z1b]*chargeb[z1b]+ycgb[z2b]*chargeb[z2b])/(chargeb[z1b]+chargeb[z2b]);
- xcgb[z1b]=(xcgb[z1b]*chargeb[z1b]+xcgb[z2b]*chargeb[z2b])/(chargeb[z1b]+chargeb[z2b]);
- column[z2b]+=2000;
- }
- else
- {
- ycgb[z2b]=(ycgb[z1b]*chargeb[z1b]+ycgb[z2b]*chargeb[z2b])/(chargeb[z1b]+chargeb[z2b]);
- xcgb[z2b]=(xcgb[z1b]*chargeb[z1b]+xcgb[z2b]*chargeb[z2b])/(chargeb[z1b]+chargeb[z2b]);
- column[z1b]+=2000;
- }
- }
- else if(column[z2b]>column[z1b]+1)
- {
- break;
- }
- else continue;
- } //2nd loop
- }//1st loop
- for(Int_t z3b=0;z3b<s1;z3b++)
- {
- if(column[z3b]<2000)
- {
- fprintf(fr,"%4.6f\t%4.2f\n",ycgb[z3b],xcgb[z3b]); // reconstructed x and
- // y positons for a cluster
- }
- }
-
- // Reconstruction of the non-bending output array-streams
-
- for(Int_t b2=0;b2<Y;b2++)
- { //Loop on Columns
- if(countery[b2]>0)
- { // Condition on Empty Columns
- Int_t iq=0;
- for(Int_t a=0;a<countery[b2];a++)
- { // A pad is obtained
- if(datan[b2][a]<18388608)
- {
- iq = (((*(*(datan+b2)+a)<<10)>>10)&4095);
- Int_t iyc = ((*(*(datan+b2)+a)>>12)&1023);
-
- // Checking Neibouring Pads of the first pad
- Int_t iq1=0,iq2=0,nyp1=0,nyn1=0,nyp2=0,nyn2=0;
- for(Int_t c=0;c<countery[b2];c++)
- {
- if(((*(*(datan+b2)+c)>>12)&1023)==iyc+1)
- {
- iq1 = (((*(*(datan+b2)+c)<<10)>>10)&4095);
- nyp1=c;
- }
- if(((*(*(datan+b2)+c)>>12)&1023)==iyc-1)
- {
- iq2 = (((*(*(datan+b2)+c)<<10)>>10)&4095);
- nyn1=c;
- }
- if(((*(*(datan+b2)+c)>>12)&1023)==iyc+2)
- {
- nyp2=c;
- }
- if(((*(*(datan+b2)+c)>>12)&1023)==iyc-2)
- {
- nyn2=c;
- }
- } // pads
-
- // Reconstruction of 3pad hits
- if(iq1>0 && iq2>0 && iq>=iq1 && iq>=iq2)
- {
- xcgnb[s2] = funcx(xdatan[b2][a],iq,iq1,iq2);
- ycgnb[s2] = ydatan[b2][a];
- row[s2]=b2;
- chargenb[s2]=iq;
- s2+=1;
- datan[b2][a]+=18388608;
- datan[b2][nyp1]+=18388608;
- datan[b2][nyn1]+=18388608;
- if(nyp2>0)
- datan[b2][nyp2]+=18388608;
- if(nyn2>0)
- datan[b2][nyn2]+=18388608;
- }// end 3pad hits
-
- // Reconstruction of 2pad hits
- else if(iq2==0 && iq1>0 && iq>=iq1)
- {
- xcgnb[s2] = funcx(xdatan[b2][a],iq,iq1,iq2);
- ycgnb[s2] = ydatan[b2][a];
- row[s2]=b2;
- chargenb[s2]=iq;
- s2+=1;
- datan[b2][a]+=18388608;
- datan[b2][nyp1]+=18388608;
- if(nyp2>0)
- datan[b2][nyp2]+=18388608;
- }//end 2pad hits1
- else if(iq1==0 && iq2>0 && iq>=iq2)
- {
- xcgnb[s2] = funcx(xdatan[b2][a],iq,iq1,iq2);
- ycgnb[s2] = ydatan[b2][a];
- row[s2]=b2;
- chargenb[s2]=iq;
- s2+=1;
- datan[b2][a]+=18388608;
- datan[b2][nyn1]+=18388608;
- if(nyn2>0)
- datan[b2][nyn2]+=18388608;
- }//end 2pad hits2
- else continue;
- }//dccut condition
- else continue;
- }//loop 2
- }//end condition1
- }// end loop on rows
- //Marging two consecutive rows in the non bending plane
-
- for(Int_t z1n=0;z1n<s2;z1n++)
- {
- for(Int_t z2n=z1n+1;z2n<s2;z2n++)
- {
- if((row[z2n]==row[z1n]+1)&&(abs(xcgnb[z1n]-xcgnb[z2n])<=DIFFX))
- {
- if(chargenb[z1n]>chargenb[z2n])
- {
- ycgnb[z1n]=(ycgnb[z1n]+ycgnb[z2n])/2;
- xcgnb[z1n]=(xcgnb[z1n]+xcgnb[z2n])/2;
- row[z2n]+=2000;
- }
- else
- {
- ycgnb[z2n]=(ycgnb[z1n]+ycgnb[z2n])/2;
- xcgnb[z2n]=(xcgnb[z1n]+xcgnb[z2n])/2;
- row[z1n]+=2000;
- }
- }
- else if(row[z2n]>row[z1n]+1)
- {
- break;
- }
- else continue;
- } //2nd loop
- }//1st loop
- for(Int_t z3n=0;z3n<s2;z3n++)
- {
- if(row[z3n]<2000)
- {
- fprintf(fs,"%4.2f\t%4.6f\n",ycgnb[z3n],xcgnb[z3n]); // reconstructed x and
- // y positons for a cluster
- }
- }
-
- timer.Stop();
- Double_t tch=timer.CpuTime();
- tch=timer.CpuTime();
- MeanTime += tch;
-
- }//event
-
- for(n3=0;n3<X;n3++)
- {
- free(datab[n3]);
- }
- free(datab);
-
- for(n4=0;n4<Y;n4++)
- {
- free(datan[n4]);
- }
- free(datan);
-
- free(counterx);
- free(countery);
-
- printf("Average processing time per event = %f s\n",(MeanTime/(evNumber2+1)));
-
- fclose(fp);
- fclose(fr);
- fclose(fs);
-}// end of funtion
-
-// Functions Calculating YCG from iy values
-
-Float_t funcy(Float_t y,Int_t q,Int_t q1,Int_t q2)
-{
- Float_t ycg1;
- Int_t tot_ch;
- ycg1 = q*y + q1*(y+0.5) + q2*(y-0.5);
- tot_ch = q + q1 + q2;
- ycg1 /= tot_ch;
- Float_t ycorr = 0.0255*(TMath::Sin(12.56637*((ycg1-y)-0.25)));
- ycg1-=ycorr;
- return(ycg1);
-}
-
-// Functions Calculating XCG from ix values
-
-Float_t funcx(Float_t x,Int_t q,Int_t q1,Int_t q2)
-{
- Float_t xcg1;
- Int_t tot_ch;
- xcg1 = q*x + q1*(x+0.5) + q2*(x-0.5);
- tot_ch = q + q1 + q2;
- xcg1 /= tot_ch;
- return(xcg1);
-}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff